Recurrence Plots and Cross Recurrence Plots

Introduction To Cross and Joint Recurrence Plots


Cross recurrence plot - a cross recurrence plot (CRP) is a graph which shows all those times at which a state in one dynamical system occurs simultaneously in a second dynamical system. With other words, the CRP reveals all the times when the phase space trajectory of the first system visits roughly the same area in the phase space where the phase space trajectory of the second system is. The data length of both systems can differ, leading to a non-square CRP matrix.

Joint recurrence plot - a joint recurrence plot (JRP) is a graph which shows all those times at which a recurrence in one dynamical system occurs simultaneously with a recurrence in a second dynamical system. With other words, the JRP is the Hadamard product of the recurrence plot of the first system and the recurrence plot of the second system. JRPs can be computed from more than two systems. The data length of the considered systems has to be the same.

Note: The following text is from 2000 and might be slightly outdated.

Following Takens' embedding theorem (1981) we reconstruct the phase space from a time series \(u_k\) by using an embedding dimension \(m\) and a time delay \(\tau\) $$ \vec{x}(t) = \vec{x}_i = (u_i,u_{i+\tau}, \ldots, u_{i+(m-1)\tau}), \qquad t = i \Delta t $$ whereas \(\vec{x}(t)\) is the vector of reconstructed states in the phase-space at the time \(t\). The choice of \(m\) and \(\tau\) should base on usually methods for detecting these parameters like method of false nearest neighbours (for \(m\)) and mutual information (for \(\tau\)), which ensures the entire covering of all free parameters and avoiding autocorrelated effects (Kantz and Schreiber, 1997). Eckmann et al. (1987) introduced the recurrence plot \(R_{i,j}\) (Fig. 1). This is an square array, that visualizes when a state in the phase space will nearly recur. Every point of the phase space trajectory \(\vec{x}_i\) (\(i = 1, \ldots, N; N = n-(m-1)\tau\)) is tested whether it is close to another point of the trajectory \(\vec{x}_j\), i.e. the distance between these two points is less than a specified threshold ε. In this case the value one (a black dot, recurrence point) in a \(N \times N\)-array $$ R_{i,j} = \Theta(\varepsilon_i - \|\vec{x}_i-\vec{x}_j\|) $$ (\(\Theta(\cdot)\) is the Heaviside-function) is placed otherwise the value zero (a white dot) is placed. Such recurrence plot represents the behaviour of a single trajectory in the phase space with typical long-scale and small-scale structures. A recurrence point contains no information about the special state.


recurrence plot
Figure 1: Recurrence plot of the Southern Oscillation Index (SOI). The diagonal structures represent the deterministic dynamics of the El Ni no/ Southern Oscillation (ENSO). From the distance between diagonal structures one can reveal cyclicities of recurrence of a especially state of ENSO.

As an extension of these recurrence plots, the recurrence quantification analysis (RQA), was introduced by Zbilut and Webber (1992). It defines measures for diagonal segments in a recurrence plot, recurrence rate, determinism, averaged length of diagonal structures, entropy and trend. The recurrence rate is the ratio of all recurrent states (recurrence points) to all possible states and is the probability of recurrence of a special state. Stochastic behaviour causes very short diagonals, whereas deterministic behaviour causes longer diagonals. Therefore, the ratio of recurrence points forming diagonal structures to all recurrence points is a measure for the determinism in the system. Diagonal structures show the range in which a piece of the trajectory is rather close to another piece of the trajectory at different time. The diagonal length is the time how long they will be close to each other and can be interpreted as the mean prediction time. The inverse of this measure is correlated with the Lyapunov-exponent (Eckmann, 1987). The entropy means the Shannon entropy of the histogram of line segment lengths and reflects the complexity of the deterministic structure in the system. Stationary systems will deliver homogeneous recurrence plots, instationary systems causes changes in the distribution of recurrence points in the plot which is visible by brightened areas. For example, a simple drift causes a brightening of the recurrence plot away from the main diagonal to edges. The value trend measures this effect by diagonalwise computation of diagonal recurrence density and its linear relation to the time distance of these diagonals to the main diagonal.

