Causality indices for bivariate time series data: A comparative review of performance Open Access

No general method exists to identify causal structures within complex systems, and there is no single consistent and unifying notion of quantitative causality estimation for time series data. Published methods can be broadly categorized into the following groups:

1. regression-based indices that use “recent history” vectors as predictors in a model (e.g., Granger causality),

2. information-theoretic indices that build upon ideas of conditional mutual information (e.g., transfer entropy),

3. indices based on state space dynamics, such as local neighborhoods and trajectories (e.g., convergent cross mapping), and

4. graphical models that scale causal inference estimation to high-dimensional multivariate systems for causal identification.

There exist common themes between these methods, and membership within these groups is sometimes somewhat blurred. Figure 1 sets out key properties and similarities between methods from groups 1–3. Previous reviews of the literature^14–17 typically focus on a subset of methods from one of these groups. The suitability, interchangeability, and performance of published methods, particularly where they span more than one of these groups, has received relatively little attention. In this work, we identified and assessed a widely used subset of indices for directed bivariate causality inference, concentrating on methods involving univariate embeddings to describe the recent history of the system (Sec. II). A review of such methods has been published previously by Lungarella et al.¹⁸ We reproduce these results for the original and newer methods. We also extend this work by proposing a set of modifications that can be made to simulated data prior to causal estimation in order to investigate sensitivity of each method to data availability, scaling, missing data, rounding, and Gaussian noise (Sec. III). Each of these reproduces phenomena that often occur in real-world data, such as when instruments have a fixed measurement precision and data are reported with rounding errors. We believe that these tests should provide in-depth benchmarking criteria for new proposed methodologies. Finally, we summarize the strengths and weaknesses of these approaches and identify key areas of further research (Sec. IV).

We observe a complex system as a set of variables within a multivariate time series. The time series $X=(x1,…,xT)$ and $Y=(y1,…,yT)$ describe a bivariate system with state $st=(xt,yt)$ at time $t$⁠. A critical implicit assumption is that the series has unit time or equivalently that the data are observed at a fixed constant frequency. Underpinning all methods is the key assumption that causes directly precedes effect. As a result, a preliminary step is the construction of time-delay embedding vectors $xtm,τ$ in the $m$-dimensional state space $X≅Rm$⁠. An important distinction in the defined methods is whether their theoretical basis is stochastic or deterministic. Both employ the same time-delay vectors, though only in deterministic systems are the vectors $xtm,τ$ considered elements within an $m$-dimensional state space $X≅Rm$⁠. We construct equivalent embedding vectors $yt$ (and state space $Y≅Rm$⁠) for $Y$ and joint embedding vectors $zt$ (and state space $Z≅R2m$⁠),

In practice, many real-world systems are stochastic, with some level of noise or randomness in at least part of the system. A further assumption for stochastic causality estimation is that of separability, which states that there is unique information about the effect variable that is contained only within the causal variable. The standard approach here is to describe or model the current value $xt$ of $X$ as conditional upon the “recent history” joint embedding vector $zt−1$ (full model). Separability means that removing the causal variable $Y$ eliminates the information it contains about the effect $X$⁠, which we observe either by identifying non-zero coefficients in the full model or constructing a reduced model, conditioned only upon $xt−1$⁠. These methods are generally described with time index shifted $t↦t+1$⁠, though the interpretation (“current” and “recent history”) remains the same.

Granger causality¹ (GC), a notable and popular method for causality estimation in time series, fits autoregressive models on the time series to this end. Extensions of GC to nonlinear systems include a locally linear version called extended Granger causality¹⁹ (EGC) and nonlinear Granger causality²⁰ (NLGC), which perform a “global” nonlinear autoregression using radial basis functions (RBFs). Predictability improvement²¹ (PI) is another locally constant linear regression of “recent history” embeddings, which measures a reduction in the mean squared error when $zt$ is used for predicting a “horizon” value $xt+h$ instead of $xt$ alone.

