A simple test for causality in complex systemsPreprint: APS/123-QED
Abstract
We provide a new solution to the long-standing problem of inferring causality from observations without modeling the unknown mechanisms. We show that the evolution of any dynamical system is related to a predictive asymmetry that quantifies causal connections from limited observations. A built-in significance criterion obviates surrogate testing and drastically improves computational efficiency. We validate our test on numerous synthetic systems exhibiting behavior commonly occurring in nature, from linear and nonlinear stochastic processes to systems exhibiting non-linear deterministic chaos, and on real-world data with known ground truths. Applied to the controversial problem of glacial-interglacial sea level and CO2 evolving in lock-step, our test uncovers empirical evidence for CO2 as a driver of sea level over the last 800 thousand years. Our findings are relevant to any discipline where time series are used to study natural systems.
Keywords:
Predictive asymmetry, causality, time seriesI Introduction
Natural systems, such as ecosystems, Earth systems, or the human brain, pose formidable challenges to causal inference. Complex underlying dynamics may render traditional statistical means of correlation powerless to resolve causal relationships in real-world data. If process modeling is impractical or if model parameters cannot be constrained, then the prospect of non-parametric detection of causality directly from observations becomes tantalizing. Hence, several methods aimed at quantifying causal interactions from observed time series have been proposed Granger 1969; Chen et al. 2004; Marinazzo et al. 2008; Schreiber 2000; Paluš and Vejmelka 2007; Rulkov et al. 1995; Schiff et al. 1996; Arnhold et al. 1999; Quiroga et al. 2002; Chicharro and Andrzejak 2009; Sugihara et al. 2012; Wiesenfeldt et al. 2001; Feldmann and Bhattacharya 2004; Krakovská and Hanzely 2016; Liang 2013; McCracken and Weigel 2016.
Widely used methods include Granger causality Granger 1969, its information-theoretic cousin, transfer entropy Schreiber 2000, and geometrically inspired methods such as convergent cross mapping Sugihara et al. 2012; Ye et al. 2015. Despite their promise, attempts to quantify the strength and directionality of causal interactions from observed time series, without recourse to modeling, remain controversial. For example, the applicability of these causality detection methods is system-dependent, can suffer from biases and low statistical power Krakovská et al. 2018; Smirnov 2013, and usually requires additional statistical testing against surrogate data Theiler et al. 1992; Lancaster et al. 2018, which inevitably introduces subjectivity when it comes to surrogate data design. Moreover, numerical estimates of information-theoretic quantities such as transfer entropy may not converge to zero for uncoupled systems, may overestimate or fail to quantify information flow, or underestimate dynamical influence Smirnov 2013; James et al. 2016. In appendix A, we demonstrate some of these issues using transfer entropy (TE) as an exemplar.
Here, we present the predictive asymmetry — a simple and robust causality test that by construction overcomes many of these issues. Our test is based on a difference between two information-theoretic functionals computed directly from observed time series. We prove that this simple difference is directly related to the flow of the underlying dynamics. By quantifying the difference between forwards-in-time and backwards-in-time prediction, our predictive asymmetry test unequivocally determine the correct causation in cases where TE alone fails (Figs. S.A1, S.A2, S.A3, S.A4). Our test provides a statistic that is zero when there is no coupling, positive in the causal direction (driver response) when directional coupling exists, and negative in the non-causal direction (response driver) when directional coupling exists. Simultaneously, by using an intrinsic, dynamically informed significance test, our method alleviates computational demands associated with surrogate testing, raising the prospect of fast quantification of causal networks from large datasets.
In the following, we formally derive the predictive asymmetry statistic, and present analytical and numerical results demonstrating its robustness as a quantifier of directional causality. We explore test performance both in the small-data frontier, as well as its asymptotic behavior for time series with more observations. Then we verify the method on multiple real-world datasets with known ground truths, showcasing our recommended workflow for data with uncertainties and a limited number of observations. Finally, we apply the method to paleoclimate time series, identifying atmospheric CO2 as a key driver of global sea level on glacial-interglacial time scales.
II A predictive test for causality
II.1 A causality statistic linked to the flow of dynamical systems
Consider a generic coupled dynamical system generated by the vector field
| (1a) | ||||
| (1b) | ||||
and denote by its evolution operator. Without loss of generality, assume that , so that there is coupling between the variables, and define the time series and . We introduce the following difference between information-theoretic functionals as a causality statistic:
| (2) |
for some prediction lag , where is the TE corresponding to the prediction lag (Appendix E). This quantity, which is a difference in predictability forwards and backwards in time using TE, measures the evolution of the system through the equality . Similar identities are also found using mutual information (Appendix B). We show that
| (3) |
where is a generalized delay reconstruction of the dynamics from and , is the invariant distribution over and is a quantity closely related to the flow of the system as
| (4) |
We show that if there is no direct coupling , then when using an appropriate delay reconstruction (Appendix B). The positivity or negativity of in the general case is not obvious, but for a widely used family of stochastic systems, we can show that if a coupling exists, then the sign of (while ), and that when no coupling exists.
II.2 Sign and magnitude of predictive asymmetry reflects underlying coupling
For systems of random variables, marginal entropies can be computed directly from the covariance matrix of the system Hahs and Pethel 2013. This allows us to obtain exact expressions for the predictive asymmetries (eq. 2) for stochastic processes with known parameters. Here, we demonstrate predictive asymmetries on the following unidirectionally coupled, stationary autoregressive system with , , and innovations and (Appendix C):
| (5a) | ||||
| (5b) | ||||
When the dynamical variables are decoupled, predictive asymmetries are zero in both directions (Fig. 1). When coupling exists, is positive and increases monotonically with . Conversely, is negative and decreases with increasing . Contributions to are most pronounced at low , and diminish for higher (Fig. 1D). The predictive asymmetry therefore plateaus at some system-specific threshold value of , which is the time horizon beyond which information about the forcing is no longer detectable in the response. Lagged information beyond this threshold does not contribute substantially to the predictive asymmetry, because beyond some system-specific and data resolution specific time lag, the extra information is irrelevant to the interaction. In other words, the influence of on fades due to vanishing covariance between states that are further apart in time.
For dynamical systems in general, this convergence can be understood as follows (Appendix D). If there is a dynamical link in the direction , then the influence that current values of have on future values of (forwards-in-time prediction) is expected to be stronger than the influence that current values of have on past values of (backwards-in-time prediction), i.e. for every . Hence, we expect if influences . In the case of bidirectional influence , we expect that both . Moreover, the relative magnitudes of and preserve the rank order of the underlying coupling strength (Fig. 1E-H).
II.3 An intrinsic significance criterion: the normalized predictive asymmetry test
In practice, a strict criterion of is not ideal for testing the hypothesis of directional causality. Governing equations of observed processes are rarely available in practice, so exact values for the statistic are unobtainable. Predictive asymmetries must therefore be approximated from phase space reconstructions from observed time series data (Appendix E). In the limit of few observations and low coupling strength, statistical fluctuations will cause numerical estimates of to deviate from the true value.
As demonstrated in appendix A, the value of TE estimates can vary greatly depending on the system. Here we leverage this system-specific TE as a dynamically informed, intrinsic significance criterion. We thus introduce a normalized causality statistic by dividing the predictive asymmetry on the TE of the observed time series integrated across the spectrum of prediction lags:
| (6) |
This normalization to the system-intrinsic TE allows the relative magnitude of predictive asymmetry to be compared across systems. Furthermore, can be treated as a binary classifier of directional causality, so that indicates a positive detection of an influence from to . The stringency of the test is then determined by the constant . We use the criterion with (i.e. normalizing to the mean TE; appendix G) to determine the statistical robustness of the test, defining as a detection of directional coupling from to (i.e. a ”positive”), and as a detection of a non-interaction (i.e. a ”negative”).
Non-parametric approaches to detecting causality from time series typically rely on the method of surrogates (Theiler et al. 1992; Lancaster et al. 2018) to avoid spurious results. A surrogate time series is a randomized or modeled version of the original time series, designed to establish a baseline for significance testing. The value of a causality statistic is deemed significant if it exceeds some threshold value obtained in a large ensemble of surrogates. An ideal surrogate for causality testing preserves all statistical properties of the original signal, except the property of being causally connected. The design of appropriate surrogate data remains a thorny problem.
Our predictive asymmetry approach solves this problem by design. The statistic compares the magnitude of forward-in-time TE to its complementary backward-in-time TE, computed on the same time series. The backward-in-time TE can we viewed as a system-intrinsic and estimator-specific ”reversed-time surrogate”, and is part of the test by construction. This time reversal preserves all properties of the signal, but explicitly breaks causality by reversing the arrow of time.
Estimator specific bias, which arises due to the fact that the asymptotic distribution of the sample statistic for TE is not known analytically Barnett et al. 2009, and due to disparate frequencies in the time series (which are known to plague TE; Paluš and Vejmelka 2007; Hannisdal 2011), are thus equally encoded in both the backward-in-time predictions and the forward-in-time predictions for stationary systems. Any asymmetry that remains is due to forwards-in-time information flow in the presence of directional coupling. With time series recorded from variables that are not connected, time-reversal has no effect, so forwards-in-time and backwards-in-time predictions are balanced (Fig. 1; formally proved in appendices B and C). Predictive asymmetry thus arises intrinsically from causal connectivity in the underlying system.
Because of its built-in significance test, predictive asymmetry does not require explicit surrogate testing, which drastically reduces its computational demands. A potentially powerful application of the method would thus be to initially screen for the presence of causal relationships in large time series ensembles. Our method could also be applied in continuous monitoring of real systems, to detect time-variable dynamical interactions.
III Identifying directional causation from time series
Here, we apply the method to multiple synthetic coupled stochastic and dynamical systems with known governing equations, and show that the predictive asymmetry yields a stand-alone criterion for the detection of directional causality from time series.
III.1 Two example systems
Consider a chaotic interaction model for two species, and Diego et al. 2018, where coupling can be absent, unidirectional, or bidirectional (Fig. 2), and a common-cause model where two non-interacting variables with nonlinear deterministic dynamics, and , respond to the same external forcing (Fig. 3). These examples serve to illustrate two important hurdles in causality testing that our statistic should reliably overcome: (1) distinguishing between uncoupled, unidirectional, and bidirectional relationships, and (2) distinguishing correlation from causation in systems with uncoupled variables responding to a common external driver, which may introduce strong correlation between the uncoupled variables. The common-cause model specifically also targets the issue of identifying causation in time series with relatively strong periodicities, which are pervasive in paleoclimate time series like the ones we analyze below.
Characteristic asymmetries for uncoupled variables. Absence of coupling consistently yields predictive asymmetry distributions centered around zero. Hence, when applied to the common-cause scenario, the normalized test reveals no evidence of directional coupling between the non-interacting variables, despite the common external forcing (Fig. 3A-B). The same holds for two-species chaotic model when there is no underlying coupling (Fig. 2). If two time series and are recorded from independent systems, then predicting future values of from present values of , and predicting past values of from present values of , are equally (un)informative.
As we expand the prediction window, statistical noise is introduced by the inclusion of more non-informative history, which results in increasing variability of for increasing . Access to more observations counteracts this effect, reducing the dispersion of , and yielding more narrow-tailed, zero-centered distributions (Appendix G.1). As expected, having more information about two unrelated variables increases our ability to reject a coupling between them. In summary, when no underlying coupling exists, predictions backwards and forwards in time are on average of similar magnitude, and is thus centered around zero across a range of prediction lags.
Characteristic asymmetries for unidirectionally coupled variables. In contrast, unidirectional coupling manifests as positive predictive asymmetry in the causal direction (driver response), and negative predictive asymmetry in the non-causal direction (response driver) (Appendix H.2).
If we reverse the direction of coupling, then the signs of will follow suit (Fig. 2).
Why does this happen? If a unidirectional coupling from to exists, forwards-in-time prediction (values of predict future values of ) is stronger than backwards-in-time prediction (values of predicts past values of ). The opposite happens in the non-causal direction (): backwards-in-time prediction (future values of predicts past values of ) becomes more successful than forwards-in-time prediction (past values of predicts future values of ). In uncoupled systems, on the other hand, there is no shared information that improves prediction neither forwards in time nor backwards in time. Time-asymmetric predictability is thus characteristic of systems with directional coupling.
Characteristic asymmetries for bidirectionally coupled variables.
Bidirectional coupling between variables yields predictive asymmetries that are on average positive in both directions (Fig. 2; Appendix H.3), and with relative magnitudes reflecting the underlying coupling strengths. If the difference between the underlying relative coupling strengths for a bidirectional system is substantial, then values of in the direction of the weaker forcing may approach zero. Thus, bidirectional coupling is most likely detected if coupling strengths in both directions are roughly equal, whereas if coupling strengths are significantly different, then the system may appear unidirectional in the eyes of the test (Appendix G.3). Unlike unidirectionally coupled stochastic AR systems, for which the predictive asymmetry works remarkably well, the results of the predictive asymmetry test for bidirectionally coupled AR systems are sensitive to model parameters (as is TE alone). We stress, however, that our approach rests on the existence of a fundamental connection between the predictive asymmetry based on attractor reconstruction and the flow of the underlying dynamical system. Whether an equivalent fundamental connection exists for stochastic systems is a topic for future research.
Characteristic asymmetries for causal chains. Relative positive/negative magnitudes of for the model systems depend on the underlying coupling strength. Pairwise application to variables of multidimensional systems with chained unidirectional coupling shows that magnitude of also depends on the number of intermediate variables: The magnitude of is greater for adjacent nodes in the interaction network and decreases with an increasing number of intermediate, indirect links (Appendix H.1).
III.2 Statistical robustness
How robust is the predictive asymmetry as a causality detection criterion? For time series generated from synthetic systems, we compare the results of the normalized predictive asymmetry test as a binary classifier (eq. 6) with the known ground truths. We repeat this process for a large number of parameterisations of different systems with varying coupling strengths, each parameterisation yielding a distinct dynamical system and a unique set of corresponding time series with varying statistical properties. Thus, we obtain counts of false positive, false negative, true positive and true negative detections. Their corresponding rates are summarized in confusion matrices, from which we compute Matthews’ correlation coefficient (MCC) Matthews 1975; Chicco 2017 and other test performance indicators. The MCC takes on values on , where MCC = 1 indicates perfect agreement between actual values and predictions, and MCC = 0 indicates no correlation between predictions and actual values.
In our sensitivity tests, we find that for sufficient coupling strength and time series length, MCC converges to high values (¿ 0.8) for all test systems, including stochastic, periodic and nonlinear dynamics with different types of coupling, and chaotic dynamics where coupling strengths are below synchronization thresholds (Fig. 4). For the common-cause model, which is strongly periodic, MCC converges to values above for time series with 500 observations or more. We emphasize that these results are obtained by applying the normalized predictive asymmetry test alone, without any surrogate testing.
IV Application to real data
An ensemble approach is needed to statistically characterize the predictive asymmetry in a system, as demonstrated above for synthetic systems. In real-world applications, time series typically represent a single realization of the system observed over a limited time window, and the governing equations are not known. Nonetheless, we can estimate ”ensemble statistics” for empirical data in two ways: (1) We generate a distribution of values computed on random segments of the time series, varying the length and position of each segment. (2) If we have information on the uncertainty of the observed data (e.g. standard errors), then we generate a distribution of values by random resampling within the uncertainty bounds of the data. In certain situations we need to resample also within the uncertainties in the time index (e.g. age estimates in paleoclimate proxy records; see example below). Note that these two resampling approaches can be combined Haaga 2019. A sliding-window approach is also possible, but we limit the present study to the analysis of ensemble predictive asymmetry averaged over the total window of observation.
In appendix I we characterize causal interactions from multiple real-world data sets for which the true causality is known. In each case where we know the ground truth, the ensemble predictive asymmetry correctly determines the underlying directional coupling. Here, we address the causal interactions between the key climate system parameters of atmospheric CO2, global sea level and summer insolation at 65∘N during the last 800 thousand years.
Since the discovery of the ice ages in the early 19th century Esmark 1824, many hypotheses have been put forward to explain the recurrent waxing and waning of Pleistocene ice sheets. During glacial intervals, these ice sheets sequestered huge volumes of fresh water, thus controlling global mean sea level, which has varied by up to m during the past 800 kyr (Fig. 5). Chief among the causal explanations is variability in Earth’s orbit Esmark 1824; Croll 1875; Milankovitch 1941; Hays et al. 1976, which affects the seasonal distribution of solar energy reaching the Earth. Difficulties remain, however, in explaining the 100 kyr saw-tooth pattern characteristic of the major late Pleistocene ice ages. Despite a similar periodicity, the energy forcing associated with changes in orbital eccentricity is negligible, hence several hypotheses have been proposed to explain the deep glacial maxima and their abrupt terminations Pisias and Moore Jr 1981; Raymo 1997; Tzedakis et al. 2017; Denton et al. 2010; Wolff et al. 2009.
When ancient air bubbles trapped in Antarctic ice cores revealed that fluctuations in atmospheric CO2 were tightly linked to ice volume changes throughout the glacial-interglacial cycles (Fig. 5C), changes in radiative forcing due to greenhouse gases were implicated in the dynamics of ice ages Petit et al. 1999; Shackleton 2000. Proxy reconstructions and transient modelling point to CO2 as a forcing of global temperature rise during the last glacial termination Shakun et al. 2012. However, the drive-response relationship between CO2 and global ice volume remains controversial. On the one hand, coupled ice sheet and general circulation models are able to recreate the saw-tooth pattern by internal feedbacks without CO2 forcing Abe-Ouchi et al. 2013. On the other hand, it has been proposed that glacial terminations were a response to CO2 release from warming southern oceans and associated changes in atmospheric and oceanic circulation (Wolff et al. 2009; Denton et al. 2010). The impetus for this mechanism is thought to be the long-lasting impact of meltwater from the massive circum-North Atlantic ice sheets that formed during glacial maxima, which created an ice sheet – CO2 feedback loop mediated by ocean circulation. In this view, the orbital variability acts as a ”pacemaker” for the ice sheet – CO2 system rather than being the primary driver of the 100 kyr ice age cycles Raymo 1997.
Here we test the causal pathways among key climate system variables in the late Pleistocene: insolation, atmospheric CO2 concentration, and global sea level (ice volume). As an external variable we use the canonical June 21 insolation at 65∘N Laskar et al. 2004 (Fig. 5A), which is typically used as a proxy for the astronomical forcing linked to the growth and decay of large ice sheets in the Northern Hemisphere through the Pleistocene epoch Milankovitch 1941; Hays et al. 1976; Raymo 1997; Huybers and Wunsch 2005; Haaga et al. 2018. We use a composite ice core record of atmospheric CO2 Bereiter et al. 2015 with reported means and standard errors for the CO2 measurements, and the AICC2012 ice core chronology with associated age uncertainties Bazin et al. 2013; Veres et al. 2013 (Fig. 5B). A global sea level stack Spratt and Lisiecki 2016, reported as first principal component scores in 1-kyr bins, with a 95 % confidence envelope accommodating uncertainty in both sea level estimates and ages, serves as a proxy for ice volume (Fig. 5C). By combining the random segment and uncertainty resampling, we obtain ensembles of predictive asymmetries for the three pairwise comparisons (Fig. 5D-F).
The dynamical evidence in the data shows that climate-intrinsic radiative forcing has a significant influence on the long-term evolution of global sea level (Fig. 5D). Changes in the planet’s energy budget and seasonal energy distribution caused by oscillations of Earth’s orbit also seem to be a strong direct driver of sea level as represented by the global stack (Fig. 5F). Relatively speaking, the influence of CO2-driven radiative forcing on sea level is greater than that of northern summer insolation. Furthermore, insolation is not a significant driver of atmospheric CO2 (Fig. 5E), indicating that the CO2 forcing of sea level is independent from the insolation forcing of sea level and not a mutual response to orbital variability.
Our analysis takes into account the reported uncertainty in the CO2 and sea level estimates as well as uncertainty in the associated ages. One important caveat, however, is that the age model of the global sea-level stack inherently assumes a lagged response to orbital forcing Spratt and Lisiecki 2016. To assess the impact of this assumption on the dynamical information in the paleoclimate records, we repeated our analysis on a 500-kyr record of sea level from Grant et al. Grant et al. 2014, which is chronologically independent of orbital parameters (Appendix J). In the orbitally independent sea level record, the evidence for insolation forcing of global ice volume is drastically reduced (Fig. S.J33), highlighting the importance of age model assumptions in determining dynamical information in geological records. Nevertheless, the results clearly confirm the strong influence of CO2 on global ice volume (Fig. S.J33). We thus conclude that state-of-the-art paleoclimate records, despite uncertainty in estimates and ages, strongly suggest that CO2 was a major dynamical driver of glacial-interglacial ice volume variability.
V Concluding Remarks
Inferring the strength and directionality of causal interactions from observed time series, without recourse to mechanistic modeling, is a subject of controversy. Transfer entropy, or generalized Granger causality, has seen widespread application across many disciplines Bossomaier et al. 2016. Related approaches to predictive causality based on dynamical systems reconstruction have also garnered considerable attention, including the geometric prediction method of convergent cross mapping Sugihara et al. 2012. On a practical level, however, both transfer entropy and cross mapping have serious limitations. Both require ad hoc interpretation of prediction skill to distinguish non-causal from causal coupling and both rely on the method of surrogate testing as a bulwark against false positives. On a more fundamental level, and to the best of our knowledge, neither transfer entropy nor cross mapping have an explicit, precise relation to the flow of a dynamical system. The novelty of our contribution lies in showing that a simple difference between transfer entropy for forwards and backwards prediction lags quantifies the flow of the underlying dynamical system. The resulting predictive asymmetry provides a theoretically founded causality statistic, with robust numerical performance for deterministic and stochastic systems. Hence, our test can unambiguously resolve causal relationships in nonlinear systems where transfer entropy alone fails (Appendix A). Moreover, our test easily characterizes the type of systems originally used to demonstrate cross mapping (Fig. 2), where the latter requires additional hypothesis testing. With its built-in significance criterion, our method eliminates costly surrogate testing, which raises the prospect of causal network reconstruction in large data sets and monitoring applications. By linking predictive asymmetry to dynamical causality, our work represents a major advance in the causal analysis of observations, with implications for any field where time series are used to study the dynamics of natural systems.
Supplementary materials
K. A. Haaga et al.
A simple test for causal asymmetry in complex systems
Appendix A Comparison of transfer entropy and predictive asymmetries
Here, we demonstrate some common issues associated with time series causality methods (exemplified using TE) and how our predictive asymmetry method overcomes these issues.
Time series generated from unrelated variables may give nonzero TE, which can be of equal magnitude in both directions (Fig. S.A1). Even for systems with unidirectional causation between variables, TE may be of similar average magnitude both for the causal and for the non-causal direction (Fig. S.A2). Without additional information, both these cases could be erroneously interpreted as a bidirectional, equal-strength influence between the variables. In a more favorable scenario, unidirectional coupling yields forwards-in-time TE that is higher in the causal direction than in the non-causal direction. One could misinterpret this result as bidirectional coupling with dominant control from one variable to the other, but at least the preferred direction of information flow is correctly detected. More disturbingly, TE values can be greater in the non-causal than in the causal direction in a unidirectionally coupled system at certain prediction lags (Fig. S.A2), which at face value would seem to imply bidirectional interaction with dominant control in the non-causal direction. Another challenge emerges in strongly periodic systems, both coupled and decoupled, which also yield oscillating TE values that are hard to interpret and could lead to the erroneous inference of two-way coupling (Figs. S.A3, S.A4). Similar false causalities also obtain with other causality statistics, and are difficult to remedy, even with external hypothesis testing using surrogate data (Krakovská et al. 2018, e.g. ).
Appendix B Formal proofs
B.1 Relation between predictive asymmetry and the underlying dynamics of the system
In this section we provide a formal derivation of the relation between the causality asymmetry and the underlying dynamics. The result applies for both discrete and continuous systems. A generic continuous system is generated by a vector field of the form
| (S.B1) | ||||
| (S.B2) |
while a generic discrete system is generated by the map
| (S.B3) | ||||
| (S.B4) |
where is smoothly invertible. Let , generically denote the evolution operator of the system. In the discrete case, and and , that is: is the -fold composition of . In the continuous case, is the solution to the differential equation with the initial condition , for all , that is: is the flow of the vector field. Let and denote the function components of corresponding to the and axes, respectively, and w.l.o.g., assume that . For fixed , define the variables
| (S.B5) | ||||
| (S.B6) | ||||
| (S.B7) | ||||
| (S.B8) |
and assume that the maps
| (S.B9) | ||||
| (S.B10) |
are diffeomorphisms over and call and . Notice that and thus the map
| (S.B12) |
generates the desired change of coordinates. Equivalently, . Its inverse map is of the same form, that is:
| (S.B13) |
Let be an invariant distribution of the system and let and be the expression of on the coordinates corresponding to and , respectively. Therefore, for any measurable set it must hold that
| (S.B14) |
where , and we have used that , with . In our case this implies
| (S.B15) | ||||
| (S.B16) |
The transfer entropies and are given by
| (S.B17) | ||||
| (S.B18) |
It is convenient to re express the above TE in terms of mutual information. In particular, defining the quantities
| (S.B19) | ||||
| (S.B20) |
and using the identities
| (S.B21) | ||||
| (S.B22) |
one easily checks that
| (S.B23) | ||||
| (S.B24) |
and hence their difference is expressed as
| (S.B25) |
We claim that . To see this, we first show that the marginals and coincide:
| (S.B26) |
where and analogously for . From this it follows that
| (S.B27) | ||||
| (S.B28) | ||||
| (S.B29) |
In addition, for a fixed , the straight line corresponds to the curve given by . Notice that neither the set nor the curve do really depend on . The presence of in their definitions is a pure formality so that the pair can be assigned a unique pair and vice versa. More specifically, and . It then holds that
| (S.B30) | ||||
| (S.B31) |
Defining the quantity
| (S.B32) |
where denotes the marginal of on the plane , we have that
| (S.B33) |
From this it follows that
| (S.B34) |
Thus the difference between the TE is given by
| (S.B35) |
Now we may use the identity
| (S.B36) | ||||
to deduce that and from here
| (S.B37) |
Finally, let denote the original phase space of the system with denoting the Lebesgue volume element in . Then the above expression can be written as
| (S.B38) | ||||
| (S.B39) |
Notice that the function carries the dependence on the embedding as .
B.2 for unidirectionally coupled systems
Suppose the dynamical system is generated by the vector field
| (S.B40) | ||||
| (S.B41) |
respectively, by the map
| (S.B42) | ||||
| (S.B43) |
| (S.B44) | ||||
| (S.B45) | ||||
| (S.B46) | ||||
| (S.B47) |
Assume that the map , for some , is (smoothly) invertible, then we can write , where . The invertibility of is equivalent to saying that the delay reconstruction reproduces the dynamics generated by (respectively, by ). In this case it follows that which we can generically denoted as and hence the map generating the change of variables (equation [S.B12]) reduces in this case to
| (S.B48) |
The key difference now is that the variable is decoupled from the pair with respect to the map connecting the embeddings and . In other words, the lack of coupling is transmitted into the embeddings and as a decoupling between and . Such a decoupling implies that will be a function of alone. In addition, the set is mapped by into the set , where we have used that . Therefore equation [S.B31] reduces to
| (S.B49) |
and from this it follows that
| (S.B50) |
and hence
| (S.B51) |
This result tells us that for unidirectionally coupled systems in the direction, say , the difference yields 0.
B.3 Predictive asymmetry based on mutual information
In this section we show that a similar relation with the underlying flow is found for a predictive asymmetry based on mutual information. More precisely, we consider again a generic dynamical system as in equations [S.B1] and [S.B2], also assuming that there is a direct coupling , and we use the same definitions as in equations [S.B5]-[S.B8]. We may consider, for instance, the asymmetry
| (S.B52) |
In this case, the equality in equation [S.B34] also implies that
| (S.B53) |
A less repetitive example could be the asymmetry
| (S.B54) |
In this case, for a fixed , one finds
| (S.B55) | ||||
| (S.B56) |
where
| (S.B57) | ||||
| (S.B58) |
and
| (S.B59) |
Therefore,
| (S.B60) |
and from this it follows
| (S.B61) |
Appendix C Exact expressions for the predictive asymmetry for autoregressive systems
We will expand on the approach from Hahs and Pethel 2013 to compute the predictive asymmetry for arbitrary prediction lags for a coupled bivariate autoregressive system. For such systems, marginal entropies may be computed analytically.
C.1 Covariance matrix for unidirectionally coupled AR1 system
Consider a simple unidirectionally coupled bivariate AR system with coefficients chosen such that the system is stationary (this is the same system as the main text’s eq. 5). Let and be the standard deviations of two independent normal distributions, where the noise draws and are independent at each time step.
| (S.C62) | ||||
| (S.C63) |
C.1.1 Variances for and
Due to stationarity, which we have by definition, and . We also find
and
C.1.2 Auto-covariances for and
The covariances between observations of separated by time steps are therefore given by the following expectations
with
thus we find
To evaluate the autocovariance of , we use that
and also that
Then,
If , then so, altogether we find that
where if and otherwise.
C.1.3 Cross-covariances for and
Using that we find that
This expression is symmetric in the time labels, so it does not matter which one of and is the smallest. With no loss of generality we may thus suppose that . Defining , we have that and therefore . Accordingly,
C.2 Filling the covariance matrix
Now that we have established the dependence of the covariance between time steps spaced arbitrary far from each other on time on the coefficients and , we can proceed with predictive asymmetry computations. Let be the maximum prediction lag, and let
and let denote the covariance matrix for . Equivalently, if choosing odd, we may shift the time indices and consider
Now, provides sufficient information to compute transfer entropy for maximum prediction lag , alternatively, the predictive asymmetry for maximum prediction lag , assuming the lags included in the history for the target variable does not exceed .
Computing the predictive asymmetry for a maximum prediction lag , the covariance matrix will have thus dimensions -by-, where , accounting for lags for and lags for , plus the zero lag cases. If dealing with a random system such as an autoregressive order-1 system (AR1), then this is all we need for transfer entropy computations. We just need to subset the relevant portions of the covariance matrix, compute relevant entropies, and from that compute the predictive asymmetry.
C.2.1 Computing the predictive asymmetry
Let and denote two generic source and target process. For convenience of notation, let and denote the present and past of the source and target variables (time series). During computation, and are kept fixed. Next, let denote the time series , but lagged time steps into the future () or into the past (). Say we want to compute the predictive asymmetry from to . We then have
In terms of entropies, the conditional mutual information (CMI) terms are
and these entropies can be computed exactly from the covariance matrix, following the approach of Hahs and Pethel 2013.
C.3 Covariance matrix for bidirectionally coupled AR1 systems
C.3.1 General setting
Consider the general linear AR system
| (S.C64) | ||||
| (S.C65) |
where and . Expressed more compactly, the above system reads
| (S.C66) |
with and . The matrix , contains the coefficients of the model. We will assume that all the eigenvalues of have norm strictly less than 1. In addition and for the system to be stationary the column vector of expectation values of the variables must verify and therefore, from equation (S.C66), . Since, in particular, can not be an eigenvalue of , it must follow that .
C.3.2 Change of variables
By redefining linearly the variables as
| (S.C67) |
with being an invertible matrix (independent of time), the system in equation (S.C66) is transformed as
| (S.C68) | ||||
| (S.C69) | ||||
The variances and covariances between the new noise terms verify that
| (S.C74) | ||||
| (S.C75) | ||||
| (S.C76) |
where we have used that and . In addition, it is clear that
| (S.C77) |
and all the covariances between the variables and the noise terms do vanish.
Suppose that in the variables , the matrix adopts a particularly simple form so that the covariances
| (S.C78) | |||||
| (S.C79) |
are easily computed. In addition, the entries of are
| (S.C80) |
and hence
| (S.C81) | ||||
| (S.C82) |
Therefore we find
| (S.C83) | ||||
| (S.C84) | ||||
| (S.C85) | ||||
| (S.C86) |
In the following, we will apply these formulas to two inequivalent examples of bivariate AR models.
C.3.3 Example 1. The matrix has two distinct real eigenvalues.
In that case there is an invertible matrix such that and the system reduces to
| (S.C87) | ||||
| (S.C88) |
From equations (S.C87) and (S.C88) we deduce that and similarly for . In addition, . From these equalities we find
| (S.C89) | ||||
| (S.C90) | ||||
| (S.C91) |
where we have used the results found in equations (S.C74)-(S.C76). It also clearly holds that
| (S.C92) | ||||
| (S.C93) |
and therefore
| (S.C94) | ||||
| (S.C95) | ||||
| (S.C96) | ||||
| (S.C97) |
As an example of such a case, we will consider the system
| (S.C98) | ||||
| (S.C99) |
with and . In this case, the matrix has eigenvalues , and hence we demand that both and , be strictly less than 1. It is easy to check that the non singular matrix
| (S.C100) |
verifies that . Therefore, by substituting the corresponding and parameters in equations (S.C94)-(S.C97), we find
| (S.C101) | ||||
| (S.C102) | ||||
| (S.C103) | ||||
| (S.C104) |
and the results on the covariances between the original variables are obtained from equations (S.C83)-(S.C86) with , , and .
C.3.4 Example 2. The matrix has only one eigenvalue (and therefore real).
The trivial case in which is proportional to the identity is already contained in the previous case example. The non trivial instance of is therefore not diagonalizable. In that case there is an invertible matrix such that
| (S.C105) |
which corresponds to the Jordan normal form. Hence with the redefinitions
| (S.C106) | ||||
| (S.C107) | ||||
| (S.C108) | ||||
| (S.C109) |
the system in equation (S.C66) reduces to
| (S.C110) | ||||
| (S.C111) |
and still it holds that
| (S.C112) | ||||
| (S.C113) | ||||
| (S.C114) |
From equations (S.C110)-(S.C111) we find,
| (S.C115) |
and hence
| (S.C116) |
Also
and hence
| (S.C117) |
Finally,
| (S.C118) |
and therefore
| (S.C119) |
To compute the -th time step covariances, it is more convenient to express the system in matrix form as
| (S.C120) |
From here we easily deduce that
| (S.C121) |
In addition, one easily checks that , hence we have
| (S.C122) | ||||
| (S.C123) |
From here we deduce that
| (S.C124) | ||||
| (S.C125) | ||||
| (S.C126) | ||||
| (S.C127) |
In summary,
| (S.C128) | ||||
| (S.C129) | ||||
| (S.C130) | ||||
| (S.C131) | ||||
A generic case of this class of systems is given by
| (S.C132) | |||
| (S.C133) |
with , and . One checks that the matrix is brought into with the linear transformation
| (S.C134) |
The covariances can thus be found by using equations (S.C128)-(S.C131) and (S.C83)-(S.C86) with
| (S.C135) |
Appendix D Heuristic explanation for the sign of the predictive asymmetry in the general case
For the unidirectional AR systems we show analytically that the predictive asymmetry is negative in the non-causal direction. In the general case, this behavior may be heuristically understood as follows.
Recall the definition of the predictive asymmetry from variable to :
| (S.D136) |
D.0.1 Unidirectional coupling
Consider first the case of unidirectional coupling . For the simplest possible TE (three-dimensional) analysis, the forward-prediction term in the expression for (eq. S.D136) becomes
| (S.D137) |
which measures how much, on average, knowing something about the present of improves our ability to predict the future of . If does actually have an influence on , then we expect .
What about the second term? It is not very intuitive to think about backwards prediction. However, after some algebraic manipulation, the backward-prediction term reads:
The backwards-lag prediction thus quantifies how well the knowledge about the past of improves our prediction of the present of , given the present of . One may erroneously conclude that this term should be trivially zero — that if the dynamical influence is , then neither the past nor present of should not have a measureable effect on the present of . But if there is dynamical influence , then information about the past of is encoded in and therefore can be considered a proxy for past values of . Including information about may thus improve our prediction of the outcome of . If , then we expect to statistically quantify some influence due to being a proxy for 11 1 measures a TE-like quantity, except the causal direction flips and the conditioning in the argument of the logarithm occurs only on one variable, not on a mixture of the two variables, as for regular TE (eqs. S.D137 and S.D138a)..
In summary, compares the direct influence resulting from the forcing with the indirect influence arising through the interaction of with . The key concept of the asymmetry test is that when an underlying coupling exists, then the latter may be statistically detectable. A reliable causality estimator should be better at detecting direct influences than indirect influences, so we expect , and hence .
In the opposite direction , we are numerically estimating the following integrals
| (S.D138a) | ||||
| (S.D138b) | ||||
The forwards-prediction here quantifies the extent to which having information about the past of improves our knowledge about the future of . This term also measures an indirect effect of on it own future through its interaction with . The backwards-prediction represents the statistical measure of a direct influence (analogous but not equal to eq. S.D137). Hence, we expect and hence .
D.0.2 Bidirectional coupling
For systems that are bidirectionally coupled , the situation is a bit more complicated, but the same basic argument applies. Recall that
and
What happens if there is a difference in the coupling strengths? The only term that would pick up any change in the coupling strength is in the argument of the logarithm. If there is bidirectional coupling with coupling strengths , then we expect , or
which is the expected from the usual TE. What about the backwards-prediction terms? If , then
Thus, for , one subtracts — relatively speaking — a smaller indirect effect from the direct effect of . For , one subtracts — relatively speaking — a larger indirect effect from the direct effect of . The effect of having is therefore that .
Next, assume the coupling strengths and are of similar magnitude. Then, and , so that and . Hence, we expect and the magnitudes of and to reflect the coupling strengths and .
D.0.3 No coupling
If has no influence on , then we expect no detectable influence neither directly from present values to future values, nor (because there is no interaction) from the present of to its own future through its interaction with . Therefore, . Simply put, if there is not coupling, then on average none of the predictions involving the other variable will be improved.
Appendix E Experimental setup
E.1 Generalized embedding of time series for TE analysis
Generically denote the time series for the source process as , and the time series for the target process as , and as the time series for any conditional processes that might act in tandem with to influence . To compute (conditional) TE, we need a Generalized embedding Sauer et al. 1991; Deyle and Sugihara 2011 incorporating all of these processes.
For convenience, define the state vectors
| (S.E139) | ||||
| (S.E140) | ||||
| (S.E141) | ||||
| (S.E142) |
where the state vectors contain future values of the target variable, contain present and past values of the target variable, contain present and past values of the source variable, contain a total of present and past values of any conditional variable(s). Here, indicates the embedding lag. In real systems, the strategy for choosing depends on the temporal resolution and the auto-correlation function of the time series data. indicates the prediction lag (the lag of the influence the source has on the target). Combining all variables, we have the Generalized embedding
| (S.E143) |
with a total embedding dimension of . Here, only depends on the prediction lag , which is to be determined by the analyst; we use multiple negative and positives s for computing . The remaining variables depend on , which may be determined from, for example, the minima of the auto-correlation or lagged mutual information function of the time series. For the synthetic examples in this paper, we push the lower limits of time series lengths, so we use to not exclude too many data points. Another reason for choosing is that theoretically it should worsen the performance of the TE method, because it leads to strongly auto-correlated reconstructed states when there is auto-correlation in the time series Kantz and Schreiber 2004. As we shall see, however, this deliberate choice does not diminish the ability of the asymmetry criterion to distinguish directional dynamical influence.
E.2 Transfer entropy
TE (in bits) from a source variable to a target variable with conditioning on variable(s) is defined as
| (S.E144) |
Without conditioning, eq. S.E144 becomes
| (S.E145) |
E.3 Numerically estimating TE and
We have used three different TE estimators to estimate . The visitation frequency estimator () Schreiber 2000 computes TE by partitioning the reconstructed state space using a regular binning. The invariant probability over the partition is then estimated by computing how often the orbit visits each box. From this invariant joint density, we obtain marginal densities, and TE can be computed using equation S.E145. The transfer operator grid estimator () is also based on a partition of the reconstructed state space using a regular grid. However, the joint probability distribution is computed as the invariant distribution of an approximation to the transfer operator associated with the system Diego et al. 2018. Lastly, the nearest neighbour based estimator () uses mutual information (MI) Cover and Thomas 2012 to compute TE through the identity . The MI is estimated using the Kraskov estimator Kraskov et al. 2004. All estimators are available in the CausalityTools.jl Julia software package (https://github.com/kahaaga/CausalityTools.jl). Finding roughly equivalent results for all three estimators, results presented in the text are generated using the visitation frequency estimator .
For our analyzes, we implement the embedding approach appearing in Krakovská et al. 2018, in which increasing the dimension of the embedded system implies conditioning on longer sequences of the past of the target variable. Thus, we use time delay embeddings of the same structure i.e., , where , and .
The TE for the binning-based estimators, and , is obtained as follows. We consider two different partitions constructed by subdividing each coordinate axis into an integer number of equal-length interval, using two separate partitions. The number of intervals for these partitions are selected according to the following heuristic and , roughly following Krakovská et al. 2018, where is the number of points in the time series and is the embedding dimension. The corresponding absolute bin sizes, and , are then computed for each and the TE is obtained as an average over those two partitions.
For the estimator, we use the Chebyshev distance metric, and use a different numbers of nearest neighbours for the estimation of each MI term. For all examples, we let be the number of nearest neighbours used for the highest dimensional MI estimate (), and be the number of nearest neighbours for the lowest dimensional MI estimate.
Because is computed as a difference between sums of TE values at symmetric prediction lags, which should not be sensitive to absolute TE values, we do not correct for estimator-intrinsic differences in absolute TE values Bossomaier et al. 2016.
Appendix F Test systems with known ground truths
F.1 Bidirectionally coupled logistic maps
For the main text, we use logistic model for the chaotic population dynamics of two interacting species given by
| (S.F146a) | ||||
| (S.F146b) | ||||
| (S.F146c) | ||||
| (S.F146d) | ||||
where the coupling strength controls how strongly species influences species , and vice versa for . To simulate time-varying influence of unobserved processes, we use the dynamical noise terms and drawn independently at each time step. If , then the influence of on is masked by dynamical noise equivalent to at the -th iteration of the map, and vice versa for .
F.2 Non-interacting variables forced by common driver
In this system, two non-interacting variables and are affected by a common external driver . All variables have nonlinear deterministic internal dynamics, overprinted by cyclic and stochastic variability, simulating typical paleoclimate time series from the most recent Quaternary period of Earth’s history (Fig. 5). The coupling between the noninteracting variables and the external forcing is highly nonlinear.
The model time step is defined as 1 kiloyears, simulating typical Quaternary paleoclimate time series. We randomly draw the periods , which yields the typical orbital-type dominant frequencies that are pervasive in paleoclimate time series. The signals are also phase-shifted by randomly assigning values to (Fig. 5). Initial conditions are drawn from a uniform distribution over the unit interval.
| (S.F147a) | ||||
| (S.F147b) | ||||
| (S.F147c) | ||||
We generate time series ensembles by drawing parameters randomly from uniform distributions as follows: , , and . are random uniformly distributed processes drawn independently from at each time step, where control the magnitude of the dynamical noise. Additionally, observational noise equivalent to 0.5 the standard deviation of each time series is added to that time series before analyses. Interaction lags and are set to 1, and internal lags are drawn randomly from the set .
F.3 Vector autoregressive (VAR) processes
Consider a -dimensional random sequence . A -dimensional process is given by
| (S.F148) |
where is a noise sequence drawn from a zero-mean normal distribution with standard deviation and is a deterministic sequence. The coefficient matrices have dimensions -by-, where the entry is the coefficient controlling the influence of the -th variable on the -th variable at time lag . Interaction strengths between variables are thus governed by the off-diagonal terms of these coefficient matrices.
For example, consider the following -dimensional AR(2) system with no deterministic sequence:
| (S.F149) |
In matrix notation we have
| (S.F150) |
Stability of the VAR process is ensured if the roots of the -by--dimensional companion matrix , as defined below, lie inside the unit circle.
| (S.F151) |
Hence, to generate stationary time series from a -process with variables, we assign coefficients such that for all .
F.4 Normally distributed noise processes
To establish a baseline for the significance threshold for the normalized predictive asymmetry test (eq. 6), we will consider various uncoupled noise processes.
| (S.F152) | |||
| (S.F153) |
where and are independent normal distributions with zero mean and standard deviations and , and values are drawn independently at each time step.
F.5 Uniformly distributed noise processes
Next, we will consider uniformly distributed noise processes
| (S.F154) | |||
| (S.F155) |
where and are independent uniform distributions with support , and values are drawn independently at each time step.
F.6 Brownian noise based on uncorrelated uniform noise
Many observed time series are trended. We simulate this phenomenon by using brownian noise processes
| (S.F156) | |||
| (S.F157) |
where and are independent uniform distributions with support , and values are drawn independently at each time step.
F.7 Autoregressive systems with periodicity and strongly nonlinear couplings
This system of unidirectionally chained autoregressive variables was extended from a simpler version in Péguin-Feissolle and Teräsvirta 1999 and Chávez et al. 2003, introducing a periodic component and variable parameters, variable internal lags and variable interaction lags.
Interaction lags are kept constant with for each instance of the system. Internal lags , as well as and are selected randomly from the set with uniform probability. The are independent normally distributed dynamical noise processes with zero mean and standard deviations of . and control the period and phase of the periodic component of the -th variable, while scales the magnitudes of the periodic component. The regulate the magnitude of the periodic components of at each time step. The coupling strength between nodes and in the chain is controlled by the parameter . The logistic function responsible for the coupling between adjacent variables and is parameterized to simulate a wide range of couplings.
Observational noise equivalent to 20% of the standard deviation of the respective variable is added to each time series. Parameters are drawn from uniform distributions as specified in the figure texts.
| (S.F158a) | ||||
| (S.F158b) | ||||
F.8 Nonlinear systems with linear coupling over multiple forcing lags
This nonlinear system with linear coupling is modified from Chen et al. 2004, but expanding the system to a chain of unidirectionally coupled variables with variable internal lags and interaction lags.
| (S.F159a) | ||||
| (S.F159b) | ||||
F.9 Nonlinear systems with nonlinear coupling over multiple forcing lags
This nonlinear system is also modified from Chen et al. 2004, but expanding the system to a chain of unidirectionally coupled variables with variable internal lags and interaction lags. Here, the coupling is nonlinear.
| (S.F160a) | ||||
| (S.F160b) | ||||
F.10 Nonlinear systems with periodic component, linear coupling
This system is modified from Chen et al. 2004, but expanding the system to a chain of unidirectionally coupled variables with variable internal lags and interaction lags. Cyclic components have also been introducing to the signals, where and controls the period and phase of the periodic component of the -th variable. The interaction between adjacent nodes in the chain is linear, and the coupling strength is controlled by the parameter .
| (S.F161a) | ||||
| (S.F161b) | ||||
F.11 Unidirectional chain of logistic maps with variable internal lags and forcing lags
The following system of difference equations describes a -dimensional system of logistic maps. Its interaction network is characterised by unidirectional coupling between adjacent nodes. The first map is independent, while for the remaining maps, the map is affected by itself and the -th map.
Here, is the value of the -th variable at time and is the value of the -th variable at time . The strength of the unidirectional forcing from to is controlled by . However, the influence from to is also masked by dynamical noise. The average magnitude of this noise is given by (given as a percentage of the allowed range of values, so a relatively low value should be chosen to not completely obscure the signals). The noise, , is dynamical noise drawn independently from a uniform distribution independently at every time step (masking the influence of on ), and then scaled by .
For the lags, is the time lag of the influence from to , and is the time lag of the influence from on itself, chosen randomly from the allowed lags for each variable .
For every realization of the system, initial conditions and parameters are randomized over uniform distributions and , respectively. Observational noise equivalent to 20% of the standard deviation of the respective variable is added to each time series. The dynamical noise is set to for all interactions.
F.12 Unidirectional chain of Henon maps
Consider a -dimensional system consisting of unidirectionally coupled Henon maps given by
| (S.F163) |
where is the -th map, and and is the coupling strength from variable to . For every realization, initial conditions are drawn from uniform distributions over . Observational noise equivalent to 20% of the standard deviation of the respective variable is added to each time series.
F.13 Rössler-Lorenz system
Here, we use two coupled Rössler and Lorenz systems, where the Rössler subsystem unidirectionally drives the Lorenz subsystem.
| (S.F164a) | ||||
| (S.F164b) | ||||
| (S.F164c) | ||||
| (S.F164d) | ||||
| (S.F164e) | ||||
| (S.F164f) |
with the coupling constant .
F.14 Bidirectional nonlinear system
This system is also modified from Chen et al. 2004, but noise and cyclic components have also been introducing to the signals, where and controls the period and phase of the periodic component of the -th variable. The interaction between adjacent nodes in the chain is linear, and the coupling strength is controlled by the parameter .
| (S.F165a) | ||||
| (S.F165b) | ||||
We generate time series ensembles by drawing parameters randomly from uniform distributions as follows: ,, , and . are random uniformly distributed processes drawn independently from at each time step, where control the magnitude of the dynamical noise. Additionally, observational noise equivalent to 0.5 the standard deviation of each time series is added to that time series before analyses. Interaction lags and are set to 1, while internal lags are drawn randomly from the set and internal lags are set to 1.
F.15 Nonlinear system without dynamical noise and periodicity
This system is also modified from Chen et al. 2004, but contains no dynamical noise or periodicity.
| (S.F166a) | ||||
| (S.F166b) | ||||
Appendix G System-specific and estimator-specific significance test
G.1 Statistical robustness when no coupling is present
By relating the to some fraction of the system-specific and estimator-specific empirical TE, can be used as a criterion for statistical significance: indicates the presence of directional coupling, while rejects coupling. First, we demonstrate the difference between the predictive asymmetry (eq. 2) and its normalized counterpart (eq. 6), and how the value of affects the ability of the test to reject coupling for uncoupled systems (Figs. S.G5, S.G6, S.G7, S.G8, S.G9). Because there are no true positives when there is no coupling, only the TNR and FPR are meaningful summary statistics to use in this context. To determine the ability of the test to correctly classify absence of coupling, we hence compute TNR and FPR for multiple realizations of the following systems.
- •
Noise time series with uniformly distributed noise (Fig. S.G5)
- •
Noise time series with normally distributed noise (Fig. S.G6)
- •
Brownian noise time series (Fig. S.G7).
- •
Chain of periodic autoregressive variables with strongly nonlinear coupling (Fig. S.G8).
- •
Chain of nonlinear variables (Fig. S.G8).
- •
Another chain of nonlinear variables (Fig. S.G8).
- •
Chain of nonlinear, periodic variables with linear coupling (Fig. S.G8).
- •
Chain of logistic maps with dynamical noise (Fig. S.G8).
- •
Chain of Henon maps (Fig. S.G8).
- •
Non-coupled autoregressive systems of maximum order for very short time series (Fig. S.G9, upper panel).
- •
Non-coupled autoregressive systems of maximum order for longer time series (Fig. S.G9, lower panel).
We find that a normalization factor (i.e. normalizing to the mean TE) is a good-trade off between statistical robustness (here: the ability to reject coupling when there is none) and sensitivity to time series length. For a given time series length, higher reduces the number of false positives. Moreover, for all the tested systems, the ability of the test to reject coupling when there is none approaches perfect for sufficient time series length.
Upper panel: TNR and FPR in the limit of very short times. The improved performance around time series length 90 correspond to when the number of subdivisions along each axis changes due to the partition heuristic. For the particular case of no coupling, we find that better TNR and FPR are achieved by finer partitions, but when coupling exists, sticking to the partition heuristic yields better performance. For each time series length, is computed on 500 unique -dimensional systems with no coupling, where the order is randomly chosen from the set (eq. S.F149) for each system. For each system, the standard deviations for the error terms are drawn from uniform distributions on . Each variable affects itself at exactly one time lag, which is also chosen randomly from . The diagonal terms of the relevant coefficient matrices are independently drawn from a uniform distribution on , while off-diagonal terms of the are set to zero, yielding no coupling.
Lower panel: TNR and FPR for longer time series. For each time series length, is computed on 1000 unique -dimensional systems with no coupling and maximum order , randomly chosen from the set (eq. S.F149) for each system. The coefficient matrices are generated as for the short time series.
G.2 Statistical robustness for systems with unidirectional coupling
For systems with unidirectional coupling, we use the following single-valued statistics to evaluate the performance of the test: accuracy, sensitivity, TPR, TNR, FPR and FNR, and for some systems PPV (positive predictive value), NPP (negative predictive value) and the F1 score. We computed these statistics as a function of coupling strength and time series length for the systems listed below. We find that, provided sufficient coupling strength and long enough time series, the performance approaches perfect across all statistical performance measures for all the tested systems.
- •
Chain of periodic autoregressive variables with strongly nonlinear coupling (Fig. S.G10).
- •
Chain of nonlinear variables with linear coupling (Fig. S.G11).
- •
Chain of nonlinear variables with nonlinear coupling (Fig. S.G12).
- •
Chain of nonlinear, periodic variables with linear coupling (Fig. S.G13).
- •
Chain of logistic maps with dynamical noise (Fig. S.G14).
- •
Chain of Henon maps (Fig. S.G15).
- •
Unidirectionally coupled autoregressive systems of maximum order for longer time series (Fig. S.G16).
Diagonal terms of the relevant coefficient matrices are independently drawn from a uniform distribution on , and each variable affects itself at exactly one lag (i.e. its coefficient appears in only one of the coefficient matrices). For a particular coupling strength interval , off-diagonal terms giving rise to the coupling generated randomly from a uniform distribution on , i.e. the uppermost row in the heatmaps are generated with coupling strengths ranging from to . Coupling terms are generated such that there is only unidirectional coupling (i.e. for variables, one off-diagonal term is zero and the other is nonzero in the coefficient matrix). Standard deviations for the error terms are drawn from uniform distributions on , and are drawn independently for each variable.
Generalized embeddings were constructed with , and was computed at with . The color scheme is such that light gray corresponds to a rate of 0.8 and black to a rate of 0.5.
G.3 Statistical robustness for systems with bidirectional coupling
Here, we use sensitivity (TPR) and FNR to characterise the statistical robustness of the normalized predictive asymmetry test for bidirectional systems.
For an ensemble of realizations of the bidirectional logistic map system from the main text, we computed TPR and FNR as a function of coupling strengths in both directions for fixed time series length (Fig. S.G17). For this system, even for short time series (here 300 observations), the test consistently detects the bidirectional coupling across all but the lowest coupling strengths (for the most part, TPR ¿ 0.8). For higher coupling strengths, the variables become partially synchronized (but not completely due to the dynamical noise), so the detection rates weaken slightly.
We also explored another bidirectional system similar to the common-cause scenario in the main text. This system consists of two bidirectionally coupled nonlinear variables with periodicity and dynamical noise, where the coupling is also nonlinear (eq. S.F165). For this system, we find that longer time series are needed to consistently detect the bidirectional relationship between the variables (Fig. S.G18). If coupling strengths in both directions are non-vanishing and roughly equal, then the bidirectional relationship is detected most of the time (TPR ¿ 0.9). If the coupling in one direction is much stronger in one direction than in the other direction, then the system will appear unidirectional in the eyes of the predictive asymmetry test (S.H28; black heat map cells away from the diagonal in Fig. S.G18). This happens because the predictive asymmetry in the direction of the strongest forcing becomes positive, while in the direction of the weakest direction, the predictive asymmetry becomes negative (S.H28), essentially rendering the detectable relationship unidirectional. Stronger mutual coupling increases the deviation between relative coupling strengths that can be tolerated before the bidirectional relationship starts to appear unidirectional (wider blue areas around the diagonal for higher coupling strengths in Fig. S.G18).
Appendix H Characteristic predictive asymmetries
Here, we further demonstrate the typical asymmetries obtained for the cases unidirectional coupling and bidirectional coupling discussed in the main text. We consider the cases of unidirectional coupling and bidirectional coupling separately, visualizing the asymmetries using two types of plots: (1) Line plots with error bars of versus at a fixed time series length and coupling strength , and (2) Heat maps with average over multiple configurations of time series lengths and coupling strengths.
For all the tested systems, the ensemble median of the predictive asymmetry converges to values around zero with increasing variability for higher prediction lags (not shown here).
H.1 Causal chains
For unidirectional causal chains, the normalized predictive asymmetry is best at detecting adjacent links. Indirect links are also detectable, but predictive asymmetries decrease in absolute magnitude with an increasing number of intermediate links (Figs. S.H19, S.H20, S.H21, S.H22, S.H23, S.H24). The exact number of intermediate links that are detectable varies between systems. For our example systems, the predictive asymmetry gets indistinguishable from non-coupled systems after after two-to-three intermediate links.
- •
Chain of periodic autoregressive variables with strongly nonlinear coupling (Fig. S.H19).
- •
Chain of nonlinear variables with linear coupling (Fig. S.H20).
- •
Chain of nonlinear variables with nonlinear coupling (Fig. S.H21).
- •
Chain of nonlinear, periodic variables with linear coupling (Fig. S.H22).
- •
Chain of logistic maps with dynamical noise (Fig. S.H23).
- •
Chain of Henon maps (Fig. S.H24).
H.2 Average magnitude of for systems with unidirectional coupling
Here, we corroborate the statement that on average, systems with unidirectional coupling in the direction where dynamical coupling exists, and that in the direction where the is dynamical no coupling. We illustrate this by heat maps of average across ensembles of realizations of different systems, as a function of time series length and coupling strength (Figs. S.G10, S.G11, S.G12, S.G13, S.G14, S.G15).
- •
Chain of periodic autoregressive variables with strongly nonlinear coupling (Fig. S.G10).
- •
Chain of nonlinear variables with linear coupling (Fig. S.G11).
- •
Chain of nonlinear variables with nonlinear coupling (Fig. S.G12).
- •
Chain of nonlinear, periodic variables with linear coupling (Fig. S.G13).
- •
Chain of logistic maps with dynamical noise (Fig. S.G14).
- •
Chain of Henon maps (Fig. S.G15).
H.3 Average magnitude of for systems with bidirectional coupling
Here, we demonstrate the median magnitude of over varying and , for fixed time series length for the bidirectional logistic map system from the main text (Fig. S.H27) and a bidirectional nonlinear system with periodicity and dynamical noise (Fig. S.H28). We find that for the bidirectional logistic maps, the ensemble median is positive in both directions, thus capturing the underlying bidirectional coupling, even for short time series (here 300 observations). For the second bidirectional nonlinear system, the system appears bidirectional if coupling strengths are roughly equal. However, if the relative coupling strengths are different, then the system appears unidirectional (positive in the direction of the strongest forcing, and negative in the direction of the weakest forcing).
Appendix I Application to real datasets
In this appendix, we apply the ensemble sub-sampling approach described in section IV in the main manuscript to infer interaction networks from time series where ground truths are known. We analyze the following pairs of time series from the cause-effect pair database of Mooji et al.’s cause-effect pair database Mooij et al. 2016 (https://webdav.tuebingen.mpg.de/cause-effect/, accessed January 15th, 2020).
- •
Dataset 1: Solar radiation ( vs average air temperature () at the same location in Furtwangen, Black Forest, Germany between January 1, 1985 and December 31, 2008 (Fig. S.I29). This is pair 0077 of the cause-effect pair database Mooij et al. 2016. Ground truth: .
- •
Dataset 2: Inside room temperature () vs. outside temperature () (Fig. S.I30). This is pair 0069 of the cause-effect pair database Mooij et al. 2016. Ground truth: Outside temperature inside temperature.
- •
Dataset 3: Average precipitation () vs. average runoff () for 438 river catchments in the US. This is pair 0093 of the cause-effect pair database Mooij et al. 2016. Ground truth: Precipitation run-off.
- •
Dataset 4: Sunspot area vs. global temperature anomalies (deviations from 1961-1990) (Fig. S.I32). This is pair 0072 of the cause-effect pair database Mooij et al. 2016. Ground truth: unclear. If there is any coupling, it is from sunspot area to global temperature.
Our test correct infers a statistically significant causal influence in the correct direction for all these datasets. One exception occurs for the sunspot-temperature data, where neither direction is significant. This may be because the putative coupling is too weak to detect with so little data, or because there is no coupling at all. Nevertheless, the dynamical coupling between sunspots and temperature on Earth is disputed.
Appendix J Predictive asymmetry analysis of Late Pleistocene paleoclimate records
Here we repeat the analysis of the causal interactions among key climate system components in the Late Pleistocene using a different sea level (ice volume) record that is chronologically independent of orbital parameters (Grant et al. 2014).
References
- Granger (1969) C. W. Granger, Investigating causal relations by econometric models and cross-spectral methods, Econometrica: Journal of the Econometric Society , 424 (1969).
- Chen et al. (2004) Y. Chen, G. Rangarajan, J. Feng, and M. Ding, Analyzing multiple nonlinear time series with extended granger causality, Physics Letters A 324, 26 (2004).
- Marinazzo et al. (2008) D. Marinazzo, M. Pellicoro, and S. Stramaglia, Kernel method for nonlinear granger causality, Physical Review Letters 100, 144103 (2008).
- Schreiber (2000) T. Schreiber, Measuring information transfer, Physical Review Letters 10.1103/PhysRevLett.85.461 (2000), arXiv:0001042v1 [nlin] .
- Paluš and Vejmelka (2007) M. Paluš and M. Vejmelka, Directionality of coupling from bivariate time series: How to avoid false causalities and missed connections, Physical Review E 75, 056211 (2007).
- Rulkov et al. (1995) N. F. Rulkov, M. M. Sushchik, L. S. Tsimring, and H. D. Abarbanel, Generalized synchronization of chaos in directionally coupled chaotic systems, Physical Review E 51, 980 (1995).
- Schiff et al. (1996) S. J. Schiff, P. So, T. Chang, R. E. Burke, and T. Sauer, Detecting dynamical interdependence and generalized synchrony through mutual prediction in a neural ensemble, Physical Review E 54, 6708 (1996).
- Arnhold et al. (1999) J. Arnhold, P. Grassberger, K. Lehnertz, and C. E. Elger, A robust method for detecting interdependences: application to intracranially recorded eeg, Physica D: Nonlinear Phenomena 134, 419 (1999).
- Quiroga et al. (2002) R. Q. Quiroga, A. Kraskov, T. Kreuz, and P. Grassberger, Performance of different synchronization measures in real data: a case study on electroencephalographic signals, Physical Review E 65, 041903 (2002).
- Chicharro and Andrzejak (2009) D. Chicharro and R. G. Andrzejak, Reliable detection of directional couplings using rank statistics, Physical Review E 80, 026217 (2009).
- Sugihara et al. (2012) G. Sugihara, R. May, H. Ye, C. H. Hsieh, E. Deyle, M. Fogarty, and S. Munch, Detecting causality in complex ecosystems, Science 10.1126/science.1227079 (2012), arXiv:0711.2729 .
- Wiesenfeldt et al. (2001) M. Wiesenfeldt, U. Parlitz, and W. Lauterborn, Mixed state analysis of multivariate time series, International Journal of Bifurcation and Chaos 11, 2217 (2001).
- Feldmann and Bhattacharya (2004) U. Feldmann and J. Bhattacharya, Predictability improvement as an asymmetrical measure of interdependence in bivariate time series, International Journal of Bifurcation and Chaos 14, 505 (2004).
- Krakovská and Hanzely (2016) A. Krakovská and F. Hanzely, Testing for causality in reconstructed state spaces by an optimized mixed prediction method, Physical Review E 94, 052203 (2016).
- Liang (2013) X. Liang, The liang-kleeman information flow: Theory and applications, Entropy 15, 327 (2013).
- McCracken and Weigel (2016) J. M. McCracken and R. S. Weigel, Nonparametric causal inference for bivariate time series, Physical Review E 93, 022207 (2016).
- Ye et al. (2015) H. Ye, E. R. Deyle, L. J. Gilarranz, and G. Sugihara, Distinguishing time-delayed causal interactions using convergent cross mapping, Scientific Reports 10.1038/srep14750 (2015).
- Krakovská et al. (2018) A. Krakovská, J. Jakubík, M. Chvosteková, D. Coufal, N. Jajcay, and M. Paluš, Comparison of six methods for the detection of causality in a bivariate time series, Physical Review E 97, 042207 (2018).
- Smirnov (2013) D. A. Smirnov, Spurious causalities with transfer entropy, Physical Review E 87, 042917 (2013).
- Theiler et al. (1992) J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, and J. Doyne Farmer, Testing for nonlinearity in time series: the method of surrogate data, Physica D: Nonlinear Phenomena 10.1016/0167-2789(92)90102-S (1992), arXiv:9909037 [chao-dyn] .
- Lancaster et al. (2018) G. Lancaster, D. Iatsenko, A. Pidde, V. Ticcinelli, and A. Stefanovska, Surrogate data for hypothesis testing of physical systems, Physics Reports 10.1016/j.physrep.2018.06.001 (2018).
- James et al. (2016) R. G. James, N. Barnett, and J. P. Crutchfield, Information flows? a critique of transfer entropies, Physical Review Letters 116, 238701 (2016).
- Hahs and Pethel (2013) D. W. Hahs and S. D. Pethel, Transfer entropy for coupled autoregressive processes, Entropy 15, 767 (2013).
- Barnett et al. (2009) L. Barnett, A. B. Barrett, and A. K. Seth, Granger causality and transfer entropy are equivalent for gaussian variables, Phys. Rev. Lett. 103, 238701 (2009).
- Hannisdal (2011) B. Hannisdal, Non-parametric inference of causal interactions from geological records, American Journal of Science 311, 315 (2011).
- Diego et al. (2018) D. Diego, K. A. Haaga, and B. Hannisdal, Transfer entropy computation using the Perron-Frobenius operator, (2018), arXiv:1811.01677 .
- Matthews (1975) B. W. Matthews, Comparison of the predicted and observed secondary structure of t4 phage lysozyme, Biochimica et Biophysica Acta (BBA)-Protein Structure 405, 442 (1975).
- Chicco (2017) D. Chicco, Ten quick tips for machine learning in computational biology, BioData mining 10, 35 (2017).
- Haaga (2019) K. Haaga, Uncertaindata. jl: a julia package for working with measurements and datasets with uncertainties., Journal of Open Source Software 4, 1666 (2019).
- Esmark (1824) J. Esmark, Bidrag til vor jordklodes historie, Magazin for Naturvidenskaberne 2, 28 (1824).
- Croll (1875) J. Croll, Climate and time, Nature 12, 329 (1875).
- Milankovitch (1941) M. Milankovitch, Kanon der Erdebestrahlung und seine Anwendung auf das Eiszeitenproblem (Königlich Serbische Akademie, 1941).
- Hays et al. (1976) J. D. Hays, J. Imbrie, and N. J. Shackleton, Variations in the Earth’s orbit: pacemaker of the ice ages, Science 194, 1121 (1976).
- Pisias and Moore Jr (1981) N. G. Pisias and T. Moore Jr, The evolution of pleistocene climate: a time series approach, Earth and Planetary Science Letters 52, 450 (1981).
- Raymo (1997) M. Raymo, The timing of major climate terminations, Paleoceanography 12, 577 (1997).
- Tzedakis et al. (2017) P. Tzedakis, M. Crucifix, T. Mitsui, and E. W. Wolff, A simple rule to determine which insolation cycles lead to interglacials, Nature 542, 427 (2017).
- Denton et al. (2010) G. H. Denton, R. F. Anderson, J. Toggweiler, R. Edwards, J. Schaefer, and A. Putnam, The last glacial termination, science 328, 1652 (2010).
- Wolff et al. (2009) E. Wolff, H. Fischer, and R. Röthlisberger, Glacial terminations as southern warmings without northern control, Nature Geoscience 2, 206 (2009).
- Petit et al. (1999) J.-R. Petit, J. Jouzel, D. Raynaud, N. I. Barkov, J.-M. Barnola, I. Basile, M. Bender, J. Chappellaz, M. Davis, G. Delaygue, et al., Climate and atmospheric history of the past 420,000 years from the vostok ice core, antarctica, Nature 399, 429 (1999).
- Shackleton (2000) N. J. Shackleton, The 100,000-year ice-age cycle identified and found to lag temperature, carbon dioxide, and orbital eccentricity, Science 289, 1897 (2000).
- Shakun et al. (2012) J. D. Shakun, P. U. Clark, F. He, S. A. Marcott, A. C. Mix, Z. Liu, B. Otto-Bliesner, A. Schmittner, and E. Bard, Global warming preceded by increasing carbon dioxide concentrations during the last deglaciation, Nature 484, 49 (2012).
- Abe-Ouchi et al. (2013) A. Abe-Ouchi, F. Saito, K. Kawamura, M. E. Raymo, J. Okuno, K. Takahashi, and H. Blatter, Insolation-driven 100,000-year glacial cycles and hysteresis of ice-sheet volume, Nature 500, 190 (2013).
- Laskar et al. (2004) J. Laskar, P. Robutel, F. Joutel, M. Gastineau, A. Correia, and B. Levrard, A long-term numerical solution for the insolation quantities of the earth, Astronomy & Astrophysics 428, 261 (2004).
- Huybers and Wunsch (2005) P. Huybers and C. Wunsch, Obliquity pacing of the late Pleistocene glacial terminations, Nature 434, 491 (2005).
- Haaga et al. (2018) K. A. Haaga, J. Brendryen, D. Diego, and B. Hannisdal, Forcing of late Pleistocene ice volume by spatially variable summer energy, Scientific Reports 10.1038/s41598-018-29916-3 (2018).
- Bereiter et al. (2015) B. Bereiter, S. Eggleston, J. Schmitt, C. Nehrbass-Ahles, T. F. Stocker, H. Fischer, S. Kipfstuhl, and J. Chappellaz, Revision of the epica dome c co2 record from 800 to 600 kyr before present, Geophysical Research Letters 42, 542 (2015).
- Bazin et al. (2013) L. Bazin, A. Landais, B. Lemieux-Dudon, H. Toyé Mahamadou Kele, D. Veres, F. Parrenin, P. Martinerie, C. Ritz, E. Capron, V. Lipenkov, et al., An optimized multi-proxy, multi-site antarctic ice and gas orbital chronology (aicc2012): 120-800 ka, Climate of the Past 9, 1715– (2013).
- Veres et al. (2013) D. Veres, L. Bazin, A. Landais, H. Toyé Mahamadou Kele, B. Lemieux-Dudon, F. Parrenin, P. Martinerie, E. Blayo, T. Blunier, E. Capron, et al., The antarctic ice core chronology (aicc2012): an optimized multi-parameter and multi-site dating approach for the last 120 thousand years, Climate of the Past 9, 1733 (2013).
- Spratt and Lisiecki (2016) R. M. Spratt and L. E. Lisiecki, A Late Pleistocene sea level stack, Climate of the Past 12, 1079 (2016).
- Grant et al. (2014) K. Grant, E. Rohling, C. B. Ramsey, H. Cheng, R. Edwards, F. Florindo, D. Heslop, F. Marra, A. Roberts, M. E. Tamisiea, et al., Sea-level variability over five glacial cycles, Nature communications 5, 1 (2014).
- Bossomaier et al. (2016) T. Bossomaier, L. Barnett, M. Harré, and J. T. Lizier, An introduction to transfer entropy, Cham, Germany: Springer International Publishing. Crossref (2016).
- (52) measures a TE-like quantity, except the causal direction flips and the conditioning in the argument of the logarithm occurs only on one variable, not on a mixture of the two variables, as for regular TE (eqs. S.D137 and S.D138a).
- Sauer et al. (1991) T. Sauer, J. A. Yorke, and M. Casdagli, Embedology, Journal of Statistical Physics 65, 579 (1991).
- Deyle and Sugihara (2011) E. R. Deyle and G. Sugihara, Generalized theorems for nonlinear state space reconstruction, PLoS One 6, e18295 (2011).
- Kantz and Schreiber (2004) H. Kantz and T. Schreiber, Nonlinear time series analysis, Vol. 7 (Cambridge university press, 2004).
- Cover and Thomas (2012) T. M. Cover and J. A. Thomas, Elements of information theory (John Wiley & Sons, 2012).
- Kraskov et al. (2004) A. Kraskov, H. Stögbauer, and P. Grassberger, Estimating mutual information, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 10.1103/PhysRevE.69.066138 (2004), arXiv:0305641 [cond-mat] .
- Péguin-Feissolle and Teräsvirta (1999) A. Péguin-Feissolle and T. Teräsvirta, A General Framework for Testing the Granger Noncausality Hypothesis (Universites d’Aix-Marseille II et III, 1999).
- Chávez et al. (2003) M. Chávez, J. Martinerie, and M. Le Van Quyen, Statistical assessment of nonlinear causality: application to epileptic EEG signals, Journal of Neuroscience Methods 124, 113 (2003).
- Krakovská et al. (2015) A. Krakovská, K. Mezeiová, and H. Budáčová, Use of false nearest neighbours for selecting variables and embedding parameters for state space reconstruction, Journal of Complex Systems 2015 (2015).
- Mooij et al. (2016) J. M. Mooij, J. Peters, D. Janzing, J. Zscheischler, and B. Schölkopf, Distinguishing cause from effect using observational data: methods and benchmarks, The Journal of Machine Learning Research 17, 1103 (2016).