Full text
Discrete Dynamics in Nature and Society, 2002 VoL. 7 (1), pp. 41-52 Taylor & Francis Taylor & Francis Group Effect of Parameter Calculation in Direct Estimation of the Lyapunov Exponent in Short Time Series A.M. L(3PEZ JIMINEZ a’*, C. CAMACHO MARTI’NEZ VARA DE REY" and A.R. GARC[A TORRES b aDepartment of Experimental Psychology, University of Seville, Avda. Camilo Jos Cela s/n. 41005 Seville, Spain; bphysics Seminar, LE.S. Los Viveros, Avda. Bias Infante, s/n. Seville, Spain (Received 21 April 2001) The literature about non-linear dynamics offers a few recommendations, which sometimes are divergent, about the criteria to be used in order to select the optimal calculus parameters in the estimation of Lyapunov exponents by direct methods. These few recommendations are circumscribed to the analysis of chaotic systems. We have found no recommendation for the estimation of A starting from the time series of classic systems. The reason for this is the interest in distinguishing variability due to a chaotic behavior of determinist dynamic systems of variability caused by white noise or linear stochastic processes, and less in the identification of non-linear terms from the analysis of time series. In this study we have centered in the dependence of the Lyapunov exponent, obtained by means of direct estimation, of the initial distance and the time evolution. We have used generated series of chaotic systems and generated series of classic systems with varying complexity. To generate the series we have used the logistic map. Keywords: Non-linear dynamic; Attractor; Chaos; Lyapunov exponent INTRODUCTION The discovery of chaotic behavior in deterministic dynamical systems has changed some philosophical aspects in the prevailing scientific paradigm and has opened new perspectives for the design and analysis of time series (Barnett and Choi, 1989; Casdagli, 1991; Casdagli et al., 1991; Sayers, 1991; Berliner, 1992; McCaffrey et al., 1992; Nychka et al., 1992; Gerr and Allen, 1993; Takens, 1993). In the 1980s, the breakthroughs in the analysis of time series based on the Qualitative Theory of Dynamical Systems have yielded a set of indexes. These, in theory, should allow us to determine if the apparently random time sequence observations of a system state, can or cannot be due to chaotic behavior generated by a system of nonlinear deterministic equations (Ashley et al., 1986; Broomhead and King, 1986; Ashley and Patterson, 1989; Brown et al., 1991; Grassberger et al., 1991; McCaffrey et al., 1992; Abarbanel et al., 1993; Palu et al., 1993; Takens, 1993). As a sub-product, it is possible to determine the number of variables which this set of unknown equations would bring into play, as well as to classify systems into universal classes (linear-non-linear, stochastic-deterministic) and relate the changes in the behavior quantifiers with changes occurred in the dynamical behavior of the system (bifurcations) (Sugihara and May, 1990; Montero and Morin, 1992). Although characterizing dynamical systems using the analysis of uni-dimensional time series has been a method widely developed since the 1980s, there are several questions that need some consideration, at least in the cases when the indicators are obtained from time series resulting from behavioral investigation. In this type of investigation, like in most situations in real life, the data combine deterministic dynamics with noise of different nature and magnitude; in addition to this, in Psychology it is difficult to maintain the same observation situation for a long time and this leads to a reduction in the length of the series, therefore the reliability of such indexes can be questionable. In this study, we have tried to provide answers to the questions arising from the calculation of dominant Lyapunov exponent using direct methods. We have *Corresponding author. Tel.: +34-954557812. Fax: +34-954551784. E-mail: [email protected] ISSN 1026-0226 (C) 2002 Taylor & Francis Ltd
42 A.M.L. JIMNEZ et al. centered on analyzing the strength of the exponent for different values of evolution times and initial distances in short time series. LYAPUNOV EXPONENTS: DIRECT ESTIMATION The dominant Lyapunov exponent is one of the most widely used indicators to describe the qualitative behavior in a dynamical system using the analysis of unidimensional time series. To define what is understood by Lyapunov exponent (A) we start from an initial condition Yo, of a discrete dynamical system and we consider a very close point, where the initial distance (do) is extremely small. Let dt be the distance after iterations. If we assume that in Eq. (3) when do draws to 0, the term within the logarithm is the derivative of the iterate off evaluated in y0((ft)(y0)). Applying the chain rule of differentiation, the derivative off can be written as a product of derivatives of y) evaluated at the successive trajectory points yo,yl,y2,.., and so on. We can then define Lyapunov exponent in a more intuitive way with the following 1 A lim- In t-1 Hft(Yl) t=0 liml (lnlf’(y0)l + lnlf(y)l +... t--,eo + lnlf’(y,-l)l) Id, Id01 exp (tA) (1) then A is what we call Lyapunov exponent (Packard et al., 1980; Schuster, 1984; Montero and Moran, 1992; Nychka et al., 1992; Abarbanel et al., 1993; Simmons, 1993; Strogatz, 1994; Hilborn, 1994; Martin et al., 1995). That is, the average exponential rate of divergence or convergence of trajectories which are very close in phase space (Wolf et al., 1985; DeSouza-Machado et al., 1990; Zeng et al., 1991). The number of Lyapunov exponents in a dynamical equals the number of state variables considered. A unidimensional system is characterized by one single exponent. The signs of Lyapunov exponents provide qualitative information about the dynamics of a system. If the sign is positive, this is an indication of chaos. If it is negative, there is convergence between close trajectories and therefore classic attractors exist. If the behavior of a dynamical system represented by functionfconverges to a fixed point (y*) or to a limit cycle of a period p which contains a y* point, then it is easy to prove that Lyapunov exponent A < 0. When Lyapunov exponent Z 0 the initial perturbation will remain with t, i.e. trajectories neither diverge nor converge, their initial distance remains constant. This kind of behavior is typical of a constant periodical orbit (Sano and Sawada, 1985; Wolf et al., 1985; McCaffrey et al., 1992). If the system is three-dimensional (i.e. it contains three state variables) the possible combination of signs and the attractors they describe are: (+, 0, -), for a strange attractor; (0, 0,-), a quasi-periodical attractor known as torus; (0,-, -), limit cycle and ), a fixed point. If in Eq. (1) we take logarithms and we replace d, by its mathematical expression we obtain A -ln In (2) Y00 7 If the expression (2) has a limit as oo we define that limit to be the Lyapunov exponent: A= t-.oo lim-llnl t lim l ln ft(yo + do) ft(yo) do (3) lim 1 lnlf’(yt)l t---,oo t=0 (4) Equation (4) tells us that the Lyapunov exponent is the average of the natural logarithm of the absolute value of the derivatives of the discrete dynamical system evaluated in the trajectory points or time series considered. If the application of the discrete dynamical system for two close trajectories eventually leads us to separate points, then the absolute value of the derivative off is greater than 1 when we evaluated at those trajectory points. If the absolute value is greater than 1, its corresponding logarithm is also positive. If the points of the trajectory continue to diverge, then the average of the logarithms of the derivatives is positive. If we calculate Lyapunov exponent for a sample of initial points and we average the results, we can define the average Lyapunov exponent (X) for the system. An unidimensional discrete system has chaotic trajectories, for certain parameter values on which its behavior depend, if the average of Lyapunov exponents (X) is positive. The expression (4) for calculating A requires that the shape of the discrete dynamical system be known. But, what happens when we do not know the system but know one time series of a relevant state variable? Studies concerning non-linear dynamics have suggested two approximations for the estimation of Lyapunov exponent: direct methods or direct estimation and Jacobian methods (Guckenheimer, 1982; Eckmann and Ruelle, 1985; Wolf et al., 1985; McCaffrey et al., 1992; Nychka et al., 1992; Abarbanel et al., 1993; Damming and Mitschke, 1993). Direct methods calculate Lyapunov exponent directly from the time series, without any additional assumptions or approximations about the subjacent dynamical system. Since the subjacent system and its dimensions are unknown, we calculate, using the reconstruction vector method and for different dimensions, exponents until that ceases to vary significantly as the dimension of reconstructed space increases. Although the disadvantage of this procedure is that it only provides the largest Lyapunov exponent, we will discuss this method below.
LYAPUNOV EXPONENT IN SHORT TIME SERIES 43 ’J’o 20 FIGURE Scale region for a series of distances obtained from yt+l 4yt(1 yt). Let us define yl,Ym,...,Yr, as the elements of a time series and dmax as the maximum distance between two points of the series required to be considered as "infinitesimally" close. If the system behaves chaotically for points Yi and yj; with a distance do--< dmax, the divergence of close trajectories will be shown in the sequence differences do lyj Yil d lYj+ Yi+ll d2 [Yj+2 Yi+m[ dr ly+r yi+rl This will show an exponential increase, of at least the mean, T. With this method for calculating , we find two close trajectories in the state space and we calculate the series of distances, which derive from these two initial conditions. Although in general, calculating Lyapunov’s maximum exponent is easy, we believe that it is worth considering a few aspects. We assume a separation rate or exponential approximation between two close trajectories. For a given time series it is necessary to demonstrate this assumption. One way of doing this is by plotting the natural logarithm of the differences (ln dr) as a function of index T (see Fig. 1). If the divergence is exponential, the points (ln dr, T) will approximate a line given by the following expression in dr In do + AT The slope of the fitted line--generally by least squares--is then the value of Lyapunov exponent The dispersion diagram is never exactly a line because the exponential divergence varies along the attractor and reaches a maximum when it is comparable to the attractor’s diameter, which is defined as the maximum distance between the points of the trajectory evolved over the attractor. Two regions are usually distinguished: one that is called scale region to which the line is fitted and another one in which ln dr remains more or less constant as T increases. Adjusting a line with least squares to the scale region provides the measurement of the dominant Lyapunov exponent and the accuracy of the adjustment On the other hand, the values of )t can, and generally depend on Yi values chosen to be the initial conditions, or rather, the values of the initial distance (Id01) between them. To characterize the subjacent attractor of a given time series we have to calculate the average corresponding value for the set of Lyapunov exponents obtained from a number of trajectories which follow the condition (5). do lyi yjl <- dmax At a practical level, several questions arise concerning the time span required between points Yi and yj so that they can be considered as initial conditions of two trajectories, the length of the series, the number of initial conditions or distances, the number of iterations or optimal evolution time required and the initial distance in order to consider points as infinitesimally close in the phase space. We do not have many answers to the matters mentioned in the last paragraph. However, we have found some suggestions, which have resulted from some simulation experiments made with specific dynamical systems in which the theoretical value of Lyapunov’s maximum exponent is known, in this study, we have gathered some of the most general suggestions made in the research material that we have revised. The initial separation required between two points so that these can be considered as initial conditions of two different trajectories in the reconstructed phase space tends to be related to what is called orbital period (Wolf et al., 1985; Theiler, 1986). This is the time, which a system takes to cover an orbit. It is recommended that the initial separation between two points should be at least one orbital period. The difficulty for applying this recommendation lies in that the shape of the dynamical system generating the data must be known since it is on system itself that the calculation of the orbital period is based. In the cases for which the shape of the dynamical system is unknown, we can follow the recommendations given by Hilborn (1994) and Theiler (1986). They maintain that the initial separation required between two points so that they can be considered as initial conditions of two different trajectories, must be greater than what is known as auto-correlation time (-) which is given by the expression (6). 1 - (6) ln(1/p) In Eq. (6) p is the autocorrelation coefficient of lag 1. When p approaches 1, Eq. (6) changes and becomes Eq. (7). 1 - (7) As far as the number of initial conditions required is concerned, we have found several recommendations. Sano
44 A.M.L. JIMINEZ et al. o. 8i o o 4 o Io 20 30 4O so FIGURE 2 Evolution of the distances between two series generated from the logistic map for a 4 and do 0.001. and Sawada (1985) establish a lower limit which will depend on the dimension of the reconstructed space; for them, the number of initial conditions required for the estimation of Lyapunov’s maximum exponent must be higher or equal to the dimension of the reconstructed space (N >-- de). Hilborn (1994) recommends from 30 to 40 initial conditions distributed over the attractor. Other authors establish a dependence on the dimension of the reconstructed phase space and on the initial tolerance in order to consider points as infinitesimally close. The expression (8), by Dimmig and Mitschke (1993), provides the number of initial conditions required for characterizing the attractor according to Lyapunov’s maximum exponent N \dmax/ (8) where d, is the dimension of the reconstructed space and alma is the maximum initial distance required to consider two trajectories as close. Wolf et al. (1985) recommend from l0 de to 30 de For each of the series of distances we calculate the local Lyapunov exponent (Ado). The average Lyapunov exponent () is given by (9) where N is the number of series of distances considered. Another aspect to take into account is the number of iterations T, optimal evolution time or length of the series of distance that are suitable for calculating the Lyapunov exponent. In the case of chaotic systems, we know that convergence to theoretic values depends, among others variables, on the use of an evolution time (T) that is not too large. Thus, the exponent brings the sensitivity of the initial conditions and not the convergence produced by the boundary of the attractor values to a region of the phases space and of the initial distance (do) to consider two neighboring trajectories. We know that there is divergence in chaotic dynamics but at the same time, due to the folding mechanism, the values go through points that are infinitesimally close to previous values. In Fig. 2 we have shown the evolution of the distances between two series generated from the logistic map, for a 4 and for the initial conditions y01 0.1 and y02 0.101. A large T value may produce an underestimation of the Lyapunov exponent. Three criteria have been proposed to fix the number of iterations or evolution times (T): (a) To establish, a priori, a fixed evolution time (i.e. T 10), (b) To establish a final distance that can be when the attractor diameter or a percentage is reached. In this case, the time of evolution would be variable and (c) To identify in the graphic In dr before T in the scale region. Finally, the distance dmax, for the limited systems that interest us, cannot be too large. Due to the fact that the values Yi are constrained in size, the initial distances cannot be larger than the difference between the maximum value, Ymix, and the minimum value, Ymin. Moreover, there are practical limits when determining the initial distance for the finite precision of the data. The number of decimals is an inferior limit for the initial distance. For example, if the data is registered with three decimals, it would be senseless to question a difference lesser than 0.001. Another effect of the finite precision is that we can encounter repeated data. To summarize, the literature about non-linear dynamics only offers a few recommendations, which sometimes are divergent, about the criteria to be used in order to select the optimal calculus parameters in the estimation of Lyapunov exponents by direct methods. These few recommendations are circumscribed to the analysis of chaotic systems. We have found no recommendation for the estimation of A starting from the time series of classic systems. The reason for this is the interest in distinguishing variability due to a chaotic behavior of determinist dynamic systems of variability caused by white noise or linear stochastic processes, and less in the identification of nonlinear terms from the analysis of time series. In this study, we have centered in the dependence of the Lyapunov exponent, obtained by means of direct estimation, of the initial distance and the time evolution. We have used generated series of chaotic systems and generated series of classic systems with varying complexity. To generate the series we have used the logistic map (10). Yt+l ayt(1 Yt) (10) We know that the logistic equation is a structurally unstable system. That is, its behavior depends on the value of the parameter a. In Table I we have specified the values of a used in this research, its behavior and the diameter of the corresponding attractor (Hilborn, 1994; Strogatz, 1994). METHOD In this section, we describe the dependence of the Lyapunov exponent both on the initial distance and the
LYAPUNOV EXPONENT IN SHORT TIME SERIES 45 TABLE Value of a Behavior Diameter of attractor Lyapunov exponent 3.2 Limit cycle (p 2) 3.52 Limit cycle (p 2) 3.55 Limit cycle (p 8) 3.58 Chaos (low extent) 4 Chaos (high extent) 0.2910.51,0.81] -0.91 0.5110.37,0.88] -0.19 0.53[0.36,0.89] -0.1 0.56[0.34,0.90] 0.1 1[0,1] 0.69 evolution time in series generated by the logistic map. We wrote a program in the Mathematica programming language (v. 2.1 for Windows) and we generated series with close initial values. The proximity criteria used was the one recommended by Sano and Sawada (1985) according to which two generated series of a dynamic system are considered to be close if the initial distance between them (do) is between 1 and 5% of the attractor diameter. With the program made to generate the data and given that sensitivity or insensitivity to initial conditions is a characteristic of a group of trajectories with close initial conditions, for the values a 3.2, 3.52, 3.55, 3.58, 4} We initially generated 50 series of N-- 100 data with initial conditions (Yo) whose distance was 1% of the attractor diameter (see Table I). From the 50 series generated we obtained: 49 series of distances with do 1%, 48 with do 2%, 47 with do 3%, 46 with do 4% and 45 with do 5%. The initial value for the first series generated (Yo) coincided with the minimum value of the attractor amplitude (see Table I) with the exception of the series generated for a 4, where we started from y0-- 0.1. For this last value, we used three additional distances among the initial conditions 0.1, 0.25 and 0.5% and 50 distance series were generated. The Lyapunov exponent for 1 <-- T --< 99 and the initial distance (do) were obtained by means of the expression A=ln In (11) Yio YjO -o We calculated the average Lyapunov exponent (X) as estimator of de A of the coefficients obtained with expression (11) for each do and T. We graphically represented X as a function of index T for each value of do in order to study the incidence of these parameters. RESULTS We organized the results in two different sections. In the first part, we show those results corresponding to the T T -,2 e)e 3.2, do=5 files 0,04 FIGURE 3 Behavior of X according to T when increasing the do for a 3.2.
46 A.M.L. JIMNEZ et al. FIGURE 4 Behavior of X when increasing do to 10 and 20%. values of a whose behavior converges on a classic attractor. In the second part, we offer the results for the values of a with a chaotic behavior. We have preferred to show the graphics corresponding to the evolution of X facing T instead of the value tables because we see them as more illustrative and easier to interpret. Classic Attractor a 3.2 In Fig. 3, we represented the average Lyapunov exponent (X) as a function of the evolution time (T) for the set of distances that were considered. In the different graphics of Fig. 3, we have drawn a broken line through the value of the theoretic Lyapunov exponent. The convergence to this value can be seen with the increase of the evolution time T. The convergence form is independent from the initial distance (do) considered. In all cases there is an oscillatory decrease of X. For the set do {1%,2%,3%} the convergence stops in T 14, obtaining the value of X closest to the theoretic when T 13. -,z, - ba=3.52, d=4%, Bias=O.01, T,,=17 ,2 0,0 -,2 T c)==3.SZ = =0.00 ’.=7 d) 3.52, d 4%, T e) 3.52, d= 5, Bias 0.05, T,, 17 FIGURE 5 Behavior of X according to the evolution time (T) for a 3.52.
LYAPUNOV EXPONENT IN SHORT TIME SERIES 47 T b) 3,55, do 2%, 8as O.01, T c) 3,55, d 3%, Odas <0.0I T d) SS, d,=,m, nab 002 Bias 0,01, T.-7 FIGURE 6 Behavior of X according to the evolution time (T) for a 3.55. For the set do {4%, 5% the convergence stops in T 16. The closest mean Lyapunov exponent (A) is obtained for T 15 and T 13 (the bias is 0.04). To see if increasing the initial distance would decrease the bias (X-A), we analyzed in 16 and 15 series of distances the behavior of A according to the evolution time (T) when increasing the initial distance (do) between the trajectories to 10 and 20% of the attractor diameter. We have presented the results in Fig. 4. It can be observed (graphic (a)), how the value of the bias and the shape of the convergence are similar to those obtained for the initial distances analyzed previously. On the contrary, when the value for the initial distance increases to 20% the bias increases dramatically to 0.51. a 3.52 Figure 5 shows the behavior of A according to the evolution time for the unit of initial distances considered: do {1%,2%,3%,4%,5%} As with the series represented in Fig. 3, for a 3.52 the mean Lyapunov exponent declines oscillatory when T increases. Moreover, some differences can be observed in relation to do. In the graphics of Fig. 5, we have drawn a continuous line. This is perpendicular to the axis of the intersection of the value of T, which provides the best estimation of A. As is customary, a broken line represents the theoretic value of the Lyapunov exponent. We can see (graphics a, b and c in Fig. 5) that X tends to the theoretic value when the initial distances are 1-3%. When increasing the initial distance to 4 and 5% (graphics d and e in Fig. 5), the limit of X is not the theoretic value but a larger one. For the values of distances and evolution times considered, the direct estimation provides values that are biased positively of A. a 3.55 Figure 6 shows the behavior of X according to T for the distances considered. In the graphics in Fig. 6, we have drawn an ordinate line A--0 (the gray line) with the aim of assessing possible qualitative errors when identifying the subjacent attractor to the data generating system. The To (optimal) is the values of T for which the bias is less. In graphs (a) and (f) in Fig. 6, we can see an initial period with great variability in which A oscillates between
48 A.M.L. JIMINEZ et al. T 3.56, dn 1% slow and practically stabilizes itself oscillating between -0.04 and -0.02 when the initial distance is 1%, as can be seen in Fig. 7. For the rest of the do with the evolution time, oscillation stabilizes between 0.03 and 0.001. In any case, we can deduce from the Fig. 6, that using short and even evolution times (first section of the dispersion diagrams) can be more problematic than using longer evolution times as, although the bias increases with T, there are no qualitative errors noted. FIGURE 7 Oscillations of X in the interval: T [60, 100]. the zone of chaotic behavior (A > 0) for the even T values, and the zone of recurrent behavior (X > 0) for the odd values. After this first span, when T grows to the value of A that is higher than the theoretical value (X -0.1), the direct estimations for converge slowly and in an oscillating way. The increase in the initial distance process brings the boundary of the convergence process of X to 0. In the interval 5 -< T --< 11 we can find the values that provide the best estimations of . A common value for To that provides reliable estimations (bias -< 0.02) for the set of distances under consideration is To 7. From this initial period, both the decline of ,( for even values of T, as well as the growth for odd values, is very Chaotic Behavior a 3.58 Amongst the many values of a whose behavior is chaotic in the interval [3.58, 4], we used precisely the extremes, which correspond to the strange attractor of the smallest and largest diameters, respectively. Figure 8 shows the results relative to the form of the dependence of X as opposed to the evolution time. In the set of graphs in Fig. 8, we can distinguish three sections: in the first section for approximately 1 <_ T <-- 4, we can see a sharp increase of X towards the theoretical value. T 0 20 40 0 80 00 b)a=3.S& do=2% To={4,8}.Bla$=O .58. =3 %, T,., 14, 8) I 0.01 of)a=3.58, do=4%, To= [4. 8), Bi= e)==3.,.’,8, do=5%. T.,={4.S}. IIZI .c0.01 do T 7,00 8,00 FIGURE 8 Behavior of X according to the evolution time (T) when increasing the initial distance for a 3.58.
LYAPUNOV EXPONENT IN SHORT TIME SERIES 49 ,10 ,05 ,05 .04 .02 .01, -.01 -,0 -,0 T c) 3.58 cJo -.3% ,05 ,03 01o -.02 -,03 -,04 " -,0 d)a 3.5 do 4% ,04 ,02 ,01, ,00 -,01 -,0. -,0 -,04 ,0 ::N FIGURE 9 Behavior of X for T > 8 for a 3.58. We can consider the second section as stabilization, in which A oscillates around a constant value very near to the theoretical value for the even values of T. This second span becomes progressively shorter as the initial distances grow. Thus, for do--1%, the values of X remain practically constant until T 25, whilst for do 5%, this section finishes in T 8. During this time of stabilization, as the initial distance increases, error is more likely in identifying the type of attractor present in the data, especially for the odd values of T. The third section, which is clearly identifiable in the figures, is the longest. In this section, A oscillates with variations that are practically constant in size, around 0, independent of the initial distance. The increase in evolution time to higher than optimal T underestimates the value of X It seems that with T, the average exponent X would be centered on 0 for all the initial distances. We have widened this section specifically in order to observe whether the behavior in this area is really independent of the initial distance or not. The results are shown in Fig. 9. Taking the values of X for T > 8, we can see that, with the initial distances do increases the possibility of qualitative error when using longer evolution times. In graph (a) in Fig. 9, practically 100% of the A values are greater than 0. Here, although the bias is sharp, nevertheless the qualitative conclusion regarding the nature of the attractor contained in the series would be correct. Graphs (b)-(d) in Fig. 9 shows how the distribution of values around A--0 inverts with the increase in the initial distance and the number of negative exponents increases progressively. The behavior of X in relation to T is similar to that observed for the values of a with recurrent behavior at intervals of 2, 4 and 8. In graph (f) in Fig. 8 we can see how, for the values of T that constitute the first and second section of the evolution of X for the set of initial distances, the bias is independent of these. a--4 Graphs (a) and (h) in Fig. 10 show the dependence of X with respect to the evolution time for the set of do that was applied. In graphs (b)-(g) (see Fig. 10) we can appreciate a very rapid initial increase which is more or less lineal, towards the theoretical value of X. This growth is interrupted abruptly in A values below the theoretical value. For growing values of the set of initial distances do--= 1%, 2%, 3%, 4%, 5% }, the maximum X values reached are: 0.59, 0.58, 0.55, 0.55 and 0.53. It appears that we can confirm a direct relationship between the bias