Information theory is a natural framework for describing causal relationships. Transfer entropy⁷ (TE) measures deviation from the generalized Markov property $p(xt+1|xt)=p(xt+1|xt,yt)$ as a conditional mutual information. With weak coupling and limited data, transfer entropy can suffer from finite sample effects, and effective transfer entropy²² (ETE) corrects for this using shuffled versions of the causal variable. TE reduces to vanilla GC under the assumption of Gaussian variables²³ (⁠$⋆$⁠, Fig. 1), and non-zero GC implies violation of the generalized Markov property and non-zero TE.²⁴ A coarse-grained transinformation rate²⁵ (CTIR) is based upon “coarse-grained entropy rates” and measures the rate of net information flow, averaged over different lags $τ$⁠. Often, the difficulty in information-theoretic methods (described in depth in Ref. 14) is the robust estimation of joint probabilities or entropy values, which in turn form building blocks for these methods. We use a histogram binning partition (H) and the (hypercube) Kraskov–Stögbauer–Grassberger (KSG) estimate,²⁶ which is a technique involving $k$-nearest neighbor statistics. All information-theoretic computation here is in “nats” (logarithm base $e$⁠).

Fully deterministic dynamical systems, which evolve according to a differential equation or a difference equation, do not necessarily satisfy the separability condition. In these systems, $xt$ can often be reformulated as a function of only past values of $X$⁠, which makes the potential causal role of $Y$ in the coupled system less clear, as highlighted by Granger.¹ Causal relationships in a coupled deterministic system are instead observed via the event that each variable belongs to a shared attractor manifold $A⊂Z$⁠. A consequence of Takens’ embedding theorem²⁷ is that the “library of historical behavior” of $X$ preserves the topology of $A$ and, by transitivity, local neighborhoods in $X$ those in $Y$ and vice versa.⁵ It is possible to detect unidirectional causal influence, where only the dynamics of a causal variable propagate to the response variable in this way. Sugihara et al.⁵ argue that the inferred direction of unidirectional causal influence is counter-intuitively reversed (i.e., cross mapping from $X$ to $Y$ reveals causal influence from $Y$ to $X$⁠).

The key assumption of cross mapped indices is that causal relationships are observed in the similarity between sets of (subscript) indices denoting the nearest neighbors for each set of embedding vectors, which can be mapped from one variable to the other to reveal interdependency. This is the idea behind the similarity indices: two similarity indices we test here, denoted $SIY→X(1)$ and $SIY→X(2)$⁠, are $H(X|Y)$ in Ref. 28 and $E(X|Y)$ in Ref. 29, respectively. Convergent cross mapping⁵ (CCM) computes the correlation $ρ$ between the cross mapped estimate and the true value, with convergence in $ρ$ as $T$ increases “a key property that distinguishes causation from simple correlation.”⁵ 

Our results are split into two parts. First, we reproduce the results from Ref. 18, evaluating the performance of all methods including the additional CTIR and CCM, plus ETE using histogram binning and TE using KSG. In these simulations, we choose the same simulation model parameters and causality index parameters as in Ref. 18 (Table I). In the second part, we investigate sensitivity to common issues relevant to real-world data using the Ulam lattice system to illustrate these.

We investigate performance on four simulated model systems (Table II). In each simulation, we assess the causality estimates of each method by varying the coupling strength $λ$⁠. These simulated systems are widely studied in chaos theory, e.g., Ref. 37, and also appear elsewhere in the literature, e.g., Ulam lattice in Ref. 7.

Linear process:

Ulam lattice:

Hénon unidirectional map:

Hénon bidirectional map:

We reproduce figures for all simulations and methods in Figs. 3–5 and summarize our results in Fig. 2, which shows correlations between each pair of indices. For linear process and Ulam lattice simulations, we report causality estimates in both directions, i.e., $iX→Y$ and $iY→X$ (where $i$ denotes any of the causality indices). For Hénon maps, we instead use the directed index $DX→Y=iX→Y−iY→X$⁠, following Ref. 18. In general, all measures exhibit a small standard deviation relative to the absolute value of the index, indicating that random initial conditions during data simulation have at most a very small influence on the causal structure when initial transients are discarded. Though we are able to replicate the findings in Ref. 18 in most cases, we occasionally find minor differences between their results and ours. In particular, we sometimes find results of a similar profile but a different magnitude, as $λ$ varies. We observe this for EGC and linear processes, PI and all simulations, SI$(1)$ and Ulam lattice, and TE and Hénon unidirectional maps. Though we handle numerical outliers differently in our visualization of results for Hénon bidirectional maps, our results for these simulations appear largely comparable in magnitude and profile. There is no mention of a data standardization step in Ref. 18 and the results we report do not involve any pre-processing, though this did not appear to rectify these differences. In one notable difference between identical Hénon bidirectional map results, Lungarella et al.¹⁸ find a region in the $λ$-space (namely, ${(λxy,λyx):λxy>0.1,λyx>0.1,λxy+λyx<0.35}$⁠) in which they identify general synchronization between $X$ and $Y$ and have difficulty estimating indices due to numerical instabilities; yet, we do not observe this. We have followed the implementation in Ref. 18 as closely as possible, and it is unclear why these differences exist.

