{'fraction_non_alphanumeric': 0.02659113243540089, 'fraction_numerical': 0.005260419677438002, 'mean_word_length': 4.595084087968952, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On Surrogate Data Testing for Linearity based on the Periodogram

The method of surrogate data is a tool to test whether data were generated by some class of model. Tests based on the periodogram have been proposed to decide if linear systems driven by Gaussian noise could have generated a sample time series. We show that this procedure based on the periodogram, in general, misspecifies the test statistic. Based on the theory of linear systems we suggest an alternative procedure to obtain the correct distribution of the test statistic and discuss problems of this approach.


Introduction
The method of surrogate data was proposed by Theiler et al. [1] as a general procedure to test whether data are consistent with some class of models. It has been used to test a proposed nonlinear dynamics in geophysics [2], certain properties of neural spike dynamics [3], but the most common application is testing for linear models [4,5,6,7,8,9,10]. In order to test the hypothesis that the data are consistent with being generated by a linear system, the FT algorithm is applied. Based on an example and the theory of linear stochastic systems we will show that this algorithm does not produce the correct distribution of time series and therefore might not generally yield the correct distribution of the test statistic.
The surrogate data method is similar to a Monte Carlo implementation of a standard hypothesis test in which one specifies a procedure for generating an ensemble of new time series with a distribution that matches a model to be tested. One generates several sample time series, the surrogate data. Next, a selected feature function f is evaluated for the original time series x and for each of the surrogates. If the distribution of the feature values is well approximated by a Gaussian it can be characterized by its mean µ surr and standard deviation σ surr . If f (x), the feature value calculated for the measured time series, is more than a few standard deviations away from µ surr , the null hypothesis that the original data were generated by the tested model is rejected, and the significance is reported as the distance z = |µ surr −f (x)| σ surr . If the distribution of the feature is not well approximated by a Gaussian distribution, nonparametric rank statistics should be used.
The surrogate data method differs from a simple Monte Carlo implementation of a hypothesis test in that it tests not against a single model, but a class of models, i.e. linear systems driven by Gaussian noise. The idea is to select a single model from the class on the basis of the measured data x, and then do a 2 Monte Carlo hypothesis test for the selected model and the original data.
The statistical properties of a time series generated by a linear process are specified by the autocovariance function (ACF) or equivalently by its Fourier transform, the power spectrum S(ω). Theiler et al. use the periodogram Per(ω), i.e., the squared modulus of the Fourier transform of the data, to specify the linear model to be tested 1 . They generate surrogates (FT algorithm) by drawing complex numbers with amplitude Per(ω) and random phases and Fourier transforming back into the time domain. There are two problems in this procedure: First, the periodogram Per(ω) is not a consistent estimator of the power spectrum S(ω), and second, the distribution of surrogates generated by randomizing the phases but not the amplitudes does not correspond to any linear stochastic system at all. The consequences of these two points on the estimated distribution of the feature, in case of Gaussianity µ surr and σ surr , are the primary focus of this paper.
In the next section we briefly summarize the main ideas of surrogate data testing for linearity based on the periodogram and demonstrate by a simple example that this procedure might misspecificy the distribution of the test statistic.
To explain this result we refer to some results of the theory of linear systems in section 3 and discuss the consequences for surrogate data testing for linearity in section 4.

Periodogram Based Surrogates
The key to a hypothesis test is the distribution of a test statistic f when the null hypothesis H 0 is true, i.e. p(f |H 0 ). Given a measured sample x, one calculates f (x) and decides if it is reasonable to suppose that f (x) was drawn from p(f |H 0 ), in general by an a priori choosen significance level. 1 An invertible nonlinearity in the observation function should be rectified by a change of coordinates that forces the distribution of the data to be Gaussian [11].

3
The surrogate data method was proposed by Theiler et al. primarily to test the null hypothesis that the data come from a linear Gaussian process. Once again, the key to designing a hypothesis test is the identification of the distribution of a test statistic p(f |H 0 ). For a class as large and complex as all linear processes, it seems impossible. The idea of Theiler et al. is to use the observation itself to limit the size of the whole class to a single member.
Surrogate data testing needs two ingredients: A procedure to generate data from the class of models to be tested and a feature or statistic to perform the testing. If a single model is given, the first step is simple. One generates a large number of time series and approximates the distribution p(f |H 0 ). For instance, see the paper by Witt et al. [2] in which a nonlinear difference equation proposed to describe the dynamics of the reversals of the earth's magnetic field was rejected using a feature that describes the transition probabilities.
Theiler et al. proposed their FT algorithm to determine if a given time series could have been produced by a linear system driven by Gaussian noise. The basic idea is that the statistical properties of a linear process are captured completely by its spectral properties. Therefore, they first calculate the periodogram Per(ω), i.e.the squared modulus of the Fourier transform, for the positive frequencies. In order to generate the surrogate data they draw complex numbers with amplitude Per(ω) and random phases, symmetrize the negative frequencies with x(−ω) = x(ω) * to obtain a real time series by Fourier transforming back into the time domain.
We now give a simple example that demonstrates that the FT algorithm, might not reproduce the parent distribution of a feature for a linear process. We used a simple Gaussian white noise process with mean zero and unit variance and evaluated the simple feature: The Kolmogorov-Smirnov-test for Gaussianity does not reject the consistency with a normal distribution at the 99% level of confidence for the shown distributions. Thus, they may be described by mean and standard deviation.
The standard deviation of the true distribution is underestimated by a factor of two. The Kolmogorov-Smirnov-test for consistency of two distributions applied to each surrogate distributions and the distribution of new realizations of the process rejected the null hypothesis in every case at the 99% level of confidence.
The choosen test statistic is indeed able to detect nonlinearities since its expectation value for a Gaussian white noise process with unit variance is 2/ √ π = 1.12... and for "gaussified" data of the logistic map, r = 4, which shows also a flat spectrum, it is 1.09... .

Linear Processes
The misallocation of µ surr and the underestimation of the true standard deviation by σ surr that appear in fig. 1 can be explained in terms of the theory of linear processes. 5

Time Domain
In order to model the fluctuations in the periodicity of sunspot data Yule [12] introduced linear stochastic systems: These are now called autoregressive (AR) processes can be generalized to autoregressive moving average (ARMA) processes by including not only past values of the observations but also past values of the noise. The general ARMA[p,q] process reads: State space models are a further generalization. They consist of a vector valued AR [1] process and a linear observation function with observation noise: where Q is the covariance matrix of the dynamical noise, and R the variance of the observation noise. State space models model observation noise explicitly.
For AR and ARMA models, observation noise causes the estimated parameters to be biased towards zero.
Once the type of model has been chosen, one must select the model order.
Generally, one fits a sequence of higher and higher order models and tests the prediction errors for consistency with white noise, e.g. by a Kolmogorov -Smirnov -test on the periodogram [13].
If a particular class and order of model is selected, when the parameters of the model are estimated using data from a random process, the estimated parameters are random variables. By restricting the number of parameters to be fit, one insures that as the number of observed data increases, the variance of the estimates goes to zero. In the next subsection we emphasize that periodograms do not have this property.

Frequency Domain
We now examine the spectrum of linear processes and the relationship between the spectrum and the periodogram. The autocovariance function R(t) is given by and the spectrum is defined by: The spectrum of a state space model is given by : Thus the spectrum S(ω) of this process is a smooth function of ω. Spectra of AR or ARMA processes are special cases of eq. (5).
We now turn to the distribution of the periodogram. For details see [13,15], we only summarize the results. Asymptotically, the Fourier transform of a linear stochastic process is a complex Gaussian random variable where S(ω) is independent of N. Furthermore, Fourier components are asymptotically uncorrelated for different frequencies. In terms of modulus and phase instead of real and imaginary part, this means that not only is the phase a random variable but so is the modulus.
The distribution of the periodogram follows from eq. (7) to be a χ 2 distritubion with two degrees of freedom χ 2 2 : again independent of N.
Since the mean and variance of χ 2 2 are 2 and 4 respectively, the standard deviation of the periodogram is equal to its mean, i.e.

Per(ω) = S(ω) ± S(ω).
Thus, the periodogram fluctuates wildly regardless of the number of data points used to calculate it. It is an unbiased but not a consistent estimator of the spectrum since its variance does not decrease with N. This is true not only for linear stochastic processes but also for nonlinear stochastic and even for chaotic processes if only the autocovariance function decays fast enough.

Consequences for the Surrogate Data Method
The characteristics of the spectrum and periodogram described in the previous section lead to the following two points concerning the FT algorithm: 8 1. The FT algorithm does not use its estimate of the spectrum correctly. Both phases and amplitudes should be randomized.

2.
A better estimator of the spectrum than the periodogram is required.
Below we elaborate on these points.

Randomizing Phases and Amplitudes
Eq. (7) means that at each frequency the Fourier transform of a time series from a linear stochastic system is a two dimensional Gaussian random variable. The Note, that this process depends on the length N of the time series.
If the true spectrum of the process were known the following algorithm would yield the correct distribution of surrogates. The algorithm is based on eq. (7) that connects the spectrum S(ω) with the variance of the complex random variable x(ω).
• For each frequency ω i = i N ∆t , i = 1, ..., N 2 , draw two Gaussian distributed random numbers, multiply them by 1 2 S(ω i ) and take the result as the real and imaginary parts of the Fourier transform of the desired data.
• For reasons of symmetry the last frequency is always real for even number of data points. Thus, in this case, only one Gaussian distributed random number has to be drawn.
• To obtain a real valued time series, chose the Fourier components for the negative frequencies as the complex conjugates of the components at positive frequencies, i.e., x(−ω) = x * (ω).

Estimating the Spectrum
The critical issue for the procedure described in the previous section is spectral estimation. There are two general approaches commonly used for this purpose: Fit an AR or ARMA model to the data and calculate the spectrum from the resulting parameters or smooth the periodogram. For more than the brief discussion given here, see [13,15].

Spectral estimation by fitting AR or ARMA models
This method cannot be recommended for two reasons. On the one hand, as discussed in section 3, in presence of observation noise these models do not give correct results even if the underlying dynamic is linear. On the other hand, if the process is nonlinear they are, in general, not consistent estimators of the spectrum. Of course, the parameters of linear models determine the spectrum. But not every spectrum, e.g. of a chaotic process, is well described by a linear model, since the prediction errors will not be consistent with white noise. Increasing the order of the model does not help, because the "spectrum" calculated from the parameters will begin to interpolate the periodogram. For the same reason, surrogate data generated by realizations from a fitted ARMA model test, in general, for some linear model, but not for that specified by the spectrum.

Estimation by smoothing the periodogram
Estimating the spectrum by smoothing the periodogram rests on three results of the theory of linear processes: the expectation value of the periodogram is the spectrum, the spectrum is a continuous function, and the random variables Per(ω) and Per(ω ′ ) are asymptotically uncorrelated. The idea is to smooth the periodogram to reduce the variance (or equivalently split the time series into segments and average the periodograms). Denoting the number of data points by N, the number of frequency bins included in the smoothing by n s , one can show that in the limit of N → ∞ and n s → ∞, with n s constrained to increase slowly enough to insure ns N → 0, this procedure leads to a consistent estimator for the spectrum.
For finite data sequences, the smoothing will introduce a bias and there is a trade-off between bias and variance of the estimator which is discussed at length in [13,15].

Surrogates based on an estimated spectrum
The bias and the variance of the spectral estimate lead to a misspecification of the surrogate data distribution that vanishes asymptotically. Please locate fig.3 near here

Conclusion
The appeal of the general idea of surrogate data testing is clarified by a comparison with the tests for linearity proposed in the stochastic systems literature.
These tests rely on asymptotic results for linear processes, mainly concerning third order properties. There are two classes of these tests. The first class is based on fitting an AR model to the data and investigating third order properties of the prediction errors. Examples of these parametric tests are given by [16,17]. The other class of tests is based on the asympotitic property that the three-point-correlations: are zero for any linear process. For examples of these nonparametric tests see [18,19]. In order to perform the test, c(τ 1 , τ 2 ) is Fourier transformed and the normalized bispectrum is estimated. For linear Gaussian processes it should be zero and the consistency with zero can be tested. Both classes of tests test data of a given process by a property that is fulfilled by all possible linear processes.
The intuition behind the surrogate data method is that only that linear process that, jugded by its spectral properties, could have generated the data needs to be tested.
The sampling properties of a test statistic of a linear process depend on the the details of the process and the statistic. In particular, the variance will often be larger for processes with sharper peaks in the spectrum or equivalently longer coherence times. For example, the variance of the sample meanx for an iid Gaussian process N (0, σ 2 ) is σ 2 N , but if the process has correlations, the variance of the sample mean will decrease more slowly as a function of N. By construction, the surrogate data method produces the process specific sampling properties.
The power of any test, i.e. the ability to detect a violation of the null hypothesis, depends on the statistic. The bispectral tests considering the three-pointcorrelations are sensitive to quadratic nonlinearities but can be fooled completely by processes showing only cubic nonlinearities like the Duffing -oscillator. The general scheme of surrogate data testing can use any statistic, and the power of the test will depend on the statistic chosen. Thus one can detect a larger class of violations of the null hypothesis.
We showed by an example that the FT algorithm proposed by Theiler et al.
[1] might specify the distribution of the test statistic incorrectly, independent of the number of data points used. Based on the theory of linear stochastic systems, we proposed an alternative that involves spectral estimation. Therefore, testing for linearity by the idea of surrogate data is, in general, also a procedure that is only valid asymptotically. For a finite number of data, simulation studies should be performed to evaluate the estimation error.