Causal network reconstruction from time series: From theoretical assumptions to practical estimation

Reconstructing the causal relations behind the phenomena we observe is a fundamental problem in all fields of science. The traditional approach is to conduct active experiments, but in many fields such as Earth system science or neuroscience, manipulations of the complex system under study are either impossible, unethical, or very expensive. On the other hand, modern science generates an ever-growing amount of data from these systems, in particular time series data. Concurrently, novel computing hardware today allows efficient processing of massive amounts of data. These developments have led to emerging interest in the problem of reconstructing causal networks or causal discovery from observational time series.

In the past few decades, a number of original causality concepts have been developed, such as Granger causality (Granger, 1969) or transfer entropy (Schreiber, 2000). Since the 1990s, computer scientists, statisticians, and philosophers have grounded causal reasoning and inference in a robust mathematical framework (Pearl, 2000; Spirtes et al., 2000). The (quite natural) definition of causality underlying this framework is that $X→Y$ if and only if an intervention or manipulation in $X$ has an effect on $Y$ (Pearl, 2000; Spirtes et al., 2000). This effect may be in changing $Y$’s mean or any change in its post-interventional distribution denoted $P[Y|do(X=x)]$ which is different from the conditional distribution $P(Y|X=x)$⁠. Unfortunately, all we can measure from observational data are statistical dependencies. These can be visualized in a graphical model (Lauritzen, 1996) or time series graph (Eichler, 2011) that represents the conditional independence relations among the variables and their time lags (Fig. 1). The theory of causal discovery lays out the assumptions under which the underlying causal dependencies can be inferred from observational data.

There are different sets of assumptions that allow us to identify a causal graph. Here, we focus on time-lagged causal discovery in the framework of conditional independence testing using the assumptions of time-order, Causal Sufficiency, the Causal Markov Condition, and Faithfulness, among others, which are all discussed thoroughly in this paper. But some of these assumptions can be replaced. Recent work (Peters et al., 2017) shows ways to use assumptions on the noise structure and dependency types in the framework of structural causal models which can complement the approach studied here and we will include references to recent work from this framework throughout the sections.

The paper is organized as follows: In Sec. II, we relate Granger causality and similar concepts to the conditional independence-framework (Spirtes et al., 2000). Section III provides the necessary definitions and notation and in Sec. IV we recapitulate the assumptions underlying time-lagged causal discovery from time series alongside illustrative examples. The practical estimation aspect from introducing some causal discovery algorithms to significance testing is presented in Sec. V, while Sec. VI discusses suggestions for performance evaluation. In Sec. VII, we present some comparison studies of common causal methods and conclude the paper with a brief discussion (Sec. VIII). The paper is accompanied by a python jupyter notebook on https://github.com/jakobrunge/tigramite to reproduce some of the examples.

Granger (1969), based on work by Wiener (1956), was the first to propose a practical, operational definition of causality based on prediction improvement. The underlying idea of measuring whether $X$Granger causes $Y$ is that there is some unique information in $X$ relevant for $Y$ that is not contained in $Y$’s past as well as the past of “all the information in the universe” (Granger, 1969). In practice, typically only $Y$’s past is used (bivariate Granger causality). Measuring prediction improvement can be operationalized in different ways. The most common framework are vector autoregressive models (VAR),

where $Xt=(Xt1,…,XtN)$⁠, $Φ(τ)$ is the $N×N$ coefficient matrix at lag $τ$⁠, $τmax$ some maximum time lag, and $η$ denotes an independent noise term. Here, $Xi$Granger-causes$Xj$ if any of the coefficients $Φji(τ)$ at lag $τ$ is non-zero. A non-zero $Φji(τ)$ can then be denoted as a causal link $Xt−τi→Xtj$ at lag $τ$⁠. Another option is to compare the residual variances of the VAR fitted with and without including the variable $Xi$⁠. The use of VARs restricts this notion of causality to a causality in mean (Granger, 1969). A more general definition is that of (bivariate) transfer entropy (Schreiber, 2000; Barnett et al., 2009)

where $I(X;Y∣Z)$ denotes the conditional mutual information (CMI). Bivariate TE is a common term, another naming option would be bivariable TE since $X$ and $Y$ could also be multivariate variables. Transfer entropy can also be phrased in a multivariate (or multi-variable) lag-specific version (Runge et al., 2012a). Many current methods are advancements of the concept of transfer entropy (Wibral et al., 2013; Staniek and Lehnertz, 2008; Vejmelka and Palus, 2008), in particular in its multivariate version (Sun and Bollt, 2014; Sun et al., 2015; Runge et al., 2012b; 2012a; Runge, 2015).