We knowingly deviated from the implementations in Ref. 18 only in the case of NLGC, in which we preferred to use $k$-means rather than fuzzy $c$-means clustering to determine RBF centers, after finding similar or improved results but with a much reduced computational cost. Lungarella et al.¹⁸ note that NLGC is numerically unstable for “small” $T$ and computationally expensive for “large” $T$⁠, which we suggest may be partly due to their use of fuzzy $c$-means. Little detail is provided about their implementation of this, but it may perhaps be that an early stopping criteria sometimes forces a “poor quality” clustering. Furthermore, we found that the performance of NLGC in Hénon bidirectional map simulations improved significantly with a different set of NLGC parameters (e.g., $P=50$ instead of $P=10$⁠), though we do not present these alternate results.

For the linear process (LP), the simplest simulation model, all indices show a very strong positive correlation in the $Y→X$ direction (Fig. 3). In the reverse $X→Y$ direction, TE and CTIR both decrease with increasing $λ$⁠, the cross mapped indices all show a marked increase, and the remaining indices are approximately zero for all $λ$⁠. This gives rise to patterns of positive and negative correlations between pairs of methods. As each $xt$ or $yt$ is a sum of Gaussian variables, we can derive theoretical values for Shannon entropy and, consequently, for the information theory methods (see the supplementary material). Figure 3 shows that TE (KSG) reliably estimates the ‘true” transfer entropy, but TE (H) significantly underestimates the theoretical values. This is a fundamental flaw that undermines any other advantageous properties of TE (H). Though computed CTIR values match the theoretical values, it is clearly negative in $X→Y$ and as such does not reflect the causal structure of the system. Increasing the size of the data $T$ alters the value of TE (H) here, but TE (KSG) remains accurate as $T$ increases. However, in all other simulations, TE (H) is more robust to increasing data size.

The Ulam lattice (UL) chains together unidirectional coupled chaotic Ulam maps. For large $NL$⁠, the causal influence from $Y$ to $X$ is negligible. UL exhibits synchronization for $λ≈0.18,0.82$⁠, where cause and effect variables are indistinguishable from each other; e.g., the system converges to a two state attractor. As a result, most indices either have values approximately equal to zero or suffer from high variance numerical instabilities. Outside of these regions of synchronization, the information-theoretic methods and regression-based indices show reasonable consistency (Figs. 2 and 4). The exception is CTIR, which slowly decreases as $λ$ increases for $T=105$⁠, albeit still correctly identifying the direction of information flow. ETE (H) successfully corrects for the small sample effects that give rise to these spurious positive TE (H) results in the $Y→X$ direction when $T=103$⁠. Both similarity indices (SIs) fail to identify any causal structure in the UL. For CCM, while the net directed index $DX→Y$ increases with $λ$⁠, it is negative for $λ<0.5$ and so misidentifies the direction of causality. The positive correlations between methods in $iY→X$ for $T=105$ occur due to a very slight peak in the value at $λ≈0.5$ for nearly all methods (except CTIR and EGC).

Synchronization occurs in the range $λ∈[0.7,1]$ for Hénon unidirectional maps (Fig. 4). All indices are consistent and perform reasonably well, and we do not observe the noisy fluctuations seen in Ref. 18 apart from for EGC when $λ>0.7$⁠. It is notable that about half the indices [TE (H), ETE (H), EGC, NLGC, PI] are fairly consistent even as the order of magnitude of the data size $T$ is increased, while the values of the other indices [TE (KSG), CTIR, SIs, CCM] increase, in some cases tripling in value. There is a strong degree of similarity between all methods for the Hénon bidirectional maps (Fig. 5). In HB (I) simulations, the exceptions to this are NLGC and CCM, though we found that setting a larger number of RBFs resulted in better performance for NLGC. For CCM, the direction of causality is sometimes incorrect and the reason for this is unclear, though this may also be a result of poor parameter choices. We observe expected symmetry in the values of $λ$⁠, and synchronization in the region is approximately equal to ${(λxy,λyx):λxy+λyx>0.28}$⁠. There are a small number of points in which numerical instabilities are present in all indices, but the consistency across all indices suggests that these are isolated points in which the system converges to some limit cycle or attractor. There are more differences between methods in HB (NI) results. Significant numerical instabilities occur in EGC, particularly when the system is in a state of synchrony: the region approximately equal to ${(λxy,λyx):0.05<λxy<0.15,0.1<λyx<0.28}$⁠. Outside of this region, EGC is broadly similar to the information-theoretic indices, which are highly correlated, and to a slightly lesser degree with the SI$(1)$ and PI. In contrast, NLGC, SI$(2)$⁠, and CCM have unusual results. The first of these is negative almost everywhere (even in a repeat analysis with more RBF kernels) and the latter two are mostly non-negative; moreover, the regions with the most extreme values occur in quite different places in all three.

