Abstract

A simple and effective algorithm for the identification of optimal time delays based on the geometrical properties of the embedded attractor is presented in this paper. A time series synchronization measure based on optimal time delays is derived. The approach is based on the comparison of optimal time delay sequences that are computed for segments of the considered time series. The proposed technique is validated using coupled chaotic Rössler systems.

1. Introduction

The reconstruction of geometry from a time series is one of the paradigmatic algorithms used for the computational analysis of nonlinear dynamical systems. There is a number of algorithms for the analysis of time series that consider phase spaces, that is, the set of all possible system states, of dynamical systems. For a deterministic dynamical system, the phase space is known from the mathematical model (equations of motion). Thus, with the knowledge of the system state at a particular time, it is always possible to determine the future states. However, real-world or experimental dynamical systems are usually too complex, and phase spaces of such systems are unknown due to the chaotic nature of the analyzed data. Therefore, phase space reconstruction methods had to be developed. A branch of applied dynamical system theory, analyzing such algorithms, that is, extracting the information about geometrical and topological properties of the phase space of time series obtained by performing measurements on an evolving system, is called embedology [1]. One of the first methods for reconstructing the phase space system from a time series by using a time delay variable was proposed several decades ago [2]. A complex time series has an inherent geometry, and Packard et al. [2] were first to show that a representative geometry of a dynamical system can be obtained by using a time series of one of its observables.

The much celebrated Takens embedding theorem proved that it is possible to reconstruct the attractor of a dynamical system using only a single time sequence of scalar measurements [3]. The vectors accomplishing this reconstruction have the form , where are the observations on the system, the integer is the embedding dimension, and is the time difference between consecutive components (also called the time lag or the delay time) and is a multiple of the sampling time step. Takens assumed an infinite sequence of noise-free measurements and proved the existence of a diffeomorphism between the original and the reconstructed attractors for almost any choice of positive delay times and a sufficiently big dimension .

Although the application of the Takens time delay technique is straightforward for almost any time series, the selection of optimal reconstruction parameters and is a nontrivial problem due to their dependence on the nature of the analyzed data. Thus, numerous studies have been made in this area. The selection of optimal time delay value is usually based on the optimization of a particular target function that corresponds to some measures of the quality of phase space reconstruction. Classical statistical approaches for the construction of the target function include the autocorrelation function method and the mutual information method [4–6]. Such techniques are based on maximizing the independence between the coordinates of the delay vector. This approach is heavily reliant on the geometry of the phase space reconstruction, which is appropriate for the study of parameters such as fractal dimension and Lyapunov exponents. However, the selection of such target function does not always yield best results for nonlinear time series prediction problems [7].

Another class of methods for the selection of an optimal is based on the concept of phase space expansion. The “fill factor” as a spatial measure of the phase space was first introduced by Buzug and Pfister in [8, 9]. The fill factor quantifies an attractor’s utilization of the embedding space as a function of . While this method does not possess the drawbacks of statistical techniques with respect to nonlinear data, it is computationally intensive (for ) and is prone to overfolding of the attractor [10]. A more computationally effective algorithm based on the expansion of the area for all pairwise planar projections of the embedded attractors is presented in [11]. Similar approaches based on the spatial properties of phase spaces include the average displacement method [10], the SVF method [12], and the wavering product [13] method.

A different viewpoint considers both phase space reconstruction parameters simultaneously since some experiments indicate that an irrelevant relationship between and may impact the congruence between the original system and the reconstructed phase space. This class of methods includes the small-window solution [14] and the C-C method [15].

Novel approaches to the selection of embedding parameters have been developed in recent years. It is shown in [16] that an improved phase space reconstruction method based on manifold embedding and Laplacian eigenmaps can reveal the structure of the chaotic attractor. A new algorithm for the estimation of the dimension of chaotic dynamical systems using neural networks and robust location estimate is proposed in [17]. Statistical analysis of the nearest neighbors is exploited for the reliable estimation of minimum embedding dimension for noisy time series in [18].

It has been demonstrated in [11, 19–21] that nonuniform embedding (when all time delays are not equal) performs better than uniform embedding in a variety of applications—typical examples are causality and coupling detection and time series prediction. However, the procedure for the selection of time delays is usually related to a particular application and is implemented by introducing a specific target function which determines the utility generated by a concrete nonuniform embedding.

This paper has two main objectives. The first objective is to introduce a simple and effective algorithm for the identification of optimal time delays which is based on the geometrical properties of the embedded attractor and is applicable for short time series.

The second objective is to introduce a time series synchronization measure based on optimal time delays. A short review of commonly applied methods for the detection of synchronization is given below.

The first linear synchronization measures based on correlation analysis, namely, cross correlation and coherence functions, are widely applied due to their computational effectiveness [22]. However, such measures can only detect the most straightforward regimes of synchronization.

