Nonlinearity in Data with Gaps: Application to Ecological and Meteorological Datasets

TL;DR

This paper investigates how data gaps affect the measurement of multifractality in ecological and meteorological datasets, demonstrating that certain gap ranges do not distort nonlinearity quantification without interpolation.

Contribution

It quantifies the impact of data gaps on multifractal spectra and identifies gap ranges that preserve nonlinear measures without interpolation, applied to real-world datasets.

Findings

Identified tolerable gap ranges for multifractal analysis

Established multifractality due to nonlinearity in datasets with gaps

Demonstrated robustness of multifractal measures in ecological and meteorological data

Abstract

Datagaps are ubiquitous in real world observational data. Quantifying nonlinearity in data having gaps can be challenging. Reported research points out that interpolation can affect nonlinear quantifiers adversely, artificially introducing signatures of nonlinearity where none exist. In this paper we attempt to quantify the effect that datagaps have on the multifractal spectrum ($f(\alpha)$), in the absence of interpolation. We identify tolerable gap ranges, where the measures defining the $f(\alpha)$ curve do not show considerable deviation from the evenly sampled case. We apply this to the multifractal spectra of multiple data-sets with missing data from the SMEAR database. The datasets we consider include ecological datasets from SMEAR I, namely CO$_2$ exchange variation, photosynthetically active radiation levels and soil moisture variation time series, and meteorological datasets from SMEAR II, namely dew point variation and air temperature variation. We could establish multifractality due to deterministic nonlinearity in two of these datasets, even in the presence of gaps.

Figures

Tables

Table 1. Table 1 : α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} values for Lorenz, Rössler, Red noise, White noise and their surrogates. α 1 d superscript subscript 𝛼 1 𝑑 \alpha_{1}^{d} and α 2 d superscript subscript 𝛼 2 𝑑 \alpha_{2}^{d} represent the values for data while α 1 s superscript subscript 𝛼 1 𝑠 \alpha_{1}^{s} and α 2 s superscript subscript 𝛼 2 𝑠 \alpha_{2}^{s} represent the corresponding values for surrogates.

Data	${\alpha }_1^d$	${\alpha }_2^d$	${\alpha }_1^s$	${\alpha }_2^s$
Lorenz	1.54	2.84	2.55$\pm$0.04	3.58$\pm$0.04
Rössler	1.60	2.99	1.95$\pm$0.03	3.44$\pm$0.19
Red Noise	2.43	3.68	2.41$\pm$0.04	3.81$\pm$0.07
White Noise	2.93	3.09	2.90$\pm$0.03	3.13$\pm$0.03

Table 2. Table 2 : α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} values for Rössler and its surrogates for two different values of m s subscript 𝑚 𝑠 m_{s} and m p subscript 𝑚 𝑝 m_{p} . We notice that in both cases, α 1 d superscript subscript 𝛼 1 𝑑 \alpha_{1}^{d} and α 1 s superscript subscript 𝛼 1 𝑠 \alpha_{1}^{s} remains different from each other, whereas α 2 d superscript subscript 𝛼 2 𝑑 \alpha_{2}^{d} merges with α 2 s superscript subscript 𝛼 2 𝑠 \alpha_{2}^{s} .

$m_s$	$m_p$	${\alpha }_1^d$	${\alpha }_2^d$	${\alpha }_1^s$	${\alpha }_2^s$
8$\tau$	7$\tau$	1.73	2.98	2.02$\pm$0.03	3.06$\pm$0.18
1$\tau$	7$\tau$	1.81	4.04	1.95$\pm$0.07	3.82$\pm$0.43

Table 3. Table 3 : τ 𝜏 \tau , m s subscript 𝑚 𝑠 m_{s} , m p subscript 𝑚 𝑝 m_{p} and D 2 subscript 𝐷 2 D_{2} for time series of CO 2 exchange, Photosynthetically Active Radiation (PAR) and Soil Moisture (SM).

Data	$\tau \left(hours\right)$	M${}_s{}^{}\left(hours\right)$	M${}_p{}^{}\left(hours\right)$	$D_2$
CO2 Ex.	5	119.6	157.3	3.18
PAR	7	119.6	157.3	4.08
SM	970	499.5	971.1	1.15