An important consideration in selecting a suitable method is any trade-off between performance and computational efficiency. The most significant factor in this is often how the algorithmic cost of each method scales with increasing data size $T$⁠. Table S.II in the supplementary material shows the mean and standard deviations of the time taken to compute each index and simulation. TE (H)/ETE (H) is the fastest in almost all cases, even though this calculation also includes ten reshuffles and recomputations for ETE (H). CCM, EGC, and TE (KSG) are similarly among methods with a smaller computational cost. Several methods have extreme values in UL simulations with $T=105$⁠, particularly CTIR and PI, but this is distorted by difficulties in computation when the system is in a synchronized state. We observed a marked difference in the computational cost for NLGC when using $k$-means for clustering instead of fuzzy $c$-means, and we suggest that $k$-means is more suitable here. Our computation was done in a high performance CPU computing cluster using SkyLake 6140 with 18 core 2.3 GHz processors and 384 GB of RAM. Although the computational times we report may be slightly faster than on a laptop computer with less processing power, we did not observe any substantial difference when internally comparing run times.

Next, we investigated the sensitivity of all causality indices to a number of modifications mirroring issues that often arise in real-world data. We choose UL with $T=103$ to illustrate the effects of these transformations, as LP is too simplistic a model to give sufficient insight. We keep all simulation parameters and causality index parameters the same. In Table III, we summarize the means $μ^$ and standard deviations $σ^$ of the directed indices $DX→Y$ for each $λ$⁠, which are normalized by their deviation from the “base” UL results and averaged over all $λ$⁠. This normalization allows us to compare across methods that take values in different ranges.

Many of the indices, with the exception of TE (KSG) and CTIR, remain consistent with increasing data size $T$⁠, while at the same time exhibit decreasing variance. These results reinforce the similar observations in HU maps. The large increases in the value of the two methods mentioned are concerning and represent a drawback of both methods that should be acknowledged in applications of these approaches. It is unclear whether there is convergence to some “correct” value as the amount of data increases or whether both are unbounded as $T→∞$⁠, but initial computations do not support the former (not shown). Though the $DX→Y$ values from both transfer entropy methods are highly correlated, they are both estimates of the same quantity, and it is difficult to reconcile their different magnitudes, particularly as we have already seen significant underestimation in TE (H) for LP simulations.

In the second set of experiments, we perform three tests: standardizing both series by their sample mean and standard deviation in the first and separately scaling each unstandardized time series by a factor of $10$ (Fig. S2 in the supplementary material). For the Ulam lattice system, sample means for both $X$ and $Y$ are typically between 0.4 and 0.7, and standard deviations are both approximately equal to 1.2 (except when the system is in synchrony). Several methods are invariant under linear scaling or shifting of the original series $X$ and $Y$⁠, including cross mapping approaches. Information-theoretic measures are also invariant in theory, but the KSG algorithm, based on $k$-nearest neighbors, does not retain this property. Similarly, EGC relies on a neighborhood size parameter, and scaling the data without changing this parameter accordingly can result in insufficient points available for the locally linear regressions, as is observed when either $X$ or $Y$ is scaled by 10. The directed index for both NLGC and PI has a vastly inflated magnitude when $Y$ is scaled by 10. With this in mind, we recommend standardization or normalization of the data before employing these methods.

We perform three experiments to investigate rounding errors, first rounding each time series separately to one decimal place and then rounding both to two decimal places (Fig. S3 in the supplementary material). TE and both GC extensions have similar performance to the baseline in all cases, while CCM suffers the most. In two experiments with missing data of 10% and 20%, all methods appear robust to this.