Tests for causality are then based on testing whether a particular CMI is greater than zero. Looking at the definition of CMI,

TEbiv and its advancements essentially test for conditional independence of $X$ and $Y$ given $Z$⁠, denoted $X⊥⊥Y|Z$ since

$Z$ then represents $Y$’s past and other included variables. The lag-specific generalization of the VAR model (1) then is the full conditional independence (FullCI) approach

where $Xt(t−1,…,t−τmax)=(Xt−1,…,Xt−τmax)$⁠. $I$ can be CMI or any other conditional dependence measure. In the case of partial correlation, a non-zero entry $Φji(τ)$ corresponds to a non-zero $Ii→jFullCI(τ)$⁠. A general concept to represent conditional independence relations among multiple variables and their time lags is that of time series graphical models (Eichler, 2011).

Consider a multivariate process $X$ of dimension $N$⁠. We define the time series graph $G=(V×Z,E)$ of $X$ as follows [Fig. 1(c)]: The set of nodes in that graph consists of the set of components $V$ at each time $t∈Z$⁠. That is, the graph is actually infinite, but in practice defined up to some maximum time lag $τmax$⁠. Compared to the general concept of graphical models (Lauritzen, 1996) for data without time-ordering, for time series graphs the time-dependence is explicitly used to define directional links in the set of edges $E$ (Eichler, 2011). For convenience, we treat $X$⁠, $Xt$⁠, and $Xt−$ as sets of random variables here and use the difference symbol “$∖$” for sets.

For $i=j$⁠, we call links autodependencies. For $τ=0$⁠, the same definition can also be used to define undirected contemporaneous links (see Eichler, 2011; Runge, 2015), which would lead to the time series graph being a mixed graph instead of a directed acyclic graph. The arrow of time is a convenient way to disambiguate independence relationships: If we do not have access to the time ordering of the variables (or there is none) and observe as the only conditional independence relation $X⊥⊥Y∣Z$ while all other relations are dependent, then this relation can be generated by any of the three causal motifs $X→Z→Y$⁠, $Y→Z→X$⁠, or $X←Z→Y$ which form a Markov equivalence class.

Here, we define links with lags $t−τ$ relative to some time point $t$⁠, but throughout this paper we assume stationarity (discussed in Sec. IV E). Then a link is repeated for every $t′<t$ if a link exists at time $t$⁠. Alternatively, the links can be estimated from different realizations at time $t$⁠. In practice, however, the links will mostly be estimated from single time series realizations requiring stationarity.

The parents of a node $Xtj$ are defined as

In the following, $A,B,S$ denote nodes or sets of nodes in the graph and $Xt−τ,Yt,Zt−τ,Ut−τ,Z$ random variables of the process $X$⁠, sometimes dropping the subscript. We will denote a general conditional dependence measure as $I(X;Y|Z)$ which can be CMI or also some other measure such as partial correlation, depending on the context.

When considering the dependency between two variables $X,Y$ given a set of conditions $Z$ as in $I(X;Y|Z)$⁠, the idea of open and blocked paths or separation between the corresponding nodes in the time series graph $G$ [Fig. 1(c)] is important. A directed path is a sequence of linked nodes containing only motifs $→∙→$⁠. But there are also other paths on which information is shared even though no causal interventions could “travel” along these. In general (Eichler, 2011), in the above defined time series graph, a path between two single nodes $A$ and $B$ is called open if it contains only the motifs $→∙→$ or $←∙→$⁠. For notational convenience, we will sometimes use left-pointing arrows, while still in the time series graph all directed links are forward in time. For example, $∙←∙→∙→∙$ is an open path. On the other hand, if any motif on a path is $→∙←$⁠, the path is blocked. Nodes in such motifs are also called colliders. If we now consider a separating or conditioning set $S$, openness and blockedness of motifs reverse, i.e., denoting a conditioned node by $◼$⁠, the motifs $→◼→$ and $←◼→$⁠, are blocked and the motif $→◼←$ becomes open. For example, the path $∙←◼→∙→∙$ is blocked, while $∙←∙→◼←∙$ is open. Note that paths can also traverse links repeatedly, e.g., forward and backward.

