theorem proposition lemma corollary definition remark

# Doubly Robust Structure Identification
from Temporal Data

Emmanouil Angelis† Helmholtz AI, Helmholtz Center Munich  Techincal University
of Munich  Francesco Quinzan† Department of Computer Science, University of
Oxford  Ashkan Soleymani  Laboratory for Information & Decision Systems
(LIDS), Massachusetts Institute of Technology
Patrick Jaillet  Laboratory for Information & Decision Systems (LIDS),
Massachusetts Institute of Technology  Stefan Bauer  Helmholtz AI, Helmholtz
Center Munich  Techincal University of Munich

###### Abstract

Learning the causes of time-series data is a fundamental task in many
applications, spanning from finance to earth sciences or bio-medical
applications. Common approaches for this task are based on vector auto-
regression, and they do not take into account unknown confounding between
potential causes. However, in settings with many potential causes and noisy
data, these approaches may be substantially biased. Furthermore, potential
causes may be correlated in practical applications. Moreover, existing
algorithms often do not work with cyclic data. To address these challenges, we
propose a new doubly robust method for Structure Identification from Temporal
Data (SITD). We provide theoretical guarantees, showing that our method
asymptotically recovers the true underlying causal structure. Our analysis
extends to cases where the potential causes have cycles and they may be
confounded. We further perform extensive experiments to showcase the superior
performance of our method.

##  1 Introduction

One of the primary objectives when working with time series data is to uncover
the causal structures between different variables over time. Learning these
causal relations and their interactions is of critical importance in many
disciplines, e.g., healthcare [Anwar et al., 2014], climate studies [Stips et
al., 2016, Runge et al., 2019a], epidemiology [Hernán et al., 2000, Robins et
al., 2000], finance [Hiemstra and Jones, 1994], ecosystems [Sugihara et al.,
2012], and many more. Interventional data is not often accessible in many of
these applications. For instance, in healthcare scenarios, conducting trials
on patients may raise ethical concerns, or in the realm of earth and climate
studies, randomized controlled trials are not feasible.

In general, understanding the underlying causal graph using only observational
data is a cumbersome task due to many reasons: i) observational data, as
opposed to interventional data, capture correlation-type relations111which
capture how changes in one variable are associated with changes in another
variable. instead of cause-effect ones. ii) unobserved confounders introduce
biases and deceive the algorithms to falsely infer causal relations instead of
the true structure, e.g., the existence of a hidden common confounder iii) the
number of possible underlying structures grows super exponentially with the
number of variables creating major statistical and computation barriers iv)
the identifiability problem, since multiple causal models can result in the
same observational distribution, thus making it impossible to uniquely
determine the true structure. To overcome these problems and determine the
underlying structure, additional assumptions are imposed. Typical assumptions
include faithfulness, linearity of relations, or even noise-free settings,
which limit the types of causal relationships that can be discovered [Pearl,
2000, Peters et al., 2017, Spirtes et al., 2000, Glymour et al., 2019]. Almost
all of these challenges 222The arrow of time assumption can slightly simplify
a very limited number of cases. See Axiom (B) for definition. extend to the
problem of identifying the underlying causal structure from observational
time-series datasets.

Subsequently, in many instances, the emphasis is placed on particular target
variables of interest and their causal features. Causal features of a target
are defined as the set of variables that conditioned on them, the target
variable is independent of the rest variables. Causal feature selection
enables training models which are much simpler, more interpretable, and more
robust [Aliferis et al., 2010a, Janzing et al., 2020]. However, learning the
causal structures between variables and a specific target is still a demanding
task, and many current approaches for causal feature selection face
limitations by making unrealistic simplifying assumptions about the data-
generating process or by lacking computational and/or statistical scalability
[Yu et al., 2021, 2020]. As an important example, under the faithfulness
assumption, causal features are equivalent to the Markov Blanket (MB) of the
target variable [Pearl, 2000]. Consequently, numerous research efforts have
been dedicated to recovering the MB of the target variable [Tsamardinos et
al., 2003b, a, Pena, 2008, Aliferis et al., 2010a, Yaramakala and Margaritis,
2005, Aliferis et al., 2010b, Margaritis and Thrun, 1999, Aliferis et al.,
2003, Borboudakis and Tsamardinos, 2019, Tsamardinos et al., 2019].
Nonetheless, this assumption is confining and these algorithms generally face
the curse of dimensionality statistically and/or computationally. These
challenges become particularly amplified in the context of time series data,
where the number of variables grows linearly with the length of the data
trajectories. As a result, to mitigate these problems, additional assumptions,
e.g., stationarity or no hidden confounders are included [Moraffah et al.,
2021, Bussmann et al., 2021, Runge, 2018] and/or weaker notions of
causality333weaker than Pearl’s structural equation model [Pearl, 2000]. such
as Granger causality have been studied extensively [Granger, 1988, 1969,
Marinazzo et al., 2011, Tank et al., 2018, Bussmann et al., 2021, Khanna and
Tan, 2019, Runge, 2018, Hasan et al., 2023].

The existing causal feature selection from time-series data algorithms often
assume some level of faithfulness or causal sufficiency [Runge et al., 2019b,
Runge, 2018, 2020]. Oftentimes, they overlook the presence of unknown
confounding factors among potential causes [Moraffah et al., 2021]. Moreover,
most cannot adapt to cyclic settings [Entner and Hoyer, 2010], which is
relatively ubiquitous in many domains [Bollen, 1989]. Furthermore, many
algorithms employ the popular vector auto-regression framework to model time-
dependence structures [Bussmann et al., 2021, Lu et al., 2016, Chen et al.,
2009, Weichwald et al., 2020, Hyvärinen et al., 2010], which again is
restrictive. To overcome these problems and inspired by the debiased machine
learning literature [Chernozhukov et al., 2018], we propose an efficient
algorithm for doubly robust structure identification from temporal data.

#### Our contribution.

  1. $\bullet$

We provide an efficient and easy-to-implement doubly robust structure
identification from temporal data algorithm (SITD) with theoretical guarantees
enjoying $\sqrt{n}$-consistency.

  2. $\bullet$

We provide an extensive technical discussion relating Granger causality to
Pearl’s framework for time series and show under which assumptions our
approach can be used for feature selection or full causal discovery.

  3. $\bullet$

We provide theoretical insights demonstrating that our algorithm doesn’t need
faithfulness or causal sufficiency while allowing for general non-linear
cyclic structures and also the presence of hidden confounders among the
covariates for the task of finding the direct causes of the target variable.
An important problem in many applications.

  4. $\bullet$

In extensive experiments we illustrate that our approach is significantly more
robust, significantly faster, and more performative than state-of-the-art
baselines.

##  2 Related Work

A longstanding line of work intends to tailor the existing causal structure
learning and Markov blanket discovery for i.i.d. data to the temporal setting.
To name a few, Entner and Hoyer [2010] adapted the Fast Causal Inference
algorithm [Spirtes et al., 2000] to time-series data. While the approach
shares the benefit of being able to deal with hidden confounders, it is not
possible to account for cyclic structures. Runge et al. [2019b] introduced
PCMCI, as an adjusted version of PC [Spirtes et al., 2000] with an additional
false positive control phase which is able to recover time-lagged causal
dependencies. PCMCI+ modified the approach further to additionally be able to
find contemporaneous causal edges [Runge, 2020]. LPCMCI extends the scope by
catering to the case of hidden confounders [Gerhardus and Runge, 2020]. Even
though methods in this category are able to provide theoretical guarantees for
learning the causal structure, SITD has several advantages over them: i) For
all of these methods, the faithfulness assumption is a key ingredient while
SITD does not need it. ii) These methods are based on conditional independence
testing which is widely recognized to be a cumbersome statistical problem
[Bergsma, 2004, Kim et al., 2022]. Shah and Peters [2020] have established
that no conditional independence (CI) test can effectively control the Type-I
error for all CI settings. Moreover in practice, conducting many conditional
independence tests from lengthy time-series is burdensome. iii) Even having
access to a perfect conditional independence test oracle, severe computational
challenges exist [Chickering, 1996, Chickering et al., 2004]444Learning
Bayesian Networks with conditional independence tests is NP-Hard. (please
refer to Section I for a detailed comparison).

Another line of work relies heavily on the vector auto-regression framework.
VarLiNGAM [Hyvärinen et al., 2010] generalizes LiNGAM [Shimizu et al., 2006]
to time-series and similar to that it assumes a linear non-Gaussian acyclic
model for the data. In the work of Huang et al. [2019], a time-varying linear
causal model is assumed, allowing for causal discovery even in the presence of
hidden confounders. More recently, deep neural networks are used to train
vector auto-regression. Tank et al. [2018] adapted neural networks (named cMLP
and cLSTM) for Granger causality by imposing group-sparsity regularizers. In a
similar fashion, Khanna and Tan [2019] used recurrent neural networks. In
another work, Bussmann et al. [2021] designed neural additive vector auto-
regression (NAVAR), a neural network approach to model non-linearities. In
contrast to previous works, they extract Granger-type causal relations by
injecting the necessary sparsity directly into the architecture of the neural
networks. This line of work is quite limited to ours as they consider
confining structural assumptions over the underlying causal structural
equations; details on these assumptions are discussed next to Axiom (A) in the
following section.

The use of double robustness in causality problems has a long history mainly
concentrated on estimating average treatment effect [Robins et al., 1994, Funk
et al., 2011, Benkeser et al., 2017, Bang and Robins, 2005, Słoczyński and
Wooldridge, 2018]. Chernozhukov et al. [2018] introduced the DoubleML
framework to achieve double robustness for structural parameters. Upon that,
Soleymani et al. [2022] defined an orthogonalized score to infer the direct
causes of partially linear models. Their approach is fast and allows for the
assumption of complicated underlying structures but unfortunately, it is
limited to only linear direct causal effects. While this was later extended to
the non-linear case [Quinzan et al., 2023], we propose a doubly robust
approach for identifying causal structures from temporal data under the
general assumptions discussed below.

$\scriptstyle X_{0}^{1}$$\scriptstyle X_{1}^{1}$$\scriptstyle
X_{2}^{1}$$\scriptstyle X_{3}^{1}$$\scriptstyle X_{4}^{1}$$\scriptstyle
X_{5}^{1}$$\scriptstyle Y_{0}$$\scriptstyle Y_{1}$$\scriptstyle
Y_{2}$$\scriptstyle Y_{3}$$\scriptstyle Y_{4}$$\scriptstyle
Y_{5}$$\scriptstyle X_{0}^{2}$$\scriptstyle X_{1}^{2}$$\scriptstyle
X_{2}^{2}$$\scriptstyle X_{3}^{2}$$\scriptstyle X_{4}^{2}$$\scriptstyle
X_{5}^{2}$$\scriptstyle U_{0}$$\scriptstyle U_{1}$$\scriptstyle
U_{2}$$\scriptstyle U_{3}$$\scriptstyle U_{4}$$\scriptstyle
U_{5}$$\scriptstyle X^{1}$$\scriptstyle Y$ $\scriptstyle X^{2}$$\scriptstyle
U$

Figure 1:  Example of a full-time graph (left), and the corresponding summary
graph (right). Note that time series $\boldsymbol{X}^{1}$ and
$\boldsymbol{X}^{2}$ causally influence the outcome $\boldsymbol{Y}$ with
different lags. The time series $\boldsymbol{U}$ is unobserved, and it acts as
a confounder for $\boldsymbol{X}^{1}$ and $\boldsymbol{X}^{2}$. Moreover,
there is a cycle between the confoundings variables $\boldsymbol{X}^{1}$ and
$\boldsymbol{X}^{2}$ of the outcome variable $\boldsymbol{Y}$.

##  3 Framework

###  3.1 Problem Description

We are given i.i.d. observations from a univariate target time series
$\boldsymbol{Y}\coloneqq\\{Y_{t}\\}_{t\in\mathbb{Z}}$ and i.i.d. observations
from a multivariate time series
$\boldsymbol{X}\coloneqq\\{X_{t}^{1},\dots,X_{t}^{m}\\}_{t\in\mathbb{Z}}$ of
potential causes. We assume that the outcome time series is specified by some
of the potential causes by a deterministic function with posterior additive
noise, and no instantaneous effects. We can formalize this model as

  1. Axiom (A)

$Y_{T}=f(\mathsf{pa}_{T}(\boldsymbol{Y}),T)+\varepsilon_{T}$,

with $\varepsilon_{T}$ exogenous independent noise and
$\mathsf{pa}_{T}(\boldsymbol{Y})\subseteq\\{X_{t}^{1},\dots,X_{t}^{m}\\}_{t\in\mathbb{Z}}$
is a subset of random variables of the multivariate time series
$\boldsymbol{X}$. Note that the independence of the additive noise is
important for identifiability. In fact, if there are dependencies between the
noise and the history, then one might run into identifiability problems. We
refer the reader to Appendix C for a counterexample. We are interested in
identifying the time series $\boldsymbol{X}^{i}$ that directly affects the
outcome $Y$. That is, we wish to identify time series $\boldsymbol{X}^{i}$
such that it holds $X^{i}_{t}\subseteq\mathsf{pa}_{T}(\mathbf{Y})$ for some
time steps $t,T$.