Nonlinear synchronization measures were introduced in order to quantify more complex synchronization effects. The mutual information measure is based on Shannon entropy and takes into account not only linear but also nonlinear dependencies [23]. Phase synchronization measure is used to quantify similarity between cyclic signals and time series. Two approaches—based on the Hilbert or wavelet transform—can be used to implement the phase synchronization detection algorithm [22].

A new form of synchronization between coupled chaotic oscillators called amplitude envelope synchronization has been discovered in [24]. Generalized synchronization of dynamical systems occurs if dynamical variables from one subsystem are a function of the variable of another subsystem [25]. A nonlinear interdependence of dynamical systems based on state space reconstruction is a similar approach to generalized synchronization but does not require a functional relationship between the dynamics of the underlying systems [26]. The function projective synchronization in relay coupled systems is studied in [27]. A novel sort of synchronization called complex antilag synchronization is introduced in [28]. A new approach for the investigation of hybrid chaos synchronization in discrete-time hyperchaotic dynamical systems based on stability theory of linear discrete-time systems and Lyapunov stability theory was proposed in [29].

This paper presents a novel time series synchronization measure based on the geometrical approach towards optimal time delays. This approach is based on the determination of time lags that maximize the volume of the state space occupied by the embedded segments of the time series. The sequences of obtained time lags reveal the level of synchronization between 2 time series. The presented algorithm for the identification of synchronization between two time series is validated using coupled chaotic Rössler systems.

2. Preliminaries

2.1. Two-Dimensional Delay Coordinate Space

Let us consider a harmonic function , where is the angular frequency and is the phase of harmonic oscillations, and the amplitude is set to 1. A pair of function values and , where is the time lag, is mapped into a point in the phase plane , where and are the coordinates of the delay coordinate plane. Two-dimensional time delay embedding maps the harmonic function into an ellipse; the equation of that ellipse reads

The geometrical shape of the ellipse can be exploited for the quantification of the quality of the embedding. Such geometric approach (for a two-dimensional phase plane) was firstly proposed by Buzug and Pfister in 1992 [9]. The radiuses of the embedded ellipse and can be directly expressed from (1):

Thus, the area of the ellipse reads

The maximum area of the ellipse is (when the ellipse becomes a circle). Now, the function representing the quality of embedding into a two-dimensional delay coordinate space can be defined as a ratio between the area of the ellipse and the area of the circle:

2.2. -Dimensional Delay Coordinate Space

2.2.1. The Embedding Quality Function

The geometrical approach for embedding a harmonic function into a two-dimensional phase plane is generalized for the -dimensional delay coordinate space in [11]. The coordinates of the reconstructed point in the -dimensional delay coordinate space read where is the time delay vector; is the total embedding window. It can be observed that there exist planar projections of the embedded ellipse. The function representing the quality of embedding into the -dimensional delay coordinate space is constructed as an arithmetic average of all quality functions is all possible planar projections [30]:

In the case of the uniform embedding (when ), (6) reads

The equality implies that the harmonic function is compressed into line segments in all possible planar projections. The maximization of in respect of time lags results in a maximum average area of ellipses in all possible planar projections.

2.2.2. A Nonharmonic Function

Let us consider a function defined in the time interval , where and are the limits of the observation window. Every harmonic component of is affected by the quality function when the appropriate harmonic signal is embedded into the -dimensional delay coordinate space. Harmonic components with frequencies where is small will be suppressed (in average) in all possible planar projections more than those harmonic components where is large. Such motivation did help to construct the optimization problem for the identification of an optimal set of time delays which do result in the richest representation of the attractor in the delay coordinate space: where

is the optimal set of time delays and is the Fourier amplitude spectrum of the signal in the observation window [11]. The integral in the denominator is used in order to normalize the target function in respect of the signal (this integral can be computed once at the beginning of the optimization process). The multiplier is used to normalize the target function in respect of the white noise (the value of the target function for the white noise signal now becomes equal to 1 for any time delays not equal to 0) [11].

3. The Proposed Algorithm

The technique proposed in [11] enables a fast and effective determination of the optimal set of time delays by assessing the geometrical properties of the embedded attractor. However, the optimization problem in [11] does not asses the phase spectrum of the Fourier transform. In other words, [11] is an approximate algorithm for the determination of optimal time delays when every discrete harmonic component of the Fourier amplitude spectrum is treated as a separate individual harmonic component. Nevertheless, such an approximating approach appears to be useful in practical applications—especially in time series forecasting algorithms based on fuzzy neural networks [31–33]. However, it is clear that the Fourier phase spectrum does have an impact on the geometrical representation of the reconstructed attractor in the phase space. In other words, a more accurate assessment of the geometrical properties of the reconstructed attractor is required. Let us consider a continuous function in the observation window and the -dimensional delay coordinate space with the set of time lags . Then, instead of computing any projections of the embedded attractor into a 2-dimensional phase planes, we compute the average distance of the attractor points to the origin of the embedding frame:

3.1. Properties of : Embedding a Harmonic Function

Let us consider that the embedding is uniform and the embedded function is a harmonic function , where , , and are the amplitude, the circular frequency, and the phase, accordingly. Then,

The change of variables yields

Now,

Let us denote . Then,

Note that when . Thus, is the maximal value for the harmonic function. On the opposite, the pole is a removable singularity at and is the minimal value for the harmonic function.

3.2. The Comparison between and

We continue with the harmonic function and uniform embedding. Then, and (9) yields , where is the Dirac delta function. Without the loss of generality, we assume that and . Then, ; . Graphs of and are illustrated in Figure 1.

(a)

(b)

(c)

(d)

(a)
(b)
(c)
(d)

Figure 1 

The comparison of functions (solid line) and (dashed line). The parameter is set to 2, 3, 4, and 20 in (a), (b), (c), and (d), respectively.

The first observation is that graphs of and are similar at (Figure 1(a)). The best embedding of the harmonic signal into the 2-dimensional delay coordinate space is produced at , (then the harmonic function is embedded into a circle). The harmonic function is embedded into a line interval on one of the diagonals of the frame at , . The difference is that now , but (note that the scale of the y-axis for is shown on the left, but the scale for is shown on the right side of Figure 1(a)).

The differences between and are more clear at (Figure 1(b)). The harmonic function cannot be mapped into circles in all 3 projections at a fixed ; therefore, the maximum value of at is lower than that at (note that the scales of the axes are the same for different in Figure 1). But remains the same and does not depend on the dimension of the delay coordinate space. Both and yield the best embedding at , where , .

The maximum value of at becomes smaller than the one at (Figure 1(c)). But now a major difference between and can be observed in terms of the best embedding. yields the best embedding at , but yields the best embedding at , where , . Clearly, misses the time delay . That can be explained by the fact that measures the average area of the mapped attractor in all possible plane projections of the -dimensional delay coordinate space. On the contrary, directly measures the geometrical volume occupied by the attractor without considering any projections of that attractor.

The maximum value of tends to 1 when (Figure 1(d)). Note that when the embedded time series is the white noise [11]. However, as mentioned previously, is an invariant value in respect of . What is more important, yields the best embedding at , where , , which is in stark contrast to (Figure 1(d)).

3.3. Embedding a Discrete Time Series

Let us consider a finite scalar discrete time series , where is the number of observation points uniformly spread in the time axis. Let the sampling rate of the time series (the time interval between adjacent observations) be . Then, all time delays are multiple values of . Let us denote . Then, (10) reads

Note that does not play any role on the right-hand side of (17)—because the definite integral in (10) is replaced by the limit sum in (17). Now, the problem of the identification of the optimal set of time delays is formulated as follows: where is the upper limit for every individual time delay. The total embedding window (the time delay between the first and the last time delay coordinates) is , and it could be also used for the estimation of . It is clear that must be greater than . However, must be at least 10 times greater than in order to produce a meaningful optimization result (18).

3.4. Properties of : Embedding the Random Noise

It is well known that embedding of the random noise is invariant to time delays. In other words, the geometric shape of the reconstructed attractor does not depend on time delays (if only the embedded time series is an ergodic random process and the realization of that process is long enough). This geometric feature is exploited in one of the tests used to measure the quality of random number generators. For example, if a random number generator yields a sequence of numbers distributed uniformly in the interval , then the density of the embedded points should be the same in any location of a two-dimensional square for a sufficiently long sequence and any time delay greater than 1 [34]. It is clear that must be invariant to time delays if only the embedded time series is random. We demonstrate the properties of for the random time series distributed uniformly in the interval and for the Gaussian time series. Of course, the invariance in respect to time delays does hold, but some features of are interesting and are presented below.

3.4.1. Embedding the Uniform Noise

Let us consider a discrete random variable distributed uniformly in the interval . Initially, let us consider and . In fact, this situation cannot be considered an embedding because all values of the time series are mapped on a single one-dimensional axis without any time delays. However, the formal result reads

Now, let us consider the situation when (for any ):

Thus, finally, if the embedded random time series is distributed uniformly in the interval .

3.4.2. Embedding the Gaussian Noise

Now, let us assume that a discrete random variable is the Gaussian noise, with the mean value and standard deviation . Firstly, let us consider and :

Now, let us consider the case when (for any ):

Thus, finally, if the embedded random time series is the Gaussian noise, with the mean and the standard deviation .