Table 4. Table 4 : α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} values for CO 2 exchange, photosynthetically active radiation, soil moisture and their surrogates. α 1 d superscript subscript 𝛼 1 𝑑 \alpha_{1}^{d} and α 2 d superscript subscript 𝛼 2 𝑑 \alpha_{2}^{d} are distinctly different from α 1 s superscript subscript 𝛼 1 𝑠 \alpha_{1}^{s} and α 2 s superscript subscript 𝛼 2 𝑠 \alpha_{2}^{s} for CO 2 exchange and PAR, while they are within errors of each other for SM. This could be either because soil moisture has no underlying deterministic nonlinearity, or because of the value of m s subscript 𝑚 𝑠 m_{s} and m p subscript 𝑚 𝑝 m_{p} lying in the identified critical region.

Data	${\alpha }_1^d$	${\alpha }_2^d$	${\alpha }_1^s$	${\alpha }_2^s$
CO2 Ex.	1.87	2.26	2.55$\pm$0.01	4.75$\pm$0.01
PAR	2.10	3.08	3.39$\pm$0.07	3.88$\pm$0.09
SM	1.45	2.91	1.35$\pm$0.02	2.88$\pm$0.03

Table 5. Table 5 : τ 𝜏 \tau , m s subscript 𝑚 𝑠 m_{s} , m p subscript 𝑚 𝑝 m_{p} and D 2 subscript 𝐷 2 D_{2} for time series of dew point and air temperature variation.

Data	$\tau \left(hours\right)$	M${}_s{}^{}\left(hours\right)$	M${}_p{}^{}\left(hours\right)$	$D_2$
Dew Pt.	1000	219.2	252.7	2.63
Air T.	667	483.2	535.1	3.28

Table 6. Table 6 : α 1 subscript 𝛼 1 \alpha_{1} , α 2 subscript 𝛼 2 \alpha_{2} values for dew point and air temperature variation time series. α 1 d superscript subscript 𝛼 1 𝑑 \alpha_{1}^{d} and α 2 d superscript subscript 𝛼 2 𝑑 \alpha_{2}^{d} are close to α 1 s superscript subscript 𝛼 1 𝑠 \alpha_{1}^{s} and α 2 s superscript subscript 𝛼 2 𝑠 \alpha_{2}^{s} for both time series. In case of dew point, we suspect this could be as a result of no underlying deterministic nonlinearity. However for air temperature the value of m s subscript 𝑚 𝑠 m_{s} and m p subscript 𝑚 𝑝 m_{p} lying in the identified critical region, may also be responsible.

Data	${\alpha }_1^d$	${\alpha }_2^d$	${\alpha }_1^s$	${\alpha }_2^s$
Dew Pt.	2.19	4.25	2.31$\pm$0.06	4.25$\pm$0.05
Air T.	2.29	3.79	1.95$\pm$0.01	4.16$\pm$0.04

Equations

$$G_S(s;m_s,{\omega }_s)=\frac{1}{{\omega }_s\sqrt{2\pi }}{e}^{-{\left(s-m_s\right)}^2/2{\omega }_s^2}$$

$$G_P(p;m_p,{\omega }_p)=\frac{1}{{\omega }_p\sqrt{2\pi }}{e}^{-{\left(p-m_p\right)}^2/2{\omega }_p^2}$$

$$p_i(R)=R^{{\alpha }_i(R)}$$

$$n(\alpha ,R)\propto R^{-f(\alpha )}$$

$$f(\alpha )=A(\alpha -{\alpha }_1{)}^{{\gamma }_1}({\alpha }_2-\alpha {)}^{{\gamma }_2}$$

Full text

\corres

g.ambika@iiserpune.ac.in

2017

\lp8

1 Introduction

