Abstract
In the present work, we present a new algorithm for assessing causality in unidirectionally coupled chaotic oscillators with small frequency mismatch embedded in heavy white Gaussian noise. This method is based on the correlation between changes in the phase dynamics of the slave oscillator and the dynamics of the phase difference between the oscillators. To recover the phase at low signaltonoise ratio, a nonlinear adaptive denoising algorithm based on finding sinusoidal fits to the local neighbourhood of the reconstructed phase space is used. Application of the proposed approach to masterslave Rössler systems showed that the new algorithm is wellsuited for assessing the presence and direction of coupling in highly noisy unidirectionally coupled chaotic oscillators, especially in the case of weak and moderate coupling.
Highlights
 The identification of directional couplings provides important information about the dynamics of the interacting nonlinear systems.
 Inferring causality from high noisy oscillators is highly nontrivial.
 The efficient approach for assessing causality from highly noisy unidirectionally coupled chaotic oscillators with small frequency mismatch is proposed.
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 causality test is based on autoregressive models and is particularly useful in stochastic linearly interconnected systems [3]. However, the 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 causality [4], transfer entropy [7], conditional mutual information [8], evolution map approach (EMA) [9, 10], convergent crossmapping (CCM) method [5, 11, 12], and measures evaluating distances of conditioned neighbours in reconstructed state spaces [1317], to name a few. A comprehensive review can be found in Focus Issue [18]. Specifically, CCM 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}$. 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. Most of the existing tools for detecting causality can make determinations of directionality in idealized conditions of sufficiently long, noisefree time series, but those determinations are relatively fragile in the presence of noise. However, the realworld datasets possess observational and process noise and may lack sufficient data. Taking the Rössler system as an example, it is shown in [19] that in the case of a small mismatch between the mutually incommensurate frequencies of the oscillators, changes in the phase dynamics of the slave oscillator are highly correlated with the dynamics of the phase difference between the oscillators. On the contrary, the phase dynamics of the master oscillator are not disturbed and do not correlate with the dynamics of the phase difference between the oscillators. This principle of asymmetries in correlations allows for efficiently assessing the presence and direction of coupling, especially in the case of weak coupling.
In this paper, we show that this method for estimating the causal effect in unidirectionally coupled chaotic oscillators can be successfully applied to a time series in highnoise environments due to its robustness against distortions caused by the noise reduction algorithms.
2. Description of the algorithm
We first give a brief introduction of the method [19] for assessing causality in unidirectionally coupled chaotic oscillators with small frequency mismatch. The approach is based on the correlation between changes in the phase dynamics of the slave oscillator and the dynamics of the phase difference between the oscillators. The instantaneous phase $\varphi \left(t\right)$ of a signal $s\left(t\right)$ is estimated using the Hilbert transform ${H}_{t}\left\{\bullet \right\}$:
where the analytic signal $z\left(t\right)$ can be understood as an embedding of the onedimensional timeseries in the twodimensional complex plane and the function ${H}_{t}\left\{s\left(t\right)\right\}$ is the Hilbert transform of $s\left(t\right)$:
and P.V. means that the integral is taken in the sense of the Cauchy principal value. In the case that $z\left(t\right)$ evolves around a central region, the cyclic phase is computed by the following expression:
and $\u2206\varphi \left(t\right)={\varphi}_{1}\left(t\right){\varphi}_{2}\left(t\right)$ is the instantaneous phase difference between the corresponding signals ${s}_{1}\left(t\right)$ and ${s}_{2}\left(t\right)$ at time $t$. Furthermore, the phase difference is determined by $\mathrm{a}\mathrm{b}\mathrm{s}\left(modulo\pi \right)$, to specify the closeness of the phases between two oscillators without indicating the direction of the difference. As a result, the instantaneous phase difference oscillates between 0 and $\pi $ at a low instantaneous frequency defined as the time rate of change of $\mathrm{\Delta}\varphi \left(t\right)$ and determined by the difference in natural frequencies of the signals and by the coupling strength (more details in [20, 21]). In order to obtain the continuous dynamics of the phases ${\varphi}_{1}\left(t\right)$ and ${\varphi}_{2}\left(t\right)$, the phases are unwrapped by tracing the ≈2$\pi $ jumps in the time course of ${\varphi}_{1}\left(t\right)$ and ${\varphi}_{2}\left(t\right)$:
where ${\widehat{\varphi}}_{i}={\varphi}_{i}+\pi $; $i=$ 1, 2; ${t}_{m}$ are time points at which ${\widehat{\varphi}}_{i}$ reaches its maximum, i.e., ${\widehat{\varphi}}_{i}\approx \text{2}\pi $, $k=$ 1, …, $N$. For practical use, there is a convenient Matlab function unwrap. The unwrapped phases are detrended to detect subtle changes in their dynamics. Furthermore, the highfrequency smallamplitude irregular fluctuations of the ${\widehat{\varphi}}_{i}\left(t\right)$ and $\mathrm{\Delta}\varphi \left(t\right)$ are neglected, and only the oscillations at a low averaged frequency (${\mathrm{\Omega}}_{D}=\u2329d\mathrm{\Delta}\varphi \left(t\right)/dt\u232a$ for the instantaneous phase difference and ${\mathrm{\Omega}}_{i}=\u2329d{\widehat{\varphi}}_{i}/dt\u232a$ for the detrended ${\widehat{\varphi}}_{i}\left(t\right)$) are considered. The method is grounded on the fact that the unwrapped and detrended phase $\widehat{\varphi}\left(t\right)$ of the slave oscillator exhibits significant lowfrequency oscillations synchronously with $\mathrm{\Delta}\varphi \left(t\right)$, i.e., is modulated according to the distance between phases and changes significantly with transition from inphase to antiphase state between coupled oscillators already at a weak coupling strength. On the contrary, no changes are observed in the unwrapped and detrended phase $\widehat{\varphi}\left(t\right)$ of the master oscillator synchronously with the phase difference. Therefore, the direction of coupling is inferred based on the asymmetries emerging from the correlations between $\mathrm{\Delta}\varphi \left(t\right)$ and $\widehat{\varphi}\left(t\right)$ for the slave oscillator and for the master oscillator. It is worth emphasizing that a high value of the correlation coefficient between $\mathrm{\Delta}\varphi \left(t\right)$ and $\widehat{\varphi}\left(t\right)$ for system $Y$ compared to the value of the correlation coefficient between $\mathrm{\Delta}\varphi \left(t\right)$ and $\widehat{\varphi}\left(t\right)$ for system $X$ indicates that system $X$ drives system $Y$. This runs in accordance with Granger’s method and counter to the CCM method. To recover the phase of the oscillating time series under high noise levels, we propose a nonlinear adaptive algorithm to recover continuoustime chaotic signals in heavy noise environments, based on a leastsquares fit of the sum of sines and cosines, each with a timevarying phase to the state vectors of the neighbourhood of the reconstructed phase space. 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 N$ contaminated by additive white Gaussian noise with zero mean and signaltonoise ratio (SNR) of 05 dB. As a first step, $D$dimensional attractor manifolds ${\mathbf{M}}_{X}$ and ${\mathbf{M}}_{Y}$ are reconstructed from lags of observable $X$ and $Y$ so that the state of the system in time $t$ is ${\mathbf{x}}_{t}={\left[{x}_{t},{x}_{t+\tau},{x}_{t+2\tau},\cdots ,{x}_{t+\left(D1\right)\tau}\right]}^{T}$ and ${\mathbf{y}}_{t}={\left[{y}_{t},{y}_{t+\tau},{y}_{t+2\tau},\cdots ,{y}_{t+\left(D1\right)\tau}\right]}^{T}$ respectively for $t=1$ to $t=L\left(D1\right)$, i.e., ${\mathbf{M}}_{X}={\left\{{\mathbf{x}}_{t}\right\}}_{t=1}^{L\left(D1\right)\tau}$ and ${\mathbf{M}}_{Y}={\left\{{\mathbf{y}}_{t}\right\}}_{t=1}^{L\left(D1\right)\tau}$, where $\tau $ is the time delay, and ${\left(\bullet \right)}^{T}$ denotes the transpose of a real matrix. One unit time delay and an embedding dimension approximately equal to one and a half of the fundamental period $\stackrel{~}{T}$ of the oscillations are used. The overembedding technique contributes to the determination of valid neighbours at strong noise and to obtain a tradeoff between noise suppression and a distortion of the local chaotic dynamics. Further, the vector time series ${\mathbf{x}}_{t}$ and ${\mathbf{y}}_{t}$ are subjected to an identical preprocessing procedure as shown for vector time series ${\mathbf{x}}_{t}$. For every $i$, $i=\text{1},W,\text{2}W,\cdots ,L\left(D1\right)\tau $, $W=D/\text{2}$, point ${\mathbf{x}}_{t}$ in the reconstructed $D$dimensional phase space, $r$, temporarily uncorrelated nearest neighbour points ${\left\{{\mathbf{x}}_{i}^{j}\right\}}_{j=1}^{r}$ are determined and mean value ${\widehat{x}}_{i}\left({t}_{k}\right)$, ${t}_{k}=\text{1},\cdots ,D$, for this cloud of nearest neighbour points is calculated. We write:
$+\cdots {\beta}_{12q}\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega {t}_{k}{\varphi}_{q}\left({t}_{k}\right)\right)+{\beta}_{22q}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega {t}_{k}{\varphi}_{q}\left({t}_{k}\right)\right)+\left.{\beta}_{0}\right),$
as an approximation to ${\widehat{x}}_{i}\left({t}_{k}\right)$, where $\omega =\text{2}\pi /\stackrel{~}{T}$ is the fundamental frequency of the sinusoids, ${\varphi}_{q}\left({t}_{k}\right)=b\left(q1\right){t}_{k}$, and ${\beta}_{11q}$, ${\beta}_{21q}$, ${\beta}_{12q}$, ${\beta}_{22q}$, and ${\beta}_{0}$ are parameters to be estimated. The parameters $Q$ and $b$ are chosen experimentally. The vector parameter $\mathbf{\beta}$ is estimated using the known leastsquares method:
where:
and ${\mathbf{M}}_{s+}$, ${\mathbf{M}}_{c+}$ are submatrices of size $D$×$Q$ encompassing $Q$ sines and $Q$ cosines respectively with a linear growth of the instantaneous frequency, ${\mathbf{M}}_{s}$, ${\mathbf{M}}_{c}$ are submatrices of size $D$×$Q$ encompassing $Q$ sines and $Q$ cosines respectively with a linear decrease of the instantaneous frequency, ${\mathbf{M}}_{I}$ is column vector $D$×1 of ones, and ${\left(\bullet \right)}^{\u2020}$ denotes the MoorePenrose pseudoinverse of a given matrix. Knowing parameter $\mathbf{\beta}$, we can rewrite Eq. (5) as:
$+\cdots {\mathbf{\beta}\left(2Q+1:3Q\right)}^{T}*{\mathbf{M}}_{s}^{T}+{\mathbf{\beta}\left(3Q+1:4Q\right)}^{T}*{\mathbf{M}}_{c}^{T}+\mathbf{\beta}\left(4Q+1\right).$
Estimated segments ${\stackrel{~}{x}}_{i}\left({t}_{k}\right)$ overlap by $W$ points. To ensure that the fitting is continuous everywhere and to eliminate the jumps or discontinuities around the boundaries of neighbouring segments, we used the fitting for the overlapped region as described in [22]. Denote the fitted curve for the $i$th and ($i+1$)th segments by ${\stackrel{~}{x}}_{i}\left({t}_{k}\right)$ and ${\stackrel{~}{x}}_{i+1}\left({t}_{k}\right)$, respectively. Then the fitting for the overlapped region can be defined as:
where ${w}_{1}=\text{1}\left({t}_{l}\text{1}\right)/W$, and ${w}_{2}=\left({t}_{l}\text{1}\right)/W$. As follows from Eq. (9), the weights decrease linearly with the distance between the point and the center of the segment. Such a weighting ensures symmetry and effectively eliminates any jumps or discontinuities around the boundaries of neighbouring segments. Numerical experiments on the Rössler system show that the proposed algorithm provides about 13 dB, on average, improvement in SNRs at initial SNR = 0 dB. The efficacy of this preprocessing algorithm enables to properly recover the phase of the noisy chaotic oscillator and, in consequence ensures the effectiveness of the proposed algorithm for inferring causality from highly noisy unidirectionally coupled chaotic oscillators.
3. Simulation results
To demonstrate the performance of the method we considered two examples of unidirectionally connected chaotic Rössler systems [12] with small frequency mismatch and having variable coupling strengths. The oscillators were coupled via a oneway driving relationship between variable ${x}_{1}$ of the driving system and variable ${y}_{1}$ of the responsive system:
where ${\omega}_{1}\left(\text{1}\right)=$ 1.015, ${\omega}_{2}\left(\text{1}\right)=$ 0.985, ${\omega}_{1}\left(\text{2}\right)=$ 0.980, and ${\omega}_{2}\left(\text{2}\right)=$ 1.020. The coupling strength c was chosen from 0 to 0.14 with the step 0.01 for the first case and from 0 to 0.045 with the step 0.005 for the second case. 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. Approximation parameters $Q$ and $b$ were chosen as $Q=$ 10, $b=$ 0.0012, i.e., the total number of sinusoids and cosines were 40, and the maximum deviation of the fundamental frequency amounted to approximately 11 %. The performance of the algorithm was evaluated by means of a 30trial Monte Carlo simulation whereby the signals are contaminated by various realizations of additive zeromean white Gaussian noise at a SNR = 0 dB. All data processing and analyses were performed using Matlab software (MathWorks, Natick, MA). Consequently, we denoted the direction from system $X$ to system $Y$ as$YX$ and the direction from $Y$ to $X$ as $XY$. Taking, for example, the measure $\rho $ (correlation coefficient): if $X$ drives $Y$, the measure $\rho \left(YX\right)$ is expected to be higher than $\rho \left(XY\right)$. As we can see from Fig. 1 and Fig. 2, the measure $\rho $ shows that system $X$ drives system $Y$ until the onset of synchronization, i.e., when their phases are locked, and when this method becomes meaningless. The mean gap between the measures $\rho \left(YX\right)$ and $\rho \left(XY\right)$ for noisefree signals and for signals contaminated with white Gaussian noise becomes roughly the same. For comparison two stateofart algorithms were chosen: CCM and EMA. In the latter case, the directionality index $d$ was calculated using the Matlab code freely available on Ref. [23]. As we can see from Fig. 3, both methods fails to detect the causality at the conditions under study. CCM method shows false, i.e., weak opposite or bidirectional coupling direction, and the directionality index $d$ randomly oscillates between coupling directions. It is worth noting that the sign of correlation coefficient $\rho $ between $\mathrm{\Delta}\varphi \left(t\right)$ and $\widehat{\varphi}\left(t\right)$ for the responsive system can be positive or negative. In our example, $\rho $ is positive when the frequency of the driving system is higher than the frequency of the responsive system and negative when the frequency of the driving system is lower than the frequency of the responsive system.
Fig. 1Correlation coefficient ρ between Δϕt and ϕ^t computed for two unidirectionally coupled Rössler oscillators with natural frequencies ω1= 1.015, ω2= 0.985: a) for noisefree signals, b) the mean values and standard deviations (with error bars) from 30 independent Monte Carlo trials consisting of the noiseadded signal at a SNR = 0 dB
a)
b)
Fig. 2Correlation coefficient ρ between Δϕt and ϕ^t computed for two unidirectionally coupled Rössler oscillators with natural frequencies ω1= 0.980, ω2= 1.020: a) for noisefree signals, b) the mean values and standard deviations (with error bars) from 30 independent Monte Carlo trials consisting of the noiseadded signal at a SNR = 0 dB
a)
b)
Fig. 3Measures ρ obtained by CCM method and directionality index d obtained by EMA method for two unidirectionally coupled Rössler oscillators with natural frequencies ω1= 1.015, ω2= 0.985 and contaminated with noise at SNR = 0 dB
4. Conclusions
A simulation involving examples of unidirectionally connected chaotic Rösslertype systems shows that the algorithm for assessing causality in unidirectionally coupled chaotic oscillators based on the correlation between changes in the phase dynamics of the slave oscillator and the dynamics of the phase difference between the oscillators in combination with the proposed nonlinear adaptive denoising algorithm is an effective method to reveal the presence and the direction of coupling in unidirectionally coupled chaotic oscillators with a small frequency mismatch and embedded in heavy white Gaussian noise. The success of this approach is due to the robustness of the phasebased algorithm for assessing causality against unavoidable distortions of the local chaotic dynamics during the denoising process, and the comparison with the stateoftheart methods shows that the proposed method is superior to the other approaches. Although in this work we explored a classical example of coupled chaotic Rössler oscillators frequently used for studying phase synchronization and causation, we believe that similar results can be obtained in other systems and future work should confirm this.
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.