As discussed in Section 2, vector auto-regression methods enforce heavy
structural assumptions on the underlying causal structural equations. NAVAR
[Bussmann et al., 2021] assumes that each variable is influenced by its causal
parents exclusively in an additive way and higher-order interactions among
them are precluded. In mathematical terms, they follow

$Y_{T}=\beta+\sum_{X\in\mathsf{pa}_{T}(\boldsymbol{Y})}^{N}f_{X}(X_{t-\kappa:t-1})+\eta_{t},$

where $\beta$ is a bias term, $\kappa$ is the time lag, and $\eta_{t}$ is an
independent noise variable. VarLiNGAM [Hyvärinen et al., 2010] takes a step
further and uses much more restrictive assumptions on the functional form of
the structural equations, confining to $f_{X}$ being a linear transformation
of $X_{t-\kappa:t-1}$. These structural assumptions are quite limited as they
do not include even simple functionals, e.g., $Y_{t}=X_{t-1}^{1}\times
X_{t-1}^{2}$ or $Y_{t}=\log(X_{t-1}^{1}+X_{t-1}^{2})$, whereas Axiom (A) is
able to capture any interaction between the covariates and timestep $T$ which
is a significant improvement over existing works.

#### Additional assumptions.

We further use the following assumptions:

  1. Axiom (B)

there are no causal effects backward in time;

  2. Axiom (C)

there are no instantaneous causal effect between $\boldsymbol{Y}$ and any of
the
potential causes $\boldsymbol{X}^{i}$.

Both Axiom (B) and Axiom (C) appear in previous related work (see, e.g.,
Peters et al. [2013], Mastakouri et al. [2021], Löwe et al. [2022]). Note that
these axioms allow for cycles and hidden common confounders between the
potential causes. Axiom (B) is a natural assumption as a system is called
causal when the output of the system only depends on the past, not the future
[Peters et al., 2017, Pearl, 2000]. Axiom (C) poses additional restrictions on
the class of models that we consider, since instantaneous effects may be
relevant in some cases and applications [Lippe et al., 2022a]555An example of
cases where time series exhibit instantaneous causal effects is given by
dynamical systems [Mogensen et al., 2018, Rubenstein et al., 2016]. In
dynamical systems, a variable may instantaneously affect another variable of
the model. Instantaneous effects have been studied in previous related work
[Gong et al., 2022] but due to identifiability issues, they rely on stronger
assumptions such as faithfulness.. However, without this assumption, it is
impossible to learn causes _from observational data_. The necessity of Axiom
(C) for causal discovery is a well-known fact [Peters et al., 2013]666In
general, Peters et al. [2013] show that causal discovery is impossible with
instantaneous effects. However, Peters et al. [2013] also provide a special
case in which the causal structure is identifiable with instantaneous effects.
This special case occurs when the random variables of the model are jointly
Gaussian, and the instantaneous effects are linear. Our framework extends to
this special case.. For the sake of completeness, we have furnished an example
of this identifiability problem tailored to time-series data in Section D. For
simplicity, we also consider the following assumption:

  1. Axiom (D)

the causal relationships are stationary, i.e., they are invariant under a
joint shift of all the variables.

Axiom (D) again follows previous related work [Peters et al., 2013, Mastakouri
et al., 2021, Löwe et al., 2022]. This assumption is not strictly necessary
for our analysis, but it is a general assumption that simplifies modeling. Our
results can be extended to models that do not fulfill this axiom.

Our framework retains some degree of generality over previous related work. In
fact, Axiom (A)-Axiom (D) allow for hidden confounding and cycles between the
potential causes (see Figure 1), providing a more general framework than the
full autoregressive model studied, e.g., by Hyvärinen et al. [2010], Peters et
al. [2013], Löwe et al. [2022], Wu et al. [2020]. Furthermore, in contrast to
several previous works [Khemakhem et al., 2020, Gresele et al., 2021,
Lachapelle et al., 2022, Lippe et al., 2022b, Yao et al., 2021], we do not
assume that the variables $Y_{T},X_{T}^{1},\dots,X_{T}^{m}$ are independent
conditioned on the observed variables at previous time steps. Importantly, we
also do not assume causal faithfulness in the full time graph, or weaker
notions such as causal minimality.777Recall that a distribution is faithful to
a causal diagram if no conditional independence relations are present, other
than the ones entailed by the Markov property. This is a major improvement
over some previous works, e.g., Mastakouri et al. [2021], Gong et al. [2022],
since there is no reason to assume that faithfulness or causal minimality hold
in practice. Theorem 1 by Gong et al. [2022] provides an identifiability
result for a model with history-dependent noise and instantaneous effects.
This result, however, requires causal minimality.

We remark that much of the previous related work considers the general problem
of full causal discovery, whereas we use Axiom (A)-Axiom (C) for the more
specific task of causal feature selection. However, our method can also be
used for full causal discovery assuming a fully-observed acyclic system (see
Section 5.2 for details).

###  3.2 Causal Structure

We are interested in direct causal effects, which are defined by distribution
changes due to interventions on the DGP. An intervention amounts to actively
manipulating the generative process of a potential cause $\boldsymbol{X}^{i}$
at some time step $t$, without altering the remaining components of the DGP.
Then, a time series $\boldsymbol{X}^{i}$ has a direct effect on
$\boldsymbol{Y}$ if performing an intervention on some temporal variable
$X^{i}_{t}$ will alter the distribution of $Y_{T}$, for some time steps $t$,
$T$. For a precise definition of these concepts, we refer the reader to
Appendix B.

Following, e.g., Mastakouri et al. [2021], we define the _full time_ graph
$\mathcal{G}$ as a directed graph whose edges represent all direct causal
effects among the variables at all time steps. Given the outcome $Y_{T}$ at a
given time step, we refer to the parent nodes in the full time graph as its
_causal parents_. We further define the _summary_ graph whose nodes are
$\boldsymbol{X}^{i}$ and $\boldsymbol{Y}$, and with directed edges
representing causal effects between the time series. We refer the reader to
Figure 1 for a visualization of these graphs. Note that the causes of
$\boldsymbol{Y}$ correspond to the parent nodes of $\boldsymbol{Y}$ in the
summary graph.

###  3.3 Granger Causality

Granger causality [Granger, 1988] is one of the most commonly used approaches
to infer causal relations from observational time-series data. Its central
assumption is that cause-effect relationships cannot work against time.
Informally, if the prediction of the future of a target time-series
$\boldsymbol{Y}$ can be improved by knowing past elements of another time-
series $\boldsymbol{X}^{i}$, then $\boldsymbol{X}^{i}$ Granger causes
$\boldsymbol{Y}$. Formally, we say that $\boldsymbol{X}^{i}$ Granger causes
$\boldsymbol{Y}$ if it holds

$\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{X}^{i}_{T}=\boldsymbol{x},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)\neq\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)$ (1)

for a non-zero probability event
$\\{\boldsymbol{X}^{i}_{T}=\boldsymbol{x},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\\}$, where $\boldsymbol{I}_{T}$ stands for the set
$\\{\boldsymbol{X}^{1}_{T},\boldsymbol{X}^{2}_{T},\dots,\boldsymbol{X}^{m}_{T}\\}$
and $\boldsymbol{I}_{T}^{\setminus i}$ represents the set
$\boldsymbol{I}_{T}\setminus\\{\boldsymbol{X}^{i}_{T}\\}$.

Granger causality is commonly used to identify causes. Assuming stationary,
multivariate Granger causality analysis usually fits a vector autoregressive
model to a given dataset. This model can be then used to determine the causes
of a target $\boldsymbol{Y}$. However, it is important to note that Granger
causality does not imply true causality in general. This limitation was
acknowledged by Granger himself in Granger [1988].

##  4 Methodology

###  4.1 Double Machine Learning (DoubleML)

DoubleML is a general framework for parameter estimation, which uses debiasing
techniques to achieve $\sqrt{n}$-consistency (see, e.g., Rotnitzky et al.
[2020], Chernozhukov et al. [2022]). In DoubleML, we consider the problem of
estimating a parameter $\theta_{0}$ as a solution of an equation of the form
$\mathbb{E}\left[\mathcal{L}(\theta_{0},\eta_{0})\right]=0$. The score
function $\mathcal{L}$ depends on two terms, the true parameter $\theta_{0}$
that we wish to estimate, and a nuisance parameter $\eta_{0}$. We do not
directly care about the correctness of our estimate of $\eta_{0}$, as long as
we get a good estimator of $\theta_{0}$. The nuisance parameter $\eta_{0}$ may
induce an unwanted bias in the estimation process, resulting in slow
convergence.

To overcome this problem, we consider score functions that fulfill the Mixed
Bias Property (MBP) [Rotnitzky et al., 2020]. Intuitively, the MBP holds if
small changes of the nuisance parameter do not significantly affect the score
function computed around the true parameters $(\theta_{0},\eta_{0})$. If the
given score function fulfills the MBP, then it is possible to partially remove
the bias induced by the nuisance parameter using (double) cross-
fitting.888Although in our analysis we focus on the MBP, strong consistency
results follow using the weaker Neyman Orthogonaloty condition. By removing
this bias, the convergence rate of the parameter estimation improves. We refer
the reader to Chernozhukov et al. [2018], Rotnitzky et al. [2020],
Chernozhukov et al. [2022] for details.

Following the notation introduced earlier, for a fixed time step $T$, let
$\boldsymbol{X}$ be a random variable in the set of random variables
$\boldsymbol{V}$, $g$ any real-valued function of $\boldsymbol{X}$ such that
$\mathbb{E}\left[g^{2}(\boldsymbol{X})\right]<\infty$. We consider parameters
of the form $\theta_{0}\coloneqq\mathbb{E}\left[m(\boldsymbol{V};g)\right]$,
where $m(\boldsymbol{V};g)$ is a linear moment functional in $g$. The
celebrated Riesz Representation Theorem ensures that, under certain
conditions, there exists a function $\alpha_{0}$ of $\boldsymbol{X}$ such that
$\mathbb{E}\left[m(\boldsymbol{V};g)\right]=\mathbb{E}\left[\alpha_{0}(\boldsymbol{X})g(\boldsymbol{X})\right]$.
The function $\alpha_{0}$ is called the _Riesz Representer_ (RR). Chernozhukov
et al. [2021] shows that the Riesz representer can be estimated from samples.
Using the RR, we can derive a score function for the parameter $\theta_{0}$
with $g_{0}(\boldsymbol{X})=\mathbb{E}\left[Y\mid\boldsymbol{X}\right]$ that
fulfills the MBP.

Chernozhukov et al. [2022], Rotnitzky et al. [2020] shows that we can learn
the parameters $\theta_{0}$ as above, with a score function of the form

$\varphi(\theta,\boldsymbol{\eta})\coloneqq
m(\boldsymbol{V};g)+\alpha(\boldsymbol{X})\cdot(Y-g(\boldsymbol{X}))-\theta.$
(2)

Here, $\boldsymbol{\eta}\coloneqq(\alpha,g)$ is a nuisance parameter
consisting of a pair of square-integrable functions. As shown by Chernozhukov
et al. [2022], the score function Eq. (2) yields
$\mathbb{E}\left[\varphi(\theta_{0},\boldsymbol{\eta})\right]=-\mathbb{E}\left[(\alpha(\boldsymbol{X})-\alpha_{0}(\boldsymbol{X}))(g(\boldsymbol{X})-g_{0}(\boldsymbol{X}))\right]$.
This equation corresponds to the MBP as in Definition 1 of Rotnitzky et al.
[2020]. Hence, if we use (double) cross-fitting, the resulting estimator
$\hat{\theta}$ has the _double robustness property_. That is, the quantity
$\sqrt{n}(\hat{\theta}-\theta_{0})$ converges in distribution to a zero-mean
Normal distribution if the product of the mean-square convergence rates of
$\alpha$ and $g$ is faster than $\sqrt{n}$.

###  4.2 Granger Causality Implies True Causation
under Axiom (A)-Axiom (C)

As discussed in Section 3.3, Granger causality does not imply true causality
in general. In our case, however, Granger causality corresponds to true
causation, as stated in the following result.

######  Theorem 4.1.

Consider a causal model as in Axiom (A)-Axiom (C). Then, a time series
$\boldsymbol{X}^{i}$ is a potential cause of $\boldsymbol{Y}$ if and only if
(iff.) $\boldsymbol{X}^{i}$ Granger causes $\boldsymbol{Y}$.

Proof of this result is given in Appendix E. Importantly, Theorem 4.1 does not
require causal faithfulness.7 We remark that Peters et al. [2013] shows that
Granger causality implies true causation for fully autoregressive models (see
also Löwe et al. [2022]). This result is based on the identifiability of
additive noise models, in which all relevant variables are observed [Peters et
al., 2011]. Our framework, however, is more general than Peters et al. [2013],
Löwe et al. [2022], since it allows for confounding among the covariates (see
Figure 1). In the special case of a fully autoregressive model, Theorem 4.1 is
equivalent to previous results [Peters et al., 2013, Löwe et al., 2022].