Intuitively, if two nodes are separated, no information is shared between the two. For example, in Fig. 1(c) $Xt1$ and $Xt−14$ are separated by $S={Xt−22,Xt−24}$ or also by $S={Xt−11,Xt−22}$⁠. Then $I(Xt−14;Xt1|S)=0$⁠. Conversely, $Xt1$ and $Xt−33$ are still connected given $S={Xt−22}$ since the “back door”-path $Xt−33←Xt−43→Xt−32→Xt−11→Xt1$ is still open. These definitions are important for the relations between the graph and the underlying process.

Causal information cannot be obtained from associations of measured variables without some assumptions. A variety of different assumptions have been shown to be sufficient to estimate the true causal graph (Spirtes et al., 2000; Peters et al., 2017). Here, we focus on three main assumptions under which the time series graph represents causal relations: Causal Sufficiency, the Causal Markov Condition, and Faithfulness. For time-lagged causal discovery from observational time series, we also need the assumptions of no instantaneous effects and stationarity. Further added to these are dependence type assumptions (e.g., linear or nonlinear) and no measurement error, and we will also assume that the joint distribution of the process has a positive density. All of these are discussed in the following.

For illustrating some of the assumptions in this paper, we will estimate the time series graphs by directly testing Definition 1 via Eq. (6) [Fig. 1(c)], mostly in the partial correlation version.

As Granger (1969) already notes, “[t]he one completely unreal aspect of the above definitions is the use of the series $Ut−$ representing all available information [in the universe].” This definition makes sure that the measured variables include all of the common causes of $X$ and $Y$⁠. However, we do not always need the whole universe. If we have available only a limited number of measured variables, we only need to assume that there exist no other unobserved (or latent) variables that directly or indirectly influence any other pair of our set of variables which is the assumption of Causal Sufficiency (Spirtes et al., 2000).

All independence-based causal discovery methods necessitate the Causal Markov Condition (Spirtes et al., 2000) which constitutes a close relationship between the process $X$ and its graph $G$⁠.

Note that if Causal Sufficiency is not fulfilled, then also generally the Markov condition will not hold (Spirtes et al., 2000). Intuitively, the Causal Markov Condition implies that once we know the values of $Yt$’s parents, all other variables in the past (⁠$t−τ$ for $τ>0$⁠) become irrelevant for predicting $Yt$⁠. Of course, $Y$’s descendants at future time points can also “predict” $Yt$⁠.

The Causal Markov Condition guarantees that separation in the graph implies independence in the process. But what can be concluded from an estimated conditional independence relation, that is, the reverse direction? Faithfulness guarantees that the graph entails all conditional independence relations that are implied by the Markov condition.

The combination of Faithfulness and the Markov property implies that $A⋈B∣S⇔XA⊥⊥XB∣XS$ and its logical contraposition $A⧸⋈B∣S⇔XA⧸⊥⊥XB∣XS$⁠. Both conditions are an important assumption for causal discovery algorithms as discussed in Spirtes et al. (2000). Intuitively, Faithfulness together with the Causal Markov Condition allow us to conclude that (in the limit of infinite sample size) a measured statistical dependency is actually due to some (not necessarily direct) causal mechanism and, conversely, a measured independence (given any set of conditions) implies that no direct causal mechanism exists (see also Remark 1 in the Discussion).

The former example only illustrated a static case of determinism. The field of nonlinear dynamics studies the properties of nonlinear and chaotic dynamical processes which has led to a plethora of nonlinear time series analysis methods (Kantz and Schreiber, 2003), often from an information-theoretic perspective (Hlavácková-Schindler et al., 2007) including transfer entropy. Many of these methods built on the assumption that no system is perfectly deterministic, for example, due to the coarse-graining of the system’s phase-space in the measurement process. In Sec. VII A, we study the effect of dynamical noise on several common time-series based causal discovery approaches for chaotic systems.

