Reliable detection of causal asymmetries in dynamical systems

Knowledge about existence, strength, and dominant direction of causal influences is of paramount importance for understanding complex systems. Current methods deduce ambiguous causal links among different observables from (complex) dynamical systems, if a limited amount of realistic data is used. It is particularly difficult to infer the dominant direction of causal influence for synchronizing systems. Missing is a statistically well defined approach that avoids false positive detection while being sensitive for weak interactions. The proposed method exploits the local inflation of manifolds to estimate upper bounds on the information loss among state reconstructions and tests for the absence of causal influences. Simulated data demonstrates that it is robust to intrinsic noise, copes with synchronization, and tolerates moderate amounts of measurement noise.

Figure 1

$x_t$ over $y_t$ and $y_{t+1}$ for the noise free unilaterally coupled logistic maps with $w_{y\to x}=0,w_{x\to y}=0.3,\text{and}{R}_x=R_y=3.82$, together with the projection of the manifold on the $y_{t+1}-y_t$ plane and the $x_t$ axis. ${10}^3$ data points are shown in blue and the 10 nearest neighbors of a reference point are shown in color in both the manifold and the projections (gray). The 10 neighbors searched in $Y$ ($X$) are shown in red (green).

Figure 2

Logarithmic neighborhood sizes $d_i^j\left(k\right)$ over $\kappa$ for noise free unilateral coupled logistic maps ($w_{x\to y}=0,w_{y\to x}=0.2$, and $R_x=R_y=3.82$) from $N={10}^3$ data points. The observables were embedded with $m=4$ and $\tau =1$ and ${10}^2$ ensembles were used for chance-level estimation. (a) $d_x^x\left(k\right)$ shown as plus signs, $d_x^y\left(k\right)$ shown as dotted line, and the respective chance-level $d_x^{y^{*}}\left(k\right)$ (dashed line). (b) $d_y^y\left(k\right)$ shown as plus signs, $d_y^x\left(k\right)$ shown as dotted line, and the respective chance-level $d_y^{x^{*}}\left(k\right)$ (dashed line). Furthermore, the fivefold standard error is shown for each $d_i^j\left(k\right)$ as a gray shade

Figure 3

Causal influence between bilaterally coupled logistic maps. $2\times {10}^4$ data points were generated with $R_x=R_y=3.92$, the observed time series were embedded with $m=4$ and $\tau =1$, and we used the $k=20$ nearest neighbors to determine the causal influence and $E={10}^3$ permutations for the chance level. Black lines show $I_{y\to x}$ corresponding to the varied coupling $w_{y\to x}$; gray lines show the reverse direction with the coupling $w_{x\to y}$ fixed at 0 (circles), 0.00316 (squares), and 0.032 (triangles).

Figure 4

Bilaterally coupled logistic maps $R_x=R_y=3.92$ embedded with $\tau =1$ and $m=4$; ${10}^3$ ensembles were generated to estimate chance level. The upper row, (a)–(c), shows the causal influence $I_{x\to y}$ for $w_{y\to x}=0.05$ and varying $w_{x\to y}$ between 0 and 0.1. The bottom row, (d)–(f), shows the same causal influence, but for a fixed $w_{x\to y}=0.05$ and varying $w_{y\to x}$ between 0 and 0.1. (a),(d) Varying amount of data between ${10}^3$ and ${10}^4$. (b),(e) Causal influence for additive internal noise. The amplitude of the noise is varied between 0 and $8\%$ of the standard deviation of the unperturbed system. (c),(f) Additive external noise is added to the observed time series. The standard deviations of noise are varied between 0 and $8\%$ of the amplitude of the unpolluted system. For all noise polluted results time series of ${10}^4$ data points were used and for better visualization contour lines mark lines of equal causal influence.

Figure 5

Asymmetry index between two Lorenz oscillators ($N={10}^4$ data points) with slightly different frequencies and coupled by their $y$ components (${\theta }_{1/2}=10,{\rho }_1=28.5,{\rho }_2=27.5,\text{and}{\theta }_{1/2}=8/3$). The time series were embedded with $m=9$ and $\tau =10$ and $E={10}^2$ were used to estimate chance level. The noise free asymmetry is shown in black. Colored circles represent internal noise with $\sigma =2$ (green) and additive external noise with $\sigma =2$ (red).

Figure 6

(a) Causal influence for the (unilaterally) coupled Lorenz $\to$ Rössler system derived from ${10}^3$ reference points chosen from $3\times {10}^4$ data points. The time series were embedded with $\tau =2$ and $m=7$ and $E={10}^3$ ensembles were used to derive a chance level. Solid lines show the noise free results and dotted lines the noise perturbed results. The causal influence in direction of the coupling $I_{x\to y}$ is shown in blue and the reverse direction in red. (b) The asymmetry $\alpha$ for the noise free (solid lines) and perturbed system (asterisks). Gray lines show the information dimension estimated from the respective time series, with solid lines being the Lorenz system and dashed line the Rössler system. (c) Causal influence (M method) $I_{x\to y}$ (blue) and $I_{y\to x}$ (red) for the noise free (solid lines) and for the systems perturbed by Gaussian white noise (asterisks). (d) Causal influence (convergent cross mapping) $I_{x\to y}$ (blue) and $I_{y\to x}$ (red) for the noise free (solid lines) and for the systems perturbed by Gaussian white noise (asterisks).