The presence of gaps in observational data is one of the major problems that affect real world datasets. Such gaps or the consequent uneven sampling in time series can be a major impediment in detecting nonlinearity in them. The most common way to deal with datagaps has been to interpolate through them to make an evenly sampled dataset. Past studies have shown that interpolation can have adverse effects on nonlinear time series quantifiers[1]. Some recent studies have been conducted on the effect datagaps have on nonlinear quantifiers when the analyzed time series is not interpolated through, which showed that there are certain regions of gap size and frequency where correlation dimension values are affected more[2]. We propose to extend these studies to the multifractal spectrum. Such an analysis helps us identify how susceptible the multifractal spectrum is to gaps in data.

Datagaps are common in real world datasets. Astrophysical data, ecological data, meteorological data etc are often plagued by gaps[3, 4]. They may occur due to multiple reasons like instrumental failure, cloud cover, adverse weather conditions etc. Most nonlinear quantifiers require large evenly sampled datasets in order to make predictions[5, 6]. One can exploit the presence of multiple smaller evenly sampled segments in order to calculate quantifiers like bicoherence[7]. However no analysis has been done so far on the effect of datagaps on the multifractal spectra. Such an analysis is extremely important in the context of quantifying the multifractal nature of the large number of datasets that come with missing data.

The multifractal spectrum, $f(\alpha)$, is a detailed characterization of the fractal structure of the system in a reconstructed phase space. Unlike an average measure like the correlation dimension, $D_{2}$, the $f(\alpha)$ spectrum takes into consideration the local contributions of different regions of the attractor. The idea of $f(\alpha)$ spectrum can be applied beyond attractors in phase space, to describe many natural objects[8]. A complete characterization of the multifractal curve was shown to be achieved with a set of four parameters[9]. In this study we aim to check how much variation these parameters show when the time series used for embedding contains missing data.

To determine the susceptibility or resilience of the $f(\alpha)$ spectrum to datagaps, we first start with an evenly sampled time series of a standard system. We subject the evenly sampled time series to multifractal analysis by using the algorithm described in [9]. We then remove data from the evenly sampled time series, progressively. The size of the data removed and frequency of the data removal are drawn from two Gaussian distributions. The time series with datagaps are also subject to multifractal analysis. The deviation of the multifractal quantifiers like the width of the spectrum, from the evenly sampled value is treated as an index of susceptibility of that quantifier to datagaps. We can hence identify the regions of gap size and gap frequency where the quantifiers can give reliable conclusions from the nature of the $f(\alpha)$ curve.

As an illustration, we proceed to analyze real world datasets with missing data. The datasets we consider are obtained from the Station for Measuring Forest

Ecosystem-Atmosphere Relation (SMEAR) in Finland [10]. We consider ecological and meteorological datasets obtained from SMEAR I and SMEAR II stations respectively. The ecological datasets we consider are related to photosynthesis levels in pine trees in the Värriö forest. We look into the multifractal spectrum of CO2 exchange in pine shoots, which acts as a proxy for photosynthesis. We also obtain the $f(\alpha)$ curves for the time series of levels of photosynthetically active radiation and soil moisture, both of which are known to influence photosynthesis rates[11, 12]. For the meteorological datasets we consider the time series of air temperature and dew point, obtained from the SMEAR II station in the Hyytiälä forest. Both the time series give a continuous $f(\alpha)$ curve, indicative of multifractality in its phase space structure.

2 Analysis of synthetic datasets

In this section we consider the effect of datagaps on the structure of the $f(\alpha)$ curves for standard nonlinear systems like Rössler and Lorenz systems. We start from large evenly sampled datasets of the $x$ variable of these systems and introduce gaps into the data. Since the gaps arise from multiple independent sources, we take that both the position and size of the gaps follow Gaussian distributions,

[TABLE]

[TABLE]

$G_{S}$ with mean, $m_{s}$ and standard deviation, $\omega_{s}$ determines the size of the gap. Gaussian, $G_{P}$ with mean $m_{p}$ and standard deviation $\omega_{p}$ determines the position of the gap.

Starting from any point in the time series, $G_{P}$ determines the position of the next gap, i.e. it determines the extend of gapless data. Hence it is an inverse measure of frequency of gaps. $G_{S}$ determines how big the gap is at a particular position.

Effect of datagaps on $f(\alpha)$ spectrum