Granger causality and the definition of time series graphs are examples for lagged definitions of causality. To guarantee that the lagged parents defined in Eq. (8) are sufficient for the Causal Markov Condition to hold, we need to assume that there are no instantaneous (contemporaneous) causal effects, i.e., $Xti→Xtj$⁠. One may argue that causality between dynamical systems cannot have instantaneous effects because the speed of light is finite and, if the process is sampled with sufficient resolution (otherwise, see the Examples of sub-sampling and aggregation), we only need to consider lagged causal effects. However, we often do not have a sufficiently sampled time series. Here, recent developments in causal inference theory (Zhang and Hyvärinen, 2009; Peters et al., 2013; Lopez-Paz et al., 2015; Spirtes and Zhang, 2016; Peters et al., 2017) address instantaneous causality within the framework of structural causal models which can be applied to determine causal directionality for contemporaneous links. These models work under assumptions on the noise in the model such as non-Gaussianity. Also, the logical causal orientation rules of the causal discovery approaches in Spirtes et al. (2000) can be used to partially orient contemporaneous links.

To estimate the time series graph defined in Definition 1 from time series data, we assume stationarity. Another option would be to utilize independent ensembles of realizations of lagged processes. Here, we define stationarity with respect to a time index set $T$⁠. For example, $T$ can contain all time indices belonging to a certain regime of a dynamical process, e.g., only winter months in the climate system.

This constitutes actually a weaker form of stationarity than the common definition of stationarity in mean, variance, spectral properties, or of the value of individual coefficients in a linear model. For example, one could require that all CMIs are stationary,

which is a much stronger statement. The strength of causal mechanisms may fluctuate over time and the causal stationarity assumption only requires conditional independence to be stationary.

To test conditional independence hypotheses $X⊥⊥Y∣Z$⁠, different test statistics can be utilized. These are typically based on making certain assumptions about the type of the underlying dependency structure. While classical statistical methods are often based on the assumption of linearity (which allows us to derive rigorous results), modern statistics, the physics community, and the recent field of machine learning have developed non-parametric or model-free methods that allow us to better capture the nonlinear reality of many dynamical complex systems—at the cost of weaker theoretical results. Conditional independence testing can be classified into regression-based and model-free approaches. Here, we only discuss tests for continuously valued variables, for discrete variables one can, for example, use methods based on contingency tables (Spirtes et al., 2000) or discrete CMI estimation (Cover and Thomas, 2006).

Regression-based conditional independence tests of $X⊥⊥Y|Z$ are based on first regressing out the influence of $Z$ from $X$ and $Y$ and then testing the dependence between the residuals. We first fit a model assuming

for centered variables $X,Y$ and independent and identically normally distributed $ϵX,Y$⁠. Now further restrictions can be laid upon the functional form of $fX,Y$⁠. For example, the partial correlation test assumes linearity, while a non-parametric regression can be based on Gaussian Process regression (Rasmussen and Williams, 2006).

Secondly, from the estimated functions $f^$⁠, the residuals are formed as

Finally, the dependence between the residuals can be tested with different pairwise association tests. For partial correlation this is a $t$-test, while the dependence between the residuals of a non-parametric regression can be tested with non-parametric tests (Gretton et al., 2008; Székely et al., 2007) such as the distance correlation coefficient$R(rX,rY)$ (Székely et al., 2007) (see also Sec. V C). Note that these models all make parametric assumptions and, thus, do not estimate conditional independence in its most general form.

The other extreme to partial correlation are model-free methods that directly test conditional independence. The most prominent test statistic is CMI as defined in Eq. (3), for which non-parametric estimators based on nearest-neighbor statistics exist (Kraskov et al., 2004; Frenzel and Pompe, 2007; Vejmelka and Palus, 2008; Póczos and Schneider, 2012) [see also Gao et al. (2015) and Lord et al. (2018) for recent progress on nearest-neighbor entropy estimators]. Other possible conditional independence tests are Kernel Conditional Independence Tests (Zhang et al., 2011; Strobl et al., 2017) which essentially test for zero Hilbert-Schmidt norm of the partial cross-covariance operator or conditional distance correlation (Wang et al., 2015). Some new recent tests are based on neural networks (Sen et al., 2017) or decision tree regression (Chalupka et al., 2018). In Runge (2018), a conditional independence test based on CMI is introduced.

Measurement error, unlike dynamical noise, contaminates the variables between which we seek to reconstruct dependencies and constitutes a difficult problem in causal network reconstruction (Scheines and Ramsey, 2016).

The previous sections concerned fundamental assumptions of time-lagged causal discovery based on conditional independence relations. We now turn to the topic of practical estimation where we introduce several common causal discovery methods and discuss their consistency, significance testing, and computational complexity. We restrict the analysis to the class of conditional independence approaches, which flexibly allows us to use different independence tests. But graphical models, in general, can also be estimated with score-based Bayesian methods, e.g., the max-min hill-climbing algorithm (Tsamardinos et al., 2006).

