On the embedding statistic

Physics Letters A 184 (1993) 83-87 North-Holland

PHYSICS LETTERS A

On the embedding statistic C.P. P r i c e a n d D. P r i c h a r d 1 Physics Department, University of Alaska Fairbanks, Fairbanks, AK 99775-1320, USA Received 2 September 1993; revised manuscript received 26 October 1993; accepted for publication 29 October 1993 Communicated by A.R. Bishop

We describe a new statistic which appears to have some utility for indicating the presence of determinacy in a time series. We motivate this statistic, describe its algorithmic construction and give some examples of its use. The statistic is robust in that it is fairly insensitive to noise.

1. Introduction A n u m b e r o f statistics have been proposed to characterize signals arising from nonlinear deterministic processes (e.g. methods o f estimating dimensions, Lyapunov exponents, nonlinear prediction error, etc. [ 1-3 ] ). However, these statistics often cannot discriminate linear stochastic processes from nonlinear discriminate ones (see, e.g., refs. [ 4 - 6 ] ). Recently Kaplan and Glass [7] have proposed a statistic to detect determinism in a time series. Their statistic (and the one we propose) use a statistical approach for detecting determinism. That is one compares the results to those for "surrogate" data sets which have the same linear properties as the original data set but are known to be stochastic [8-10,5 ]. Surrogate data sets can easily be made by taking the Fourier transform o f the signal, randomizing the phase information, and inverting the transform. It can be shown that the amplitude distribution o f a surrogate data set made via scrambling the phases o f a Fourier transform is (usually) statistically Gaussian, in the sense that it goes to a Gaussian distribution in the limit that the size o f the data set becomes large. While investigating another related matter, it was observed that this property is preserved after embedding [ 1 l - 13 ]: if the embedded vectors are projected along the principal axes o f the i Current address: Theoretical Division (T-13), Los Alamos National Laboratory, Los Alamos, USA.

distribution in embedding space, the distribution o f the projections are also statistically Gaussian. A 8192 point data set was made by integrating the Lorenz equations with parameters a = 1 0 . 0 , b = 2 8 . 0 , and c - ~ , with a time step dt = 0.0 l; the first 1000 points were discarded to insure that the sequence had converged to the attractor. F r o m that data set, a phase scrambled surrogate was constructed. Figure 1 shows the amplitude distribution of that surrogate data set, and fig. 2 shows the distributions along the principal axes for a four-dimensional embedding o f that surrogate. One sees that all five o f the distributions are essentially Gaussian. One might surmise that the principal axis distri-

< Amplitude Fig. 1. The distribution of a 8192 point surrogate of the Lorenz attractor. The dashed line is a Gaussian distribution with the same mean and variance. The surrogate was made by the phase scrambring technique; the Lorenz attractor data set was made with parameters a= 10.0, b= 28.0, and c= }, with a time step dt=0.01.

0375-9601/93/$ 06.00 © 1993 Elsevier Science Publishers B.V. All fights reserved.

83

Volume 184, number 1

Axis 1 Amplitude

3 Amplitude

PHYSICS LETTERS A

]

IA" Axis 2 Amplitude

Axis 1 .~r~lit~le

Axis 2 Amplitude

Axis 4 An~plitude

Axis 3 Amplitude

Axis 4 Amplitude

Fig, 2. The distributions along the principal axes for a fourdimensional embedding of the 8192 point surrogate of the Lorenz attractor. (a) The distribution along the most significant principal axis. (b) The distribution along the second most significant principal axis. (c) The distribution along the third most significant principal axis. (d) The distribution along the least significant principal axis. The time delay is 12 (first zero crossing of the autocorrelation function). The dashed fines are Gaussian distributions with the same mean and variance; notice that the distributions are basically Gaussian.

butions in fact mimic the data set distribution; but in fact this is not the case. The same Lorenz attractor data set from which the surrogate was made was amplitude corrected to a strict Gaussian distribution (see the appendix). Figure 3 shows the distributions along the principal axes for a four-dimensional embedding of the amplitude corrected data set. We see that the first three projection distributions are highly non-Gaussian. Embeddings in higher dimensions show that most of the projection distributions are highly non-Gaussian, not just the first three as in this example.

2. Developing a statistic

As a result of these observations, a statistic was devised based on ,the distributions in embedding space. The embedding statistic, like the statistic of Kaplan and Glass [ 7 ], is envisaged as a discriminant of the determinacy of a time series. To be precise, the embedding statistic is based on comparison of a statistic with the same statistic for surrogate data sets. 84