Analogous to Webber and Zbilut (1998) we expanded the method of recurrence plot to the method of cross recurrence plot, which compares the dynamical behaviour of two time series, which are embedded in phase space (Marwan, 1999). Here two time series are simultaneously embedded in the phase space. The test for closeness of each point of first trajecory \(\vec{x}_i\) (\(i=1,\ldots,N_x\)) with each point of second trajectory \(\vec{y}_j\) (\(j=1,\ldots,N_y\)) results in a \(N \times N\)-array $$ CR_{i,j} = \Theta(\varepsilon_i - \|\vec{x}_i - \vec{y}_j\|), \qquad \vec{x}_i, \vec{y}_j \in \mathbb{R}^m, $$

which is called cross recurrence plot. Long diagonal structures show similar phase space behaviour of both time series. We have modified some RQA measures for quantifying the phase space similarity. The recurrence rate \(RR\), determinism \(DET\) and the averaged diagonal line length \(L\) were determined as a function of the distance from the main diagonal, e.g. the diagonal-wise recurrence rate would be $$ RR_k = \frac{1}{N-k} \sum_{j-i=k}^{N-k} CR_{i,j}. $$ It is therefore possible to assess the similarity in dynamics depending on a certain delay. The second time series was embedded with positive and negative sign, what gives us the type of a possible correlation.

An example illustrates the application of our method. We use a linear correlated noise process (auto-regressive process) and couple it non-linear with the \(x\)-component of a Lorenz-System \(x_n\) (the quadratic exponent in the coupling term) to a second order autoregressive process (Fig. 2) $$ y_n = 1.095 y_{n-1} - 0.400 y_{n-2} + 0.700 \xi_n + 0.300 x_n^2, $$ (where \(\xi\) is white noise and \(x\) is normalized to standard deviation \(\sigma\)). The coupling is realized without any lag.


example system
Figure 2: Time series of a nonlinear related system consisting of a driven second-order autoregressive process, forced by the squared \(x\)-component of the Lorenz system.

The linear correlation analysis is not able to detect any significant coupling or correlation (Fig. 3b). In contrast, our new method finds a positive coupling with a lag of zero (Figs. 3c-e).

cross recurrence plot
Figure 3a: Cross recurrence plot between the driven second-order autoregressive process and the forcing function (\(x\)-component of the Lorenz system).


cross recurrence quantification analysis
Figure 3b-e: Crosscorrelation (b), \(RR\) (c), \(DET\) (d) and \(L\) (e) between the driven second-order autoregressive process and the forcing function (\(x\)-component of the Lorenz system). The solid lines show positive relation, dashed lines show negative relation. The point-dashed lines represent the significance levels. The crosscorrelation function find none significant correlation but the \(RR\), \(DET\) and \(L\) functions show a positive relation with a lag of zero.

