Abstract
In this paper, we present some improvements to the convergent cross mapping (CCM) algorithm for detecting causality in unidirectionally connected chaotic systems. The basic concept of the CCM algorithm is that the causal influence of system $X$ on system $Y$ appears as mapping of the neighbouring states in the reconstructed $d$dimensional manifold, ${M}_{y}$, to the neighbouring states in the reconstructed $d$dimensional manifold, ${M}_{x}$, and this effect is evaluated using the correlation coefficient between the estimated and observed values of ${M}_{x}$. We proposed a composite indicator of causality as the ratio between the correlation coefficient and the Shannon entropy of the distribution of the residuals between the estimated and observed values of ${M}_{x}$. Application of the proposed approach to four masterslave Rössler and Lorenz systems and realworld data showed that the new algorithm allowed a slight increase in capability to reveal the presence and direction of couplings.
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 wellknown 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 [46]. 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 [913], to name a few. In 2012, Sugihara et al. [5] introduced another method based on state space reconstruction. The method was called convergent crossmapping (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 crossmapped and observed values of the reconstructed state space, and (2) Shannon entropy of the distribution of the residuals between the crossmapped 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}_{t2\tau},\cdots ,{x}_{t\left(E1\right)\tau}\right]}^{T}$ for $t=1+\left(E1\right)\tau $ to $t=L$, i.e., ${\mathbf{M}}_{X}={\left\{{\mathbf{x}}_{t}\right\}}_{t=1+\left(E1\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(E1\right)\tau ,\cdots ,L$ by:
where ${w}_{i}={u}_{i}\u2215\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)\u2215d\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 crossmap estimates is indicated by the correlation coefficient $\left(\rho \right)$ between ${\left\{{y}_{t}\right\}}_{t=1+\left(E1\right)\tau}^{L}$ and ${\left\{{\widehat{y}}_{t}{\mathbf{M}}_{X}\right\}}_{t=1+\left(E1\right)\tau}^{L}$ (or between ${\left\{{x}_{t}\right\}}_{t=1+\left(E1\right)\tau}^{L}$ and ${\left\{{\widehat{x}}_{t}{\mathbf{M}}_{Y}\right\}}_{t=1+\left(E1\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}$”. Crossmapping 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(E1\right)\tau}^{L}$ and ${\left\{{\widehat{x}}_{t}{\mathbf{M}}_{Y}\right\}}_{t=1+\left(E1\right)\tau}^{L}$ comparing to the value of the correlation coefficient between ${\left\{{y}_{t}\right\}}_{t=1+\left(E1\right)\tau}^{L}$ and ${\left\{{\widehat{y}}_{t}{\mathbf{M}}_{X}\right\}}_{t=1+\left(E1\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 unidirectionally coupled nonidentical 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(E1\right)\tau}^{L}$ and manifold ${\mathbf{M}}_{\widehat{X}}{\mathbf{M}}_{Y}={\left\{{\widehat{\mathbf{x}}}_{t}\right\}}_{t=1+\left(E1\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(E1\right)\tau}^{L}$ and manifold ${\mathbf{M}}_{\widehat{Y}}{\mathbf{M}}_{X}={\left\{{\widehat{\mathbf{y}}}_{t}\right\}}_{t=1+\left(E1\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 crossmapped estimate ${\mathbf{M}}_{\widehat{X}}{\mathbf{M}}_{Y}$ and between vectors of the manifold ${\mathbf{M}}_{Y}$ and its crossmapped estimate ${\mathbf{M}}_{\widehat{Y}}{\mathbf{M}}_{X}$, by determining the 2D correlation coefficient, and introduce interdependence measures ${CME}_{XY}$ and ${CME}_{YX}$ as:
where ${\rho}_{{M}_{X}}$ is the 2D correlation coefficient between vectors of the manifold ${\mathbf{M}}_{X}$ and its crossmapped 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 crossmapped 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 crossmapping evaluated by the proposed interdependence measures ($CME$), we compared it with the crossmapping evaluated using the correlation coefficient ($CM$) on four test examples of unidirectionally connected chaotic Rössler and Lorenztype 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 oneway 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 RungeKutta 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 timedelayed 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 $XY$ and the direction from $Y$ to $X$ as $YX$. Taking, for example, the measure CM: if $X$ drives $Y$, the measure $CM\left(XY\right)$ is expected to be higher than $CM\left(YX\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 $YX$ direction is little, but is considerable in the $XY$ direction, especially at moderate and high coupling strengths. As a result, the gap between $CME\left(XY\right)$ and $CME\left(YX\right)$ becomes larger, compared to the gap between $CM\left(XY\right)$ and $CM\left(YX\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 independentsamples ttest 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 unidirectionally 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 oneway 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 timedelayed 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 unidirectionally 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 realworld 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 Rpeaks 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 independentsamples ttest showed a statistically significant ($p=$ 0.0097) difference between $CM\left(XY\right)$ and $CM\left(YX\right)$ and much more statistically significant ($p=$ 0.000454) difference between$CME\left(XY\right)$ and $CME\left(YX\right)$.
5. Conclusions
A simulation involving examples of unidirectionally connected chaotic of Rössler and Lorenztype systems and realworld 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 unidirectionally connected chaotic systems.
References

Papana A, Kyrtsou C., Kugiumtzis D., Diks C. Simulation study of direct causality measures in multivariate time series. Entropy, Vol. 15, 2013, p. 26352661.

Granger C. W. J. Investigating causal relations by econometric models and crossspectral methods. Econometrica, Vol. 37, 1969, p. 424438.

Krakovska A., Hanzely F. Testing for causality in reconstructed state spaces by an optimized mixed prediction method. Physical Review E, Vol. 94, 2016, p. 052203.

Chen Y., Rangarajan G., Feng J., Ding M. Analyzing multiple nonlinear time series with extended Granger causality. Physics Letters A, Vol. 324, 2004, p. 2635.

Sugihara G., May R., Ye H., Hsieh C.H., Deyle E., Fogarty M., Munch S. Detecting causality in complex ecosystems. Science, Vol. 338, 2012, p. 496500.

Lusch B., Maia P. D., Kutz J. N. Inferring connectivity in networked dynamical systems: challenges using Granger causality. Physical. Review E, Vol. 94, 2016, p. 032220.

Schreiber T. Measuring information transfer. Physical Review Letters, Vol. 85, 2000, p. 461464.

Paluš M., Vejmelka M. Directionality of coupling from bivariate time series: how to avoid false causalities and missed connections. Physical Review E, Vol. 75, 2007, p. 056211.

Arnhold J., Grassberger P., Lehnertz K., Elger C. E. A robust method for detecting interdependences: application to intracranially recorded EEG. Physica D, Vol. 134, Issue 4, 1999, p. 419430.

Quiroga R. Q., Kraskov A., Kreuz T., Grassberger P. Performance of different synchronization measures in real data: a case study on electroencephalographic signals. Physical Review E, Vol. 65, 2002, p. 041903.

Andrzejak Kraskov R. G. A., Stögbauer H., Mormann F., Kreuz T. Bivariate surrogate techniques: Necessity, strengths, and caveats. Physical Review E, Vol. 68, 2003, p. 066202.

Smirnov D. A., Andrzejak R. G. Detection of weak directional coupling: phasedynamics approach versus statespace approach. Physical Review E, Vol. 71, 2005, p. 036207.

Chicharro D., Andrzejak R. G. Reliable detection of directional couplings using rank statistics. Physical Review E, Vol. 80, 2009, p. 026217.

Coufal D., Jakubík J., Jajcay N., Hlinka J., Krakovská A., Paluš M. Detection of coupling delay: a problem not yet solved. Chaos, Vol. 27, 2017, p. 083109.

Kantz H., Schreiber T. Nonlinear Time Series Analysis. Cambridge University Press, Cambridge, 2004.

Ye H., Deyle E. R., Gilarranz L. J., Sugihara G. Distinguishing timedelayed causal interactions using convergent cross mapping. Scientific Reports, Vol. 5, 2015, p. 14750.

Krakovská A., Jakubík J., Budáčová H., Holecyová M. Causality Studied in Reconstructed State Space. Examples of unidirectionally connected chaotic systems. 2016, https://arxiv.org/abs/1511.00505.

Jakubik J. Convergent Cross Mapping. Version 1.3, http://www.mathworks.com/matlabcentral/ fileexchange/52964convergentcrossmapping.