11: split the dataset $\mathsf{D}$ into $k$ random disjoint subsets
$\mathsf{D}_{j}$ of the same size;

2: for each dataset partition $j$ do

3: perform regressions to estimate $\alpha^{0}_{j}$, $g^{0}_{j}$ on dataset
$\mathsf{D}\setminus\mathsf{D}_{j}$;

4: $\theta_{j}^{0}\leftarrow\hat{\mathbb{E}}_{\mathsf{D}_{j}}[Y_{T}\cdot
g_{j}^{0}+\alpha_{j}^{0}\cdot(Y_{T}-g_{j}^{0})]$;

5: end for

26: $\theta^{0}\leftarrow k^{-1}\sum_{j}\theta_{j}^{0}$;

7: for each potential cause $\boldsymbol{X}^{i}$ do

8: for each dataset partition $j$ do

9: apply zero-masking to $\alpha^{0}_{j}$, $g^{0}_{j}$ for the variables
$\boldsymbol{X}^{i}$, and denote with $\tilde{\alpha}^{0}_{j}$,
$\tilde{g}^{0}_{j}$ the resulting models;

10:
$\theta_{j}^{i}\leftarrow\hat{\mathbb{E}}_{\mathsf{D}_{j}}[Y_{T}\cdot\tilde{g}_{j}^{i}+\tilde{\alpha}_{j}^{i}\cdot(Y_{T}-\tilde{g}_{j}^{i})]$;

11: end for

12: $\theta^{i}\leftarrow k^{-1}\sum_{j}\theta_{j}^{i}$

13: perform a paired Student’s $t$-test to determine if
$\theta^{i}\approx\theta^{0}$, and select time series $\boldsymbol{X}^{i}$ as
a cause if the null-hypotheses is rejected;

314: end for

15: return the selected time series;

Algorithm 1 The SITD

###  4.3 Testing Granger Causality with DoubleML

Our approach to identifying potential causes consists of performing a
statistical test to determine if Eq. 1 holds. Due to Theorem 4.1, a
straightforward approach would then consist of using a conditional
independence test, to select or discard a time series $\boldsymbol{X}^{i}$ as
a cause of the outcome $\boldsymbol{Y}$. However, conditional independence
testing is challenging in high-dimensional settings. Furthermore, kernel-based
conditional independence tests [Fukumizu et al., 2007, Zhang et al., 2011,
Park and Muandet, 2020, Sheng and Sriperumbudur, 2020] are computationally
expensive. Instead, we provide a new statistical test based on DoubleML. Our
approach is based on the observation that under Axiom (A)-Axiom (C), the
condition in Eq. 1 can be written in terms of simple linear moment
functionals. The following theorem holds.

######  Theorem 4.2.

Consider the notation as in Eq. 1, and fix a time step $T$. For any square-
integrable random variable
$g^{0}(\boldsymbol{X}^{i}_{T},\boldsymbol{I}_{T}^{\setminus i})$, consider the
moment functional $m_{0}(\boldsymbol{V};g^{0})\coloneqq Y_{T}\cdot g^{0}$.
Similarly, for any square-integrable random variable
$g^{i}(\boldsymbol{I}_{T}^{\setminus i})$, consider the moment functional
$m_{i}(\boldsymbol{V};g^{i})\coloneqq Y_{T}\cdot g^{i}$. Assuming Axiom
(A)-Axiom (D), $\boldsymbol{X}^{i}$ Granger causes $\boldsymbol{Y}$ iff. it
holds

$\mathbb{E}\left[m_{0}(\boldsymbol{V};g_{0}^{0})\right]-\mathbb{E}\left[m_{i}(\boldsymbol{V};g_{0}^{i})\right]\neq
0,$ (3)

with
$g_{0}^{0}(\boldsymbol{x},\boldsymbol{i})=\mathbb{E}[Y_{T}\mid\boldsymbol{X}^{i}_{T}=\boldsymbol{x},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}]$, and
$g_{0}^{i}(\boldsymbol{i})=\mathbb{E}[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}]$.

The proof is deferred to Appendix F. By this theorem, we can identify the
causal parents of $Y$ by testing Eq. 3. This boils down to learning parameters
$\theta^{0}_{0}\coloneqq\mathbb{E}[Y_{T}\cdot g^{0}_{0}]$,
$\theta^{i}_{0}\coloneqq\mathbb{E}[Y_{T}\cdot g^{i}_{0}]$, and then performing
a sample test to verify if it holds $\theta^{0}-\theta^{i}\neq 0$.

Intuitively, the parameter
$\theta\coloneqq\mathbb{E}[Y_{T}.\mathbb{E}[Y_{T}|I]]$ quantifies the cross-
correlation between the target variable $Y_{T}$ and its conditional mean
$\mathbb{E}[Y_{T}|I]$, given the set of variables $I$. In general terms,
variations in the parameter $\theta$ due to changes in the set $I$ suggest the
presence of a causal relationship between alterations in the set $I$ and the
target variable $Y_{T}$.

#### Outline of the approach.

Given a sample dataset $\mathsf{D}$, we iterate the following steps:

  1. Step (1)

Select a potential cause $\boldsymbol{X}^{i}$ to test if $\boldsymbol{X}^{i}$
Granger causes $\boldsymbol{Y}$. Split the dataset $\mathsf{D}$ into $k$
disjoint sets $\mathsf{D}_{j}$ with $k\geq 2$.

  2. Step (2)

Learn $g^{0}_{0}$, $g^{i}_{0}$ as in Theorem 4.2 on dataset
$\mathsf{D}\setminus\mathsf{D}_{j}$. Denote with $g^{0}_{j}$, $g^{i}_{j}$ the
resulting estimates. Learn the RRs $\alpha^{0}_{0}$, $\alpha^{i}_{0}$
coresponding to $\mathbb{E}\left[Y_{T}\cdot g_{0}^{0}\right]$ and
$\mathbb{E}\left[Y_{T}\cdot g_{0}^{i}\right]$ respectively, on dataset
$\mathsf{D}\setminus\mathsf{D}_{j}$. Denote with $\alpha^{0}_{j}$,
$\alpha^{i}_{j}$ the resulting estimates.

  3. Step (3)

Learn the parameter $\theta^{0}_{0}=\mathbb{E}[Y_{T}\cdot g^{0}_{0}]$ on
dataset $\mathsf{D}_{j}$ for all $j$, by minimizing the orthogonal score
function as in Eq. 2. Learn the parameter
$\theta^{i}_{0}=\mathbb{E}[Y_{T}\cdot g^{i}_{0}]$ on dataset $\mathsf{D}_{j}$
for all $j$, by minimizing the orthogonal loss as in Eq. 2. Denote with
$\theta_{j}^{0}$, $\theta_{j}^{i}$ the resulting estimates, and define
$\theta^{0}\coloneqq k^{-1}\sum_{j}\theta_{j}^{0}$ and $\theta^{i}\coloneqq
k^{-1}\sum_{j}\theta_{j}^{i}$.

  4. Step (4)

Use a two-sample Student’s $t$-test to determine if
$\mathbb{E}\left[\theta^{0}-\theta^{i}\right]$ is approximately zero. Select
$\boldsymbol{X}^{i}$ as a cause of $\boldsymbol{Y}$ if the null hypothesis is
rejected.

This procedure can be iterated for each of the features.

#### Strong consistency guarantees.

Step (4) requires some explanation. Under mild conditions on the convergence
of $g^{0}_{j},g^{i}_{j}$ and $\alpha^{0}_{j},\alpha^{i}_{j}$, the quantity
$\theta^{0}-\theta^{i}$ has $\sqrt{n}$-consistency. As a consequence of strong
consistency, it holds

$\lim_{n\to+\infty}\sqrt{n}\mathbb{E}\left[\theta^{0}-\theta^{i}\right]=0,$
(4)

iff.
$\mathbb{E}\left[Y_{T}\cdot\eta_{0}^{0}\right]-\mathbb{E}\left[Y_{T}\cdot\eta_{0}^{i}\right]=0$.
Then, by Theorem 4.2, Eq. 4 holds iff. $\boldsymbol{X}^{i}$ does not Granger
causes $\boldsymbol{Y}$. Step (4) in the procedure outlined above uses a two-
sample Student’s $t$-test to determine if $\boldsymbol{X}^{i}$ Granger causes
$\boldsymbol{Y}$. We refer the reader to Chernozhukov et al. [2022, 2018] for
a proof of the $\sqrt{n}$-consistency for estimates as $\theta^{0}$ and
$\theta^{i}$.

##  5 The Algorithm

###  5.1 Overview

Our approach to test Granger causality essentially follows Step (1)-Step (4)
outlined in the previous section. We refer to our approach as the SITD
(Structure Identification from Temporal Data). Our method is presented in
Algorithm 1, where we use the notation $\hat{\mathbb{E}}_{\mathsf{D}}$ to
denote the empirical expected value over a given dataset $\mathsf{D}$. We note
that, however, this procedure requires training $g^{0}_{j}$, $g^{i}_{j}$ and
$\alpha^{0}_{j}$, $\alpha^{i}_{j}$ for each potential cause, and for each
partition $\mathsf{D}_{j}$. This can be time-consuming in large instances. To
overcome this bottleneck, we make the following improvements:

  1. $\bullet$

The parameters $g^{0}_{j}$, $\alpha^{0}_{j}$, $\theta^{0}_{j}$, and
$\theta^{0}$ are the same for each one of the potential causes, and we only
trained them once at the beginning of the process (Line 2-6 of Algorithm 1).

  2. $\bullet$

The RR corresponding to
$\mathbb{E}\left[m_{0}(\boldsymbol{V};g_{0}^{0})\right]$ as in Theorem 4.1 can
be written in explicit form as
$\alpha_{0}^{0}=\mathbb{E}[Y_{T}\mid\boldsymbol{X}^{i}_{T}=\boldsymbol{x},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}]$. Similarly, the RR of
$\mathbb{E}\left[m_{i}(\boldsymbol{V};g_{i}^{0})\right]$ as in Theorem 4.1 can
be written in explicit form as
$\alpha_{0}^{i}=\mathbb{E}[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}]$. Hence, we estimate $\alpha_{0}^{0}$, $\alpha_{0}^{i}$ by
performing regressions.

  3. $\bullet$

Instead of computing $g_{0}^{i}$ and $\alpha_{0}^{i}$ directly, we apply a
zero-masking layer to $g^{0}_{j}$ and $\alpha^{0}_{j}$ for the features
$\boldsymbol{X}^{i}$. This masking layer tells the sequence-processing layers
that the input values for features $\boldsymbol{X}^{i}$ should be skipped. We
then compute $\theta^{i}_{0}$ using the resulting surrogate models
$\tilde{g}^{i}_{j}$, $\tilde{\alpha}^{i}_{j}$ (Line 9-10 in Algorithm 1).
Please refer to Section G on why zero-masking is not hurting the estimations.

#### Runtime of Algorithm 1.

Much of the run time of our algorithm consists of performing a regression to
learn $\eta^{0}_{j}$ on dataset $\mathsf{D}_{j}$. Denote with $d$ the time
complexity of performing such a regression. For a given $k$-partition of the
dataset, we can upper-bound the time complexity of our algorithm as
$\mathcal{O}(dk)$. Here, $d$ depends on the specific techniques used for the
regression. Non-parametric regression can be performed efficiently in the
problem size. We analyse the computational complexity in Appendix I further
and showp runtime plots for extensive experiments in Appendix J.3.

###  5.2 Full Causal Discovery

It is possible to use Algorithm 1 for full causal discovery, for fully-
observed acyclic auto-regressive models with no instantaneous effects (see
Peters et al. [2013], Löwe et al. [2022] for a precise definition of this
restricted framework). In fact, under these more restrictive assumptions, we
can identify the causes of each random variable of the model by testing
Granger causality. For these models, we can learn the full summary graph, by
identifying the causes of each variable with Algorithm 1. The resulting run
time can be quantified as $\mathcal{O}(mdk)$, where $d$ is the time complexity
of performing a regression, as outlined above, $m$ is the number of time
series considered, and $k$ is the number of dataset partitions used for double
cross-fitting. Since $k=\mathcal{O}(1)$ w.r.t. parameters of the problem (see
Section H), the runtime of SITD is $\mathcal{O}(md)$.

##  6 Experiments

###  6.1 Semi-Synthetic Experiments

#### The Dream3 dataset.

Following previous related work [Tank et al., 2018, Khanna and Tan, 2019,
Nauta et al., 2019, Bussmann et al., 2021, Gong et al., 2022], we evaluate
performance with the Dream3 benchmark [Prill et al., 2010, Marbach et al.,
2009]. Dream3 consists of measurements of gene expression levels for five
different networks. Each network contains $100$ genes. In this dataset, time
series correspond to different perturbation trajectories. For each one of the
five networks, there are $m=46$ different trajectories.