27 December 1993

Fig. 3. The distributions along the principal axes for a fourdimensional embedding of the 8192 point set of the Lorenz attractor, after amplitude correction to a strict Gaussian distribution. (a) The distribution along the most significant principal axis. (b) The distribution along the second most significant principal axis. (c) The distribution along the third most significant principal axis. (d) The distribution along the least significant principal axis. The time delay is 12 (first zero crossing of the autocorrelation function), The dashed lines are Gaussian distributions with the same mean and variance; notice that the first three distributions are highly non-Gaussian.

One forms the f f statistic for the embedding space distribution relative to a Gaussian having the same mean and variance. Then, the X2 is compared to the average and r.m.s, deviation o f z 2 for a set often other surrogate data sets. Specifically, one forms an average (over the principal axes) Z 2 statistic (see below) for the time series after embedding in DE dimensions. Then, this (Z 2 ) is compared to the average and r.m.s, deviation of the (Z 2) for a set of ten similarly constructed surrogate data sets, (if) and o. The embedding statistic is defined to be (_)/~. Somewhat arbitrarily, we consider a set with an embedding statistic greater than 3 to be deterministic. After we describe the algorithm for calculating in detail, and then how that algorithm is used to make and o, we will calculate the embedding statistic for a number of test data sets. Algorithm for : (i) Given a series {xi}. Amplitude correct the series to have a strict Gaussian amplitude distribution of zero mean and unit variance (see the appendix);

Volume 184,number I

PHYSICS LETTERSA

call the amplitude corrected series {y~}.Note that no further reference to the original series {xt} will be made. (ii) Embed {Yi} in Ds dimensions by the method of time delays to make the vector set {y~}. (iii) For the set {y~} determine the principal axes ~t/) about the center of mass, and project out the series {z~) =y;.e u)}. Calculate the mean and the variance of {z~) }. Determine the distribution P ( z c/)) via some appropriate binning. (iv) Determine the distribution P' (again, in some appropriate binning) which is Gaussian and has the same mean and variance as the set {zy ) }. (v) Compute the chi-squared for the distributions P ( z U)) and P', Z2. --_ ~., (ek-P~,)2

J

"ff Pk +P'k

where k is the bin index. (vi) (Z 2) is the average over the DE %2's. Now, we are to compare (Z 2) to an average over surrogates. We make ten surrogates from the amplitude corrected series {yi}, call them {sy(m)i}. Amplitude correct each of those surrogates to a strict Gaussian distribution of zero mean and unit variance; call them {Sy(m);}. Use the method previously described to form (%(m) 2) from each series {Sy(rn)j}, and then find the mean and the r.m.s, deviation of the set { ( z ( m ) 2) }, (Z 2) and a.

3. Testing the statistic We have calculated the embedding statistic for a number of different data sets in order to test it. Each data set has 8192 elements. The data sets comprise three distinct classes. The first class is composed of deterministic data sets: a set made from the H~non attractor [ 14] (parameters a = 1.4 and b=0.3), a set made from the Lorenz attractor [ 15] (parameters a=10.0, b=28.0, and c=~, with a time step dt= 0.01 ), a set made from the RSssler attractor [ 16 ] (parameters a=0.15, b=0.2, and c=10.0, with a time step d t = ~ n ) , a set made from the hyperchaos attractor [17] (parameter a = 2 . 2 with a time step d t = 0 . 2 ) , and a set made from a three-torus (made of three equal amplitude sinusoids with cycle lengths 3.5/8192 times 41, 71 and 97). The second class is composed of stochastic data sets: a power law series,

27 December1993