Rosenblum M. G., Pikovsky A. S. Detecting direction of coupling in interacting oscillators. Physical Review E, Vol. 64, 2001, p. 045202.

Smirnov D. A., Bezruchko B. P. Estimation of interaction strength and direction from short and noisy time series. Physical Review E, Vol. 68, 2003, p. 046209.

Coufal D., Jakubik J., Jajcay N., Hlinka J., Krakovska A., Palus M. Detection of coupling delay: a problem not yet solved. Chaos, Vol. 27, 2017, p. 083109.

Krakovská A., Jakubík J., Chvosteková M. Comparison of six methods for the detection of causality in a bivariate time series. Physical Review E, Vol. 97, 2018, p. 042207.

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 R. G. Kraskov 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.

Bollt E. M., Sun J., Runge J. Introduction to focus issue: causation inference and information flow in dynamical systems: theory and applications. Chaos, Vol. 28, 2018, p. 075201.

Pukenas K. Detecting causality in unidirectionally coupled chaotic oscillators with small frequency mismatch, accepted for publication. Journal of Applied Nonlinear Dynamics, 2020, (in press).

Rosenblum M. G., Pikovsky A. S., Kurths J. Phase Synchronization of Chaotic Oscillators. Physical Review Letters, Vol. 76, Issue 11, 1996, p. 18041807.

Boccaletti S., Kurths J., Osipov G., Valladares D. L., Zhou C. S. The synchronization of chaotic systems. Physics Reports, Vol. 366, Issues 12, 2002, p. 1101.

Tung W. W., Gao J., Hu J., Yang L. Detecting chaos in heavynoise environments. Physical Review E, Vol. 83, 2011, p. 046210.

Kralemann B., Rosenblum M., Pikovsky A. DAMOCO: Data analysis with models of coupled oscillators Matlab Toolbox for multivariate time series analysis. Version 2.0, 2014, http://www.stat.physik.unipotsdam.de/~mros/damoco2.html.