Nonlinear dynamics, delay times, and embedding windows

Physica D 127 (1999) 48–60 Nonlinear dynamics, delay times, and embedding windows H.S. Kim1;a, R. Eykholt b;�, J.D. Salas c a Department of Civil Engineering, Colorado State University, Fort Collins, CO 80523, USA b Department of Physics, Colorado State University, Fort Collins, CO 80523, USA c Hydrologic Science and Engineering Program, Department of Civil Engineering, Colorado State University, Fort Collins, CO 80523, USA Received 11 September 1996; received in revised form 10 January 1998; accepted 12 August 1998 Communicated by A.M. Albano Abstract In order to construct an embedding of a nonlinear time series, one must choose an appropriate delay time �d. Often, �d is estimated using the autocorrelation function; however, this does not treat the nonlinearity appropriately, and it may yield an incorrect value for �d. On the other hand, the correct value of �d can be found from the mutual information, but this process is rather cumbersome computationally. Here, we suggest a simpler method for estimating �d using the correlation integral. We call this the C–C method, and we test it on several nonlinear time series, obtaining estimates of �d in agreement with those obtained using the mutual information. Furthermore, some researchers have suggested that one should not choose a fixed delay time �d, independent of the embedding dimension m, but, rather, one should choose an appropriate value for the delay time window �w D .m − 1/� , which is the total time spanned by the components of each embedded point. Unfortunately, �w cannot be estimated using the autocorrelation function or the mutual information, and no standard procedure for estimating �w has emerged. However, we show that the C–C method can also be used to estimate �w. Basically �w is the optimal time for independence of the data, while �d is the first locally optimal time. As tests, we apply the C–C method to the Lorenz system, a three-dimensional irrational torus, the Rossler system, and the Rabinovich–Fabrikant system. We also demonstrate the robustness of this method to the presence of noise. c©1999 Elsevier Science B.V. All rights reserved. Keywords: Delay time; Correlation integral; Embedding; Time series PACS: 05:45: C b; 47:52: C j 1. Introduction Analysis of chaotic time series is common in many fields of science and engineering, and the method of delays has become popular for attractor reconstruction from scalar time series. From the attractor dynamics, one can estimate the correlation dimension and other quantities to see whether the scalar time series is chaotic or stochastic. Therefore, attractor reconstruction is the first stage in chaotic time series analyses. Since the choice of the delay � Corresponding author Tel.: +1-970-491-7366; fax: +1-970-491-7947; e-mail: eykholt@lamar.colostate.edu 1 Present Address: Department of Construction Engineering, Sun Moon University, Korea 0167-2789/99/$ – see front matter c©1999 Elsevier Science B.V. All rights reserved. PII: S 0167 - 278 9 (98 )00240 -1 More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn H.S. Kim et al. / Physica D 127 (1999) 48–60 49 time �d for attractor reconstruction using the method of delays has not been fully developed, many researchers use the autocorrelation function, which is computationally convenient and does not require large data sets. However, it has been pointed out [1] that the autocorrelation function is not appropriate for nonlinear systems, and, instead, �d should be chosen as the first local minimum of the mutual information. Unfortunately, this approach is cumbersome computationally and requires large data sets [2]. According to Packard et al. [3] and Takens [4], the method of delays can be used to embed a scalar time series fxig; i D 1; 2; : : : ; into an m-dimensional space as follows: xxxi D .xi; xiC1; : : : ; xiC.m−1/t /; xxxi 2 Rm; (1) where t is the index lag. If the sampling time is �s, the delay time is �d D t�s. Takens’ theorem assumes that we have an infinite noise-free data set, in which case, we can choose the delay time almost arbitrarily. However, since real data sets are finite and noisy, the choice of the delay time plays an important role in the reconstruction of the attractor from the scalar time series. If �d is too small, the reconstructed attractor is compressed along the identity line, and this is called redundance. If �d is too large, the attractor dynamics may become causally disconnected, and this is called irrelevance [5]. In common practice, the delay time �d is chosen so as to ensure that the components of xxxi are independent, and the same delay time is used for all embedding dimensions m. However, in recent years, it has been suggested that the delay time window �w D .m − 1/� , which is the entire time spanned by the components of xxxi , should be independent of m instead [6–11]. In this case, the delay time � varies with the embedding dimension m. Rosenstein et al. [6] compared several approaches for estimating �d and �w and indicated the critical disadvantages of each, such as inconsistencies and long computation times. The geometrical concepts of redundance and irrelevance were used to evaluate the quality of an attractor’s reconstruction. It has been suggested [11] that the delay time window may be set to �w � �p, where �p is the mean orbital period, which can be approximated by examining the oscillations of the time series. Martinerie et al. [7] examined the delay time window and compared it with the delay times estimated using the autocorrelation function and the mutual information, but they concluded that �w cannot be estimaed using either of these methods. Basically, �w is the optimal time for independence of the data, but these methods estimate only the first locally optimal time, which is �d. In this work, we develop a technique for choosing either the delay time �d or the delay time window �w using the correlation integral. Although �w may be a more useful quantity for the estimation of the dimension than �d, many researchers continue to use �d for attractor reconstruction using the method of delays. Also, the estimation of both �d and �w using this technique helps to illuminate the difference between these two quantities. Brock et al. [12,13], in their development of a test for nonlinearity in a time series, used the statistic S.m; N; r/ D C.m; N; r/ − Cm.1; N; r/, where C.m; N; r/ is the correlation integral (see Eq. (2)). Here, we use the similar statistic S.m; N; r; t/ D C.m; N; r; t/ − Cm.1; N; r; t/, and we examine its dependence on the index lag t . Hence, we call our technique the C–C method. Compared with the use of the mutual information, the C–C method is easier to implement, useful for smaller data sets, and less demanding computationally. Also, its estimate of �d agrees with the first local minimum of the mutual information. 2. Correlation integral and BDS statistic The correlation dimension introduced by Grassberger and Procaccia [14,15] is widely used in many fields for the characterization of strange attractors. The correlation integral for the embedded time series is the following function: C.m; N; r; t/ D 2 M.M − 1/ X 1�i 0; (2) More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn 50 H.S. Kim et al. / Physica D 127 (1999) 48–60 where 2.a/ D 0; if a < 0; 2.a/ D 1; if a � 0; N is the size of the data set, t is the index lag, M D N − .m − 1/t is the number of embedded points in m- dimensional space, and k � � � k denotes the sup-norm. C.m; N; r; t/ measures the fraction of the pairs of points xxxi; i D 1; 2; : : : ; M; whose sup-norm separation is no greater than r . If the limit of C.m; N; r; t/ as N ! 1 exists for each r , then the fraction of all state vector points that are within r of each other is denoted by C.m; r; t/ D limN!1C.m; N; r; t/; and the correlation dimension is defined as D2.m; t/ D limr!0[logC.m; r; t/=logr]. Since N remains finite for real data sets, then we cannot let r go to 0; instead, we look for a linear region of slope D2.m; t/ in the plot of logC.m; N; r; t/ vs: log r . Brock et al. [12,13] studied the BDS statistic, which is based on the correlation integral, to test the null hypothesis that a given data set is independently and identically distributed (iid). The test has been particularly useful for chaotic systems and nonlinear stochastic systems. For completeness, we briefly review the BDS statistic. Let F be the invariant distribution of a variable X in the state space, and define the spatial correlation integral as follows: C.m; r/ D Z Z 2.r − kxxx − yyyk/dF.xxx/dF.yyy/; r > 0: (3) If X is iid, then using 2.r − kxxx − yyyk/ D QmkD12.r − jxk − ykj/ yields C.m; r/ D Cm.1; r/, where C.1; r/ D Z [F.x C r/ − F.x − r/]dF.x/ � C: (4) Denker and Keller [16] showed that C.m; N; r/ is a U-statistic estimator. Using the U-statistics theory for an absolutely regular process, Brock et al. [12,13] proved that, as N ! 1; pN [C.m; N; r/ − Cm.1; r/] approaches a normal distribution with mean zero and variance � 2.m; r/ D 4 " Km − C2m C 2 m−1X iD1 .Km−iC2i − C2m/ # ; (5) where K � Z [.F .x C r/ − F.x − r/]2dF.x/ (6) (assuming that K > C2). Thus, the BDS statistic defined by BDS.m; N; r/ D p N �.m; r/ [C.m; N; r/ − Cm.1; r/] (7) approaches a standard normal distribution. However, since the distribution F is generally unknown, then we can- not obtain the values of C and K and the variance � 2.m; r/ from the above definitions. Instead, the correlation integral C.1; r/ and the variance � 2.m; r/ must be estimated from the sample data. Thus, C.1; r/ is estimated by C.1; N; r; t/, and � 2.m; r/ is estimated by O� 2 D 4 ( m.m−1/ OC2.m−1/. OK− OC2/C OKm− OC2m C 2 m−1X iD1 [ OC2i . OKm−i− OC2.m−i//−m OC2.m−i/. OK− OC2/] ) ; (8) where OC D C.m; N; r; t/ and OK D 6 M.M − 1/.M − 2/ X 1�i OC2 and m > 1, the BDS statistic becomes BDS.m; N; r/ D p N O� [C.m; N; r; t/ − C m.1; N; r; t/]; (10) and this converges to a standard normal distribution as N ! 1. The BDS statistic originates from the statistical properties of the correlation integral, and it measures the statistical significance of calculations of the correlation dimension. Even though the BDS statistic cannot be used to distinguish between a nonlinear deterministic system and a nonlinear stochastic system, it is a powerful tool for distinguishing random time series from chaotic or nonlinear stochastic time series. Further statistical properties and proofs can be found in [12,13]. 3. Measure of nonlinear dependence 3.1. C–C method This study is concerned with the properties of S.m; N; r; t/ D C.m; N; r; t/ − Cm.1; N; r; t/. We refer to a comment by Brock et al. [12]: “If the stochastic process fxig is iid, it will be shown that C.m; r/ D Cm.l; r/ for all m and r . That is to say, the correlation integral behaves much like the characteristic function of a serial string in that the correlation integral of a serial string of independent random variables is the product of the correlation integrals of component substrings.” This leads us to interpret the statistic S.m; N; r; t/ as the serial correlation of a nonlinear time series. Therefore, it can be regarded as a dimensionless measure of nonlinear dependence. For fixed m; N , and r the plot of S.m; N; r; t/ vs. t is a nonlinear analog of the plot of the autocorrelation function vs. t . In order to study the nonlinear dependence and eliminate spurious temporal correlations, we must subdivide the time series fxig; i D 1; 2; : : : ; N , into t disjoint time series, where t is the index lag. S.m; N; r; t/ is then computed from these disjoint time series as follows: For t D 1, we have the single time series fx1; x2; : : : ; xN g, and S.m; N; r; 1/ D C.m; N; r; 1/ − Cm.1; N; r; 1/: (11) For t D 2, we have the two disjoint time series fx1; x3; : : : ; xN−1g and fx2; x4; : : : ; xN g, each of length N=2, and we average the values of S.m; N=2; r; 2/ for these two series: S.m; N; r; 2/ D 12 f[C1.m; N=2; r; 2/ − Cm1 .1; N=2; r; 2/] C [C2.m; N=2; r; 2/ − Cm2 .1; N=2; r; 2/]g: (12) For general t , this becomes S.m; N; r; t/ D 1 t tX sD1 [Cs.m; N=t; r; t/ − Cms .1; N=t; r; t/]: (13) Finally, as N ! 1, we can write S.m; r; t/ D 1 t tX sD1 [Cs.m; r; t/ − Cms .1; r; t/]; m D 2; 3; : : : (14) For fixed m and t , S.m; r; t/ will be identically equal to 0 for all r if the data are iid and N ! 1. However, real data sets are finite, and the data may be serially correlated; so, in general, we will have S.m; r; t/ 6D 0. Thus, the locally optimal times may be either the zero crossings of S.m; r; t/ or the times at which S.m; r; t/ shows the least More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn 52 H.S. Kim et al. / Physica D 127 (1999) 48–60 variation with r , since this indicates a nearly uniform distribution of points. Hence, we select several representative values rj , and we define the quantity 1S.m; t/ D maxfS.m; rj ; t/g − minfS.m; rj ; t/g; (15) which is a measure of the variation of S.m; r; t/ with r . The locally optimal times t are then the zero crossings of S.m; r; t/ and the minima of 1S.m; t/. The zero crossings of S.m; r; t/ should be nearly the same for all m and r , and the minima of 1S.m; t/ should be nearly the same for all m (otherwise, the time is not locally optimal). The delay time �d will correspond to the first of these locally optimal times. 3.2. Choosing m and r for finite sample sizes In determining the nonlinear dependence of a finite time series by using the statistic S.m; N; r; t/, one must have criteria for selecting the values of m and r . In addition, one must know the role of the sample size N . For a fixed value of N , as m becomes large, the data become very sparse, so that C.m; N; r; t/ becomes vanishingly small. Also, if r exceeds the size of the attractor, then C.m; N; r; t/ saturates, since most pairs of points are within the distance r . Thus, neither m nor r should be too large. Appropriate choices for m; N , and r may be found by examining the BDS statistic. Brock et al. [12] carried out such an investigation using time series generated from a variety of asymptotic distributions. Time series with three sample sizes, N D 100; 500 and 1000, were generated by Monte Carlo simulations from six asymptotic distributions: a standard normal distribution, a student-t distribution with three degrees of freedom, a double exponential distribution, a chi-square distribution with four degrees of freedom, a uniform distribution, and a bimodal mixture of normal distributions. These studies led to the conclusion that m should be between 2 and 5 and r should be between �=2 and 2� . In addition, the asymptotic distributions were well approximated by finite time series when N � 500. As an example, we illustrate the dependence of the statistic S.m; N; r; 1/ on m; N , and r for the Lorenz system of three coupled differential equations [17]: Px D −a.x − y/; Py D −xz C cx − y; Pz D xy − bz; (16) where a; b, and c are constants. We solve this system of equations for a D 16:0; b D 4:0, and c D 45:92 to generate a time series of the variable x with �s D 0:01. The values of S.m; N; r; 1/ for various sample sizes N and various values of m and r are shown in Fig. 1. For N > 100, S.m; N; r; 1/ is positive for all r between 0 and 4� , where � is the standard deviation of the data set, and it is nearly 0 for r > 4� . For N > 1000, the behavior of S.m; N; r; 1/ as a function of r within the range 0 < r < 4� is nearly the same as for N D 1000, and this represents the true correlation of the time series. However, as N decreases below 100, S.m; N; r; 1/ fails to represent the true correlation. For example, when N D 15 and m D 10, S.m; N; r; 1/ vanishes at r D .0:5/� , and it departs substantially from the curves for N D 1000. On the other hand, for m D 2, S.m; N; r; 1/ has roughly the correct behavior even for N as small as 15. These results are in agreement with the conclusions of Brock et al. Within the ranges 2 � m � 5 and �=2 � r � 2� , the quantity S.m; N; r; 1/ represents the serial correlation for the wide range of values of N considered in Fig. 1. Finally, note that, while the asymptotic distributions are well approximated by finite time series of size N � 500, it is not necessary for N to be this large in order to use S.m; N; r; 1/ to study the nonlinear dependence of the time series. For the range �=2 � r � 2� , we choose four representative values r1 D .0:5/�; r2 D .1:0/�; r3 D .1:5/� , and r4 D .2:0/� . We then define the following averages of the quantities given by Eqs. (14) and (15): NS.t/ D 1 16 5X mD2 4X jD1 S.m; rj ; t/; (17) More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn H.S. Kim et al. / Physica D 127 (1999) 48–60 53 Fig. 1. S.m; N; r; 1/ vs. R D r=� for various sample sizes N and for m D 2; 5, and 10 using time series of the variable x from the Lorenz system of Eq. (16) with a D 16:0; b D 4:0; c D 45:92, and �s D 0:01. 1 NS.t/ D 1 4 5X mD2 1S.m; t/; (18) and we look for the first zero crossing of NS.t/ or the first local minimum of 1 NS.t/ to find the first locally optimal time for independence of the data, which gives the delay time �d D t�s. The optimal time is the index lag t for whichNS.t/ and 1 NS.t/ are both closest to 0. If we assign equal importance to these two quantities, then we may simply More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn 54 H.S. Kim et al. / Physica D 127 (1999) 48–60 Fig. 2. S.m; r; t/ for the variable x from the Lorenz system of Eq. (16) with a D 16:0; b D 4:0; c D 45:92, and �s D 0:01 using 3000 data points. The circles indicate the vicinity of �d, where the first local minimum occurs in the variation of S.m; r; t/ with r . look for the minimum of the quantity Scor.t/ D 1 NS.t/ C j NS.t/j; (19) and this optimal time gives the delay time window �w D t�s. 4. Numerical examples The following examples show that the C–C method is relatively easy to implement, that the delay times �d agree with those obtained using mutual information, that the delay time windows �w agree with those obtained by Martinerie et al. [7], and that the method is robust to the presence of noise. 4.1. Applications to dynamical systems We solve the Lorenz system of Eq. (16) for a D 16:0; b D 4:0; and c D 45:92 to generate a time series of 3000 data points with �s D 0:01. We then compute S.m; r; t/ from Eq. (14), and the results are shown in Fig. 2. The circles in Fig. 2 indicate the index lag t where the variation of S.m; r; t/ with r is at its first local minimum, and Fig. 3(a) shows this first local minimum of 1S.m; t/ more clearly. We choose the delay time at this point, which gives More documents and datum download website Lu Zhenbo's Blog: blog.sina.com.cn/luzhenbo2 Communication & Cooperation: luzhenbo@yahoo.com.cn H.S. Kim et al. / Physica D 127 (1999) 48–60 55 Fig. 3. 1S.m; t/; 1 NS.t/; NS.t/, and Scor.t/ for the variable x from the Lorenz system of Fig. 2. The solid line locates �d D 10�s, and the minimum of Scor.t/ yields �w D 100�s. �d D 10�s D 0:10. This agrees with the delay time found by Abarbanel et al. [18] using the first local minimum of the mutual information. Als