A detailed characterization of the non uniformities of the attractor is captured by the generalized dimensions, $D_{q}$ [8, 9]. In an analogous formulation, if one covers the attractor with boxes of size R, the probability $p_{i}$ for points to fall inside the $i^{th}$ box, is found to scale as

[TABLE]

Then the number of such boxes that have $\alpha$ between $\alpha$ and $\alpha+\Delta\alpha$ is[13]

[TABLE]

$D_{q}$ and $f(\alpha$) are related through Legendre transformations. This provides a method to compute the $f(\alpha)$ spectrum numerically from the $D_{q}$ values. We can characterize the $f(\alpha)$ curve using the following function fit

[TABLE]

and $\alpha_{1}$, $\alpha_{2}$, $\gamma_{1}$ and $\gamma_{2}$ serve as unique parameters that can characterize the $f(\alpha)$ completely[9]. The $f(\alpha)$ spectrum finds widespread use across a variety of different fields[14, 15, 16, 17].

In this section, we consider the variation of the $\Delta\alpha$ = $\alpha_{2}$-$\alpha_{1}$, $\gamma_{2}$ and $\gamma_{1}$ with varying gap size, $m_{s}$ and position, $m_{p}$. As a sample case, we show the calculated $f(\alpha)$ spectrum for the evenly sampled time series and two cases with datagaps in Fig. 1, for the Rössler system. The full variation of $\Delta\alpha$ with $m_{s}$ and $m_{p}$ for the Rössler system is shown in Fig. 2, where we plot the relative variation, $\delta\Delta\alpha^{rel}$ = $\frac{|\Delta\alpha-\Delta\alpha_{E}|}{\Delta\alpha_{E}}$. For a fixed value of $m_{p}$, as we increase $m_{s}$, the value of $\Delta\alpha$ deviates heavily from the evenly sampled value and reaches a peak at around 1$\tau$. On increasing $m_{s}$ further, it falls again and saturates at a value close to the evenly sampled value. We find that in the critical region identified, the f($\alpha$) spectrum widens considerably for the Rössler system. Since the shape of the curve changes considerably when gaps are present, we find that $\gamma$ values fluctuate heavily. The variation of $\delta\gamma^{rel}$ = $\frac{|\gamma-\gamma_{E}|}{\gamma_{E}}$ values with $m_{s}$ and $m_{p}$ is shown in Fig. 2(b) and 2(c). The heavy fluctuations in the figure makes it clear that $\gamma$ values derived from a time series contaminated with gaps cannot be trusted, in general.

The analysis was conducted for the Lorenz system as well, where the variation is much smaller in comparison. Hence we conclude that while the width of the $f(\alpha)$ spectrum may remain broadly unchanged, except at small $m_{p}$ and for $m_{s}\approx 1\tau$, the overall shape of the curve changes, as signified by the change in the $\gamma$ values. We point out that widening of the $f(\alpha)$ curve can also happen due to noise contamination[9].

Surrogate Analysis

A wide $f(\alpha)$ spectrum is indicative of multifractality in the system. However this is not conclusive evidence of deterministic nonlinearity. Linear stochastic processes were shown to have saturating $D_{2}$ and wide multifractal curves in literature[18, 19]. A surrogate analysis is required to confirm the existence of multifractality due to deterministic nonlinearity in the system. For this we initially produce 5 Iterative Amplitude Adjusted Fourier Transform (IAAFT) surrogates of the evenly sampled data implemented using the TISEAN package[20, 21]. The $f(\alpha)$ curves obtained from the Lorenz system, Rössler system, red noise, white noise and their surrogates are shown in Fig. 3. The $\alpha_{1}$ and $\alpha_{2}$ values in all the cases studied are tabulated in Table 1. It is clear that for the deterministic systems, the data and surrogates have clearly differing $f(\alpha)$ curves, while for noise, they are not distinguishable.

We repeat the analysis using data with gaps. The same profile of gaps is used to remove data from the surrogates. We find that in the presence of gaps, the original data and surrogates become increasingly difficult to differentiate. In Fig. 4 we show function fits for two cases of the Rössler system from a region where the gaps are tolerable and where gaps cause high variation. We see that data merges with surrogates in the region where gaps are not tolerable. This is also evident from Table 2 which lists the values for $\alpha_{1}$ and $\alpha_{2}$ for data and the corresponding values averaged over five surrogate time series. Hence we conclude that in the regions of gap distribution we identify to be non tolerable, multifractality alone cannot confirm deterministic nonlinearity.