We prepare a test for significance with a set of \(N\) surrogate data. On the one hand these surrogates should reveal some features like our original data but on the other hand features which represent the randomness of possible correlation (stochastic processes). Linear correlated noise are such systems. We calculate a new surrogate time series for such system with the following iteration function (autoregressive process of 3rd order) $$ x_n = a_1 x_{n-1} + a_2 x_{n-2} + a_3 x_{n-3} + b \xi_n, $$ whereas \(\xi\) is white noise, \(n\) runs from 4 to the value of the desired length of the time series and \(a_k\) are coefficients which determine the behaviour of the system and allow us to adapt this stochastic system to our natural systems. After this calculation we apply our method to the natural data (e.g. SOI) and to the surrogate data and determine the introduced measures (\(RR\), \(DET\), \(L\). We repeat this procedure (calculation of surrogate data and application of the method) for \(N\)-times. The result is a distribution of our measures for the unrelated systems. For each measure we obtain the \((p-1)\)-quantiles for the error level of \(p\) (e.g. 10%), which is the one-side significance level for this error level. With this significance levels we assess the significance of the revealed measures of CRP and the interrelations of the natural systems.

An further interesting feature of cross recurrence plots is the line of synchronization (LOS). This line reveals the relationship of the both systems in the time domain.

Lines of synchronization in cross recurrence plots. Lines of synchronization in cross recurrence plots.
Figure 4: Lines of synchronization in cross recurrence plots between harmonic functions, where the second function has a modification in the time domain (black: no modification, blue: small variation, red: stronger variation in the time domain of the second function).

animated plot

A fitting of this LOS with a nonparametric function allows to re-synchronize both systems. However, both systems should have a similar dynamical evolution. Applications e. g., in adjustment of geophysical data which were gained from different boreholes or cores.

Another multivariate approach to recurrence plots was introduced by Romano et al. and is called joint recurrence plot (2004), and is simply the product of two (or more) recurrence plots of the data series $$ JR_{i,j} = \Theta(\varepsilon_x - \|\vec{x}_i - \vec{x}_j\|) \cdot \Theta(\varepsilon_y - \|\vec{y}_i - \vec{y}_j\|), \qquad \vec{x}_i \in \mathbb{R}^m, \quad \vec{y}_i \in \mathbb{R}^n, \quad i,j=1,\ldots,N. $$ It represents the simultaneous recurrences in different systems. Such joint recurrence plots have the advantage, that the data can be different observables and can have different magnitudes. They can be used for the detection of phase synchronisation.

The purpose of CRPs an JRPs is different. To illustrate this let us consider the Roessler system (Fig. 5). We will compare the original system with transformed version of itself (Fig. 6), where we apply (1) a simple rotation of the phase space and (2) a transformation of the time axis (compressing the time).

Time series of the Roessler system.
Figure 5: Time series of the Roessler system: (a) original, (b) rotated phase space trajectory, and (c) compressed time. The rotation of the phase space trajectory is almost impossible to see (b). However, the compressed time axis is clearly visible in (c).

Whereas the rotated system can be clearly identified in the phase space, the trajectory of the time transformed system looks almost identical to the original system. However, the time series give a completely opposite picture.

Phase space representation of the realisations of the Roessler system.
Figure 6: Phase space representation of the realisations of the Roessler system as given in Fig. 5: (a) original, (b) rotated phase space, and (c) compressed time. Where the rotated version is obvious (b), the time axis transformation cannot be seen from the phase space (c).

The RPs of the original and rotated Roessler systems are identical (Fig. 7a and b), because the recurrence structure remains unaltered by the rotation. Therefore, the JRP is identical to these auto-RPs (Fig. 7d). However, because the phase space trajectory of the rotated system now visits different regions of the phase space, the CRP is more sparse (Fig. 7c).

Figure 7: (a) RP of the original Roessler system, (b) RP of the rotated system, (c) CRP between the original and rotated system, and (d) JRP of the original and rotated system.

In contrast, the RP of the time transformed Roessler system already reaveals the changing time scale (Fig. 8b). In the CRP, the line of synchronization has now a parabolic shape (Fig. 8c), reflecting the function that had been applied in order to get the time transformation (which was indeed the square function). However, the recurrence structure is now completely different, causing an almost empty JRP (Fig. 8d).

Figure 8: (a) RP of the original Roessler system, (b) RP of the time transformed system, (c) CRP between the original and transformed system, and (d) JRP of the original and transformed system.

We can conclude from this example that a CRP visualizes the simultaneous occurrence of similar states whereas a JRP visualizes the simultaneous occurrence of recurrences.


Further references about RPs, RQA and their applications.

Creative Commons License © 2000-2017 SOME RIGHTS RESERVED
The material of this web site is licensed under a Creative Commons Attribution-NonCommercial-NoDerivs 2.0 Germany License.

Please respect the copyrights! The content of this web site is protected by a Creative Commons License. You may use the text or figures, but you have to cite this source ( as well as N. Marwan, M. C. Romano, M. Thiel, J. Kurths: Recurrence Plots for the Analysis of Complex Systems, Physics Reports, 438(5-6), 237-329, 2007.

Spam Harvester Protection Network
provided by Unspam