An algorithm based on the convergent cross mapping method for the detection of causality in uni-directionally connected chaotic systems

1. Introduction

Quantification of the causal effects between simultaneously observed systems from the analysis of time series recordings is essential in many scientific fields, including economics, climatology, ecosystems, electrical activity of the brain, or cardiorespiratory relations. Estimating the interdependence between the observed variables provides valuable knowledge regarding the processes that generate time series [1]. The Granger causality method [2] is the most well-known and principal method for identifying directional interactions between variables from their time series, and many modifications and extensions of the Granger causality test has been developed. The Granger test focuses on determining whether one time series is useful in forecasting another. If useful, the first system causally affects the second system. The conventional Granger’s causality test is based on autoregressive models and is particularly useful in stochastic linearly interconnected systems [3]. However, its appropriateness of direct application to nonlinear systems depends on the specific problem to be analyzed [4-6]. Therefore, new approaches have been proposed, including nonlinear extensions of Granger’s causality [4], transfer entropy [7], conditional mutual information [8] and measures evaluating distances of conditioned neighbors in reconstructed state spaces [9-13], to name a few. In 2012, Sugihara et al. [5] introduced another method based on state space reconstruction. The method was called convergent cross-mapping (CCM) and involves evaluating distances between conditioned neighbours in reconstructed state spaces and tests for causation between the driver (“master”) system $X$ and the driven (“slave”) system Y, by measuring the extent to which the historical states of reconstructed state space, $M_Y$, can reliably estimate states of reconstructed state space, $M_X$. Similarity between points, $X_n$, and estimates, ${\widehat{X}}_n$, is evaluated by assessing their correlation. To infer the direction of coupling, the opposite direction must also be assessed. Then, the direction of coupling is inferred based on asymmetries emerging from the calculations of the two possible causal directions.

In this paper, we propose a simple improvement to the CCM method by introducing a composite parameter for inferring the direction of causation composed of two components: (1) the correlation coefficient between the cross-mapped and observed values of the reconstructed state space, and (2) Shannon entropy of the distribution of the residuals between the cross-mapped and observed values of the reconstructed state space.

2. Convergent cross mapping

CCM [5, 14] is a method for detecting causality between two systems represented by time series. CCM is based on Takens’ embedding theorem, exploiting the geometry of attractors of coupled dynamical systems and constructs a map between mutual neighborhoods in state spaces of the coupled dynamical systems under study.

Let $X={\left\{x_t\right\}}_{t=1}^L$, and $Y={\left\{y_t\right\}}_{t=1}^L$ be two time series of finite length $L\in \mathbb{N}$. To establish cross mapping from $X$ to $Y$, first we reconstruct the attractor manifold ${\mathbf{M}}_X$ as a set of $E$-dimensional vectors [15] ${\mathbf{x}}_t={\left[x_t,x_{t-\tau },x_{t-2\tau },\cdots ,x_{t-\left(E-1\right)\tau }\right]}^T$ for $t=1+\left(E-1\right)\tau$ to $t=L$, i.e., ${\mathbf{M}}_X={\left\{{\mathbf{x}}_t\right\}}_{t=1+\left(E-1\right)\tau }^L$, where $E$ is the embedding dimension, $\tau$ is the time delay, and ${\left(\bullet \right)}^T$ denotes the transpose of a real matrix. We find the $E+1$ nearest neighbors of ${\mathbf{x}}_t$ in ${\mathbf{M}}_X$ and denote their time indices (from the closest to the farthest) by $t_1,\cdots ,t_{E+1}$. These indices will be used in the construction of the cross mapping of ${\widehat{y}}_t$ for $y_t$, $t=1+\left(E-1\right)\tau ,\cdots ,L$ by:

where $w_i=u_i∕\sum _{j=1}^{E+1}{u}_j$, $u_i=\mathrm{e}\mathrm{x}\mathrm{p}\left[-d\left({\mathbf{x}}_t,{\mathbf{x}}_{t_i}\right)∕d\left({\mathbf{x}}_t,{\mathbf{x}}_{t_1}\right)\right]$ and $d\left(A,B\right)$ is the Euclidean distance between vectors $A$ and $B$. The cross mapping from Y to X is defined analogously. The skill of cross-map estimates is indicated by the correlation coefficient $\left(\rho \right)$ between ${\left\{y_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and ${\left\{{\widehat{y}}_t|{\mathbf{M}}_X\right\}}_{t=1+\left(E-1\right)\tau }^L$ (or between ${\left\{x_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and ${\left\{{\widehat{x}}_t|{\mathbf{M}}_Y\right\}}_{t=1+\left(E-1\right)\tau }^L$, respectively). Sugihara et al. [5] declare in the supplement: “If $X$ and $Y$ are dynamically coupled, the nearest neighbors of ${\mathbf{M}}_X$ should identify the time indices of corresponding nearest neighbors on ${\mathbf{M}}_Y$. As $L$ increases, the attractor manifold fills in and the distances among the $E+1$ nearest neighbors shrinks. Consequently, ${\widehat{y}}_t|{\mathbf{M}}_X$ should converge to $y_t$ and ${\widehat{x}}_t|{\mathbf{M}}_Y$ should converge to $x_t$. In this way, we use convergence of the nearest neighbors to test whether there is a correspondence between states on ${\mathbf{M}}_X$ and states on ${\mathbf{M}}_Y$”. Cross-mapping that converges in only one direction is the criterion for unidirectional causality. One of the fundamental concepts of CCM is that when causation is unilateral, $X\to Y$ ($X$ drives $Y$), then it is possible to estimate $X$ from $Y$, but not $Y$ from $X$. It follows that a high value of the correlation coefficient between ${\left\{x_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and ${\left\{{\widehat{x}}_t|{\mathbf{M}}_Y\right\}}_{t=1+\left(E-1\right)\tau }^L$ comparing to the value of the correlation coefficient between ${\left\{y_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and ${\left\{{\widehat{y}}_t|{\mathbf{M}}_X\right\}}_{t=1+\left(E-1\right)\tau }^L$ indicates that system $X$ drives system $Y$. This runs counter to Granger’s intuitive scheme. For more details on the CCM method see [5, 14, 16, 17].

3. Description of the algorithm

The simulation results obtained from the CCM method applied to uni-directionally coupled non-identical oscillators ($X$ drives $Y$) shows that probability distribution of the residuals taken as the Euclidean distance $R\left({\mathbf{M}}_X,{\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y\right)$ between manifold ${\mathbf{M}}_X={\left\{{\mathbf{x}}_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and manifold ${\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y={\left\{{\widehat{\mathbf{x}}}_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ obtained by cross mapping using a manifold ${\mathbf{M}}_Y$, and the Euclidean distance $R\left({\mathbf{M}}_Y,{\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X\right)$ between manifold ${\mathbf{M}}_Y={\left\{{\mathbf{y}}_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ and manifold ${\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X={\left\{{\widehat{\mathbf{y}}}_t\right\}}_{t=1+\left(E-1\right)\tau }^L$ obtained by cross mapping using a manifold ${\mathbf{M}}_X$, differs. More specifically, the probability distribution of $R\left({\mathbf{M}}_X,{\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y\right)$ has a taller, narrower shape (Fig. 1). It is therefore reasonable to combine the Shannon entropy of the probability distributions of $R\left({\mathbf{M}}_X,{\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y\right)$ or $R\left({\mathbf{M}}_Y,{\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X\right)$, and the corresponding correlation coefficient to enhance the inference of the direction of instantaneous causality. Furthermore, we compute the relationship between vectors of the manifold ${\mathbf{M}}_X$ and its cross-mapped estimate ${\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y$ and between vectors of the manifold ${\mathbf{M}}_Y$ and its cross-mapped estimate ${\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X$, by determining the 2D correlation coefficient, and introduce interdependence measures ${CME}_{X|Y}$ and ${CME}_{Y|X}$ as:

where ${\rho }_{M_X}$ is the 2D correlation coefficient between vectors of the manifold ${\mathbf{M}}_X$ and its cross-mapped estimate ${\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y$, ${\rho }_{M_Y}$ is the 2D correlation coefficient between vectors of the manifold ${\mathbf{M}}_Y$ and its cross-mapped estimate ${\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X$, ${En}_{R_X}$ is the Shannon entropy of the probability distribution of $R\left({\mathbf{M}}_X,{\mathbf{M}}_{\widehat{X}}|{\mathbf{M}}_Y\right)$, and ${En}_{R_Y}$ is the Shannon entropy of the probability distribution of $R\left({\mathbf{M}}_Y,{\mathbf{M}}_{\widehat{Y}}|{\mathbf{M}}_X\right)$.

Fig. 1Example of the probability distribution of the Euclidean distances RMX,MX^|MY and RMY,MY^|MX

4. Results

To demonstrate the effectiveness of the cross-mapping evaluated by the proposed interdependence measures ($CME$), we compared it with the cross-mapping evaluated using the correlation coefficient ($CM$) on four test examples of uni-directionally connected chaotic Rössler and Lorenz-type systems that were coupled with variable coupling strengths. In the examples used in this paper, the possibility of correlation instead of causation was excluded, and therefore the aspect of convergence was not used. The first and second datasets originated from the coupling of two Rössler systems [17], which were coupled via a one-way driving relationship between variables $x_1$ of the driving system and variable $y_1$ of the responsive system:

where ${\omega }_1\left(1\right)=$ 1.015, ${\omega }_2\left(1\right)=$ 0.985 and ${\omega }_1\left(2\right)=$ 1.075, ${\omega }_2\left(2\right)=$ 1.0. The coupling strength $c$ was chosen from 0 to 0:25 with the step 0:01. The data were generated by Runge-Kutta integration with a step size of 0.1. The first 2000 data points were discarded. The total number of obtained data was 14000 and this resulted in around 60 samples per one average orbit around the attractor. The causal relationship between the two systems was calculated using 7000 time-delayed vectors of $x_1$ and $y_1$ with a time delay equal to 3 and embedding dimension of 7. Eight nearest neighbours were used. All data processing and analyses were performed using Matlab software (MathWorks, Natick, MA) by adopting the ideas from reference [18], among others. Consequently, we denoted the direction from $X$ to $Y$ as $X|Y$ and the direction from $Y$ to $X$ as $Y|X$. Taking, for example, the measure CM: if $X$ drives $Y$, the measure $CM\left(X|Y\right)$ is expected to be higher than $CM\left(Y|X\right)$. As we can see from Fig. 2, both CME and CM measures show that $X$ drives $Y$ until the onset of synchronization. The effect of the entropy contribution in the $Y|X$ direction is little, but is considerable in the $X|Y$ direction, especially at moderate and high coupling strengths. As a result, the gap between $CME\left(X|Y\right)$ and $CME\left(Y|X\right)$ becomes larger, compared to the gap between $CM\left(X|Y\right)$ and $CM\left(Y|X\right)$, i.e. the direction of causation can be determined more reliable. Furthermore, in the case of the coupled Rössler systems with ${\omega }_1=$ 1.075, ${\omega }_2=$ 1.0, the independent-samples t-test showed that this difference between gaps across all coupling strengths is near to statistical significance ($p=$ 0.0866).

Fig. 2Measures CME and CM computed for two uni-directionally coupled Rössler systems: a) ω1= 1.015, ω2=0 .985, b) ω1= 1.075, ω2= 1.0

a)

b)

The third and fourth datasets originated from the coupling of two Lorenz systems [17], which are coupled via one-way driving relationship between variables $x_1$ of the driving system and variable $y_1$ of the responsive system:

where ${\omega }_1\left(1\right)=$ 28,5, ${\omega }_2\left(1\right)=$ 27.5 and ${\omega }_1\left(2\right)=$ 39, ${\omega }_2\left(2\right)=$ 35. The coupling strength $c$ was chosen from 0 to 11.6 with a step size 0:4. The data were generated using the Matlab solver of ordinary differential equations ode45. The first 2000 data points were discarded. The total number of obtained data was 14000 and the causal relationship between the two systems was calculated using 7000 time-delayed vectors of $x_1$ and $y_1$, with a time delay equal to 3 and embedding dimension of 7. As we can see from Fig. 3, the results are similar to the coupled Rössler systems, although the effect of entropy is greater.

Fig. 3Measures CME and CM computed for two uni-directionally coupled Lorenz systems: a) ω1= 28.5, ω2= 27.5, b) ω1= 39, ω2= 35

a)

b)

To demonstrate the performance of the proposed algorithm on real-world data, we applied the algorithm for the quantification of coupling between respiration and heart rate variability (HRV) because is well known that HRV is strongly related to respiration. For this purpose, we used data from the Fantasia database (https://physionet.org/physiobank/database/fantasia/). The data that are free of artifacts were acquired from 12 young subjects. The R-peaks of the electrocardiogram (ECG) signal were detected by means of a threshold method and the RR time series were obtained. In this study, we removed the mean of every RR time series either as of every respiration signal so that they oscillate around zero. The independent-samples t-test showed a statistically significant ($p=$ 0.0097) difference between $CM\left(X|Y\right)$ and $CM\left(Y|X\right)$ and much more statistically significant ($p=$ 0.000454) difference between$CME\left(X|Y\right)$ and $CME\left(Y|X\right)$.

5. Conclusions

A simulation involving examples of uni-directionally connected chaotic of Rössler and Lorenz-type systems and real-world time series shows that the proposed substitutions in the CCM method reinforce the asymmetries between the calculations of the two possible causal directions and thus improve determination of the presence and direction of coupling, and also detection of the onset of full synchronization in uni-directionally connected chaotic systems.