Transfer entropy, as introduced by Schreiber (2000), is a direct information-theoretic implementation of Granger causality (Barnett et al., 2009). In a lag-specific implementation, as given by FullCI [Eq. (6)], it tests for conditional independence between each $Xt−τi$ and $Xtj$ conditioned on the entire past $Xt(t−1,…,t−τmax)$ (excluding $Xt−τi$⁠). As the experiments in Sec. VII will demonstrate, this approach strongly suffers from the curse of dimensionality.

Sun and Bollt (2014) and Sun et al. (2015) developed a discovery algorithm based on the information-theoretic optimal causation entropy principle (algorithm abbreviated as OCE) which reconstructs the lagged parents of a variable $Xtj$ by an iterative procedure alleviating the curse of dimensionality: Starting with an empty parent set $P^OCE(Xtj)=∅$⁠, first the MIs $I(Xt−τi;Xtj)$ for all $Xt−τi∈Xt−$ are estimated. As the first parent $X(1)$⁠, the one with the largest MI with $Xtj$ is selected. The next parent $X(2)$⁠, however, is chosen according to the largest CMI $I(Xt−τi;Xtj|X(1))$ among all remaining variables, the third parent is the one with largest CMI conditional on the two previously selected parents, etc. The process is continued until the CMI of a selected parent is non-significant. This forward-selection stage is followed by a backward elimination whereby the significance of each of the parents $Xt−τi∈P^OCE(Xtj)$ is tested conditional on the remaining parents:

The significance of CMIs can be tested with a nearest-neighbor CMI estimator (Kraskov et al., 2004; Frenzel and Pompe, 2007; Vejmelka and Palus, 2008) in combination with a permutation test where $Xt−τi$ is randomly shuffled. Of course, the conditional independencies in Eq. (26) can also be tested with other test statistics.

An alternative to this forward-backward scheme is the PC algorithm (named after its inventors Peter and Clark) (Spirtes and Glymour, 1991). The original PC algorithm was formulated for general random variables without assuming a time order. It consists of several phases where first, in the skeleton-discovery phase, an undirected graphical model (Lauritzen, 1996) is estimated whose links are then oriented using a set of logical rules (Spirtes and Glymour, 1991; Spirtes et al., 2000). A later improvement led to the more robust modification called PC-stable (Colombo and Maathuis, 2014).

For the case of time series, we can use the information of time order which naturally provides an orientation rule for links. The algorithm then is as follows: For every variable $Xtj∈Xt$ it starts by initializing the preliminary parents $P^(Xtj)=(Xt−1,Xt−2,…,Xt−τmax)$⁠. In the first iteration (⁠$p=0$⁠), we remove a variable $Xt−τi$ from $P^(Xtj)$ if the null hypothesis

cannot be rejected at a significance threshold $α$⁠. Then, in each next iteration, we increase $p→p+1$ and remove a variable $Xt−τi$ from $P^(Xtj)$ if any of the null hypotheses

cannot be rejected, where $S$ iterates (in an inner loop) through all combinations of subsets $S⊆P^(Xtj)∖{Xt−τi}$ with cardinality $p$⁠. The algorithm converges for a link $Xt−τi→Xtj$ once $S=P^(Xtj)∖{Xt−τi}$ and the null hypothesis $Xt−τi⊥⊥Xtj∣P^(Xtj)∖{Xt−τi}$ is rejected (if the null hypothesis cannot be rejected, the link is removed). Runge et al. (2018) provide pseudo-code for this algorithm.

The forward-backward scheme of OCE conducts conditional independence tests only using the conditions with highest CMI in the preceding stage and quickly increases the number of conditions. This can lead to wrong parents being kept in $P^OCE(Xtj)$ (see Sec. V B) which are only removed in the backward stage where the dimensionality of the set $P^OCE(Xtj)$ can be already quite high. High dimensionality, in principle, leads to lower detection power (Example VII C). The PC algorithm conducts conditional independence tests not only using the condition with highest association, but it goes through (in theory all) combinations of conditions $S$ which can help to alleviate the curse of dimensionality regarding the estimation dimension of $I(Xt−τi;Xtj∣S)$ compared to OCE. On the other hand, the PC algorithm conducts many more tests which increases other problems (see Sec. V C).