3 Analysis of datasets with gaps

The analysis conducted in the sections above are with the aim of studying variations in $f(\alpha)$ in the presence of datagaps and identifying what the tolerable regions for gap size and frequency are, for the validity of the conclusions drawn from the spectrum to hold. We now apply the analysis to real world datasets. We use datasets with missing data, from the Station for Measuring Forest Ecosystem-Atmosphere Relation (SMEAR), Finland. We look into the time series of air temperature and dew point data measured by the SMEAR II station in the Hyytiälä forest, and datasets of variation of CO2 exchange, photosynthetically active radiation and soil moisture from the SMEAR I station situated in the Värriö forest.

In all the datasets, we calculate the correlation time, $\tau$ and estimate $m_{s}$ and $m_{p}$ in terms of $\tau$. In order to find the dimension used to embed the attractor in, we calculate saturated $D_{2}$ for the dataset and pick the smallest integer greater than $D_{2}$ as the embedding dimension. The calculation of $D_{2}$ filters out all non saturating datasets (i.e. noise processes). Hence the aim of calculating the $f(\alpha)$ spectrum is not just to detect nonlinearity, but to characterize the multifractal properties of the attractor.

Datasets from SMEAR I

The SMEAR I station in Värriö, Finland measures a number of meteorological and ecological parameters. The station is situated about 200 km north of the Arctic circle. The datasets produced by the instruments are however not evenly sampled and contain datagaps[4]. Some of the primary reasons for gaps in data in this context include thunderstorms and breaks in electricity [22].

In this work we consider the time series of CO2 exchange in pine trees, which acts as a proxy for photosynthesis rate[23] and two factors that affect the photosynthesis rate, namely variation of photosynthetically active radiation on top of the measurement chamber and variation of soil moisture. The rate of photosynthesis has already been suspected of chaotic variation by multiple authors[24, 25]. However, to the best of our knowledge no study has been conducted to determine the multifractal characteristics of its embedded phase space structure. Photosynthesis is well known to be dependent on ambient temperature, concentration of atmospheric CO2, amount of photosynthetically active radiation(PAR), soil moisture (SM) etc[11, 12]. In this context the study of the nonlinear properties of these variables become important as well. We also study the multifractal properties of the embedded attractors of the time series of variation of photosynthetically active radiation and soil moisture.

All three time series are averaged over half an hour for a period of 12 years starting from 2005 to 2017 and the $m_{s}$, $m_{p}$, $\tau$ and $D_{2}$ values are shown in Table 3.

We find that barring the time series for soil moisture, most of the $m_{s}$ and $m_{p}$ are outside the identified critical region. Hence in the $f(\alpha)$ computed, $\alpha_{1}$ and $\alpha_{2}$ are likely to be largely unaffected. In the case of soil moisture it is possible that the width of the curve in the absence of datagaps, may be much smaller than the one calculated. The plots of the $f(\alpha)$ spectra are shown in Fig. 5.

We then proceed to do surrogate analysis to check for signs of deterministic nonlinearity. Since the time series has missing data we cannot directly apply the method of IAAFT mentioned above. Instead we take evenly sampled segments of the data to construct surrogates. We first take segments of data that are continuously sampled. Surrogates are constructed for each segment, and joined together to act as a surrogate for the whole data. We caution that this may not always be possible and is dependent upon the availability of gap-less chunks in the the dataset. As before we find the $\alpha_{1}$ and $\alpha_{2}$ values for data and surrogates, presented in Table 4. We find that for CO2 exchange and photosynthetically active radiation have values of $\alpha_{1}$ and $\alpha_{2}$ distinctly different from their surrogates, suggesting that multifractality due to underlying deterministic nonlinearity can be concluded. This adds weight to the claims in literature that photosynthesis may be varying chaotically[25, 24]. No such conclusion can be drawn for the variation of soil moisture, as the data and surrogates lie very close to each other. This could be as a result of the value of $m_{s}$ and $m_{p}$ lying very close to the identified critical range or due to a lack of deterministic nonlinearity or noise contamination.