In the case of the earlier LP simulations, Gaussian noise forms an integral component of the system itself, and the theoretical expression for TE shows that this depends only on the ratio of the variances $σx/σy$ (see the supplementary material). However, this noise is inherent in the simulation process (i.e., it does not arise in observation of the system). In our UL experiments, we added additional Gaussian noise after simulation. The inclusion of this “observation” noise does not alter the state of the system or the information flow between variables, but it does obscure the causal structure. In the first of these experiments (Fig. S3 in the supplementary material), in which we added small variance Gaussian noise (⁠$σG=0.1$⁠), the amplitude of this noise is an order of magnitude less than the original UL values and the inclusion of this noise has a small effect for all indices. In the latter experiments, we added Gaussian noise (with $σG=1$⁠) to each variable individually, and the effect is more pronounced. NLGC performs best in general and appears very resilient to noise added to $Y$ (effect variable), though it drops slightly in value when Gaussian noise is added to $X$ (cause variable). It is interesting that the two SIs have quite different results, with SI$(1)$ more robust to noise, though both methods are not in general able to successfully identify the direction of causality.

In-depth comparative studies of this kind are relatively rare in the mathematical literature (examples include Refs. 16 and 38), particularly in evaluating performance of methods for estimating a concept, such as causality, that does not have a consistent, fundamental mathematical definition. Even without this, causal inference has a huge importance in how we can model, predict, and exploit real-world applications from many scientific disciplines. Asymmetric bivariate causal inference is the first key step to providing this insight into interactions between components in complex networks. In the context of causality indices, review papers^14–16 have previously had a narrower focus in some manner, for example, on only one group of methods or on a few bivariate methods and their multivariate extensions. We follow the template of Ref. 18 in reviewing methods drawn from diverse mathematical foundations, but we extend this review with additional methods, and crucially, we investigate the impact of common issues that are relevant to real-world data. In reproducing and updating their work, we are also able to resolve some computational stability issues and comment on the computational costs of each method, while we also make our code publicly available for other researchers to develop further.

A primary concern in causality inference is the difficulties with model misspecification, specifically causal identification in multivariate systems. Omission of confounding variables can create spurious false-positive causal relationships. There may also be redundancy across multiple variables that provide similar information to the effect variable or sets of variables that interact synergistically such that their combined causal influence is greater than the “sum of their parts.” These are key concerns outlined in Ref. 15; consequently, results from bivariate indices cannot be definitively interpreted as the existence of a fundamental direct causal relationship between two variables.²⁸ A key avenue for further work is to advance this analysis beyond a bivariate setting by including possible confounding variables, in line with conditional extensions to Granger causality^8,19,39,40 and transfer entropy.⁴¹ Recent work with graphical models of multivariate systems² is an important step toward high-dimensional causal identification.

Separate univariate embedding is not without some limitations and is not necessarily the optimal multivariate embedding. Aside from in Ref. 42, mixed embeddings are as yet uncommon in causality estimation. There is not yet a theoretical framework for longitudinal data that is recorded non-simultaneously and irregularly. Often, a typical workflow for such data involves pre-processing to transform the data into a multivariate series with constant time intervals. However, many imputation methods result in significant and poorly quantified biases in information content and flow, which inevitably propagate through to estimates of causality, and more work is needed to explicitly factor this into a causal inference framework.

Each causality index has strengths and weaknesses, and there is no single method whose all-round performance exceeds all others. Transfer entropy and Granger causality have long been regarded as the leading methods for systems that contain a small number of variables, and these have had wide applications.^8,13 Transfer entropy has the distinct advantage that it is built upon the principles of Shannon entropy in a well-established and universal information-theoretic framework. It performs solidly throughout, though there is some tension between algorithms for TE, with the estimates rarely in complete agreement. We have shown that a histogram fixed partition approach is biased even in the simplest model despite TE (H) having general consistency and computational efficiency. Therefore, we recommend the KSG algorithm for transfer entropy computation unless perhaps data are extremely scarce (⁠$T<103$⁠). However, there are some unanswered concerns about TE (KSG), particularly that it appears to increase in magnitude as more data are available. TE (KSG) also suffers in performance when data are unequally scaled due to the resultant difficulties with identifying unique nearest neighbors. CTIR, while sometimes not wholly dissimilar in value from TE, did not seem to offer any obvious advantage to compensate for its much higher computational cost or occasional unusual behavior. Vanilla GC is widely favored but has restrictive assumptions and is ill-suited to complex nonlinear problems. Of the two nonlinear extensions to Granger causality, Lungarella et al.¹⁸ appear to prefer EGC. Some of the computational challenges and numerical instability that they experienced with NLGC may have been a result of their choice of a fuzzy $c$-means for determining RBF kernels, and alternate parameter choices appear to resolve some of their concerns. We find that NLGC is one of the most robust methods to rounding errors, missing data, and Gaussian noise. They rightly note that “If the rank of the data is small, kernel-based methods tend to overfit,”¹⁸ but we did not observe any issues with this in our simulation experiments. Predictability improvement (PI) likewise performed solidly and has a slight advantage among the regression-based indices in that it is perhaps less reliant on parameter choices. Finally, dynamical systems theory offers a different insight into causal inference that should not be readily dismissed despite our mixed results here, even though the deterministic simulation models appeared to be well-suited to the underlying theory. Convergent cross mapping is a more recent and popular method, and this offered a broad improvement on the similarity indices (SIs), which did not consistently identify the strength or direction of causality. However, CCM too did not always manage to determine the correct causal flow in our simulations.