3.5. Computational Experiments: Optimal Embedding of Mackey-Glass Time Series

Equation (18) is a nonlinear integer optimization problem. The full sort method could be used for the identification of the optimal set of time delays—only if , , and are not very large and computational resources are not strictly limited. Otherwise, soft computing techniques could be exploited for the determination of the near-optimal set time delays [30]. Let us consider Mackey-Glass time series [35] the discrete numerical solution to chaotic Mackey-Glass delay differential equation that reads where , , and ; the observation window is and the time step (Figure 2). The standard FNN algorithm [36] yields (5) as the minimum embedding dimension for the investigated chaotic Mackey-Glass series.

Figure 2 

Mackey-Glass time series at , , and .

Now, we use (18) to determine the optimal uniform time lag in the interval (Figure 3). The maximum value of is reached at (Figure 3). Note that this value of can be used as a lower bound for the integer programming problem when the near-optimal set of nonuniform time delays is sought using branch-and-bound [37] or other soft optimization algorithms.

Figure 3 

Plot of for the Mackey-Glass time series (uniform embedding in the 5-dimensional delay coordinate space). The maximum value of is attained at .

The full sort method is used to determine the optimal time delay vector for nonuniform embedding. 100⁴ different combinations are assessed in 8.14 hours on a two-kernel 2.20 GHz processor with 8 GB RAM. As mentioned previously, soft computing techniques can be effectively exploited for the determination of the near-optimal set of time delays. Finally, the optimal set of time delays reads ; . All possible planar projections of the embedded Mackey-Glass attractor are illustrated in Figure 4.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)

Figure 4 

All possible planar projections of the embedded Mackey-Glass attractor occupy the largest area (in average) at , , , and .

3.6. Synchronization Measure Based on a Geometric Approach to Attractor Embedding

The proposed technique for the determination of optimal time delays is based on a straightforward geometric approach. Such an approach can be effectively exploited for comparing geometric properties of attractors embedded from different time series. Moreover, we will demonstrate that such an approach is working well with short time series.

The basic idea of the proposed synchronization feature is based neither on the direct comparison between geometric shapes of the reconstructed attractors nor on global characteristics of those attractors (like averaged maximum Lyapunov exponent). We will use the optimal time delay as a single numerical characteristic describing the geometry of the reconstructed attractor in a 2-dimensional delay coordinate space. The reasoning for such a measure is straightforward. It follows from the properties of the reconstructed attractor (Figure 4). Such synchronization measure is much more general compared to synchronization measures used to detect identical, phase, amplitude envelope, lag, or even generalized synchronization. Two attractors are considered to be similar if their optimal time lags are the same. The higher the difference between optimal time lags, the more different these attractors are considered (from the geometrical point of view).

We will demonstrate that such an approach for detecting synchronization between two different time series is highly feasible—especially if those time series are far from being similar.

Let us consider two paradigmatic chaotic Rössler systems with diffusive coupling: where , , , and ; is the coupling parameter. The initial condition is set to , , and .

Attractors of both Rössler systems are illustrated in Figure 5. The difference between the uncoupled attractors can be seen by a naked eye (Figures 5(a) and 5(b)). However, the attractors become similar as the coupling increases (Figures 5(e) and 5(f)).

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 5 

Coupled Rössler attractors (26) in the phase space: (a, b) uncoupled systems at ; (c, d) coupled systems at ; (e, f) coupled systems at (, , , and ).

Analogously, attractors of coupled Rössler systems (26) are visualized in the embedding space (Figure 6). Note that a short observation window (400 scalar data points) is used to reconstruct the attractor in the planar delay-coordinate space. One time delay is used for every embedding; optimal values of time delays are depicted in Figure 6. Note that both the geometric shape of the embedded attractors and the optimal time delays become closer as the coupling parameter grows (Figure 6).

(a)

(b)

(c)

(d)

(e)

(f)

(a)
(b)
(c)
(d)
(e)
(f)

Figure 6 

Coupled Rössler attractors (26) in the planar delay coordinate space: (a, b) uncoupled systems at ; (c, d) coupled systems at ; (e, f) coupled systems at (, , , and ).

It is also interesting to compare the maximal Lyapunov exponents calculated for the system (26) and the one in the embedding space. Averaged maximal Lyapunov exponents for the uncoupled Rössler systems are computed separately: for the first system and for the second system. However, the coupled Rössler systems (26) are considered a single system. The 6-dimensional Jacobian matrix of the coupled system yields (after averaging) a single maximal Lyapunov exponent: at and at . It can be observed that the single maximal Lyapunov exponent of the coupled system at approaches the average of and . But the maximal Lyapunov exponent decreases (the coupled system becomes less chaotic) at .