three series made by autoregression of Gaussian random noise ("Gaussian autoregression", x~=axi_ l+Su where s~ is a Gaussian random variable and a=0.90, 0.95 and 0.99), and three series made by autoregression of a strongly non-Gaussian innovation, ("non-Gaussian autoregression", x~=ax~_ 1+ t;, where t~ randomly takes on the values 1.0 or - 1.0 and a=0.90, 0.95 and 0.99). The power law series has a power spectrum ocf -5/3 and was obtained by taking the inverse Fourier transform of a f - 5/3 power spectrum with phases uniformly distributed on the interval (0, 2n] [5]. The third class is used to test the ability of the embedding statistic to discriminate in the presence of added noise. This class is composed of four sets made by adding Gaussian random noise to the previously referenced set made from the Rfssler attractor: sets with the variance of the Gaussian random variable (noise amplitude) equal to 0.0, 1.0, 2.0, and 5.0 (roughly 0, 5, 10 and 25% noise). The embedding statistic was calculated for each of the data sets for embedding dimension DE = 4 through 10 inclusive. For the H~non set, the embedding statistic ranges from 293 to 455. The embedding statistics for the remaining members of the first class of data sets and for all the members of the second class of data sets are shown in fig. 4. The embedding statistic for the Lorenz data set ranges from 32 to 88. The embedding statistic for the R6ssler data set ranges from 8.7 to 29. The embedding statistic for the hyperchaos data set ranges from 16 to 27. The embedding statistic for the three-torus data set ranges from 15 to 35. The embedding statistic for the power law data set ranges from - 1 . 5 to 0.3. The embedding statistic ranges for the GAR sets are - 0 . 8 to 1.4 ( a = 0 . 9 0 ) , - 1.3 to 0.4 ( a = 0 . 9 5 ) and - 0 . 7 to 0.0 (a=0.95). The embedding statistic ranges for the NGAR sets are - 2 . 1 to - 0 . 4 ( a = 0 . 9 0 ) , - 2 . 4 to - 0 . 2 ( a = 0 . 9 5 ) and - 1 . 0 to 0.5 (a--0.95). The embedding statistic is clearly able to distinguish the deterministic data sets from the stochastic data sets. The embedding statistics for the members of the third class of data sets are shown in fig. 5. The embedding statistic for the Rfssler data set ranges from 8.7 to 29. The embedding statistic for the R6ssler data set with 5% noise ranges from 4.8 to 21. The embedding statistic for the R6ssler data set with 10% noise ranges from 2.3 to 16. The embedding statistic for the R6ssler data set with 25% noise ranges 85

Volume 184, number 1

PHYSICS LETTERS A

' 'T"

40