A more recent approach that addresses some of the shortcomings of the PC algorithm above is PCMCI (Runge et al., 2018). PCMCI is a two-step approach which uses a version of the PC-algorithm only as a condition-selection step (PC$1$ algorithm) to obtain $P^(Xtj)$ for all $Xtj∈Xt$⁠, followed by the momentary conditional independence (MCI) test for $Xt−τi→Xtj$ defined as

with $Xt−=(Xt−1,Xt−2,…,Xt−τmax)$⁠. The main difference between PC and PC$1$ is that PC$1$ tests only the condition subset $S$ with largest association instead of going through all possible combinations. To this end $P^(Xtj)$ is sorted after every iteration according to the test statistic value and $S$ is determined by the first $p$ variables in $P^(Xtj)$ (excluding $Xt−τi$⁠). This leads to less tests, but still provably removes incorrect links.

The MCI test is the most important difference to the PC algorithm and the approach by Sun and Bollt (2014). The additional conditioning on the parents $P^(Xt−τi)$ in MCI accounts for autocorrelation leading to well-controlled false positive rates at the expected level (Runge et al., 2018). A variant (PCMCI$0$⁠) where the condition on the parents $P^(Xt−τi)$ is dropped leads to a very similar approach to OCE. PCMCI, like FullCI, OCE, and PC, can be implemented with different conditional independence tests. For further details on PCMCI, see Runge et al. (2018). In Sec. VII we compare FullCI, OCE, and PCMCI in a number of numerical comparison studies.

Consistency is an important property of causal methods that tells us whether the method provably converges to the true causal graph in the limit of infinite sample size. Consistency concerns the conditional independence tests on the one hand, but also the causal algorithm in the case of iterative approaches such as those discussed in Sec. V A.

For example, for the consistency of the non-parametric regression independence tests in Eq. (23), we need to assume that the function estimator converges to the true function, that the noise in the model is additive and independent, and finally that we have a consistent unconditional independence test for the residuals. With a consistent test, the time series graph can be directly estimated based on Definition 1. For iterative causal algorithms, we can define universal consistency as follows.

That is, the probability of estimating the wrong graph becomes arbitrarily small if enough data is available, for any distribution $P$ (hence “universal”). Consistency has been proven for classical causal discovery algorithms such as the PC-algorithm (Spirtes et al., 2000), the optimal causation approach Sun and Bollt (2014), Sun et al. (2015) and PCMCI Runge et al. (2018), as an approach based on PC.

However, universal consistency is a weaker statement than, for example, uniform consistency which bounds the error probability as a function of the sample size $n$ giving a rate of convergence. Thus, for a non-uniform, but only universally consistent method, the sample size at which a given error can be guaranteed and can be different for every distribution $P$⁠. Robins et al. (2003) showed that no uniformly consistent causal discovery algorithm from the class of independence-based approaches (Spirtes et al., 2000) exists since the convergence can always be made arbitrarily slow by a distribution that is almost unfaithful with some dependencies made arbitrarily small. Uniform consistency for conditional-independence based algorithms can only be achieved under further assumptions such as having strong enough dependencies (Kalisch, 2008).

How can we assess the statistical significance of conditional independence tests on which the causal algorithms in Sec. V A are based, such as the tests discussed in Sec. IV F? Using a test statistic $I^n(x;y∣z)$ (⁠$I$ here stands not only for CMI, but any conditional independence test statistic) for the observed samples $(x,y,z)={xi,yi,zi}i=1n$ we wish to test the hypothesis

versus the general alternative

To assess the significance of an outcome of such a test using $p$-values, we need to know the distribution $Pr(In∣H0)$ of the test statistic under the null hypothesis. For partial correlation tests, exact analytical expressions of the null distribution exist under certain assumptions, but for nonlinear tests (such as CMI or also non-parametric regression tests) these are typically not available except for some asymptotic large sample size cases (Strobl et al., 2017). The alternative then are permutation tests as discussed in Example 12.

Given the null distribution, the $p$-value is defined as the probability—given $H0$—of observing a value of the test statistic that is the same or more extreme than what was actually observed. If our test statistic is non-negative (such as CMI), the $p$-value for an observed test statistic value $I^n$ is defined as $p=Pr(In≥I^n∣H0)$⁠. Choosing a significance level $α$⁠, we reject the null hypothesis if $p<α$⁠.