Figure 7

(a) $\mathrm{\Delta }I$ derived for a system of two coupled Lorenz oscillators [${\theta }_{1/2}=10,{\rho }_1=39,{\rho }_2=35,{\theta }_{1/2}=8/3, {\mathrm{\Omega }}_{1\to 2}=w_{1\to 2}(x_1-x_2)$, and ${\mathrm{\Omega }}_{2\to 1}=0$]. The unilateral coupling weight $w_{x\to y}$ is increasing from 0 to 17. A fourth-order Runge-Kutta algorithm with step size 0.005 and down sampled in intervals of 0.03 time units integrated the system. After removing the transient $N=2^{11}$ data points are used for further analysis and are injected with additive noise either in $x_1, x_2$ or kept free of noise. We use Gaussian white noise with amplitude of $\sigma =0.95{\sigma }_{x_i}$ and all time series are embedded with $m=8$ and $\tau =4$. The solid black line shows the noise free $\mathrm{\Delta }I$; dashed blue lines show perturbations in $x_1$ and dashed red lines in $x_2$. All results are averaged over 100 ensembles for each coupling weight and noise condition. Panel (b) shows corresponding results using a moving average filter with a window of eight time steps before further analysis. Dots mark $\mathrm{\Delta }I$ with a significant difference from 0 derived by a Wilcoxon signed rank test ($p=0.001$).

Figure 8

(a) Average rank of neighbors over the number of neighbors $k$ for different noise levels shown in shades of gray. The dashed line marks the base line of the neighbor rank of an unperturbed time series. (b) Average rank of the first $k=20$ neighbors for four different Minkowski distances (0.1 to 2). To improve the neighborhood relations, moving average (blue), shared nearest neighbors (red), and no method (black). For each curve or bar a time series of $N=2^{11}$ data points was generated by Eq. (1). The moving average filter with a window size of eight time steps was used before searching neighbors. (Blue bars) The shared nearest neighbors were chosen as a subset of the first 200 neighbors that also shared 40 neighbors with a reference point (red bars). The black bars show the rank for a perturbed time series

Figure 9

(a) $\mathrm{\Delta }TC$, (b) $\mathrm{\Delta }M$, and (c) $\mathrm{\Delta }CCM$ derived for a system of two coupled Lorenz oscillators [${\theta }_{1/2}=10,{\rho }_1=39,{\rho }_2=35,{\theta }_{1/2}=8/3, {\mathrm{\Omega }}_{1\to 2}=w_{1\to 2}(x_1-x_2)$, and ${\mathrm{\Omega }}_{2\to 1}=0$]. The unilateral coupling weight $w_{x\to y}$ is increasing from 0 to 17. A fourth-order Runge-Kutta algorithm with step size 0.005 and down sampled in intervals of 0.03 time units integrated the system. After removing the transient $N=2^{11}$ data points are used for further analysis and are injected with either noise in $x_1, x_2$ or kept free of noise. The injected noise is Gaussian white noise with amplitude of $\sigma =0.95{\sigma }_{x_i}$ and all time series are embedded with $m=8$ and $\tau =4$. The solid black line shows the noise free $\mathrm{\Delta }I$; dashed blue lines show perturbations in $x_1$ and dashed red lines in $x_2$. All results are averaged over 100 ensembles for each coupling weight and noise condition. Panels (d), (e), and (f) show corresponding results using a moving average filter with a window of eight time steps before further analysis. Dots mark a significant positive or negative $\mathrm{\Delta }X$ derived by a Wilcoxon signed rank test ($p=0.001$).

Figure 10

(a) Coupling matrix $W$ of a 10-species Ricker model. A darker color signifies a larger entry in $W$. (b) Interaction matrix $W$ derived using CPM; circles represent true positive entries and squares false positives. The interaction matrix was derived from time series of length 300 and embedded using a delay of $\tau =1$ and embedding dimension of $m=8$ using the first $k=20$ neighbors. The simulation was repeated 25 times; the radius of the circle is proportional to the proportion of detection. The color of the circle is the average value of CPM over all 25 simulations.

Figure 11

Heart rate (a) and breath rate (b) of a sleeping apnea patient. Panels (c), (d) show logarithmic neighborhood sizes for the heart and breathing rate. The data was sampled at $2Hz$ and embedded with $\tau =6$ and $m=5$; all $n=601$ data points were chosen as reference points. (c) $d_H^H\left(k\right)$ shown as solid line, $d_H^B\left(k\right)$ shown as dotted line, and the respective chance level $d_H^{B^{*}}\left(k\right)$ (dashed line). Due to the small amount of data only the threefold standard deviation was used to determine the chance level. (d) $d_B^B\left(k\right)$ shown as solid line, $d_B^H\left(k\right)$ shown as dotted line, and the respective chance level $d_B^{H^{*}}\left(k\right)$ (dashed line).