#### Description of the experiments.

Our goal is to infer the causal structure of each group of time series
(trajectories) of Dream3, as detailed above. This gives us a total of five
tasks, i.e., E.Coli 1, E.Coli 2, Yeast 1, Yeast 2, and Yeast 3. We run our
algorithm on this dataset and we use the area under the ROC curve (AUROC) as
the performance metric. Following [Gong et al., 2022], we compare against the
following baselines: cMLP [Tank et al., 2018], cLSTM [Tank et al., 2018], TCDF
[Nauta et al., 2019], SRU [Khanna and Tan, 2019], eSRU [Khanna and Tan, 2019],
Dynotears [Pamfil et al., 2020], Rhino+g [Gong et al., 2022], and Rhino [Gong
et al., 2022].

#### Results.

The results for cMLP, cLSTM, TCDF, SRU and SRU are taken directly from Khanna
and Tan [2019], Gong et al. [2022]. Regarding Rhino+g and Rhino, we partially
reproduce the experiments by Gong et al. [2022], using their source code. Our
implementation of Rhino+g and Rhino differs from Gong et al. [2022] only in
the choice for the hyper-parameters, which is the same on all five tasks. We
specifically use the following hyper-parameters for Rhino: Node Embedd. =
$16$, Instantaneous eff. = False, Node Embedd. (flow) = $16$, lag = $2$,
$\lambda_{s}$ = $19$, Auglag = $30$. And we use the following for Rhino+g:
Node Embedd. = $16$, Instantaneous eff. = False, lag = $2$, $\lambda_{s}$ =
$15$, Auglag = $60$. This is the setting that is reported for the Ecoli1
subtask and found in the corresponding code implementation. We opted for this
approach because in the experiment by Gong et al. [2022], it seems that Rhino
overfitted to the dataset. We learn $\eta^{i}_{j}$ as in Line 3 of Algorithm
1, using a simple kernel ridge regression model with polynomial kernels of
degree three combined with zero masking.

Overall, we observe that our method gives competitive results. Specifically,
SITD outperforms all the other benchmarks on two tasks (E.Coli 1, E.Coli 2),
and it obtains the same AUROC as Rhino+g on Yeast 2. Furthermore, SITD obtains
comparable performance on the remaining tasks (Yeast 1 and Yeast 3). We also
would like to emphasize that our estimator is much simpler than deep nets such
as Rhino or Rhino+g. As such, it has lower sample complexity and lower run
time than the other algorithms.

###  6.2 Additional experiments

We provide additional extensive synthetic experiments for various noise-to-
signal ratios and number of causes in Appendix J. One appealing property of
SITD is the double robustness property. Due to this property, the dependence
of SITD on the estimator is relatively lower compared to existing baselines.
Using a fast and simple estimator with low statistical complexity in Algorithm
1, we are able to get competitive results compared to richer models especially
in low sample complexity settings. The illustrated case with a simple kernel
ridge regression estimator with polynomial kernels of degree three in Section
6.1 is a compelling example of this. In Appendix J, we will thus focus on
these practically important aspects in more detail and show in particular
further experiments on low sample settings in Appendix J.2 and runtime plots
in Appendix J.3.

Table 1: AUROC score for the SITD and common algorithms for Granger causality. Models with "${}^{\textrm{\scalebox{1.0}{\char 70}}}$" sign use deep neural networks. We remark that our estimator is much simpler than deep nets, and it has low sample complexity. We observe that our method achieves state-of-the-art performance on three tasks (E.Coli 1, E.Coli 2, Yeast 2) and that it obtains comparable performance on the remaining two (Yeast 1, Yeast 3). Method | E.Coli 1 | E.Coli 2 | Yeast 1 | Yeast 2 | Yeast 3
---|---|---|---|---|---
cMLP${}^{\textrm{\scalebox{1.0}{\char 70}}}$ | 0.644 | 0.568 | 0.585 | 0.506 | 0.528
cLSTM${}^{\textrm{\scalebox{1.0}{\char 70}}}$ | 0.629 | 0.609 | 0.579 | 0.519 | 0.555
TCDF${}^{\textrm{\scalebox{1.0}{\char 70}}}$ | 0.614 | 0.647 | 0.581 | 0.556 | 0.557
SRU${}^{\textrm{\scalebox{1.0}{\char 70}}}$ | 0.657 | 0.666 | 0.617 | 0.575 | 0.55
eSRU${}^{\textrm{\scalebox{1.0}{\char 70}}}$ | 0.66 | 0.629 | 0.627 | 0.557 | 0.55
DYNO. | 0.590 | 0.547 | 0.527 | 0.526 | 0.510
PCMCI+ | 0.530 | 0.519 | 0.530 | 0.510 | 0.512
Rhino${}^{\textrm{\scalebox{1.0}{\char 70}}}$ Reprod. | 0.686 | 0.625 | 0.667 | 0.580 | 0.540
Rhino+g${}^{\textrm{\scalebox{1.0}{\char 70}}}$ Reprod. | 0.662 | 0.606 | 0.663 | 0.585 | 0.577
SITD (ours) | 0.704$\pm 0.005$ | 0.680$\pm 0.004$ |  0.653$\pm 0.001$ | 0.585$\pm 0.003$ |  0.544$\pm 0.003$

##  7 Discussion

In this work, we propose an efficient algorithm for doubly robust structure
identification from temporal data. We further provide asymptotical guarantees
that our method is able to discover the direct causes even when there are
cycles or hidden confounding and that our algorithm has
$\sqrt{n}$-consistency. We extensively discuss the relations of the approach
between the popular frameworks of Granger and Pearl’s causality as well as
relate and extend approaches from debiased machine learning to structure
discovery from temporal data. We hope to apply these approaches to important
real-world applications in bio-medicine where robustness to confounding,
sample efficiency, and ease of use are important for causal discovery from
observational time series.