There are two types of errors we can make. Rejecting $H0$ when $H0$ is true is called a type I error or false positive. Retaining $H0$ when $H1$ is true is called a type II error or false negative. We also call one minus the type II error rate of a test the true positive rate or detection rate.

If the test statistic has a continuous distribution, then under $H0$ the $p$-value has a uniform $(0,1)$ distribution (Wasserman, 2004). Therefore, if we reject $H0$ when the $p$-value is less than $α$⁠, the probability of a type I error is $α$⁠. For a well-calibrated test under $H0$⁠, we thus expect to measure on average a false positive rate of $α$⁠. If this is not the case, the test is ill-calibrated indicating that we got the null distribution wrong because certain assumptions, such as independent samples, are violated. In Examples 12 and 13, we discuss such cases.

Application areas of causal discovery methods vary in the typical numbers of variables as well as available sample sizes $n$⁠. Next to the properties discussed before, an important issue then is how a method scales with dimensionality and sample size. High-dimensionality arises from the number of included variables $N$ and the maximum time lag $τmax$ [see Fig. 1(c)] and has at least two consequences: (1) higher computational complexity leading to longer runtimes and (2) typically lower detection power. Independence tests may also become ill-calibrated in high-dimensions.

For directly testing causal links via Definition 1 (FullCI), the computational complexity depends on the complexity of a single high-dimensional conditional independence test. In the linear partial correlation case, OLS regression scales $∼O(n(Nτmax)2)$⁠. FullCI estimated using nearest-neighbor estimators of CMI (Kraskov et al., 2004; Frenzel and Pompe, 2007), on the other hand, will scale $∼O(nlog⁡n)$ regarding time complexity while the complexity in $Nτmax$ will depend on algorithmic details such as using efficient KD-tree nearest-neighbor search procedures (Maneewongvatana and Mount, 1999).

The methods PC, OCE, or PCMCI (Sec. V A) based on a condition-selection step avoid high-dimensional conditional independence estimation by conducting more tests with lower dimensional conditioning sets. Their theoretical complexities are difficult to evaluate, for numerical evaluations see Sun et al. (2015) and Runge et al. (2018), but typically they scale polynomially in time.

The other major challenge with high dimensionality is detection power as analyzed in Example VII C.

How can causal discovery methods, such as those described in Sec. V be evaluated? Typically, we want to know which method performs best on data from the kind of system we wish to study. Ideally, we would like to compare different methods on a data sample where the underlying causal truth is known or evaluate methods by experimentally manipulating a system, i.e., actually performing the do-experiment (Pearl, 2000) mentioned in the introduction which forms the theoretical basis of the present concept of causality. Since both of these options are mostly not available, an alternative is to construct synthetic model data where the underlying ground truth is known. These can then be used to study the performance of causal methods for realistic finite sample situations.

To evaluate causal methods on synthetic data, several aspects for constructing model systems are relevant:

1. Model realism: The model systems should mimic the domain-specific properties of real data in terms of nonlinearity, autocorrelation, spectral properties, noise structure (dynamical as well as observational), etc.

2. Model diversity: To avoid biased conclusions, a large number of different randomly selected connectivity structures should be tested [including link density as well as properties such as small-worldness (Watts and Strogatz, 1998)]. For example, the aforementioned forward-selection approach failed for the example shown in Fig. 11 but works for many other graphs. But also consistent methods may have biases for finite samples as studied in Runge et al. (2018).

3. Model dimensionality: As studied in Fig. 16, a method may perform well only for a small number of variables and the performance for high-dimensional settings (large networks) can quickly degrade.

4. Sample sizes: The comparative performance of different models may vary widely for different sample sizes which needs to be studied if no uniform consistency results are available.

The performance of a causal method on a single realization of a model does not allow for reliable conclusions. Therefore, each model needs to be evaluated from many realizations. Then the most straightforward evaluation metric is to measure false positive rates and true positive rates for a given $α$ as shown in the comparison studies in Sec. VII. The number of realizations should be chosen high enough since the error of a false or true positive rate $r$ for $B$ realizations is given by $σr=r(1−r)/B$⁠. Alternative evaluation metrics that do not depend on a particular significance level but directly on the $p$-values are the Kullback-Leibler divergence to evaluate whether the $p$-values are uniformly distributed (to measure how well-calibrated a test is) and the Area Under the Power Curve (AUPC) to evaluate true positives.