(Z

•

T

'

"T"

'

G

T

'

'T'

G

27 December 1993

30

'

0t

1 30

2O

Y

8 20

8

Y

y

2

Y

8

3

Y

3

I

1

4

P 2

,

4

,

l

,

I

,

,

.

6

,

I

8

,

I

,

t

,

10

,

4

D~ Fig. 4. The embedding statistic for DE=4-10 for most of the members o f the first class o f data sets (deterministic sets) and for the members of the second class o f data sets (stochastic sets). The points marked with the symbols tt, 13, T and 8 are for the Lorenz, Rtissler, hyperchaos and three-toms data sets, respectively. The points marked with the symbols a, b, e, d, e, f and g are for the power-law, GAR 0.90, GAR 0.95, GAR 0.99, NGAR 0.90, NGAR 0.95 and NGAR 0.99 data sets, respectively. The dashed line is at a value of 3.

from 1.4 to 9.1. The statistic is sufficiently robust to be able to identify the deterministic nature of the Rfssler data sets for noise up to about the 25% level. Note the alternation of values of the embedding statistic for the various Rtissler data sets; we do not have an explanation for this feature.

4. Conclusions

We are not claiming that the embedding statistic is the only one extant which is able to distinguish the underlying deterministic nature of a noisy signal. For example, using the null hypothesis that the process is distinct from a Gaussian linear stochastic with static monotonic transform, the Takens dimension estimator applied to the correlation integral provides nearly as good a test for the presence of a deterministic signal in the noisy Rfssler sets; however, even using Theiler's box assisted algorithm [ 18 ], the embedding statistic is about 50 times faster to corn86

2

2

!

b

o

1

4

,

,

I

,

,

,

.

I

6

t

,

,

8

,

I0

DB Fig. 5. The embedding statistic for DE=4-10 for the members of the third class of data sets (RSssler data sets plus added Gaussian noise). The points marked with the symbol 1 are for the set with no added noise; the points marked with the symbols 2, 3 and 4 are for the sets with 5, 10 and 25% added noise, respectively. The dashed line is at a value of 3.

pute than the correlation imegral for these data sets. We have described what we term the "embedding statistic" which appears to have some utility for indicating the degree of determinacy in a time series. It was observed that if the embedded vectors of a surrogate data set are projected along the principal axes of the distribution in embedding space, the distributions of the projections are also statistically Gaussian, while the similar projection distributions for deterministic data sets are non-Gaussian. The embedding statistic is a measure of statistical deviation of the Z 2 for a given distribution from the average X2 for a set of surrogate data sets. The embedding statistic is robust in that it is fairly insensitive to noise, as demonstrated by the addition of noise to a data set from the Rt~ssler attractor. Finally, the embedding statistic can be computed fairly rapidly.

Acknowledgement

We appreciate the comments of J. Theiler. D.P. acknowledges the support of a UAF Chancellors

Volume 184, number I

PHYSICS LETTERS A

27 December 1993

G r a d u a t e Fellowship a n d an A W U / D O E Thesis Parts Fellowship for study with the c o m p l e x systems group ( T - 1 3 ) at Los A l a m o s N a t i o n a l Laboratory. C o m p u t a t i o n a l resources were p r o v i d e d by the U A F Physics D e p a r t m e n t a n d the Arctic Regions Superc o m p u t e r Center. This w o r k was s u p p o r t e d b y N S F grant ATM-9213522.

which can be stably integrated with a simple center differenced predictor. The resulting series {Gxi} has a very s m o o t h ( " s t r i c t " ) G a u s s i a n distribution. Finally, rank the series (Gx;} to m a k e the series {y; = Gxrx,}. {y~}is an amplitude corrected series with strict G a u s s i a n d i s t r i b u t i o n o f zero m e a n a n d unit variance.

Appendix

References

The m e t h o d for a m p l i t u d e correcting a series to have a strict G a u s s i a n d i s t r i b u t i o n is similar to the first step o f the Theiler et al. A A F T algorithm o r surrogate m e t h o d I algorithm [ 9 ]. We begin with a ser e s {x~} with N elements. This series is then sorted from least to greatest to m a k e the series {sx~} where sXi+l ~ sx i. T h e ranking series {rx;} inverts the sorting operation: x = SXrx,. NOW we m a k e a sorted series with a strict G a u s sian d i s t r i b u t i o n o f zero m e a n a n d unit variance, {Gx~}. Such a series will obey the relation

[1 ] J. Theiler, J. Opt. Soc. Am. A 7 (1990) 1055. [2 ] J.-P. Eekmann and D. Ruelle, Rev. Mod. Phys, 57 (1985 ) 617. [ 3 ] J.D. Farmer and J.J. Sidorowieh, Phys. Rev. Lett. 59 (1987) 845. [4] J. Theiler, Phys. Rev. A 34 (1986) 2427. [ 5 ] A.R. Osborne and A. Provenzale, Physica D 35 (1989) 357. [6] A. Provenzale, A.R. Osborne and R. Soj, Physica D 47 (1991) 361. [7] D.T. Kaplan and L. Glass, Phys. Rev. Lett. 68 (1992) 427. [8]J. Thailer, B. Galdfikian, A. Longtin, S. Eubank and J.D. Farmer, in: SH studies in the sciences of complexity, Vol. XII. Nonlinear modeling and forecasting, eds. M. Casdagli and S. Enbank (Addison-Wesley, Reading, M_A,1992) pp. 163-188. [9]J. Theiler, B. Galdrikian, A. Longtin, S. Euhank and J.D. Farmer, in: Proe. 1st Experimental Chaos Conference, eds. S. Vohora, M. Spano, M. Sehlesinger,L. Pecora and W. Ditto (World Scientific, Singapore, 1992) pp. 47-53. [ 10] M.B. Kennel and S. Isabelle, Phys. Rev. A 46 (1992) 3111. [ 11 ] N.H. Packard, J.P. Crutchfield, J.D. Farmer and R.S. Shaw, Phys. Rev. Lett. 45 (1980) 712. [ 12 ] F. Takens, Detecting strange attractors in fluid turbulence, in: Dynamical systems and turbulence, eds. D. Rand and L.-S. Young (Springer, Berlin, 1981 ). [13]T. Sauer, J.A. Yorke and M. Casdagli, J. Slat. Phys. 65 (1991) 579. [ 14] M. H~non, Commun. Math. Phys. 50 (1976) 69. [ 15 ] E.N. Lorenz, J. Atmos. Sei. 20 (1963) 130. [ 16] O.E. RSssler, Phys. Lett. A 57 (1976) 397. [ 17 ] O.E. RSssler, Phys. Lett. A 71 (1979) 155. [ 18] J. Theiler, Phys. Rev. A 36 (1987) 4456.

Gx~

i r a J e x p ( - ½ u 2) d u , --O/

where we choose ot such that ½= A exp( - ½0t2). T h e n self-consistently we d e m a n d

A=N/f?~, e x p ( - ½u 2) du ~ N / f~_oo e x p ( - ½u2) d u = N / x / ~ for sufficiently large values o f N. Differentiate the relation with respect to the index i, replace the d e r i v a t i v e dGx/di with the ratio o f differences AGx/Ai=AGx, a n d find the spacing between elements A G x - - A - l e x P ( ½ G x 2) ,

87