## References

  * Aliferis et al. [2003] C. F. Aliferis, I. Tsamardinos, and A. Statnikov.  Hiton: a novel markov blanket algorithm for optimal variable selection.  In _AMIA annual symposium proceedings_ , volume 2003, page 21. American Medical Informatics Association, 2003.
  * Aliferis et al. [2010a] C. F. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos.  Local causal and markov blanket induction for causal discovery and feature selection for classification part i: algorithms and empirical evaluation.  _Journal of Machine Learning Research_ , 11(1), 2010a.
  * Aliferis et al. [2010b] C. F. Aliferis, A. Statnikov, I. Tsamardinos, S. Mani, and X. D. Koutsoukos.  Local causal and markov blanket induction for causal discovery and feature selection for classification part ii: analysis and extensions.  _Journal of Machine Learning Research_ , 11(1), 2010b.
  * Anwar et al. [2014] A. R. Anwar, K. G. Mideska, H. Hellriegel, N. Hoogenboom, H. Krause, A. Schnitzler, G. Deuschl, J. Raethjen, U. Heute, and M. Muthuraman.  Multi-modal causality analysis of eyes-open and eyes-closed data from simultaneously recorded eeg and meg.  In _2014 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , pages 2825–2828. IEEE, 2014.
  * Bang and Robins [2005] H. Bang and J. M. Robins.  Doubly robust estimation in missing data and causal inference models.  _Biometrics_ , 61(4):962–973, 2005.
  * Benkeser et al. [2017] D. Benkeser, M. Carone, M. V. D. Laan, and P. B. Gilbert.  Doubly robust nonparametric inference on the average treatment effect.  _Biometrika_ , 104(4):863–880, 2017.
  * Bergsma [2004] W. P. Bergsma.  Testing conditional independence for continuous random variables.  2004\.
  * Bollen [1989] K. A. Bollen.  _Structural equations with latent variables_ , volume 210.  John Wiley & Sons, 1989.
  * Borboudakis and Tsamardinos [2019] G. Borboudakis and I. Tsamardinos.  Forward-backward selection with early dropping.  _The Journal of Machine Learning Research_ , 20(1):276–314, 2019.
  * Bussmann et al. [2021] B. Bussmann, J. Nys, and S. Latré.  Neural additive vector autoregression models for causal discovery in time series.  In _Discovery Science: 24th International Conference, DS 2021, Halifax, NS, Canada, October 11–13, 2021, Proceedings 24_ , pages 446–460. Springer, 2021.
  * Chen et al. [2009] G. Chen, J. P. Hamilton, M. E. Thomason, I. H. Gotlib, Z. S. Saad, and R. W. Cox.  Granger causality via vector auto-regression tuned for fmri data analysis.  In _Proc Intl Soc Mag Reson Med_ , volume 17, page 1718, 2009.
  * Chernozhukov et al. [2018] V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and structural parameters.  _The Econometrics Journal_ , 21(1):C1–C68, 01 2018.
  * Chernozhukov et al. [2021] V. Chernozhukov, W. K. Newey, V. Quintas-Martinez, and V. Syrgkanis.  Automatic debiased machine learning via neural nets for generalized linear regression.  _arXiv:2104.14737_ , 2021.
  * Chernozhukov et al. [2022] V. Chernozhukov, W. Newey, V. Quintas-Martinez, and V. Syrgkanis.  RieszNet and ForestRiesz: Automatic debiased machine learning with neural nets and random forests.  In _Proc. of ICML_ , pages 3901–3914, 2022.
  * Chickering [1996] D. M. Chickering.  Learning bayesian networks is np-complete.  _Learning from data: Artificial intelligence and statistics V_ , pages 121–130, 1996.
  * Chickering et al. [2004] M. Chickering, D. Heckerman, and C. Meek.  Large-sample learning of bayesian networks is np-hard.  _Journal of Machine Learning Research_ , 5:1287–1330, 2004.
  * Entner and Hoyer [2010] D. Entner and P. O. Hoyer.  On causal discovery from time series data using fci.  _Probabilistic graphical models_ , pages 121–128, 2010.
  * Fukumizu et al. [2007] K. Fukumizu, A. Gretton, X. Sun, and B. Schölkopf.  Kernel measures of conditional dependence.  In _Advances in Neural Information Processing Systems_ , pages 489–496, 2007.
  * Funk et al. [2011] M. J. Funk, D. Westreich, C. Wiesen, T. Stürmer, M. A. Brookhart, and M. Davidian.  Doubly robust estimation of causal effects.  _American journal of epidemiology_ , 173(7):761–767, 2011.
  * Gerhardus and Runge [2020] A. Gerhardus and J. Runge.  High-recall causal discovery for autocorrelated time series with latent confounders.  _Advances in Neural Information Processing Systems_ , 33:12615–12625, 2020.
  * Glymour et al. [2019] C. Glymour, K. Zhang, and P. Spirtes.  Review of causal discovery methods based on graphical models.  _Frontiers in genetics_ , 10:524, 2019.
  * Gong et al. [2022] W. Gong, J. Jennings, C. Zhang, and N. Pawlowski.  Rhino: Deep causal temporal relationship learning with history-dependent noise.  _CoRR_ , abs/2210.14706, 2022.
  * Granger [1969] C. W. Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica: journal of the Econometric Society_ , pages 424–438, 1969.
  * Granger [1988] C. W. Granger.  Some recent development in a concept of causality.  _Journal of econometrics_ , 39(1-2):199–211, 1988.
  * Gresele et al. [2021] L. Gresele, J. Von Kügelgen, V. Stimper, B. Schölkopf, and M. Besserve.  Independent mechanism analysis, a new concept?  In _Proc. of NeurIPS_ , volume 34, pages 28233–28248, 2021.
  * Hasan et al. [2023] U. Hasan, E. Hossain, and M. O. Gani.  A survey on causal discovery methods for temporal and non-temporal data.  _arXiv preprint arXiv:2303.15027_ , 2023.
  * Hernán et al. [2000] M. Á. Hernán, B. Brumback, and J. M. Robins.  Marginal structural models to estimate the causal effect of zidovudine on the survival of hiv-positive men.  _Epidemiology_ , pages 561–570, 2000.
  * Hiemstra and Jones [1994] C. Hiemstra and J. D. Jones.  Testing for linear and nonlinear granger causality in the stock price-volume relation.  _The Journal of Finance_ , 49(5):1639–1664, 1994.
  * Huang et al. [2019] B. Huang, K. Zhang, M. Gong, and C. Glymour.  Causal discovery and forecasting in nonstationary environments with state-space models.  In _International conference on machine learning_ , pages 2901–2910. PMLR, 2019.
  * Hyvärinen et al. [2010] A. Hyvärinen, K. Zhang, S. Shimizu, and P. O. Hoyer.  Estimation of a structural vector autoregression model using non-gaussianity.  _Journal of Machine Learning Research_ , 11:1709–1731, 2010.
  * Hyvärinen et al. [2010] A. Hyvärinen, K. Zhang, S. Shimizu, and P. O. Hoyer.  Estimation of a structural vector autoregression model using non-gaussianity.  _Journal of Machine Learning Research_ , 11(5), 2010.
  * Janzing et al. [2020] D. Janzing, L. Minorics, and P. Blöbaum.  Feature relevance quantification in explainable ai: A causal problem.  In _International Conference on artificial intelligence and statistics_ , pages 2907–2916. PMLR, 2020.
  * Khanna and Tan [2019] S. Khanna and V. Y. Tan.  Economy statistical recurrent units for inferring nonlinear granger causality.  _arXiv preprint arXiv:1911.09879_ , 2019.
  * Khemakhem et al. [2020] I. Khemakhem, D. Kingma, R. Monti, and A. Hyvarinen.  Variational autoencoders and nonlinear ica: A unifying framework.  In _Proc. of AISTATS_ , pages 2207–2217, 2020.
  * Kim et al. [2022] I. Kim, M. Neykov, S. Balakrishnan, and L. Wasserman.  Local permutation tests for conditional independence.  _The Annals of Statistics_ , 50(6):3388–3414, 2022.
  * Lachapelle et al. [2022] S. Lachapelle, P. Rodriguez, Y. Sharma, K. E. Everett, R. Le Priol, A. Lacoste, and S. Lacoste-Julien.  Disentanglement via mechanism sparsity regularization: A new principle for nonlinear ica.  In _Proc. of CLearR_ , pages 428–484, 2022.
  * Lippe et al. [2022a] P. Lippe, S. Magliacane, S. Löwe, Y. M. Asano, T. Cohen, and E. Gavves.  icitris: Causal representation learning for instantaneous temporal effects.  _CoRR_ , abs/2206.06169, 2022a.
  * Lippe et al. [2022b] P. Lippe, S. Magliacane, S. Löwe, Y. M. Asano, T. Cohen, and S. Gavves.  Citris: Causal identifiability from temporal intervened sequences.  In _Proc. of ICML_ , pages 13557–13603, 2022b.
  * Löwe et al. [2022] S. Löwe, D. Madras, R. Z. Shilling, and M. Welling.  Amortized causal discovery: Learning to infer causal graphs from time-series data.  In _Proc. of CLeaR_ , pages 509–525, 2022.
  * Lu et al. [2016] Z. Lu, C. Zhou, J. Wu, H. Jiang, and S. Cui.  Integrating granger causality and vector auto-regression for traffic prediction of large-scale wlans.  _KSII Transactions on Internet and Information Systems (TIIS)_ , 10(1):136–151, 2016.
  * Marbach et al. [2009] D. Marbach, T. Schaffter, C. Mattiussi, and D. Floreano.  Generating realistic in silico gene networks for performance assessment of reverse engineering methods.  _Journal of Computational Biology_ , 16(2):229–239, 2009.
  * Margaritis and Thrun [1999] D. Margaritis and S. Thrun.  Bayesian network induction via local neighborhoods.  _Advances in neural information processing systems_ , 12, 1999.
  * Marinazzo et al. [2011] D. Marinazzo, W. Liao, H. Chen, and S. Stramaglia.  Nonlinear connectivity by granger causality.  _Neuroimage_ , 58(2):330–338, 2011.
  * Mastakouri et al. [2021] A. Mastakouri, B. Schölkopf, and D. Janzing.  Necessary and sufficient conditions for causal feature selection in time series with latent common causes.  In _Proceedings of ICML_ , volume 139, pages 7502–7511, 2021.
  * Mogensen et al. [2018] S. W. Mogensen, D. Malinsky, and N. R. Hansen.  Causal learning for partially observed stochastic dynamical systems.  In _UAI_ , pages 350–360, 2018.
  * Moraffah et al. [2021] R. Moraffah, P. Sheth, M. Karami, A. Bhattacharya, Q. Wang, A. Tahir, A. Raglin, and H. Liu.  Causal inference for time series analysis: Problems, methods and evaluation.  _Knowledge and Information Systems_ , 63:3041–3085, 2021.
  * Nauta et al. [2019] M. Nauta, D. Bucur, and C. Seifert.  Causal discovery with attention-based convolutional neural networks.  _Machine Learning and Knowledge Extraction_ , 1(1):312–340, 2019.
  * Pamfil et al. [2020] R. Pamfil, N. Sriwattanaworachai, S. Desai, P. Pilgerstorfer, K. Georgatzis, P. Beaumont, and B. Aragam.  Dynotears: Structure learning from time-series data.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1595–1605. PMLR, 2020.
  * Park and Muandet [2020] J. Park and K. Muandet.  A measure-theoretic approach to kernel conditional mean embeddings.  In _Advances in Neural Information Processing Systems_ , 2020.
  * Pearl [2000] J. Pearl.  _Causality: Models, Reasoning and Inference_.  Cambridge University Press, 2000.
  * Pena [2008] J. M. Pena.  Learning gaussian graphical models of gene networks with false discovery rate control.  In _European conference on evolutionary computation, machine learning and data mining in bioinformatics_ , pages 165–176. Springer, 2008.
  * Peters et al. [2011] J. Peters, J. M. Mooij, D. Janzing, and B. Schölkopf.  Identifiability of causal graphs using functional models.  In _Proc. of UAI_ , pages 589–598, 2011.
  * Peters et al. [2013] J. Peters, D. Janzing, and B. Schölkopf.  Causal inference on time series using restricted structural equation models.  In _Proc. of NIPS_ , pages 154–162, 2013.
  * Peters et al. [2017] J. Peters, D. Janzing, and B. Schölkopf.  _Elements of causal inference: foundations and learning algorithms_.  The MIT Press, 2017.
  * Prill et al. [2010] R. J. Prill, D. Marbach, J. Saez-Rodriguez, P. K. Sorger, L. G. Alexopoulos, X. Xue, N. D. Clarke, G. Altan-Bonnet, and G. Stolovitzky.  Towards a rigorous assessment of systems biology models: the dream3 challenges.  _PloS one_ , 5(2):e9202, 2010.
  * Quinzan et al. [2023] F. Quinzan, A. Soleymani, P. Jaillet, C. R. Rojas, and S. Bauer.  Drcfs: Doubly robust causal feature selection.  In _International Conference on Machine Learning_ , pages 28468–28491. PMLR, 2023.
  * Robins et al. [1994] J. M. Robins, A. Rotnitzky, and L. P. Zhao.  Estimation of regression coefficients when some regressors are not always observed.  _Journal of the American statistical Association_ , 89(427):846–866, 1994.
  * Robins et al. [2000] J. M. Robins, M. A. Hernan, and B. Brumback.  Marginal structural models and causal inference in epidemiology.  _Epidemiology_ , pages 550–560, 2000.
  * Rotnitzky et al. [2020] A. Rotnitzky, E. Smucler, and J. M. Robins.  Characterization of parameters with a mixed bias property.  _Biometrika_ , 108(1):231–238, 08 2020.
  * Rubenstein et al. [2016] P. K. Rubenstein, S. Bongers, B. Schölkopf, and J. M. Mooij.  From deterministic odes to dynamic structural causal models.  _arXiv preprint arXiv:1608.08028_ , 2016.
  * Runge [2018] J. Runge.  Causal network reconstruction from time series: From theoretical assumptions to practical estimation.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 28(7):075310, 2018.
  * Runge [2020] J. Runge.  Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets.  In _Conference on Uncertainty in Artificial Intelligence_ , pages 1388–1397. PMLR, 2020.
  * Runge et al. [2019a] J. Runge, S. Bathiany, E. Bollt, G. Camps-Valls, D. Coumou, E. Deyle, C. Glymour, M. Kretschmer, M. D. Mahecha, J. Muñoz-Marí, et al.  Inferring causation from time series in earth system sciences.  _Nature communications_ , 10(1):2553, 2019a.
  * Runge et al. [2019b] J. Runge, P. Nowack, M. Kretschmer, S. Flaxman, and D. Sejdinovic.  Detecting and quantifying causal associations in large nonlinear time series datasets.  _Science advances_ , 5(11):eaau4996, 2019b.
  * Shah and Peters [2020] R. D. Shah and J. Peters.  The hardness of conditional independence testing and the generalised covariance measure.  2020\.
  * Sheng and Sriperumbudur [2020] T. Sheng and B. K. Sriperumbudur.  On distance and kernel measures of conditional independence, 2020.
  * Shimizu et al. [2006] S. Shimizu, P. O. Hoyer, A. Hyvärinen, A. Kerminen, and M. Jordan.  A linear non-gaussian acyclic model for causal discovery.  _Journal of Machine Learning Research_ , 7(10), 2006.
  * Słoczyński and Wooldridge [2018] T. Słoczyński and J. M. Wooldridge.  A general double robustness result for estimating average treatment effects.  _Econometric Theory_ , 34(1):112–133, 2018.
  * Soleymani et al. [2022] A. Soleymani, A. Raj, S. Bauer, B. Schölkopf, and M. Besserve.  Causal feature selection via orthogonal search.  _Transactions on Machine Learning Research_ , 2022.
  * Spirtes et al. [2000] P. Spirtes, C. N. Glymour, and R. Scheines.  _Causation, prediction, and search_.  MIT press, 2000.
  * Stips et al. [2016] A. Stips, D. Macias, C. Coughlan, E. Garcia-Gorriz, and X. S. Liang.  On the causal structure between co2 and global temperature.  _Scientific reports_ , 6(1):21691, 2016.
  * 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_ , 338(6106):496–500, 2012.
  * Tank et al. [2018] A. Tank, I. Covert, N. Foti, A. Shojaie, and E. Fox.  Neural granger causality for nonlinear time series.  _stat_ , 1050:16, 2018.
  * Tsamardinos et al. [2003a] I. Tsamardinos, C. F. Aliferis, and A. Statnikov.  Time and sample efficient discovery of markov blankets and direct causal relations.  In _Proceedings of the ninth ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 673–678, 2003a.
  * Tsamardinos et al. [2003b] I. Tsamardinos, C. F. Aliferis, A. R. Statnikov, and E. Statnikov.  Algorithms for large scale markov blanket discovery.  In _FLAIRS conference_ , volume 2, pages 376–380. St. Augustine, FL, 2003b.
  * Tsamardinos et al. [2019] I. Tsamardinos, G. Borboudakis, P. Katsogridakis, P. Pratikakis, and V. Christophides.  A greedy feature selection algorithm for big data of high dimensionality.  _Machine learning_ , 108:149–202, 2019.
  * Weichwald et al. [2020] S. Weichwald, M. E. Jakobsen, P. B. Mogensen, L. Petersen, N. Thams, and G. Varando.  Causal structure learning from time series: Large regression coefficients may predict causal links better in practice than small p-values.  In _NeurIPS 2019 Competition and Demonstration Track_ , pages 27–36. PMLR, 2020.
  * Williams [1991] D. Williams.  _Probability with Martingales_.  Cambridge University Press, 1991.
  * Wu et al. [2020] T. Wu, T. M. Breuel, M. Skuhersky, and J. Kautz.  Discovering nonlinear relations with minimum predictive information regularization.  _CoRR_ , abs/2001.01885, 2020.
  * Yao et al. [2021] W. Yao, Y. Sun, A. Ho, C. Sun, and K. Zhang.  Learning temporally causal latent processes from general temporal data, 2021.
  * Yaramakala and Margaritis [2005] S. Yaramakala and D. Margaritis.  Speculative markov blanket discovery for optimal feature selection.  In _Fifth IEEE International Conference on Data Mining (ICDM’05)_ , pages 4–pp. IEEE, 2005.
  * Yu et al. [2020] K. Yu, X. Guo, L. Liu, J. Li, H. Wang, Z. Ling, and X. Wu.  Causality-based feature selection: Methods and evaluations.  _ACM Computing Surveys (CSUR)_ , 53(5):1–36, 2020.
  * Yu et al. [2021] K. Yu, L. Liu, and J. Li.  A unified view of causal and non-causal feature selection.  _ACM Transactions on Knowledge Discovery from Data (TKDD)_ , 15(4):1–46, 2021.
  * Zhang et al. [2011] K. Zhang, J. Peters, D. Janzing, and B. Schölkopf.  Kernel-based conditional independence test and application in causal discovery.  In _Conference on Uncertainty in Artificial Intelligence_ , pages 804–813, 2011.

Appendix

##  A Broader Impact and Limitations

In many safety-critical applications, one needs to infer causal structure from
observational data e.g. in healthcare, economics, or geosciences since
interventions are unethical or simply impossible. For these cases, one
requires algorithms that have theoretical guarantees, are robust to hidden
confounding, and are computationally efficient. We hope that in the long term,
our work can thus contribute to better machine learning-driven decision-making
in safety-critical environments.

While theoretical guarantees and extensive experiments are critical components
for the evaluation of algorithms, especially in safety-critical environments
this can potentially lead to a false sense of security and trust. Constant
monitoring and assessment by domain experts are needed before and during the
deployment of any machine learning algorithm, especially in safety-critical
environments.