Next to the true and false positives of a causal method for finite samples, another performance criterion is computational runtime, though this may strongly depend on a given implementation.

In this section, we compare several common causal discovery methods in three numerical comparison studies highlighting the effect of dynamical noise in deterministic chaotic systems, autocorrelation, and high dimensionality.

In Example 6, we studied a static example of determinism. Here, we evaluate the effect of dynamical noise in a system of coupled chaotic logistic maps:

with uniformly distributed independent noise $η$ and $r=4$ leading to chaotic dynamics. Here, $σ$ controls the amount of dynamical noise in the system. To evaluate true positive rates (correctly detecting $Zt−1→Xt$ and $Zt−1→Yt$⁠) and false positive rates (incorrect detections for any other variable pair, direction, or lag), $200$ realizations with time series length $n=150$ were generated.

We compare three methods from an information-theoretic framework (FullCI, OCE, PCMCI) with convergent-cross mapping (CCM, Sugihara et al., 2012) (see also Arnhold et al., 1999; Hirata et al., 2016) as a nonlinear dynamics-inspired approach. The significance of CMIs in FullCI and OCE is tested with a nearest-neighbor CMI estimator (Kraskov et al., 2004; Frenzel and Pompe, 2007; Vejmelka and Palus, 2008) in combination with a permutation test where $Xt−τi$ is randomly shuffled. PCMCI is implemented with the CMIknn independence test (Runge, 2018) as discussed in Example 12. We also evaluate a variant (PCMCI$0$⁠) where the condition on the parents $P^(Xt−τi)$ is dropped and the only difference to OCE is the condition-selection step. CCM reconstructs the variable’s state-spaces using lagged coordinate embedding and concludes on $X→Y$ if points on $X$ can be well predicted using nearest neighbors in the state space of $Y$⁠. Note that CCM and related works only use the time series of $X$ and $Y$ with the underlying assumption that the dynamics of $Z$ can be reconstructed using delay embedding. All methods were evaluated at a significance level of $0.05$⁠. For implementation details, see  Appendix B 1. In the top panel of Fig. 14, we depict the true causal graph.

Figure 14 shows that in the purely deterministic regime with $σ=0$⁠, PCMCI has almost no power, while PCMCI$0$ features a detection rate of $0.8$ and FullCI, OCE, and CCM almost always detect the couplings. On the other hand, CCM also has the highest number of false positives around $0.2$ which exceeds the significance level indicating an ill-calibrated test. FullCI, OCE, and PCMCI$0$ also do not control false positives well, while with PCMCI they are around the expected level.

For higher dynamical noise levels interesting behavior emerges: FullCI and CCM have continuously decreasing power (and slightly decreasing false positives) dropping to almost zero at $σ=0.4$⁠, while the power of OCE and both PCMCI versions steadily increases up to $σ=0.2$ after which power decreases with a more pronounced decrease for PCMCI and PCMCI$0$ having the highest power. False positives are above the $0.05$ threshold for FullCI and OCE for a wide range of noise levels, while for PCMCI false positives are always well-controlled at the expected $0.05$⁠.

How can these results be understood? CCM attempts to reconstruct the attractor manifolds underlying $X$ and $Y$⁠. With more dynamical noise, this reconstruction becomes more difficult and CCM looses power. The fact that CCM does not control false positives well, especially in the deterministic regime deserved further study. FullCI suffers from the curse of dimensionality especially for higher noise levels.

The MCI test statistic for the link $Zt−1→Xt$ estimates $IZt−1→XtMCI=I(Zt−1;Xt∣Zt−2,Xt−1)$⁠. Now since for $σ=0$ $Zt−1$ is a deterministic mapping of $Zt−2$⁠, in theory we have $IMCI=0$ for the same reasons as in the simple deterministic model (19) above (and analogously for $Zt−1→Yt$⁠). In practice, we can only measure the entropies at some coarse-grained level (here determined by the CMI nearest-neighbor parameter) and the deterministic dependency is never exactly recovered leading to a non-zero MCI. In OCE and PCMCI$0$⁠, only the parents of $Xt$ are included in the condition. The fact that FullCI, despite conditioning on the whole past, also detects the link deserves further study. Possible explanations are an ill-calibrated test or dynamical properties of the logistic-map system. In summary, purely deterministic dynamics here seem to generate too little momentary information which is necessary for information-theoretic coupling detection with MCI (see also Pompe and Runge, 2011).