We have highlighted the value of a standardization pre-processing step in avoiding algorithmic issues, which is also important in comparing results from different data for each method. Rounding errors give rise to practical issues within the implementation of several of the algorithms. For instance, in $k$-nearest neighbor approaches, it is typically assumed that the distances between pairs of points are unique and not discrete. Subsequent edge cases can be treated by adding random noise with low amplitude to the data before estimating the causal relationships,²⁶ though propagation of this noise to final estimates is something that should be analyzed. Likewise, many existing implementations of the methods are not equipped to handle missing data (e.g., Refs. 41 and 43). We believe that this is broadly straightforward to implement across all indices, as it can be handled exclusively within the time-delay embedding vector step, by performing embedding and then removing any embedding vectors missing at least one component. As Mønster et al.³⁶ put it, “Noise in real-world data is ubiquitous, the inclusion of noise in model investigations has been largely ignored.” Added Gaussian noise leads to the biggest changes in the value for most methods, particularly noise at observation in the causal variable. However, provided the magnitude of noise is small compared to the values themselves, all methods perform adequately.

On the basis of this work, we conclude that the strongest choice for identifying and quantifying bivariate causal relationships is, in our view, either transfer entropy (KSG) or nonlinear Granger causality. Predictability improvement is a reasonable alternative and perhaps the next best candidate. A more cautious approach may involve using more than one method, from different theoretical backgrounds. Where possible, it is advantageous to identify a base case for the system, where subsequent results can be reliably compared against. For new methodologies, we recommend investigation into the real-world issues that we have discussed.

See the supplementary material that contains additional tables and figures. Tables S.I and S.II show full parameter choices and computational time requirements of each method, respectively. Figures S1–S3 show the results of real-world relevant transformation experiments. The supplementary material also contains theoretical results for information-theoretic measures in the linear process simulation.

T.E. is funded by the Engineering and Physical Sciences Research Council (EPSRC) National Productivity Investment Fund (NPIF) No. EP/S515334/1, Reference No. 2089662. A CC BY or equivalent license is applied to the AAM arising from this submission. We would like to acknowledge Marcel Stimberg and Daniel Nüst for their work with the CODECHECK. We found several existing open-source code repositories, listed in the Data Availability section. It was insightful to view and test these packages, though we still decided to develop our own code for these methods. In addition, we appreciate advice from George Sugihara on the use of CCM (pyEDM) in email correspondence. We also acknowledge the Python community for core packages that this work depends upon, including IPython⁴⁷ v7.16.1, Matplotlib⁴⁸ v3.3.2, NumPy⁴⁹ v1.18.5, Palettable⁵⁰ v3.3.0, Pandas⁵¹ v1.0.5, Python⁵² v3.8.3, scikit-learn⁵³ v0.23.1, SciPy⁵⁴ v1.5.0, and statsmodels⁵⁵ v0.11.1.

The data that support the findings of this study are openly available in GitHub at https://github.com/tedinburgh/causality-review, Ref. 44. Our code is available at the same repository. A CODECHECK certificate is available confirming that the computations underlying this article could be independently executed: doi.org/10.5281/zenodo.4720843. The existing open-access code for some indices includes repositories for information theory and transfer entropy: IDTxl⁴¹ v1.1, PyIF,⁴⁵ and for convergent cross mapping: pyEDM⁴³ v1.7.4. We also adapted the fuzzy $c$-means code based on Ref. 46. We checked our results for transfer entropy and convergent cross mapping against those from IDTxl and pyEDM, respectively. All the codes in our repository and in these others are in Python.