One key limitation of our approach as well as of all causal discovery
approaches from observational data is the reliance on additional assumptions.
These assumptions are required since causal discovery from observational data
is impossible without them. We related our assumptions to previous ones in the
Framework Section 3. These assumptions of our approach need to be checked
before applying the approach but are less restrictive than for comparable
baselines. Moreover, since our approach is inspired by debiased machine
learning, the focus is on inferring the true or "debiased" underlying
structure rather than obtaining low variance estimates of the structure.

##  B Direct Causal Effects and Interventions

In this section, we clarify the notion of an intervention and direct causal
effects. We will also introduce some notation that will be used later on in
the proofs.

We consider interventions by which a random variable $X^{j}_{t}$ is set to a
constant $X^{j}_{t}\leftarrow x$. We denote with $Y_{T}\mid do(X^{i}_{t}=x)$
the outcome time series $\boldsymbol{Y}$ at time step $T$, after performing an
intervention as described above. We can likewise perform multiple joint
interventions, by setting a group of random variables $\boldsymbol{I}$ at
different time steps, to pre-determined constants specified by an array
$\boldsymbol{i}$. We use the symbol $Y_{T}\mid
do(\boldsymbol{I}=\boldsymbol{i})$ to denote the resulting post-interventional
outcome, and we denote with $\mathbb{P}\left(Y_{T}=y\mid
do(\boldsymbol{I}=\boldsymbol{i})\right)$ the probability of the event
$\\{Y_{T}\mid do(\boldsymbol{I}=\boldsymbol{i})=y\\}$.

Using this notation, a time series $\boldsymbol{X}^{i}$ has a direct effect on
the outcome $\boldsymbol{Y}$, if performing different interventions on the
variables $\boldsymbol{X}^{i}$, while keeping the remaining variables fixed,
will alter the probability distribution of the outcome $\boldsymbol{Y}$.
Formally, define the sets of random variables
$\boldsymbol{I}_{T}\coloneqq\\{X_{t}^{1},\dots,X_{t}^{n},Y_{t}\\}_{t<T}$,
which consists of all the information before time step $T$. Similarly, define
the random variable $\boldsymbol{X}^{i}_{T}\coloneqq\\{X^{i}_{t}\\}_{t<T}$,
consisting of all the information of time series $\boldsymbol{X}^{i}$ before
time step $T$. Define the variable $\boldsymbol{I}_{T}^{\setminus
i}\coloneqq\boldsymbol{I}_{T}\setminus\boldsymbol{X}^{i}_{T}$, which consists
of all the variables in $\boldsymbol{I}_{T}$ except for
$\boldsymbol{X}^{i}_{T}$. Then, a time series $\boldsymbol{X}^{i}$ has a
direct effect on the outcome $\boldsymbol{Y}$ if it holds

$\mathbb{P}\left(Y_{T}=y\mid
do\left(\boldsymbol{X}^{i}_{T}=\boldsymbol{x}^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)\right)\neq\mathbb{P}\left(Y_{T}=y\mid
do\left(\boldsymbol{X}^{i}_{T}=\boldsymbol{x}^{\prime\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)\right)$ (5)

We say that a time series $\boldsymbol{X}^{i}$ causes $\boldsymbol{Y}$, if
there is a direct effect between $\boldsymbol{X}^{i}$ and $\boldsymbol{Y}$ as
in Eq. 5, for any time step $T$.

##  C Necessity of the Statistical Independence of $\varepsilon_{t}$

We provide a counterexample, to show that if there are dependencies between
the noise and the historical data, then the causal structure may not be
identifiable from observational data. To this end, we consider a first dataset
$\\{X_{t},Y_{t}\\}_{t\in\mathbb{Z}}$, defined as

$\begin{array}[]{l}X_{t-1},Y_{t}\sim\mathcal{N}(\boldsymbol{0},\Sigma)\end{array}$

Here, $\mathcal{N}(\boldsymbol{0},\Sigma)$ is a zero-mean joint Gaussian
distribution with covariance matrix

$\Sigma=\left[\begin{array}[]{cc}1&1\\\ 1&1\end{array}\right].$

We also consider a second dataset $\\{W_{t},Z_{t}\\}_{t\in\mathbb{Z}}$,
defined as

$W_{t}\sim\mathcal{N}(0,1),\quad Z_{t}=W_{t-1}$

The parameter $\Sigma$ is defined as above. Both datasets entail the same
joint probability distribution. However, the causal diagrams change from one
dataset to the other. Hence, the causal structure cannot be recovered from
observational data, if the posterior additive noise $\varepsilon_{t}$ is
correlated with some of the covariates.

##  D Necessity of No Instantaneous Causal Effects
between $\boldsymbol{Y}$ and the Potential Causes $\boldsymbol{X}^{i}$

Here, we provide a counterexample to show that without the no instantaneous
causal effect, the causal structure may not be identifiable from observational
data. Consider the following two models:

  * •

Model 1: We consider time series $\\{X_{t}\\}$, $\\{Y_{t}\\}$ of the form
$X_{t}=\mathbb{E}[X_{t-1}]+c$ and $Y_{t}=\mathbb{E}[Y_{t-1}-X_{t-1}]+X_{t}$.
In this model, $c$ is a random variable drawn from a Gaussian distribution
with a mean of 0 and covariance of 1.

  * •

Model 2: We consider time series $\\{X_{t}\\}$, $\\{Y_{t}\\}$ of the form
$Y_{t}=\mathbb{E}[Y_{t-1}]+c$ and $X_{t}=\mathbb{E}[X_{t-1}-Y_{t-1}]+Y_{t}$.
In this model, $c$ is a random variable drawn from a Gaussian distribution
with a mean of 0 and covariance of 1.

Both models entail the same joint distribution. However, in Model 1 $X$ has a
causal effect on $Y$, whereas in Model 2 $Y$ has a causal effect on $X$.
Hence, in this example, the causal structure is not identifiable.

##  E Proof of Theorem 4.1

We prove the following result. See 4.1

###### Proof.

We first prove that it holds

$\mathbb{P}\left(Y_{T}=y\mid do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right)=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$ (6)

for any non-zero event $\\{Y_{T}=y\\}$. To this end, define the group
$\boldsymbol{P}$ consisting of all the causal parents of the outcome. Note
that $\boldsymbol{P}\subseteq\\{X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\\}$. By Axiom (A), the outcome can be described as
$Y=f(\boldsymbol{P})+\varepsilon$, where $\varepsilon$ is independent of
$\\{X^{i}_{t},\boldsymbol{I}_{T}^{\setminus i}\\}$. Hence by Rule 2 of the do-
calculus (see Pearl [2000], page 85) Eq. 6 holds, since $Y$ becomes
independent of $\\{X^{i}_{t},\boldsymbol{I}_{T}^{\setminus i}\\}$ once all
arrows from $\boldsymbol{P}$ to $Y$ are removed from the graph of the DGP.

We now prove the claim using Eq. 6. To this end, assume that Eq. 6 holds and
suppose that $X^{i}$ does not Granger causes $Y$, i.e., it holds

$\mathbb{P}\left(Y_{T}=y\mid X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)=\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right),$ (7)

for any non-zero event $\\{Y_{T}=y\\}$. Then,

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i})\right)$
$\displaystyle=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
$\displaystyle[\text{Eq. \ref{eq:2_new}}]$
$\displaystyle=\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)$ $\displaystyle[\text{Eq. \ref{eq:granger}}]$
$\displaystyle=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
$\displaystyle[\text{Eq. \ref{eq:granger}}]$
$\displaystyle=\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right).$ $\displaystyle[\text{Eq. \ref{eq:2_new}}]$

Hence, causality implies Granger causality.

We now prove that Granger causality implies causality. To this end, suppose
that $X^{i}$ is not a potential cause of $Y$. By the definition of direct
effects, it holds

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i})\right)$
$\displaystyle\qquad\qquad\qquad=\mathbb{E}\left[\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right)\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right].$ (8)

Hence,

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i})\right)$
$\displaystyle\qquad\qquad=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
$\displaystyle[\text{Eq. \ref{eq:2_new}}]$
$\displaystyle\qquad\qquad=\mathbb{E}\left[\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i})\right)\mid\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
$\displaystyle[\text{Eq. \ref{eq:lemma101}}]$
$\displaystyle\qquad\qquad=\mathbb{E}\left[\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right)\mid\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
$\displaystyle[\text{Eq. \ref{eq:2_new}}]$
$\displaystyle\qquad\qquad=\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right),$

and the claim follows. ∎

##  F Proof of Theorem 4.2

We prove the following result. See 4.2 In order to prove Theorem 4.2, we use
the following auxiliary lemma.

######  Lemma F.1.

Consider a causal model as in Axiom (A)-Axiom (C). Then, the following
conditions are equivalent:

  1. 1.

$\mathbb{E}\left[Y_{T}\mid X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right]=\mathbb{E}\left[Y_{T}\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
a.s. ;

  2. 2.

$\mathbb{P}\left(Y_{T}=y\mid X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right)=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
a.s.

###### Proof.

Clearly, Item 2 implies Item 1.

We now prove the converse, i.e., we show that Item 1 implies Item 2. To this
end, define the group $\boldsymbol{P}_{T}$ consisting of all the causal
parents of $Y_{T}$. Note that it holds
$\boldsymbol{P}_{T}\subseteq\\{\boldsymbol{I}_{T}^{\setminus
i},X^{i}_{t}\\}\subseteq\\{\boldsymbol{I}_{T}^{\setminus i},X^{i}_{t}\\}$.
Hence, the joint intervention
$\\{X^{i}_{t},\boldsymbol{i}^{i,t}_{T-1}\\}\leftarrow\\{x,\boldsymbol{i}\\}$
define an intervention on the parents
$\boldsymbol{P}_{T}\leftarrow\boldsymbol{p}$. Further, we can write the
potential outcome as

$Y_{T}\mid do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})=f(\boldsymbol{p})+\varepsilon.$ (9)

Similarly, the joint intervention
$\\{X^{i}_{t},\boldsymbol{I}^{i,t}_{T-1}\\}\leftarrow\\{x,\boldsymbol{i}\\}$,
define an intervention on the parents
$\boldsymbol{P}_{T}\leftarrow\boldsymbol{p}^{\prime}$. We can write the
potential outcome as

$Y_{T}\mid do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})=f(\boldsymbol{p}^{\prime})+\varepsilon.$ (10)

Hence, it holds

$\displaystyle f(\boldsymbol{p})+\mathbb{E}\left[\varepsilon\right]$
$\displaystyle=\mathbb{E}\left[Y_{T}\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i})\right]$
$\displaystyle[\text{Eq. \ref{eq:new_eq1}}]$
$\displaystyle=\mathbb{E}\left[Y_{T}\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
$\displaystyle[\text{Eq. \ref{eq:2_new}, Theorem \ref{thm:granger}}]$
$\displaystyle=\mathbb{E}\left[Y_{T}\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
$\displaystyle[\text{by assumption}]$ $\displaystyle=\mathbb{E}\left[Y_{T}\mid
do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right]$ $\displaystyle[\text{Eq. \ref{eq:2_new}, Theorem
\ref{thm:granger}}]$
$\displaystyle=f(\boldsymbol{p}^{\prime})+\mathbb{E}\left[\varepsilon\right].$
$\displaystyle[\text{Eq. \ref{eq:new_eq2}}]$

By Axiom (A), the variable $\varepsilon$ is exogenous independent noise. From
the chain of equations above it follows that
$f(\boldsymbol{p})=f(\boldsymbol{p}^{\prime})$. Hence,

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right)=\mathbb{P}\left(f(\boldsymbol{p})+\varepsilon=y\right)$
$\displaystyle\qquad\qquad\qquad=\mathbb{P}\left(f(\boldsymbol{p}^{\prime})+\varepsilon=y\right)=\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right)$ (11)

We conclude that it holds

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
$\displaystyle\qquad=\mathbb{P}\left(Y_{T}\mid
do(X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i})\right)$
$\displaystyle[\text{Eq. \ref{eq:2_new}, Theorem \ref{thm:granger}}]$
$\displaystyle\qquad=\mathbb{P}\left(Y_{T}=y\mid
do(X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i})\right)$ $\displaystyle[\text{Eq. \ref{eq:new_3}}]$
$\displaystyle\qquad=\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right),$
$\displaystyle[\text{Eq. \ref{eq:2_new}, Theorem \ref{thm:granger}}]$

as claimed. ∎

We can now prove the main result.

###### Proof of Theorem 4.2.

We first prove that $X^{i}$ Granger causes $Y$ iff. it holds

$\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]-\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)^{2}\right]\neq 0.$ (12)

First, suppose that Eq. 12 does not hold. Then, it holds $\mathbb{e}[Y_{T}\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}]=\mathbb{E}[Y_{T}\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}]$, a.s.
Combining this equation with Lemma F.1 yields

$\displaystyle\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right)$
$\displaystyle=\mathbb{E}\left[\mathbb{P}\left(Y_{T}=y\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right)\mid\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right]$
$\displaystyle=\mathbb{P}\left(Y_{T}=y\mid\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right),$