Datasets from SMEAR II

The SMEAR II station is established in the Hyytiälä forest in Finland and measures atmospheric aerosols, eco-physiology, soil and water measurements, solar and terrestrial radiation and meteorological measurements. In this work we primarily consider two meteorological datasets, the dew-point and air temperature. These were chosen partially due to the high instance of datagaps in them[4]. Climate time series have been often subject to studies of multifractality, in the past[14, 16]. We consider the data sampled every half hour over the period from 2008-2017. As before we first determine the $m_{s}$, $m_{p}$, $\tau$ and $D_{2}$, shown in Table 5.

We notice that $m_{s}$ for the air temperature is in the critical region identified and hence the $f(\alpha)$ curve may be significantly affected. The falpha curves for both are plotted in Fig. 6.

With surrogates constructed as before, we perform a surrogate analysis for datasets from SMEAR II. The results are presented in Table 6. For the dew point variation data, both $\alpha_{1}$ and $\alpha_{2}$ lie close to the value for the surrogate datasets. Since the $m_{s}$ and $m_{p}$ values are away from the identified critical region, it appears that dew point variation has no underlying deterministic nonlinearity and is probably a linear stochastic process. The air temperature time series has $\alpha_{1}$ and $\alpha_{2}$ values close to the surrogates, but not within error bars. This could be as a result of the value of $m_{s}$ and $m_{p}$ lying in the identified critical region.

4 Conclusion

We investigate the effect of gaps in data on the $f(\alpha)$ spectrum of the reconstructed attractor. We use data from two standard systems, the Rössler and Lorenz systems, to check the effect of datagaps on the computed multifractal spectrum. From the analysis on standard chaotic data we could arrive at tolerable limits on gap parameters, within which conclusions from nonlinear time series analysis can be meaningful. Since these limits are linked to the embedding technique itself, we take that the same can be valid for real world datasets analyzed using the same procedure of embedding. In this work, we illustrate this by analyzing 5 ecological and meteorological datasets with gaps of varying ranges. We support our conclusions by carrying out detailed surrogate analysis on all of them.

We specifically study the variations in the four quantifiers of $f(\alpha)$, viz. $\alpha_{1}$, $\alpha_{2}$, $\gamma_{1}$ and $\gamma_{2}$, in all the datasets. To the best of our knowledge, our work for the first time establishes multifractality due to deterministic nonlinearity in photosynthesis variation and factors affecting it.

The absence of large evenly sampled datasets is a major impediment for calculating measures applying nonlinear time series analysis techniques and especially deriving quantifiers from $f(\alpha)$. The dangers of interpolation makes it an undesirable solution in these cases. Hence it becomes important to consider how much deviation from the evenly sampled case, actually happens if we ignore gaps. We find that certain broad conclusions about the data can still be retained even in the presence of missing data. We suggest the following procedure while dealing with time series data contaminated with gaps. The means of the gap distributions, $m_{s}$ and $m_{p}$ should initially be quantified. If they are not within tolerable limits, binning can be suggested as a method to change the values of $m_{s}$ and $m_{p}$ and bring them into the tolerable range identified[2]. Once $m_{s}$ and $m_{p}$ are within tolerable ranges, reliable conclusions can be drawn about $\alpha_{1}$ and $\alpha_{2}$. For any conclusions about underlying deterministic nonlinearity, one needs to conduct an analysis for comparison with phase randomized surrogates. If the surrogates give very different values for quantifiers, from the values for data, one can conclude the existence of an underlying nonlinearity. If surrogates and data merge within tolerable range, the underlying process can be assumed to be purely stochastic or deterministic with high noise contamination. In cases where, even with considerable binning the $m_{s}$ and $m_{p}$ values cannot be brought to tolerable ranges, no conclusions can be drawn based on multifractal analysis with confidence.

Acknowledgements

We acknowledge with thanks the data from the SMEAR database that are used in this research. SVG thanks Ron Sunny, IISER, Pune for help with the SMEAR datasets.