a.s. Hence, if $X^{i}$ Granger causes $Y$, then Eq. 12 holds.

$\displaystyle\mathbb{E}\left[Y_{T}\mid
X^{i}_{t}=x,\boldsymbol{I}_{T}^{\setminus
i}=\boldsymbol{i}\right]\neq\mathbb{E}\left[Y_{T}\mid
X^{i}_{t}=x^{\prime},\boldsymbol{I}_{T}^{\setminus i}=\boldsymbol{i}\right],$
(13)

for a triple $\\{x,x^{\prime},\mathbf{w}\\}$. By combining Eq. 13 with Lemma
F.1 we see that Eq. 12 implies causality. However, by Theorem 4.1 Granger
causality is equivalent to causality in this case.

We now prove the claim. By the tower property of the expectation Williams
[1991] that

$\displaystyle\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]-\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]-\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]-\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)^{2}\mid\boldsymbol{I}_{T}^{\setminus i}\right]\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]^{2}-\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)\mid\boldsymbol{I}_{T}^{\setminus i}\right]\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]^{2}\right]-\mathbb{E}\left[\mathbb{E}\left[\left(\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right)\mid\boldsymbol{I}_{T}^{\setminus i}\right]\right]$
$\displaystyle=\mathbb{E}\left[\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]^{2}\right]-\mathbb{E}\left[\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]^{2}\right]$
$\displaystyle=\mathbb{E}\left[Y_{T}\mathbb{E}\left[Y_{T}\mid
X^{i}_{t},\boldsymbol{I}_{T}^{\setminus
i}\right]\right]-\mathbb{E}\left[Y_{T}\mathbb{E}\left[Y_{T}\mid\boldsymbol{I}_{T}^{\setminus
i}\right]\right],$

as claimed. ∎

##  G Intuition on Zero-Masking

We provide intuition why masking is a reasonably good idea. Assume that the
function $\widehat{f}$ is a $\epsilon$-close estimator of the true function
$f^{*}$ in the $\mathcal{L}^{2}(P_{X_{1},X_{2},\dots,X_{m}})$ norm, where the
functions $\widehat{f},f^{*}$ and the joint probability distribution
$P_{X_{2},\dots,X_{m}}$ are defined on the set
$\\{X_{1},X_{2},\dots,X_{m}\\}$,

$||\widehat{f}-f^{*}||_{\mathcal{L}^{2}(P_{X_{1},X_{2},\dots,X_{m}})}\leq\epsilon$

Now, let’s mask the random variable $X_{1}$. We are interested to see how
close is the estimator $\mathbb{E}_{X_{1}}\widehat{f}$ to the true function
$\mathbb{E}_{X_{1}}f^{*}$ in the $\mathcal{L}^{2}(P_{X_{2},\dots,X_{m}})$
norm, where the functions
$\mathbb{E}_{X_{1}}\widehat{f},\mathbb{E}_{X_{1}}f^{*}$ and the marginal
probability distribution $P_{X_{2},\dots,X_{m}}$ are defined on the rest of
variables $\\{X_{2},\dots,X_{m}\\}$, By Jensen’s inequality, we infer that for
any realization of $X_{2}=x_{2},X_{3}=x_{3},\dots,X_{m}=x_{m}$,

$\displaystyle(\mathbb{E}_{X_{1}}\widehat{f}(X_{1},x_{2},\dots,x_{m})$
$\displaystyle-\mathbb{E}_{X_{1}}f^{*}(X_{1},x_{2},\dots,x_{m}))^{2}$
$\displaystyle\leq\mathbb{E}_{X_{1}}[(\widehat{f}(X_{1},x_{2},\dots,x_{m})-f^{*}(X_{1},x_{2},\dots,x_{m}))^{2}]$

Plugging it in the $\epsilon$-closeness assumption leads to,

$||\mathbb{E}_{X_{1}}\widehat{f}-\mathbb{E}_{X_{1}}f^{*}||_{\mathcal{L}^{2}(P_{X_{2},\dots,X_{m}})}\leq||\widehat{f}-f^{*}||_{\mathcal{L}^{2}(P_{X_{1},X_{2},\dots,X_{m}})}\leq\epsilon,$

which guarantees that $\mathbb{E}_{X_{1}}\widehat{f}$ is also $\epsilon$-closs
to $\mathbb{E}_{X_{1}}f^{*}$ and hence it’s a good estimator. In the sequel, a
natural solution would be to estimate $\mathbb{E}_{X_{1}}\widehat{f}$ by
taking averages of $\widehat{f}$ over different samples of $X_{1}$. However,
for the linear regression problem that the estimator has a linear structure of
the input, it is straightforward to show that it is enough to evaluate
$\widehat{f}$ at $\mathbb{E}[X_{1}]$. And finally due to the zero-centering
step of data preprocessing, $\mathbb{E}[X_{1}]=0$. Thus, the aforementioned
procedure is equivalent to zero-masking.

##  H A Note on the Number of Partitions

The number of partitions $k$ affects the performance of our algorithm in
practice since a larger number of partitions will help in removing a bias in
the estimates. However, in our experiments, we observe that a small number of
partitions is sufficient to achieve good results. Furthermore, an excessive
number of random partitions may have a detrimental effect, due to the possible
small number of samples in each partition. Hence, we believe that the number
of partitions will not drastically affect performance in practice. Reasonable
choices of $k$ for our experiments range between 3-7, hence $k=\mathcal{O}(1)$
w.r.t. parameters of the problem. Thus, the resulting runtime can be reported
as $\mathcal{O}(md)$.

##  I Computational Complexity and Comparison

As discussed in Section 2, compared to SITD, conditional independence-based
approaches such as PCMCI [Runge et al., 2019b], PCMCI+ [Runge, 2020], and
LPCMCI [Gerhardus and Runge, 2020] face exponential computational barriers. It
is widely known that even endowed with a perfect infinite sample independence
testing oracle, learning Bayesian Networks becomes NP-Hard [Chickering et al.,
2004, Chickering, 1996]. Consequently, computational challenges arise not only
due to the nature of the conditional independence tests themselves but also
from the computational intractability of searching through the exponentially
large space of possible network structures. Hence, the runtime of the
$\mathcal{O}(m)$ number of regressions that SITD demands is negligible
compared to the exponential number of conditional independence tests from
lengthy time-series. To support this argument in practice, we provide a
runtime comparison between SITD and PCMCI+ w.r.t. the number of nodes $m$ in
Table 2.

| 10 | 20 | 30 | 40 | 50 | 100 | 200 | 400
---|---|---|---|---|---|---|---|---
SITD | $42\pm 13$ | $35\pm 6$ | $29\pm 1$ | $30\pm 4$ | $35\pm 15$ | $41\pm 6$ | $62\pm 6$ | $77\pm 10$
PCMCI+ | $4.2\pm 0.4$ | $18.6\pm 0.8$ | $48.6\pm 1.4$ | $99.2\pm 10$ | $216.4\pm 42.2$ | $1091\pm 65$ | $5678\pm 264$ |  $\approx 8$ hours

Table 2: Table with runtime means and standard deviations for SITD and PCMCI+
(in seconds).

##  J Additional Experiments

###  J.1 Synthetic Experiments

#### Dataset generation.

We first set the number of potential causes $m$, and we fix the lag $\Delta$
and the number of time steps for the dataset $T$. We create a covariate
adjacency 3-dimensional tensor $\Sigma$ of dimensions $\Delta\times m\times
m$. This tensor has $0$-$1$ coefficients, where $\Sigma_{k,i,j}=1$ if
$X^{i}_{t-k}$ has a casual effect on $X^{j}_{t}$ for all time steps $t$.
Similarly, we create an adjacency matrix $\Sigma^{\boldsymbol{Y}}$ for the
outcome $\boldsymbol{Y}$. $\Sigma^{\boldsymbol{Y}}$ is a binary $\Delta\times
m$ array, such that $\Sigma^{\boldsymbol{Y}}_{k,j}=1$ if $X_{t-k}^{j}$ has a
causal effect on $Y_{t}$ for all time steps $t$. The entries of $\Sigma$ and
$\Sigma^{\boldsymbol{Y}}$ follow the Bernouli distribution with parameter
$p=0.5$. Note that the resulting casual structure fulfills Axiom (B)-Axiom
(D).

We then create $m$ transformations that are used to produce the potential
causes $\boldsymbol{X}^{1},\dots,\boldsymbol{X}^{m}$. Each one of these
transformations is modeled by an MLP with $1$ hidden layer and $200$ hidden
units. We use Tanh nonlinearities (included also in the output layer) in order
to control the scale of the values. The final output value is further scaled
up so that all transforms generate values in the range $[-10,10]$. The input
layer of each MLP is coming from the corresponding causal parents of the
corresponding time series, as calculated from $\Sigma$.

In order to generate the potential causes $\boldsymbol{X}^{i}$, we use
$\Sigma$ and the MLP transforms. The first value of each time series is
randomly generated from a uniform distribution in $[-10,10]$. Then, each
$\boldsymbol{X}^{i}$ is produced by applying the appropriate transform to its
causal parents, as determined by $\Sigma$, and a zero-mean unit variance
Gaussian noise is added. We generate the target time series in a similar
fashion. Each variable $Y_{t}$ is produced by applying the MLP transform to
its causal parents, as determined by the target adjacency matrix
$\Sigma^{\boldsymbol{Y}}$. We then add zero-mean Gaussian noise. The scale of
the posterior additive noise for the outcome is referred to as the _noise-to-
signal ratio_ (NTS).

#### Description of the experiments.

We are given a dataset as described above with $m$ potential causes and a
fixed NTS for the generation of the outcome $\boldsymbol{Y}$. We determine
which series $\boldsymbol{X}^{1},\dots,\boldsymbol{X}^{n}$ are the causal
parents of the outcome, using Algorithm 1. For a given choice of $m$ and NTS,
we repeat the experiments five times, and we report on the sample mean for the
Accuracy (see Table 3). We additionally report on the F1 Score and CSI Score
in Table 5-4 in the Appendix. This experiment is repeated for an increasing
number of potential causes, and increasing noise-to-signal ratio, to evaluate
the performance of Algorithm 1 on challenging datasets.

In this set of experiments, we learn $\eta_{j}^{i}$ as in Line 3 of Algorithm
1, using for the regression task an MLP model. We found that this simple
approach, combined with zero-masking, dramatically reduces the false positives
of the Student’s t-test in Line 13 of Algorithm 1.

#### Results.

Table 3 shows the performance of our method with respect to the accuracy, for
an increasing number of potential causes and increasing NSR. We observe that
our method maintains a good score. These results are partially confirmed when
we look at the F1 and CSI scores (see Table 5-4 in the Appendix). In fact, we
see that the SITD is stable for an increasingly higher posterior noise.

Figure 2: AUROC metric for RHINO and SITD for various noise levels on the
synthetic dataset.

In Figure 2 we plot the AUROC performance of SITD against RHINO on the
synthetic dataset, confirming the competitive performance of STID. In order to
calculate the AUROC for SITD, we sort our predictions (for existence of an
edge) on the standard deviation of
$Z:=Y_{T}\cdot\tilde{g}_{j}^{i}+\tilde{\alpha}_{j}^{i}\cdot(Y_{T}-\tilde{g}_{j}^{i})-Y_{T}\cdot
g_{j}^{0}-\alpha_{j}^{0}\cdot(Y_{T}-g_{j}^{0})$, which is simply the
difference of the doubly robust statistics for $\theta^{i}$ and $\theta^{0}$
for a datapoint in partition $D_{j}$.

Table 3: Accuracy of our method for increasing number of potential causes $m$, and different noise-to-signal ration (NSR). We observe that our method maintains good accuracy, even in challenging settings with many potential causes and high noise. | Accuracy
---|---
$\boldsymbol{m}$ | $\boldsymbol{\textbf{{NSR}}=0}$ | $\boldsymbol{\textbf{{NSR}}=0.05}$ | $\boldsymbol{\textbf{{NSR}}=0.1}$ | $\boldsymbol{\textbf{{NSR}}=0.15}$ | $\boldsymbol{\textbf{{NSR}}=0.2}$ | $\boldsymbol{\textbf{{NSR}}=0.25}$ | $\boldsymbol{\textbf{{NSR}}=0.3}$
$\boldsymbol{5}$ |  0.60 $\pm$ 0.09 |  0.74 $\pm$ 0.25 |  0.66 $\pm$ 0.19 |  0.90 $\pm$ 0.11 |  0.80 $\pm$ 0.06 |  0.82 $\pm$ 0.10 |  0.94 $\pm$ 0.12
$\boldsymbol{10}$ |  0.99 $\pm$ 0.02 |  0.92 $\pm$ 0.08 |  0.97 $\pm$ 0.04 |  0.79 $\pm$ 0.18 |  0.89 $\pm$ 0.06 |  0.94 $\pm$ 0.04 |  0.90 $\pm$ 0.08
$\boldsymbol{15}$ |  0.89 $\pm$ 0.05 |  0.89 $\pm$ 0.05 |  0.88 $\pm$ 0.10 |  0.85 $\pm$ 0.02 |  0.79 $\pm$ 0.11 |  0.79 $\pm$ 0.09 |  0.81 $\pm$ 0.03
$\boldsymbol{20}$ |  0.83 $\pm$ 0.07 |  0.73 $\pm$ 0.05 |  0.72 $\pm$ 0.07 |  0.77 $\pm$ 0.10 |  0.73 $\pm$ 0.05 |  0.75 $\pm$ 0.05 |  0.69 $\pm$ 0.06
$\boldsymbol{25}$ |  0.76 $\pm$ 0.04 |  0.72 $\pm$ 0.10 |  0.71 $\pm$ 0.06 |  0.63 $\pm$ 0.08 |  0.68 $\pm$ 0.04 |  0.71 $\pm$ 0.07 |  0.66 $\pm$ 0.04
$\boldsymbol{30}$ |  0.76 $\pm$ 0.04 |  0.74 $\pm$ 0.04 |  0.72 $\pm$ 0.07 |  0.70 $\pm$ 0.10 |  0.66 $\pm$ 0.05 |  0.68 $\pm$ 0.04 |  0.64 $\pm$ 0.07
$\boldsymbol{35}$ |  0.68 $\pm$ 0.02 |  0.65 $\pm$ 0.07 |  0.72 $\pm$ 0.06 |  0.65 $\pm$ 0.03 |  0.66 $\pm$ 0.05 |  0.61 $\pm$ 0.07 |  0.64 $\pm$ 0.09
$\boldsymbol{40}$ |  0.64 $\pm$ 0.03 |  0.69 $\pm$ 0.05 |  0.67 $\pm$ 0.05 |  0.62 $\pm$ 0.03 |  0.63 $\pm$ 0.06 |  0.58 $\pm$ 0.07 |  0.60 $\pm$ 0.07
$\boldsymbol{45}$ |  0.65 $\pm$ 0.04 |  0.68 $\pm$ 0.08 |  0.58 $\pm$ 0.04 |  0.61 $\pm$ 0.03 |  0.64 $\pm$ 0.04 |  0.59 $\pm$ 0.04 |  0.61 $\pm$ 0.06
$\boldsymbol{50}$ |  0.68 $\pm$ 0.05 |  0.63 $\pm$ 0.05 |  0.64 $\pm$ 0.06 |  0.63 $\pm$ 0.05 |  0.66 $\pm$ 0.05 |  0.59 $\pm$ 0.08 |  0.64 $\pm$ 0.07

We provide the CSI Score in Table 4 and the F1 Score in Table 5 for the
experiments on synthetic datasets. Recall that in this set of experiments we
considered synthetic datasets with varying numbers of potential causes ($m$)
and noise-to-signal ratio (NSR). Increasing $m$ and NSR give more challenging
settings. We observe that our method maintains good CSI and F1 scores, for
increasing NSR.

Table 4: CSI Score of Algorithm 1 for increasing number of potential causes $m$, and different noise-to-signal ration (NSR). Again, we observe that our method is robust to increasing NSR. | CSI Score
---|---
$\boldsymbol{m}$ | $\boldsymbol{\textbf{{NSR}}=0}$ | $\boldsymbol{\textbf{{NSR}}=0.05}$ | $\boldsymbol{\textbf{{NSR}}=0.1}$ | $\boldsymbol{\textbf{{NSR}}=0.15}$ | $\boldsymbol{\textbf{{NSR}}=0.2}$ | $\boldsymbol{\textbf{{NSR}}=0.25}$ | $\boldsymbol{\textbf{{NSR}}=0.3}$
$\boldsymbol{5}$ |  0.57 $\pm$ 0.07 |  0.71 $\pm$ 0.28 |  0.63 $\pm$ 0.19 |  0.86 $\pm$ 0.12 |  0.69 $\pm$ 0.10 |  0.73 $\pm$ 0.12 |  0.91 $\pm$ 0.17
$\boldsymbol{10}$ |  0.98 $\pm$ 0.04 |  0.86 $\pm$ 0.14 |  0.95 $\pm$ 0.06 |  0.69 $\pm$ 0.21 |  0.80 $\pm$ 0.09 |  0.87 $\pm$ 0.07 |  0.82 $\pm$ 0.14
$\boldsymbol{15}$ |  0.80 $\pm$ 0.07 |  0.77 $\pm$ 0.11 |  0.75 $\pm$ 0.16 |  0.68 $\pm$ 0.06 |  0.59 $\pm$ 0.19 |  0.57 $\pm$ 0.15 |  0.59 $\pm$ 0.03
$\boldsymbol{20}$ |  0.66 $\pm$ 0.12 |  0.52 $\pm$ 0.12 |  0.46 $\pm$ 0.07 |  0.56 $\pm$ 0.10 |  0.46 $\pm$ 0.08 |  0.51 $\pm$ 0.07 |  0.39 $\pm$ 0.09
$\boldsymbol{25}$ |  0.51 $\pm$ 0.06 |  0.45 $\pm$ 0.12 |  0.43 $\pm$ 0.09 |  0.37 $\pm$ 0.09 |  0.38 $\pm$ 0.05 |  0.41 $\pm$ 0.08 |  0.33 $\pm$ 0.06
$\boldsymbol{30}$ |  0.47 $\pm$ 0.05 |  0.50 $\pm$ 0.07 |  0.46 $\pm$ 0.12 |  0.38 $\pm$ 0.08 |  0.35 $\pm$ 0.08 |  0.34 $\pm$ 0.05 |  0.31 $\pm$ 0.07
$\boldsymbol{35}$ |  0.42 $\pm$ 0.04 |  0.31 $\pm$ 0.06 |  0.39 $\pm$ 0.07 |  0.29 $\pm$ 0.07 |  0.30 $\pm$ 0.09 |  0.22 $\pm$ 0.09 |  0.26 $\pm$ 0.09
$\boldsymbol{40}$ |  0.32 $\pm$ 0.06 |  0.38 $\pm$ 0.05 |  0.33 $\pm$ 0.08 |  0.29 $\pm$ 0.06 |  0.24 $\pm$ 0.06 |  0.20 $\pm$ 0.07 |  0.19 $\pm$ 0.11
$\boldsymbol{45}$ |  0.33 $\pm$ 0.10 |  0.34 $\pm$ 0.07 |  0.20 $\pm$ 0.02 |  0.22 $\pm$ 0.05 |  0.25 $\pm$ 0.03 |  0.20 $\pm$ 0.05 |  0.19 $\pm$ 0.06
$\boldsymbol{50}$ |  0.33 $\pm$ 0.03 |  0.29 $\pm$ 0.06 |  0.29 $\pm$ 0.06 |  0.23 $\pm$ 0.04 |  0.26 $\pm$ 0.06 |  0.21 $\pm$ 0.08 |  0.26 $\pm$ 0.07
Table 5: F1 Score of the SITD for increasing number of potential causes $m$, and different noise-to-signal ration (NSR). Interestingly, our method maintains a good F1 score for increasing NSR. | F1 Score
---|---
$\boldsymbol{m}$ | $\boldsymbol{\textbf{{NSR}}=0}$ | $\boldsymbol{\textbf{{NSR}}=0.05}$ | $\boldsymbol{\textbf{{NSR}}=0.1}$ | $\boldsymbol{\textbf{{NSR}}=0.15}$ | $\boldsymbol{\textbf{{NSR}}=0.2}$ | $\boldsymbol{\textbf{{NSR}}=0.25}$ | $\boldsymbol{\textbf{{NSR}}=0.3}$
$\boldsymbol{5}$ |  0.73 $\pm$ 0.06 |  0.79 $\pm$ 0.20 |  0.75 $\pm$ 0.12 |  0.92 $\pm$ 0.07 |  0.81 $\pm$ 0.07 |  0.84 $\pm$ 0.08 |  0.95 $\pm$ 0.11
$\boldsymbol{10}$ |  0.99 $\pm$ 0.02 |  0.92 $\pm$ 0.09 |  0.98 $\pm$ 0.03 |  0.80 $\pm$ 0.15 |  0.89 $\pm$ 0.06 |  0.93 $\pm$ 0.04 |  0.90 $\pm$ 0.09
$\boldsymbol{15}$ |  0.89 $\pm$ 0.05 |  0.86 $\pm$ 0.08 |  0.85 $\pm$ 0.11 |  0.81 $\pm$ 0.04 |  0.73 $\pm$ 0.15 |  0.72 $\pm$ 0.12 |  0.74 $\pm$ 0.03
$\boldsymbol{20}$ |  0.79 $\pm$ 0.08 |  0.67 $\pm$ 0.10 |  0.62 $\pm$ 0.07 |  0.71 $\pm$ 0.09 |  0.63 $\pm$ 0.07 |  0.67 $\pm$ 0.06 |  0.56 $\pm$ 0.09
$\boldsymbol{25}$ |  0.68 $\pm$ 0.05 |  0.61 $\pm$ 0.11 |  0.60 $\pm$ 0.09 |  0.53 $\pm$ 0.11 |  0.55 $\pm$ 0.05 |  0.57 $\pm$ 0.08 |  0.49 $\pm$ 0.07
$\boldsymbol{30}$ |  0.64 $\pm$ 0.05 |  0.66 $\pm$ 0.06 |  0.62 $\pm$ 0.10 |  0.54 $\pm$ 0.09 |  0.51 $\pm$ 0.09 |  0.50 $\pm$ 0.06 |  0.46 $\pm$ 0.09
$\boldsymbol{35}$ |  0.59 $\pm$ 0.05 |  0.47 $\pm$ 0.07 |  0.56 $\pm$ 0.07 |  0.45 $\pm$ 0.08 |  0.46 $\pm$ 0.11 |  0.35 $\pm$ 0.11 |  0.41 $\pm$ 0.11
$\boldsymbol{40}$ |  0.49 $\pm$ 0.07 |  0.55 $\pm$ 0.05 |  0.49 $\pm$ 0.08 |  0.45 $\pm$ 0.07 |  0.38 $\pm$ 0.08 |  0.33 $\pm$ 0.09 |  0.30 $\pm$ 0.15
$\boldsymbol{45}$ |  0.49 $\pm$ 0.11 |  0.50 $\pm$ 0.08 |  0.34 $\pm$ 0.03 |  0.36 $\pm$ 0.06 |  0.40 $\pm$ 0.04 |  0.33 $\pm$ 0.06 |  0.32 $\pm$ 0.08
$\boldsymbol{50}$ |  0.50 $\pm$ 0.03 |  0.44 $\pm$ 0.07 |  0.44 $\pm$ 0.08 |  0.38 $\pm$ 0.05 |  0.41 $\pm$ 0.07 |  0.34 $\pm$ 0.11 |  0.40 $\pm$ 0.09

###  J.2 Performance in Low-Sample Regimes

Here we report the results for additional experiments in the practically
important low sample setting. In Figure 3, the plots are provided to support
Section 6.2. The double robustness property enables our algorithm to rely on
simple estimators with low statistical complexity. As a result, our method
shows more consistent performance in low-sample regimes as opposed to existing
approaches that are based on overparameterized models demanding so many data
points.

(a) Task: E.Coli 1

(b) Task: E.Coli 2

(c) Task: Yeast 1

(d) Task: Yeast 2

(e) Task: Yeast 3

Figure 3: This figure demonstrates the consistent performance of SITD w.r.t
number of observations compared to state-of-the-art methods Rhino and Rhino+g.
Note that Rhino and Rhino+g are built on neural networks. SITD significantly
outperforms Rhino and Rhino+g in E.Coli 1 and E.Coli 2 and shows competitive
results in Yeast 1. Thanks to the double robustness property of SITD, the
dependence of our algorithm on the estimator is much lower than the well-
established approaches. In this regard, SITD with a simple kernel regression
with polynomial kernels has superior performance compared to state-of-the-art
methods Rhino and Rhino+g. This superiority gets magnified in the low number
of observation regimes due to the high sample complexity required by Rhino and
Rhino+g.

###  J.3 Run Time and Hardware

Table 6: Run time and Hardware used for our method (SITD) and the state-of-the-art baseline Rhino. Category\Method | Rhino | SITD (ours)
---|---|---
Runtime |  18 mins 40 sec $\pm$ 30 sec |  57.12 sec $\pm$ 1.6 sec
Hardware | 1 NVIDIA A100 GPU + AMD EPYC 7402 24-Core CPU | 11th Gen Core i5-1140F CPU

To give a taste of the computational efficiency of our method, here we report
the average runtime of SITD and the competitive rival Rhino for experiment
Table 1 across all five tasks (E.Coli 1, E.Coli 2, Yeast 1, Yeast 2 and Yeast
3) in Table 6. Despite having access only to a single CPU in contrast to the
GPU-equipped execution of Rhino our method is almost 20x faster. This is
because of the fact that SITD algorithm will provide reasonable results even
when employing simple fast efficient estimators (in this case a kernel
regression).