Signature Kernel Conditional Independence Tests in Causal Discovery for Stochastic Processes

## Abstract

Abstract Inferring the causal structure underlying stochastic dynamical systems from observational data holds great promise in domains ranging from science and health to finance. Such processes can often be accurately modeled via stochastic differential equations (SDEs), which naturally imply causal relationships via ‘which variables enter the differential of which other variables’. In this paper, we develop conditional independence (CI) constraints on coordinate processes over selected intervals that are Markov with respect to the acyclic dependence graph (allowing self-loops) induced by a general SDE model. We then provide a sound and complete causal discovery algorithm, capable of handling both fully and partially observed data, and uniquely recovering the underlying or induced ancestral graph by exploiting time directionality assuming a CI oracle. Finally, to make our algorithm practically usable, we also propose a flexible, consistent signature kernel-based CI test to infer these constraints from data. We extensively benchmark the CI test in isolation and as part of our causal discovery algorithms, outperforming existing approaches in SDE models and beyond.

## 1 Introduction

Understanding cause-effect relationships from observational data can help identify causal drivers for disease progression in longitudinal data and aid the development of new treatments, act upon the underlying influences of stock prices to support lucrative trading strategies, or speed up scientific discovery by uncovering interactions in complex biological systems such as gene regulatory pathways.
Causal discovery (or causal structure learning) has received continued attention from the scientific community for at least two decades (Glymour et al., ; Spirtes et al., ; Vowels et al., ) with a notable uptick in previous years particularly regarding differentiable score-based methods (Zheng et al., ; Brouillard et al., ; Charpentier et al., ; Hägele et al., ; Lorch et al., ; Annadani et al., ; Zheng et al., ), building on score matching algorithms (Rolland et al., ; Montagna et al., ; ), and other deep-learning based approaches (Chen et al., ; Ke et al., ; Yu et al., ; Ke et al., ).
Many of these approaches aim at improving the scalability of causal discovery in the number of variables and observations as well as at incorporating uncertainty or efficiently making use of interventional data.

However, causal discovery from time series data has received much less attention and been mostly neglected in these recent advances.
At a fundamental level, causal effects in dynamical systems can only ‘point into the future’, which should make causal discovery in time resolved data intuitively simpler.
Nevertheless, except for restricted settings, this promise has not been realized methodologically and causal discovery in time series, especially for systems evolving in continuous time, remains a major challenge (Singer, ; Runge et al., ; Lawrence et al., ; Runge et al., ).

#### Data generating process.

We assume data to follow a stochastic process $X\coloneqq(X^{1},\ldots,X^{d})$, with $X^{k}$ taking values in $\mathbb{R}^{n_{k}}$ ($n_{k}\geq 1$), and $X$ satisfying the following system of stationary, path-dependent SDEs

$$ $\begin{cases}\mathrm{d}X^{k}_{t}=\mu^{k}(X_{[0,t]})\mathrm{d}t+\sigma^{k}(X_{[ 0,t]})\mathrm{d}W^{k}_{t},\\ X^{k}_{0}=x^{k}_{0}\qquad\text{for }k\in[d]:=\{1,\ldots,d\}\>.\end{cases}$ (1) $$

We call this the *SDE model*.
The subscript $[0,t]$ at $X$ means that $\mu^{k}$ (the ‘drift’) and $\sigma^{k}$ (the ‘diffusion’) are functions of the entire solution up to time $t$, i.e., they are defined on $C([0,+\infty),\mathbb{R}^{n})$ (or some suitable subspace thereof), where $n:=n_{1}+\ldots+n_{d}$.
Each $W^{k}$ is an $m_{k}$-dimensional Brownian motion, and $\sigma^{k}$ maps to $n_{k}\times m_{k}$-dimensional matrices.
The noises $W^{k}$ together with the (possibly random) initial conditions $x^{k}_{0}$ are jointly independent. We refer to Rogers & Williams () for details, Evans () for a concise introduction, and Section A.1 for more intuition.
Hence, the SDE together with a distribution over initial conditions defines our data generating process.
Individual observations are paths, which in practice translate to stochastic, potentially irregularly sampled time series observations for $X^{k}_{[0,T]}$ with a maximum observation time $T$, see Figure 1.

Figure: Figure 1: Illustration of causal discovery in the SDE model (A), leveraging conditional independencies in the observed samples (B) to infer the dependence graph of the SDE (C).
Refer to caption: x1.png

#### Induced causal graph.

Eq. 1 naturally implies cause-effect relationships:
we call $i\in[d]$ a parent of $j\in V:=[d]$ ($i\in\mathrm{pa}^{\mathcal{G}}_{j}$) when either $\mu^{j}$ or $\sigma^{j}$ is not constant in the $i$-th argument.
These parental relationships define a directed graph $\mathcal{G}=(V,E)$, which we call *dependence graph* of the SDE model.
The goal of causal discovery is to infer the graph $\mathcal{G}$ induced by the SDE model from a sample of observed solution paths, depicted in Figure 1.
Here, we only consider data generating processes leading to directed graphs without cycles of length greater than one (i.e., only ‘loops’ $X^{k}\to X^{k}$ are allowed).
We still call these directed acyclic graphs (DAGs) for simplicity.

#### Limitations and requirements.

The key limitation of our setting is the assumptions of acyclicity (except for self-loops).
We discuss fundamental issues arising from dropping the acyclicity constraint in Section A.2.
Further, we focus on stationary SDE models, i.e., the coefficients in eq. 1 are not explicitly time-dependent, motivated by the requirement that causal relationships do not change over time.
Despite these limitations, our setting far exceeds the scope of existing causal discovery methods from time series data:
(a) It does not rely on the ‘discrete-time’ assumption, i.e., that all variables are observed at the same, typically homogeneously spaced, time points (neither in the underlying model nor the actual observations). Instead, different variables can be observed at different, irregularly spaced time points in our setting.
(b) We can handle partially observed systems with unobserved confounders.
(c) Path-dependence ($\mu^{k},\sigma^{k}$ may depend on the entire previous history $X_{[0,t]}$ of paths), including delayed SDEs, is typically not captured in existing approaches (Peters et al., ).
(d) Beyond ‘additive observation noise’ it incorporates ‘driving noise’, where causal dependence may stem from the diffusion, not captured by the means of observed paths.

#### Contributions.

Our main contributions include (a) developing CI constraints on coordinate processes over selected intervals that we prove to be Markov with respect to the acyclic dependence graph (with self-loops) induced by a general SDE model.
(b) We then propose an efficient constraint based causal discovery algorithm and prove it to be sound and complete, assuming a CI-oracle, for both fully and partially observed data. Unlike existing constraint based methods in static settings, our algorithm uniquely recovers the full underlying graph or induced ancestral graph by exploiting the direction of time.
(c) Finally, we propose a flexible, practical CI test on path space, rendering our method highly practicable.
(d) We extensively demonstrate the test’s efficacy and achieve superior performance in causal discovery compared to existing methods.

## 2 Background and Related Work

#### Causal discovery on time series.

While virtually all causal discovery methods on time series exploit that the direction of time constrains edges to point forward (except for Laumann et al. ()), most existing work assumes observations to be sampled on a regular time grid with a discrete (auto-regressive) law, where past observations $(X_{t-\tau},\ldots,X_{t-1})$ determine the present $X_{t}=f(X_{t-\tau},\ldots,X_{t-1},\varepsilon_{t})$ for some fixed lag $\tau$, which can be seen as a static structural causal model (SCM) over an ‘unrolled graph’ with variables at different time steps considered as distinct nodes (Assaad et al., ; Hasan et al., ; Peters et al., ).
Based on the seminal work (Granger, ), various Granger-type approaches leverage the assumption that past values of a variable can help predict future values to detect lagged causal influences both for linear (Diks & Panchenko, ; Granger, ) and non-linear (Marinazzo et al., ; Shojaie & Fox, ; Runge, ; Pamfil et al., ) functional relationships.
While non-parametric constraint-based methods for such discrete-time models (Runge et al., ; ; ) can handle partial observations, they are also fundamentally limited by the ‘discrete-time’ assumption making them unfit for irregularly sampled observations of continuous time systems, path-dependence, or diffusion dependence.
Furthermore, they require numerous tests across all time points and rely on the correct estimation of a fixed ‘lookback window’ $\tau$.
The fundamental problems of the ‘discrete-time’ assumption have also been discussed in Runge ().

Other work focuses on the (static) equilibrium behavior of differential equations (Mooij et al., ; Bongers et al., ; Bongers & Mooij, ), exploiting invariance in mass-action kinetics (Peters et al., ), or attempts to directly learn the full dynamical law via non-convex optimization of heavily overparameterized neural network models (Aliee et al., ; ; Bellot et al., ; Wang et al., ).
The latter are fundamentally limited by the fixed functional modeling structures and therefore also unable to model partially observed settings.
Albeit its limited applicability to fully-observed, Markovian SDE models, SCOTCH (Wang et al., ), which uses a variational formulation to infer posterior distributions over possible graphs, is the current state-of-the-art baseline.

Constraint-based methods for continuous-time stochastic processes mostly rely on (conditional) local independence, an asymmetric independence relation (Schweder, ; Mogensen et al., ), which can be used to infer a (partial) causal graph in, e.g., point process models (Didelez, ; Meek, ; Mogensen et al., ; Mogensen & Hansen, ).
However, besides the limitation to drift-dependencies, there exists no practical test of local independence for diffusions.
The only work that actually tests local independence from data is (Christgau et al., ), which is heavily restricted to counting processes.

In this work, we overcome the challenges of irregular and partial observations, path- and diffusion dependence, and practical applicability by transforming CI statements of the form

$$ $X^{I}_{\mathcal{I}}\operatorname{\perp\perp}X^{J}_{\mathcal{J}}\mid X^{K}_{ \mathcal{K}},$ (2) $$

for disjoint $I,J,K\subseteq[d]$ between subsets of coordinate-processes $X^{I}$, $X^{J}$, $X^{K}$ over intervals $\mathcal{I},\mathcal{J},\mathcal{K}\subseteq[0,T]$ into a practical hypothesis test using the signature kernel and the arrow of time.

#### Conditional independence tests.

As a key building block of graphical models, there is a large literature on conditional independence tests, i.e., testing the null hypothesis $H_{0}:X\operatorname{\perp\perp}Y\mid Z$.
When $Z$ is discrete, CI testing reduces to a series of unconditional tests for each value of $Z$ (Tsamardinos & Borboudakis, ).
For continuous $Z$, nonparametric kernel-based methods often either check independence of residuals after kernel (ridge) regressing $X$ and $Y$ on $Z$ (Shah & Peters, ; Lundborg et al., ; Zhang et al., ; Daudin, ; Strobl et al., ), or measure the distance between joint and marginal kernel mean embeddings (Muandet et al., ; Park & Muandet, ).
For example, KCIPT (Doran et al., ) phrases it as a two-sample test of the null $H_{0}:P(X,Y,Z)=P(X\mid Z)P(Y\mid Z)P(Z)$ and simulates samples from $P(X\mid Z)P(Y\mid Z)P(Z)$ via a permutation of the data that approximately preserves $Z$. Lee & Honavar () improve on KCIPT with an unbiased estimate of Maximum Mean Discrepancy (MMD) (Gretton et al., ).
Other permutation-based approaches require large sample sizes (Sen et al., ) or rely on densities (Kim et al., ), unfit for our setting.

Laumann et al. () develop a CI test for functional data using the Hilbert-Schmidt conditional independence criterion (HSCIC) (Park & Muandet, ) and a permutation test (Berrett et al., ), which they use as part of the PC-algorithm (Glymour et al., ) for causal discovery.
Even though their approach assumes $P_{X\mid Z}$ to be known (not the case in our setting), is inapplicable under partial observations, and does not exploit the arrow of time (to identify graphs beyond the Markov equivalence class), we benchmark the work by Laumann et al. () against our method in the experiments.
CI testing faces fundamental practical challenges, such as the exponential growth in the number of values to condition on as dimensionality of $Z$ increases, and theoretical limits, as any CI test’s power is bounded by its size (Shah & Peters, ; Lundborg et al., ).
This ‘no free lunch’ statement underscores the need for carefully selecting the right test.
Most existing approaches assume Euclidean spaces and do not generalize to path-valued random variables without densities, a gap we fill by combining kernel-based permutation tests (KCIPT, SDCIT) with the signature kernel for a tailored CI test on path space and proving its consistency.

#### Signature kernels.

Signature kernels (Király & Oberhauser, ; Salvi et al., ), a universal class for sequential data, have received attention recently for their efficiency in handling path-dependent problems (Lemercier et al., ; Salvi et al., ; Cochrane et al., ; Salvi et al., ; Cirone et al., ; Issa et al., ; Pannier & Salvi, ).
The definition of the signature kernel requires an initial algebraic setup, which we keep as concise as possible – yet self-contained – here and refer the reader to Cass & Salvi (, Chapter 2) for more details.
Let $\langle\cdot,\cdot\rangle_{1}$ be the Euclidean inner product on $\mathbb{R}^{d}$.
Denote by $\otimes$ the standard outer product of vector spaces.
For any $n\in\mathbb{N}$, we denote by $\langle\cdot,\cdot\rangle_{n}$ on $(\mathbb{R}^{d})^{\otimes n}$ the canonical Hilbert-Schmidt inner product defined for any $a=(a_{1},\ldots,a_{n})$ and $b=(b_{1},\ldots,b_{n})$ in $(\mathbb{R}^{d})^{\otimes n}$ as $\left\langle a,b\right\rangle_{n}=\prod_{i=1}^{n}\left\langle a_{i},b_{i}
\right\rangle_{1}$.
Define the direct sum of vector spaces $T(\mathbb{R}^{d}):=\bigoplus_{n=1}^{\infty}(\mathbb{R}^{d})^{\otimes n}$, where it is understood that the direct sum runs over finitely many non-zero levels.
The inner product $\langle\cdot,\cdot\rangle_{n}$ on $(\mathbb{R}^{d})^{\otimes n}$ can then be extended by linearity to an inner product $\langle\cdot,\cdot\rangle$ on $T(\mathbb{R}^{d})$ defined for any $a=(a_{0},a_{1},\ldots)$ and $b=(b_{0},b_{1},\ldots)$ in $T(\mathbb{R}^{d})$ as $\langle a,b\rangle=\sum_{n=0}^{\infty}\left\langle a_{n},b_{n}\right\rangle_{n}$.
Other choices of linear extensions have been studied in Cass et al. ().
We denote by $\mathcal{T}(\mathbb{R}^{d})$ the Hilbert space obtained by completing $T(\mathbb{R}^{d})$ with respect to $\langle\cdot,\cdot\rangle$.

The *signature transform* is a classical path-transform from stochastic analysis. For any sub-interval $[s,t]\subset[0,T]$ and any continuous path $X\in C_{p}([0,T],\mathbb{R}^{d})$ of finite $p$-variation, with $1\leq p<3$, it is (canonically) defined as
$S(X)_{s,t}:=\bigl{(}1,S(X)^{(1)}_{s,t},\ldots,S(X)^{(n)}_{s,t},\ldots\bigr{)}
\in\mathcal{T}(\mathbb{R}^{d})$,
where $S(X)^{(n)}_{s,t}\in(\mathbb{R}^{d})^{\otimes n}$ is the $n$-fold iterated integral
$S(X)^{(n)}_{s,t}=\int_{s<u_{1}<\ldots<u_{n}<t}\mathrm{d}X_{u_{1}}\otimes
\mathrm{d}X_{u_{2}}\otimes\ldots\otimes\mathrm{d}X_{u_{n}}$. Given two arbitrary sub-intervals $[a,b],[c,d]\subset[0,T]$, the *signature kernel* $K_{S}:C_{p}([s,t],\mathbb{R}^{d})\times C_{p}([s^{\prime},t^{\prime}],\mathbb{
R}^{d})\to\mathbb{R}$ is a positive definite kernel on continuous paths of bounded variation defined as

$$ $K_{S}(X,Y)=\left\langle S(X)_{s,t},S(Y)_{s^{\prime},t^{\prime}}\right\rangle\>.$ (3) $$

Salvi et al. (, Thm. 2.5) shows that $K_{S}(X,Y)=f(t,t^{\prime})$, where $f:[s,t]\times[s^{\prime},t^{\prime}]\to\mathbb{R}$ is the solution of the following path-dependent integral equation:

$$ $f(t,t^{\prime})=1+\int_{s}^{t}\int_{s^{\prime}}^{t^{\prime}}f(u,v)\langle dX_{ u},dY_{v}\rangle_{1}\>,\quad\text{with }f(0,\cdot)=f(\cdot,0)=1\>.$ (4) $$

This ‘kernel trick’ allows us to evaluate the signature kernel without explicit computation of the signature transform by solving the partial differential equation (PDE) in eq. 4.
We refer to Salvi et al. () for a numerical approximation scheme to solve this hyperbolic PDE and its error rates and to Section A.3 for more mathematical details on the applicability in our setting.
In our experiments, we use the JAX library [sigkerax](https://github.com/crispitagorico/sigkerax) to efficiently solve eq. 4.

In summary, the signature kernel is a universal kernel that takes (multi-variate) continuous paths – potentially on different intervals – as input, is efficiently computable, and performs well in capturing characteristics of sequential data (Lee & Oberhauser, ).

## 3 Methodology

We start with our constraint-based causal discovery algorithms, where we assume an oracle for CI testing on path space for eq. 2 in Section 3.1.
In Section 3.2 we then introduce the signature kernel CI test and state its consistency, which requires a new proof that does not rely on densities.
Our test, applicable to arbitrary stochastic processes, time series, or functional data, performs well empirically in our extensive experiments in Section 4 making it a contribution of independent interest.

### 3.1 Causal Discovery in the Acyclic SDE Model

Throughout this section, we assume access to a CI oracle for eq. 2, where $X^{H}_{[a,b]}$ denotes the $C([a,b],\mathbb{R}^{|H|})$-valued random variable $\omega\mapsto([a,b]\ni t\mapsto X^{h}_{t}(\omega)^{h\in H})$.
Asymptotically, such an oracle can be replaced in practice by a consistent finite-sample CI test for path-valued random variables (see Section 3.2).
In Section A.2, we prove that when allowing for cycles in the SDE model, constraint-based causal discovery of the full graph is impossible using arbitrary expressions of the form in eq. 2 despite the flexibility of these conditional independencies.
This impossibility result is a key motivation to study the acyclic setting still allowing for loops, which are crucial in dynamic settings as variables should be allowed to depend on themselves infinitesimally into the past.

We will use the oracle in the following ways for intervals $[0,s]$ and $[s,s+h]$ with $h>0$:

- •
$X^{i}$ is *symmetrically CI* of $X^{j}$ given $X^{K}$ on $[0,T]$ if
$X^{i}_{[0,T]}\operatorname{\perp\perp}X^{j}_{[0,T]}\mid X^{K}_{[0,T]}$;
we then write $X^{i}\operatorname{\perp\perp}_{\mathrm{sym}}X^{j}\mid X^{K}$
- •
$X^{i}$ is *future-extended $h$-locally CI* of $X^{j}$ given $X^{K}$ at $s$ if
$X^{i}_{[0,s]}\operatorname{\perp\perp}X^{j}_{[s,s+h]}\mid X^{j}_{[0,s]},X^{K}_
{[0,s+h]}$;
(^††\dagger†^††\dagger†$\dagger$The ‘$+$’ denotes conditioning on the future of $X^{K}$, not common in the literature (Florens & Fougere, ).)
we then write $X^{i}\operatorname{\perp\perp}_{s,h}^{+}X^{j}\mid X^{K}$
- •
$X^{i}$ is *conditionally $h$-locally self-independent* given $X^{K}$ at $s$ if $X^{i}_{[0,s]}\operatorname{\perp\perp}X^{i}_{[s,s+h]}\mid X^{K}_{[0,s+h]}$;
we then write $X^{i}\operatorname{\perp\perp}^{\circlearrowleft}_{s,h}\mid X^{K}$

Figure: Figure 2: Illustration of the lifted dependence graph $\tilde{\mathcal{G}}$ for DAG $\mathcal{G}$ with colors highlighting the correspondence of some selected edges.
Refer to caption: x2.png

The key idea to leverage the unidirectional flow of time is to split the observed paths into ‘past’ and ‘future’ at some time $t\in[0,T]$.
Therefore, we define the *lifted dependence graph* $\tilde{\mathcal{G}}=(\tilde{V},\tilde{E})$ by setting $\tilde{V}\coloneqq V_{0}\sqcup V_{1}$, where $\sqcup$ denotes disjoint union and $V_{0},V_{1}$ are two copies of $V$, whose elements we subscript with 0 0 and $1$; we include an edge $(i_{0}\to i_{1})\in\tilde{E}$ if and only if there is a loop $(i\to i)\in E$, and for each edge $(i\to j)\in E$ with $i\neq j$, include edges $(i_{0}\to j_{0}),(i_{1}\to j_{1}),(i_{0}\to j_{1})\in\tilde{E}$, see Figure 2.
Conversely, we say that the graph with edges $E$ is obtained by *collapsing* the graph with edges $\tilde{E}$.
Intuitively, $V_{0}$ contains all ‘past’ and $V_{1}$ all ‘future’ variables with the same edges among them (causal relations) with additional edges from ‘past’ to ‘future’.
We can now establish the following Markov property.

###### Proposition 3.1 (Markov property) .

Assume the dependence graph $\mathcal{G}$ of the SDE model is acyclic except for loops and associate $X^{i}_{[0,s]}$ with $i_{0}\in V_{0}$ and $X^{i}_{[s,s+h]}$ with $i_{1}\in V_{1}$. Then $\tilde{\mathcal{G}}$ is acyclic and the independence relation $\operatorname{\perp\perp}_{s,h}^{+}$ satisfies the global Markov property w.r.t. d-separation in $\tilde{\mathcal{G}}$.

Our proof (see Section A.5) leverages the structure of the SDE-generating mechanism, which factorizes according to the dependence graph, thereby allowing us to formally establish the global Markov property using arguments from stochastic analysis.
With an oracle for eq. 2 and a faithfulness type assumption (for more details see Section A.7), we now develop new constraint-based causal discovery algorithms that infer the full dependence graph instead of just equivalence classes.

Figure: Algorithm 1 Causal discovery for acyclic SDEs.

###### Proposition 3.1 (Markov property) .

#### Discovering the full DAG 𝒢 𝒢 \mathcal{G} caligraphic_G .

To discover the full graph including loops we propose Algorithm 1 that makes use of the time-ordering via our $h$-local independence models (consider $s,h>0$ arbitrary but fixed).

###### Theorem 3.2 .

Algorithm 1 is sound and complete for the SDE model, assuming acyclicity except for loops and faithfulness:
its output is the true dependence graph $\mathcal{G}$.

The proof (see Section A.6) establishes completeness by identifying a separating set and leveraging the global Markov property in Proposition 3.1; soundness is shown via the faithfulness assumption.
In particular, our proof does not require the ‘strong faithfulness’ assumption, but a weaker version, ‘parent faithfulness’, suffices—see Section A.7 for a detailed discussion of faithfulness in our setting.

###### Theorem 3.2 .

#### Discovering and post-processing the CPDAG.

In Algorithm 1, we leverage our Markov property with respect to the lifted graph that also incorporates time directionality to infer the full graph.
We now show that we can also use the symmetric criterion $\operatorname{\perp\perp}_{\mathrm{sym}}$ to recover the ‘Markov equivalence class’—all graphs that are Markov equivalent to $\mathcal{G}$—via the completed partially directed acyclic graph (CPDAG) of $\mathcal{G}$ (Peters et al., , Def. 6.24).
Loosely speaking, the CPDAG has the same adjacencies as $\mathcal{G}$, but some edges may remain undirected.
Due to space limitations, we present the required global Markov property with respect to $\operatorname{\perp\perp}_{\mathrm{sym}}$ in Section A.8.
This Markov property allows us to apply the sound and complete PC algorithm to infer the CPDAG assuming an CI oracle for $\operatorname{\perp\perp}_{\mathrm{sym}}$ and faithfulness (Spirtes et al., ).
While Algorithm 1 is strictly more informative (returns $\mathcal{G}$ instead of just its CPDAG), the reliability of $\operatorname{\perp\perp}_{s,h}^{+}$ in practice might be negatively affected compared to $\operatorname{\perp\perp}_{\mathrm{sym}}$ by (a) conditioning on larger sets, (b) shorter time segments (potentially losing information), and (c) the additional choice of parameters $s,h$.
Hence, in real-world applications inferring the CPDAG using $\operatorname{\perp\perp}_{\mathrm{sym}}$ may be more robust than inferring $\mathcal{G}$ fully using $\operatorname{\perp\perp}_{s,h}^{+}$.
Building on these potential benefits, we also develop a post-processing procedure for the CPDAG that again leverages the directionality of time to provably also orient all remaining unoriented edges.

###### Corollary 3.3 (post-processing) .

Let $\mathcal{G}=(V,E)$ be the dependence graph of the acyclic (except for loops) SDE model and $\tilde{\mathcal{G}}=(V,\tilde{E})$ its CPDAG.
Then we have $X^{j}_{[0,T]}\perp\not\hskip 1.00006pt\perp X^{i}_{0}$ but $X^{i}_{[0,T]}\operatorname{\perp\perp}X^{j}_{0}$ for all $(i,j)\in E\subset\tilde{E}$ with $(j,i)\in\tilde{E}$ and $i\neq j$.

The proof in Section A.9 relies on the joint independence of Brownian motions $\mathrm{d}W^{k}_{t}$ and initial conditions $X_{0}^{k}$.
Corollary 3.3 directly motivates an alternative algorithm for causal discovery of $\mathcal{G}$, written out as Algorithm 2 in Section A.9, by first constructing the CPDAG (PC algorithm with $\operatorname{\perp\perp}_{\mathrm{sym}}$) followed by testing unconditionally the yet unoriented edges as in Corollary 3.3.

###### Corollary 3.3 (post-processing) .

#### Partially observed setting.

Finally, we consider the partially observed setting, where only $\{X^{v_{k}}_{t}\}_{v_{k}\in V_{obs},t\in[0,1]}$ of the process $\{X^{v_{i}}_{t}\}_{v_{i}\in V_{obs}\sqcup V_{L},t\in[0,1]}$ from the model in eq. 1 with DAG $\mathcal{G}=(V_{obs}\sqcup V_{L},E)$ are observed. Asymmetric local independence based methods for causal discovery are challenged by requiring exponentially many oracle calls in the partially observed setting (Mogensen & Hansen, ; Mogensen, ).
However, since our symmetric conditional independence constraint satisfies the global Markov property with respect to the induced acyclic directed mixed graph (ADMG) $\mathcal{G}_{obs}=(V_{obs},E^{\prime})$ (also known as latent projection, with potentially arbitrary unobserved processes) over the observed variables, we can take inspiration from the Fast Causal Inference (FCI) algorithm (Zhang, ) by first running its skeleton-discovery-part and then again leveraging the direction of time to uniquely infer the edge-type $\{\leftarrow,\rightarrow,\leftrightarrow\}$ for found adjacencies.
This reduces the typical partial ancestral graph output of FCI (representing an equivalence class of ancestral graphs) to a single, maximally informative graph (as informative as an ancestral graph can be).
We provide all details in Appendix D.

### 3.2 Signature Kernel Conditional Independence Test

Our theoretical results in Section 3.1 assume a CI oracle.
To make our algorithms usable in practice, we now propose a flexible CI test for path-valued random variables on different time intervals for expressions of the form eq. 2.
Recent work using the signature kernel is limited to unconditional hypothesis tests (Chevyrev & Oberhauser, ) or only use it as a heuristic measure of conditional independence (not developing a hypothesis test) (Salvi et al., ), and both do not use the time order.
Using terms of the signature to detect causality (not specifically with hypothesis testing) was also proposed by Giusti & Lee (); Glad & Woolf ().
Hence, we are the first to combine ideas from kernel-based CI tests and the signature kernel for a practically usable, consistent CI test on path space.
While the existing consistency proof for KCIPT (also applying to SDCIT) (Doran et al., ) relies on the existence of densities, we prove consistency for testing on path-valued random variables in Section A.11 requiring novel arguments.
Given $n$ samples of path segments $X^{(i)}_{\mathcal{I}}$, $Y^{(i)}_{\mathcal{J}},Z^{(i)}_{\mathcal{K}}$ on intervals $\mathcal{I},\mathcal{J},\mathcal{K}\subset[0,T]$ for $i\in[n]$, we compute the Gram matrices $k_{XX},k_{YY},k_{ZZ}\in\mathbb{R}^{n\times n}$ using the signature kernel $K_{S}$ and run a kernel-based permutation CI test like KCIPT (Doran et al., ) or SDCIT (Lee & Honavar, ) to test $X_{\mathcal{I}}\operatorname{\perp\perp}Y_{\mathcal{J}}\mid Z_{\mathcal{K}}$.
The theoretical foundation combined with the strong empirical performance of our CI test suggest that it may be of independent interest outside of causal discovery applications.

## 4 Experiments

In this section, after introducing the baselines and implementation details, we empirically evaluate our CI test in a variety of settings.
We then demonstrate strong performance in various causal discovery tasks and also outperform traditional methods in a real-world finance case study.

#### Baselines.

We compare against CCM (Sugihara et al., ), PCMCI (Runge et al., ), Granger causality (Granger, ), a kernel-based approach (Lau) (Laumann et al., ), and SCOTCH (Wang et al., ), a variational neural SDE method. CCM, Granger, and PCMCI are well-established methods for time series, while Lau and SCOTCH represent the latest advances in causal discovery for functional data and SDE models, respectively. CCM and Granger are limited to bivariate cases, while Lau is only applicable to unconditional tests here. Details are in Appendix C.

#### Implementation details and metrics.

We use [sigkerax](https://github.com/crispitagorico/sigkerax) for the signature kernel with an RBF kernel with length scale selected via a median heuristic (see Section B.1). For bivariate cases, the causal structure is $X^{1}\rightarrow X^{2}$, for $d>2$, we draw Erdös–Rényi DAGs. Performance is measured by Structural Hamming Distance (SHD, Section B.8). For $\operatorname{\perp\perp}_{s,h}^{+}$, $s=0.1\cdot T$ (and a fixed $T=1$) performed best (Table 4). The SDE model is

$$ $\displaystyle\mathrm{d}X_{t}=$ $\displaystyle(AX_{t}+c)\mathrm{d}t+\operatorname{\mathrm{Diag}}(BX_{t}+d) \mathrm{d}W_{t},\quad\text{with }A,B\in\mathbb{R}^{d\times d},c,d\in\mathbb{R} ^{d}\>.$ (5) $$

except for the 2-dimensional non-linear SDE, which is given by

$$ $\displaystyle\mathrm{d}\begin{pmatrix}X^{1}_{t}\\ X^{2}_{t}\end{pmatrix}=\begin{pmatrix}-r\omega\sin(\omega t)\\ r\omega\tanh(X^{2}_{t})\end{pmatrix}\mathrm{d}t+\operatorname{\mathrm{Diag}}(d ^{\top})\mathrm{d}W_{t}\>.$ (6) $$

#### Choosing (conditional) independence tests.

In initial experiments on two- and three-dimensional structures (e.g., forks, colliders, chains) detailed in Section B.2, we found bootstrapped SDCIT and HSIC to be the best performing variants of kernel-based (C)I tests (see Figures 6, 5, 7 and 6 for details) and use them for all subsequent experiments.

#### Power analysis of the unconditional test.

We first demonstrate consistently superior performance of our unconditional independence test (HSIC) SigKer (ours) over the only existing baseline for this setting from Laumann et al. () (Lau) in the linear SDE-setting eq. 5 (with $B\equiv 0$, $d_{i}=0.4$), reaching test power near 1 already for $n\geq 40$ when the causal interaction $a_{21}$ (strength of $X^{1}\to X^{2}$) is comparable to $X^{2}$’s self-dependence ($a_{22}$).
Figure 7 further highlights SigKer’s robustness even with *high data missingness*, and Figure 8 demonstrates its effectiveness on *fractional Brownian motions* outside the semi-martingale framework.
We are not aware of other CI tests that apply to such settings—establishing it as state-of-the-art CI test for stochastic processes.

Figure: Figure 3: Test power for $X^{1}_{[0,t]}\operatorname{\perp\perp}X^{2}_{[0,t]}$ over $\tfrac{a_{21}}{a_{22}}$. Lines (shades) are means (standard errors) over $1000$ SDE instances.
Refer to caption: x3.png

#### Leveraging the direction of time.

Next, we leverage time to orient the edge in the much-studied bivariate causal-discovery setting ($X^{1}\to X^{2}$) via our $\operatorname{\perp\perp}_{s,h}^{+}$ constraint with $K=\emptyset$ for a variety of settings: (i) linear SDEs with dependence in the drift, (ii) diffusion, (iii) or path-dependent, as well as (iv) nonlinear interactions.
Table 1 shows that SigKer decisively dominates all baselines in all settings for different sample sizes except diffusion dependence.
In this setting, SCOTCH is better *for specific hyperparameter settings*, while performing worse for others.
Crucially, one cannot ‘cross-validate’ or otherwise select such hyperparameters for the causal-discovery task in a data-driven fashion without ground-truth knowledge.
The lack of robustness (different optimal hyperparameters for different sample sizes) renders SCOTCH unreliable in practice even for diffusion dependence.

Specifically, we draw the settings for these experiments as follows:
for linear drift interactions we sample $a_{21}\sim\mathcal{U}([1,2.5])$, $a_{11},a_{22}\sim\mathcal{U}([-0.5,0.5])$, $a_{12}=0$, $B\equiv 0$, $d_{i}\sim\mathcal{U}([0.1,0.2])$ in a two-dimensional linear SDE model eq. 5.
For linear diffusion interactions we draw $a_{11},a_{22}\sim\mathcal{U}([0.5,1])$ and $b_{21}\sim\mathcal{U}([1,4.5])$), the rest set to zero in a two-dimensional linear SDE model eq. 5.
For path-dependence, we simulate $\mathrm{d}X_{t}^{2}=\mathbf{\mu}^{2}(X_{[0,t]}^{1})\mathrm{d}t+d_{2}\mathrm{d}
W_{t}^{1}$ via a three-dimensional SDE eq. 5 with $B\equiv 0$, $c\equiv 0$, $a_{23},a_{31}\sim\mathcal{U}([-3.5,-1]\cup[1,3.5])$ and $d=(d_{1},d_{2},0)^{\top}$, $d_{i}\sim\mathcal{U}([0.1,0.2])$, the rest set to zero. For the nonlinear SDEs we use eq. 6 with $\omega\sim\mathcal{U}([6\pi,8\pi])$, $r\sim\mathcal{U}([0.5,1])$, $d_{i}\sim\mathcal{U}([2,2.5])$.

**Table 1: SHD ($\times 10^{2}$) comparison of SigKer to the baselines in four bivariate SDE settings: linear, path-dependence, non-linear, and diffusion dependence. We ran SCOTCH for different sparsity parameters $\lambda$ and numbers of epochs $n_{e}$; there is no clear trend in either parameter, but $\lambda=100,n_{e}=2000$ performed best overall. Different learning rates performed worse across the board.**
|  |  | linear | path-dependence | non-linear | diffusion dependence |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $\times 10^{2}$ | $(\lambda,n_{e})$ | $n=200$ | $n=400$ | $n=200$ | $n=400$ | $n=200$ | $n=400$ | $n=200$ | $n=400$ |
| CCM |  | $176\pm 7$ | $166\pm 8$ | $186\pm 5$ | $194\pm 3$ | $120\pm 10$ | $106\pm 10$ | $100\pm 10$ | $84\pm 10$ |
| Granger |  | $92\pm 5$ | $87\pm 5$ | $85\pm 5$ | $95\pm 6$ | $103\pm 5$ | $103\pm 5$ | $105\pm 5$ | $111\pm 4$ |
| PCMCI |  | $89\pm 9$ | $100\pm 10$ | $91\pm 4$ | $125\pm 15$ | $41\pm 8$ | $55\pm 15$ | $98\pm 13$ | $75\pm 23$ |
| SCOTCH | 100, 1k | $74\pm 14$ | $77\pm 17$ | $50\pm 12$ | $46\pm 12$ | $64\pm 15$ | $80\pm 10$ | $75\pm 14$ | $62\pm 10$ |
| SCOTCH | 50, 2k | $79\pm 18$ | $44\pm 17$ | $23\pm 12$ | $73\pm 23$ | $83\pm 11$ | $73\pm 13$ | $\mathbf{25\pm 13}$ | $44\pm 28$ |
| SCOTCH | 100, 2k | $50\pm 17$ | $36\pm 19$ | $64\pm 19$ | $55\pm 15$ | $85\pm 7$ | $91\pm 9$ | $33\pm 16$ | $\mathbf{9\pm 8}$ |
| SigKer |  | $\mathbf{14\pm 4}$ | $\mathbf{7\pm 3}$ | $\mathbf{5\pm 2}$ | $\mathbf{6\pm 2}$ | $\mathbf{28\pm 5}$ | $\mathbf{5\pm 2}$ | $72\pm 6$ | $63\pm 5$ |

**Table 2: SHD ($\times 10^{2}$) comparison of SigKer to PCMCI and SCOTCH (different $\lambda$ and $n_{e}$) in causal discovery. Means and standard errors are over 40 SDE instances.**
|  | $d=3$ | $d=5$ | $d=10$ | $d=20$ | $d=50$ |
| --- | --- | --- | --- | --- | --- |
| PCMCI | $29\pm 16$ | $243\pm 37$ | $793\pm 84$ | $3530\pm 159$ | $19.6k\pm 857$ |
| SCOTCH 100, 2k | $188\pm 28$ | $417\pm 86$ | $250\pm 61$ | $1525\pm 1160$ | $10275\pm 6176$ |
| SCOTCH 200, 2k | $110\pm 21$ | $270\pm 48$ | $530\pm 223$ | $\mathbf{370\pm 174}$ | $\mathbf{538\pm 70}$ |
| SCOTCH 200, 1k | $400\pm 17$ | $1370\pm 29$ | $6425\pm 80$ | $705\pm 77$ | $7863\pm 2670$ |
| $\operatorname{\perp\perp}_{s,h}^{+}$ | $26\pm 5$ | $80\pm 8$ | $284\pm 19$ | $1026\pm 40$ | $4946\pm 113$ |
| $\operatorname{\perp\perp}_{\mathrm{sym}}$+p.p. | $\mathbf{13\pm 4}$ | $\mathbf{38\pm 7}$ | $\mathbf{157\pm 15}$ | $725\pm 439$ | $4593\pm 93$ |
| $\operatorname{\perp\perp}_{\mathrm{sym}}$ only | $57\pm 84$ | $147\pm 11$ | $436\pm 19$ | $1294\pm 48$ | $6005\pm 98$ |

#### Causal discovery.

For causal discovery, we sample $40$ linear SDEs from eq. 5 (with $a_{ij}\sim\mathcal{U}([-2,-1]\cup[1,2])$ for $j\neq i$ and $a_{ii}\sim\mathcal{U}([-0.5,0.5])$) with random DAG adjacency structures for $d\in\{3,5,10,20,50\}$ and use $n=200$ sample paths from each SDE as input for the algorithms.
Table 2 shows that our Algorithm 1 (denoted by $\operatorname{\perp\perp}_{s,h}^{+}$) and Algorithm 2 ($\operatorname{\perp\perp}_{\mathrm{sym}}$+pp and $\operatorname{\perp\perp}_{\mathrm{sym}}$ the result before post-processing) clearly outperform PCMCI and SCOTCH up to $d=10$.
SCOTCH’s heavy dependence on hyperparameter choices for $\lambda$ (graph sparsity) and $n_{e}$ (the number of epochs), can sometimes render it superior when charitably picking the best setting.
However, since good values particularly for $\lambda$ cannot be known up front nor selected in a data-driven fashion, one must interpret SCOTCH’s performance more conservatively—arguably in terms of its worst case performance over a set of reasonable hyperparameter settings.
Hence, SigKer – free of any hyperparameter choices – broadly outperforms the state-of-the-art even in the setting it was tailored to (SDE models).

#### Summary and discussion of results.

Our CI test for $\operatorname{\perp\perp}_{s,h}^{+}$ reliably detects the direction of time and outperforms strong baselines (CCM, Granger, PCMCI), consistently improving with larger sample sizes (Table 1).
SigKer also dominates PCMCI in causal discovery across dimensions and consistently beats most hyperparameter settings of SCOTCH.
While isolated hyperparameter settings for SCOTCH perform better than SigKer in diffusion dependence and high-dimensional causal discovery, its strong dependence on hyperparameters, particularly the sparsity parameter $\lambda$, renders it unreliable in practice.
Finally, SCOTCH and PCMCI are tailored to Markovian SDEs (or SDEs with a fixed lag).
Instead, our CI test is broadly applicable and effective across data modalities which we further demonstrate empirically in a functional data example where both PCMCI and SCOTCH fail to detect path-dependence, see Table 9.

Figure: Figure 4: An example graph demonstrating the advantage of our approach under partial observations.
Refer to caption: x4.png

#### The partially observed setting.

One of the main benefits of constraint-based causal discovery (especially with a non-parametric test) is the ability to handle partial observations. Score-based methods like SCOTCH are fundamentally challenged in this setting, since there could in principle be infinitely many unobserved variables, impossible to model, e.g., with neural-network based approaches.
We showcase this advantage on the challenging example graph in Figure 4, where SCOTCH falsely infers an adjacency $A\operatorname{\text{---}}D$ in 88 out of 100 runs, whereas our algorithm correctly handles the unobserved confounder and only (falsely) predicts this adjacency in 8 instances.

**Table 3: SigKer outperforms baselines in most pairs trading performance metrics. $\uparrow/\downarrow$ indicates ‘higher/lower is better’.**
|  | return $\uparrow$ | APR $\uparrow$ | Sharpe $\uparrow$ | maxDD $\downarrow$ | maxDDD $\downarrow$ |
| --- | --- | --- | --- | --- | --- |
| ADF | 0.004 | 0.004 | 0.090 | 0.087 | 230 |
| Granger | -0.010 | -0.011 | -0.230 | 0.056 | 219 |
| ADF $\&$ Granger | 0.008 | 0.008 | 0.242 | 0.022 | 153 |
| SigKer | 0.076 | 0.077 | 1.500 | 0.027 | 21 |

#### Real-world pairs trading example.

To demonstrate the applicability of our developed methods on real-data, we evaluate pairs trading strategies on ten stocks from the VBR Small-Cap ETF over a three-year period (2010/01/01–2012/12/31).
We only provide a concise summary of our results here and refer to Section B.7 for more details on this proof-of-concept study.
We assess our method’s effectiveness by quantifying the profit-and-loss (P&L) profile of the generated pairs trading strategy as a substitution for ground truth.
Pairs were selected based on pairwise p-values from different hypothesis tests: 1) cointegration via the Augmented Dickey-Fuller (ADF) test, 2) Granger causality, and 3) our method.
Table 3 shows that our strategy’s P&L substantially outperforms the baselines in total return, APR, Sharpe ratio, and maxDDD while being almost on par with ADF & Granger in the fifth relevant metric maxDD.
This case-study highlights the broad applicability and potential downstream impact of the findings in this work.

## 5 Conclusion and Future Work

We introduce constraint based causal discovery algorithms for both fully and partially observed data in the form of stochastic processes, proving them to be sound and complete assuming a CI oracle.
These are based on novel Markov properties corresponding to the developed CI constraints with respect to the induced dependence graph.
Our algorithms critically leverage the directionality of time to uniquely identifying the full underlying graph or induced ancestral graph, critically improving over existing constraint-based algorithms in the static case that only output equivalence classes.
Our framework also efficiently captures path- or diffusion-dependence.
Finally, we propose a practical and consistent kernel-based CI test on path-space leveraging recent advances in signature kernels that empirically outperforms existing alternatives across a wide range of settings also beyond the SDE model such as functional data and fractional Brownian motions.

Due to the identified fundamental limitations in the cyclic setting, we limit ourselves to acyclic (except for self-loops) dependence graphs and assume causal relationships to not change over time.
Relaxing both of these assumptions and assessing how much of the causal structure can in principle be learned in the cyclic settings are immediate interesting directions for future work.
We highlight that the field of causality represents just one of the many potential applications of our conditional independence test for path-valued random variables.
Exploring applications in different domains is a worthwhile direction for future work that could also catalyze new methodological developments.

## Appendix

## Appendix A Proofs and Theoretical Digressions

### A.1 Intuition and Details for the SDE Model

Making the actual dependence of $\mu^{k}$ and $\sigma^{k}$ on their arguments in eq. 1 explicit, we can rewrite it as

$$ $\begin{cases}\mathrm{d}X^{k}_{t}=\mu^{k}(X_{[0,t]}^{\mathrm{pa}_{k}})\mathrm{d }t+\sigma^{k}(X_{[0,t]}^{\mathrm{pa}_{k}})\mathrm{d}W^{k}_{t},\\ X^{k}_{0}=x^{k}_{0}\qquad k\in[d].\end{cases}$ (7) $$

Then $\mu^{k}$, $\sigma^{k}$ are functions defined on $C([0,+\infty),\mathbb{R}^{\mathrm{dim}(\mathrm{pa}_{k})})$ (or some suitable subspace thereof): Lipschitz conditions on the coefficients that guarantee existence and uniqueness for this type of SDE can be found in Rogers & Williams (), which we assume to hold throughout.
Each $W^{k}$ is an $m_{k}$-dimensional Brownian motion (a collection of $m_{k}$ independent $1$-dimensional Brownian motions, that is), $\sigma^{k}$ maps to the space of $n_{k}\times m_{k}$-dimensional matrices, and the noises $W^{k}$ together with the (possibly random) initial conditions $x^{k}_{0}$ are jointly independent.
In other words, the system can be written as an $n\coloneqq n_{1}+\ldots+n_{d}$-dimensional SDE driven by an $m\coloneqq m_{1}+\ldots+m_{d}$-dimensional Brownian motion, with a block-diagonal diffusion coefficient $\sigma$ (since the noise is unobserved, this structure has to be imposed if we wish not to deal with unobserved confounding).
The superscript $\mathrm{pa}_{k}$ refers to the parents of the $k^{\text{th}}$ node in $\mathcal{G}$, which may (and most often does) contain $k$ itself.
Intuitively, all of this means that, for all times $t$ and small increments $\Delta t$, the increment of the solution $X^{k}_{t+\Delta t}-X^{k}_{t}$ is random with distribution that is well-approximated by a multivariate normal with mean $\mu^{k}(X_{[0,t]}^{\mathrm{pa}_{k}})\Delta t$ and covariance function $\sigma^{k}(X_{[0,t]}^{\mathrm{pa}_{k}})\sigma^{k}(X_{[0,t]}^{\mathrm{pa}_{k}})
^{\intercal}\Delta t$, and independent of the history of the system up to time $t$.
This interpretation can actually be made precise by showing that these piecewise constant paths converge in law to the true solution (usually known as the weak solution, when viewed in this way).

We will refer to $\mathcal{G}$ as the *dependence graph* of the eq. 1, which we refer to as the *SDE model* for brevity. Compared to that of Peters et al. (), this model is slightly more general in that (i) it allows for path-dependence and (ii) for each node to represent a multidimensional process; the special case of state-dependent SDE—i.e., in which $\mu$ and $\sigma$ only depend on $X_{t}$, the value of $X$ at time $t$—continues to be an important special case, although our model also accommodates delayed SDEs (the coefficients depend on the value of $X$ at a prior instant in time, e.g., $X_{t-\tau}$ for fixed or possibly random/time-dependent $\tau>0$), and SDEs that depend on quantities involving the whole past of $X$, such as the average $t^{-1}\int_{0}^{t}X_{s}\mathrm{d}s$. We note that the choice not to make the coefficients explicitly time-dependent is deliberate and motivated by the requirement that the system be *causally stationary* (the causal relations between the variables do not change over time).

Given that we are considering the dynamic setting, it is generally natural to allow for $\mathcal{G}$ to have cycles. This comes at no additional requirement of consistency constraints as it does in the static case (see for example Bongers et al. ()), since the causal arrows in the model eq. 1 should be thought of as ‘pointing towards the infinitesimal future, with infinitesimal magnitude’, integrated over the whole time interval considered. Indeed, the system of SDEs does not require any global consistency to be well-posed, other than the regularity and growth conditions that guarantee global-in-time existence and uniqueness. On the other hand, we will be interested in the potential for constraint-based causal discovery of such systems, qualified by the following assumption on the types of conditional independencies that we allow to be tested in continuous time:

### A.2 Counterexample for Cyclic Causal Discovery in the SDE Model

The following example demonstrates that constraint-based causal discovery of the full graph is impossible using expressions of the form eq. 2 when allowing for cycles.

###### Example A.1 .

An oracle for eq. 2 is not powerful enough to rule out a directed edge $X^{1}\to X^{3}$ for an SDE model with the following dependence graph:

We now go through this example and describe a different type of oracle that would be required for causal discovery in cyclic SDE models.
Concretely, we claim that an oracle for eq. 2 is not powerful enough to rule out a directed edge $X^{1}\to X^{3}$ for an SDE model with the graph in Example A.1.
Testing $X^{1}_{[0,s]}\operatorname{\perp\perp}X^{3}_{[s,s+h]}\mid X^{2,3}_{[0,s]}$, will not remove the edge, due to the open path

$$ $X^{1}_{[0,s]}\to X^{1}_{[s,s+h]}\to X^{2}_{[s,s+h]}\to X^{3}_{[s,s+h]}.$ $$

Testing $X^{1}_{[0,s]}\operatorname{\perp\perp}X^{3}_{[s,s+h]}\mid X^{3}_{[0,s]},X^{2}_
{[0,s+h]}$, on the other hand, runs into the collider

$$ $X^{1}_{[0,s]}\to X^{1}_{[s,s+h]}\to X^{2}_{[s,s+h]}\leftarrow X^{3}_{[s,s+h]}.$ $$

There is no other way of splitting up the time interval (even allowing for more than 2 subintervals) that overcomes both these problems: when testing $X^{1}_{[0,s]}\operatorname{\perp\perp}X^{3}_{[s,s+h]}$ (as is necessary in order to rely on time to obtain direction of the arrow), if $X^{2}$ is not conditioned on over $[s,s+h]$ it will be a mediator, and if it is, it will be a collider.

Figure: Figure 5: The blue highlighter is for intervals over which the path is being tested for independence, while orange is for conditioning. The first figure illustrates the path opened by the unconditioned-on mediator $X^{2}_{[s,s+h]}$, while the second illustrates the path opened by the same variable acting as a conditioned-on collider. While not quite possible to draw extra/figures this way, it is helpful to think of the arrows as pointing from values to increments $X_{t}\to(Y_{t+\mathrm{d}t}-Y_{t})$, infinitesimally, with the total causal effect accrued over time.

What would be needed to detect the edge in the above example (and to perform causal discovery more generally in the cyclic case), is the availability of tests of *strong instantaneous non-causality in the Granger sense* as defined in Florens & Fougere ().
This reflects the fact that the SDE model is ‘acyclic at infinitesimal scales’.
Unfortunately, such a property—which has to do with the Doob-Meyer decompositions of functions of the process w.r.t. two different filtrations—is much more difficult to test for, as it is infinitesimal in nature.
Recent progress in this direction was made in Christgau et al. () for the case of SDEs driven by jump processes; this local independence criterion, which goes back to Schweder (), is however *weak* in the sense that it only takes into account the Doob-Meyer decomposition of $X$ and not functions of it, and is not therefore able, for example, to detect dependence in the diffusion coefficient (Gégout-Petit & Commenges, ).
Here we take the view that a continuously-indexed path, queried over an interval, can be considered as an observational distribution (cf. the discussion in Sokol & Hansen () on whether the infinitesimal generator can be considered ‘observational’).

###### Example A.1 .

### A.3 Signature Kernel Details

For any $p\in[1,2)$, let $\mathcal{K}\subset C_{p}([0,T],\mathbb{R}^{d})$ be a compact subset such that the signature transform is a continuous injection from $\mathcal{K}$ to $\mathcal{T}(\mathbb{R}^{d})$. An example of compact set $\mathcal{K}$ satisfying such conditions is the set of time-augmented paths started at a common origin, where $C_{p}([0,T],\mathbb{R}^{d})$ is endowed with the $p$-variation topology. Further examples of compact sets and choices of topologies on path-space are discussed by Cass & Turner (). Then, the signature kernel is *universal* over $\mathcal{K}$ in the sense that the associated RKHS, whose elements are restricted to act on $\mathcal{K}$, is dense in $C(\mathcal{K})$, the space of continuous real-valued functions on $\mathcal{K}$, endowed in the topology of uniform convergence.
See Cass et al. (, Proposition 3.3) for a proof. The compactness assumption can be relaxed by either renormalizing the signature transform so that its range is always guaranteed to be bounded (Chevyrev & Oberhauser, ) or by operating on weighted path-spaces with a different choice of topology Cuchiero et al. ().

The signature kernel is well-defined and can be shown, under similar compactness conditions, to be universal also when $p\geq 2$ (Salvi et al., ; Lemercier & Lyons, ). In the specific setting of stochastic differential equations (SDEs) considered in this paper, i.e. when $p\in[2,3)$, one can use a classical limiting procedure from rough analysis to show that the integral equation in eq. 4 takes the Stratonovich form $f(t,t^{\prime})=1+\int_{s}^{t}\int_{s^{\prime}}^{t^{\prime}}f(u,v)\langle\circ
dX
_{u},\circ dY_{v}\rangle_{1}$. We refer the interested reader to Salvi et al. (, Section 4.3) for further details.

Another important property that a kernel should possess in order to ensure theoretical guarantees of most kernel methods, including tests for conditional independence (Muandet et al., ), is that of being *characteristic*. Loosely speaking, a kernel $k$ is said characteristic to a topological space if for any (Borel) probability measure $\mu$ on that space, the $\mu$-expectation of the feature map $k(x,\cdot)$ uniquely characterises $\mu$. An important result from the theory of kernels relates to the equivalence between the notions of universality and characteristicness. In particular, these three notions can be shown to be equivalent for kernels defined over general locally convex topological vector spaces (Simon-Gabriel & Schölkopf, , Theorem 6).

### A.4 Testing on Time Points

It is also possible to replace any of the path-valued random variables in eq. 2 by a random variable at a single time-point (as for example done in the post-processing by testing on $X_{0}$) using a regular rbf-kernel.
Note, that in addition, if we knew in advance that there is no path-dependence, the unconditional independence $X^{i}\operatorname{\perp\perp}_{s,h}^{+}X^{j}$ is equivalent to the conditional independence $X^{i}_{s}\operatorname{\perp\perp}X^{j}_{s+h}-X^{j}_{s}\mid X^{j}_{s}$ (at fixed times $s,s+h$): this is because without path-dependence $X^{j}_{u}$ for $u\in[s,t]$ only depends on $X_{[s,t]}^{j}$ via $X^{j}_{s}$.
Testing on entire (segments of) paths is always required when allowing for path-dependence.

### A.5 Proof of Proposition 3.1

###### Proof of Proposition 3.1 .

Since $\mathcal{G}$ is acyclic, a directed cycle in $\tilde{\mathcal{G}}$ that is not a loop must contain nodes both in $V_{0},V_{1}$. But there can be no such directed cycles, since edges only travel in the direction $V_{0}\to V_{1}$, by construction of $\tilde{\mathcal{G}}$. Since $\tilde{\mathcal{G}}$ also free of loops, it is a DAG. We now show that there exists a structural causal model (SCM) over $\tilde{\mathcal{G}}$ with the path-valued random variables of the statement: Markovianity will then follow from Peters et al. (, Proposition 6.31), original to Verma & Pearl (). In other words, we must show there exist Borel-measurable functions (setting $t\coloneqq s+h$)

$$ $\displaystyle F^{k}_{0}\colon\mathbb{R}^{\mathrm{dim}(\mathrm{pa}_{k})}\times C ([0,s],\mathbb{R}^{m_{k}})\times C([0,s],\mathbb{R}^{\mathrm{dim}(\mathrm{pa}_ {k}\setminus\{k\})})$ $\displaystyle\rightharpoonup C([0,s],\mathbb{R}^{n_{k}})$ $\displaystyle(x_{0}^{\mathrm{pa}(k)},W_{[0,s]}^{k};X_{[0,s]}^{\mathrm{pa}_{k} \setminus\{k\}})$ $\displaystyle\mapsto X^{k}_{[0,s]}$ $$

and

$$ $\displaystyle F^{k}_{1}\colon C([s,t],\mathbb{R}^{m_{k}})\times C([0,s], \mathbb{R}^{\mathrm{dim}(\mathrm{pa}_{k})})\times C([s,t],\mathbb{R}^{\mathrm{ dim}(\mathrm{pa}_{k}\setminus\{k\})})$ $\displaystyle\rightharpoonup C([0,s],\mathbb{R}^{n_{k}})$ $\displaystyle(W_{[s,t]}^{k};X_{[0,s]}^{\mathrm{pa}_{k}},X_{[s,t]}^{\mathrm{pa} _{k}\setminus\{k\}})$ $\displaystyle\mapsto X^{k}_{[0,s]}$ $$

with the initial conditions $x^{0}$ and Brownian path segments $W^{k}_{[0,s]}$, $W^{k}_{[s,t]}$ collectively independent.
This independence follows from the definition of the SDE model and from independence of Brownian increments over disjoint or consecutive intervals.
These functions are partial in that they are not defined on the whole space of continuous functions on the interval: we only need to show them to be defined on a measurable set of paths that contains all solutions of SDEs on the required interval.
$F^{k}_{1}$ is defined as follows, with care to make explicit the dependence on the solution on the two intervals via the operation of path-concatenation $*$:

$$ $\displaystyle F^{k}_{1}(X^{k}_{[s,t]})=\text{solution of}\quad X^{k}_{u}=X^{k} _{s}+\int_{0}^{u}$ $\displaystyle\mu^{k}(X_{[0,s]}^{\mathrm{pa}_{k}\setminus\{k\}}*X_{[s,u]}^{ \mathrm{pa}_{k}\setminus\{k\}},X_{[0,s]}^{k}*X_{[s,u]}^{k})\mathrm{d}t$ $\displaystyle+\sigma^{k}(X_{[0,s]}^{\mathrm{pa}_{k}\setminus\{k\}}*X_{[s,u]}^{ \mathrm{pa}_{k}\setminus\{k\}},X_{[0,s]}^{k}*X_{[s,u]}^{k})\mathrm{d}W^{k}_{t}$ $$

We have also separated out the $k$-th component of the solution from the rest of the arguments on which $\sigma^{k},\mu^{k}$ are dependent, which we consider part of the measurable adapted dependence on $W^{\mathrm{pa}\setminus\{k\}}$ needed in the existence and uniqueness theorem.
$F^{k}_{0}$ is defined similarly (with no dependence on past path-segments and with $x_{0}^{k}$ replacing $X_{s}^{k}$).
Such functions are well-known to be well-defined and measurable (Rogers & Williams, , Theorem 10.4).
∎

###### Remark A.2 .

In testing, using the symmetric criterion $\operatorname{\perp\perp}_{\mathrm{sym}}$ has proven to be more reliable than the asymmetric one used in the test above. A procedure, which is a little redundant, but which makes good use of the superior performance of $\operatorname{\perp\perp}_{\mathrm{sym}}$ would work as follows:

- 1.
Run the PC algorithm to discover the CPDAG corresponding to the variables $X_{[0,T]}^{k}$.
- 2.
Pass this CPDAG as input to Algorithm 1: by this we mean that $\tilde{V}$ is the same and the input $\tilde{E}$ is given by
$\{i_{0}\to j_{0},\ i_{1}\to j_{1},i_{0}\to j_{1}\mid i,j\in V,\ (i\to j)\text{
or }(i\relbar\mkern-9.0mu\relbar j)\in E^{\prime}\}\cup\{k_{0}\to k_{1}\mid k
\in V\}$
where $(i\to j)$ denotes an undirected arrow and $E^{\prime}$ is the CPDAG output by the classical PC algorithm.
- 3.
The loop-recovery step is identical as before.

The algorithm would run very similarly as before, but performing fewer CI tests in the first phase, since it can make use of the information contained in the partially-directed skeleton.

###### Remark A.3 .

If we knew in advance that there is no path-dependence, the test $X^{i}\operatorname{\perp\perp}^{+}_{s,h}X^{j}$ (with no additional conditioning set) is equivalent to the test on values of the process $X^{i}_{s}\operatorname{\perp\perp}X^{j}_{s+h}\mid X^{j}_{s}$: this is because $X^{j}_{u}$, for $u\in[s,t]$ only depends on $X_{[s,t]}^{j}$ via $X^{j}_{s}$.
However, this kind of statement no longer works in the conditional case (by arguments similar to Example A.1), and testing on whole paths is strictly necessary when allowing for path-dependence.

###### Remark A.4 .

CI testing of SDEs is conceptually similar to static CI on discretely sampled data, when the sampling rate is lower than the actual frequency at which causal effects propagate.
This is because, in continuous time and for any $h>0$, there are going to be causal effects that occur at time scales less than $h$.
Whenever we just discretely sample time series with no instantaneous effects, Peters et al. (, Thm. 10.3) provides full causal discovery, even in the presence of cycles in the summary graph, by performing $d^{2}$ tests $X^{i}_{s}\operatorname{\perp\perp}X^{j}_{s+\Delta s}\mid X^{[d]\setminus\{i,j
\}}_{[0,s]}$.
Potentially absent loops could then be removed as in Algorithm 1.
Ignoring that this might suffer from the conditioning set being too large, our signature method can be of help in this setting too, if there is path-dependence: the conditioning variable can be taken to be $S(X^{[d]\setminus\{i,j\}})_{0,s}$.

###### Proof of Proposition 3.1 .

###### Remark A.2 .

###### Remark A.3 .

###### Remark A.4 .

### A.6 Proof of Theorem 3.2

###### Proof of Theorem 3.2 .

We have to show that $E^{Alg}=E$ (where $E^{Alg}$ the edges detected by Algorithm 1).
"$\subseteq$:" This direction can be shown by finding an appropriate separation set in $\tilde{\mathcal{G}}$ by exploiting the global Markov property Proposition 3.1. We begin by proving correctness of the recovery of the skeleton modulo loops. This will follow if we show that for $(i\to j)\not\in E$ for $i,j\in V$ distinct

$$ $i_{0}\operatorname{\perp\perp}^{\tilde{\mathcal{G}}}_{d-sep}j_{1}\mid\{j_{0}\} \cup\{k_{0},k_{1}\mid k\in\mathrm{pa}^{\mathcal{G}}_{j}\setminus\{j\}\},$ $$

where $\operatorname{\perp\perp}^{\tilde{\mathcal{G}}}_{d-sep}$ denotes $d$-separation in the graph $\tilde{\mathcal{G}}$. Indeed, eventually, the set $K$ in the algorithm will take the value $\{k_{0},k_{1}\mid k\in\mathrm{pa}^{\mathcal{G}}_{j}\setminus\{j\}\}$, at which point the three edges $(i_{0},j_{0})$, $(i_{1},j_{1})$, $(i_{0},j_{1})$ (which are either all present or all absent, inductively, by construction of $\tilde{E}$) are deleted. (Of course it may be that these edges are deleted before this, if a smaller or different $d$-separating set is found, but note that edges are never added.) All undirected paths between $i_{0}$ and $j_{1}$ factor as $i_{0}\cdots h_{0}\to l_{1}\cdots j_{1}$, where the dots stand for a possibly empty undirected path. All such paths with $l_{1}=j$ are blocked by $\{j_{0}\}\cup\{k_{0}\mid k\in\mathrm{pa}^{\mathcal{G}}_{j}\setminus\{j\}\}$, so we focus our attention on the case in which $l_{1}\cdots j_{1}$ is non-empty. Here we follow a similar argument to that of Verma & Pearl (, Lemma 1). If $l_{1}\cdots j_{1}$ is of the form $l_{1}^{1}\cdots l_{1}^{r}\to j_{1}$ (with $l^{1}_{1}=l_{1}$) it is blocked by $\{k_{1}\mid k\in\mathrm{pa}^{\mathcal{G}}_{j}\setminus\{j\}\}$. Assume instead it is of the form $l_{1}^{1}\cdots l_{1}^{r}\leftarrow j_{1}$ and let $1\leq q\leq r$ be such that $l^{q}_{1}$ is the first collider on the entire path $i_{0}\cdots j_{1}$, starting from $j_{1}$ and travelling back: certainly this exists (and belongs to $V_{1}$), thanks to the arrows $h_{0}\to l_{1}^{1}$ and $l^{r}_{1}\leftarrow j_{1}$. For the path $i_{0}\cdots j_{1}$ to be open given the conditioning set, it must be the case that $l^{p}_{1}$ must be an ancestor (in $\tilde{\mathcal{G}}$) or member of $\{k_{1}\mid k\in\mathrm{pa}_{j}^{\mathcal{G}}\setminus\{j\}\}$ ($l^{p}_{1}$ cannot be an ancestor of nodes in $V_{0}$): this would yield a directed cycle, which is impossible.

We argue similarly for the loop-removal phase. If $(i\in i)\not\in E\iff(i_{0},i_{1})\not\in\tilde{E}$, we must show

$$ $i_{0}\operatorname{\perp\perp}^{\tilde{\mathcal{G}}}_{d-sep}i_{1}\mid\{k_{0},k _{1}\mid k\in\mathrm{pa}^{\mathcal{G}}_{i}\setminus\{i\}\}$ $$

Consider a path $i_{0}\cdots h_{0}\to l_{1}\cdots i_{1}$ in $\tilde{\mathcal{G}}$. If the segment $l_{1}\cdots i_{1}$ is non-empty we conclude that the path is blocked by the same argument as above (with $i_{1}$ replacing $j_{1}$). Assume that the path is of the form $i_{0}\cdots h_{0}\to i_{1}$ (with $i_{0}\neq h_{0}$ since $(i_{0},i_{1})\not\in\tilde{E}$): in this case $h\in\mathrm{pa}^{\mathcal{G}}_{i}\setminus\{i\}$ and the path is again blocked. Thus all paths are blocked and the implication follows.
$\supseteq:$ Assume $(i,j)\in E$, thus $(i_{0},j_{1})\in\tilde{E}$. Hence as $\tilde{\mathcal{G}}$ a DAG by construction, $\nexists$ $S\subseteq V\setminus\{i,j\}$ s.t. $\{i_{0}\}\operatorname{\perp\perp}_{d}^{\tilde{\mathcal{G}}}\{j_{1}\}\mid S_{0
}\cup\{j_{0}\}\cup S_{1}$. By the faithfulness assumption Definition A.5 one therefore has $X_{[0,s]}^{i}\perp\not\hskip 1.00006pt\perp X_{[s,s+h]}^{j}\mid X_{[0,s]}^{S},
X_{[0,s]}^{j},X_{[s,s+h]}^{S}$ $\forall$ $S\subseteq V\setminus\{i,j\}$. Hence $(i,j)\in E^{Alg}$.
∎

###### Proof of Theorem 3.2 .

### A.7 Faithfulness and Causal Minimality

In the above proof, we have used the following notion of *faithfulness*:

###### Definition A.5 ((Global) Faithfulness of $\operatorname{\perp\perp}_{s,h}^{+}$ in $\tilde{\mathcal{G}}$ ) .

Let $\{X^{i}_{t}\}_{i\in[d]}$ be the coordinate processes of a solution of eq. 1 for $t\in[0,1]$, $\mathcal{G}=(V\cong[d],E)$ the dependence graph defined by eq. 1. The independence relation $\operatorname{\perp\perp}_{s,h}^{+}$ is called *(globally) faithful* w.r.t. the DAG $\tilde{\mathcal{G}}$, if

$$ $\displaystyle X^{A}\operatorname{\perp\perp}_{s,h}^{+}X^{B}\mid X^{C}\> \Rightarrow\>A_{0}\operatorname{\perp\perp}_{d}^{\tilde{\mathcal{G}}}B_{1}\mid B _{0},C_{0},C_{1}$ (8) $$

for all $A,B,C\subseteq V$ pairwise disjoint.

Faithfulness is a standard assumption in the causal discovery literature. Some type of faithfulness-like assumption is strictly necessary to do causal discovery since, otherwise, it would be impossible to make any inference about adjacencies in the underlying graph from the observed distribution. While we used the strong faithfulness assumption Definition A.5 for the proof above, various weaker faithfulness-like assumptions in the case of DAG-based causal discovery could be used, e.g., the strictly weaker assumption known as ‘parent faithfulness’ (Mogensen () gives a definition of parent faithfulness in the context of stochastic processes). This assumption suffices in our algorithm without making any changes to it (or the respective proofs) in the acyclic case.

The strength of the faithfulness assumption in real-world data is heavily debated.
Besides Weinberger (), who states that ‘proposed counterexamples to the causal faithfulness condition are not genuine’, Spirtes et al. () argues that for parametric models, ‘the measure of the set of free parameter values for any DAG that lead to such cancellations is zero for any smooth prior probability density, such as the Gaussian or exponential one, over the free parameters’. What is meant in our case is, that when allowing the diffusion coefficients $\sigma$ in eq. 1 to range in some ‘reasonable’ finite-dimensional family (such as linear or affine functions) and drawing the parameters of such a family from a distribution which has positive density w.r.t. the Lebesgue measure (and independently of the noise), the resulting distribution on path-space, split over $[0,s]$ and $[s,t]$, will almost surely be faithful for any $0<s<t$ w.r.t. the lifted dependence graph. We do not attempt a proof of this statement here, as it may require considerable effort and technique.
On the other side, Andersen () or Uhler et al. () (among others) caution against making the strong faithfulness assumption easily. Many of their concerns are ameliorated by only assuming the weaker parent faithfulness.

Causal minimality, on the other hand, is generally easier to understand and can be expected to hold for the arrow $X^{i}_{[0,s]}\to X^{j}_{[s,t]}$ whenever the following condition is satisfied: the conditional distribution $X^{j}_{u}\mid X^{\mathrm{pa}^{\mathcal{G}}_{j}\setminus\{i\}}=x^{j}$ admits a positive density on some submanifold $M_{u}^{x}$ of $\mathbb{R}^{n}$ and for all $u$ there exist $x^{i}_{1},x^{i}_{2}\in M_{u}^{x}$, $x^{i}_{1}\neq x^{i}_{2}$, s.t. $(\mu^{j},\sigma^{j})(x,x^{i}_{1})\neq(\mu^{j},\sigma^{j})(x,x^{i}_{2})$.

###### Example A.6 (Causal minimality) .

The fact that we are allowing for variables to have more than one dimension means the causal minimality condition is not a trivial one, as the following examples demonstrate. Let $X^{1}=(W,0)$ where $W$ is a $1$-dimensional Brownian motion (so that $n_{1}=2$), and let $\mu^{1}$ and/or $\sigma^{1}$ depend on $X^{1}$ only through its second coordinate. Then, even though this dependence may be non-trivial, causal minimality does not hold. This may lead one to believe that one can generically only expect causal minimality to hold if $m=n$, but this is not the case. Take, for example, $d=4=n$ and assume $\sigma^{3}\equiv 0$, i.e., $X^{3}$ has no driving Brownian motion, so that $m$ is only $3$. Assume, furthermore, that $\mathrm{pa}_{4}=\{4\}$. Then, by Hörmander’s Lie bracket-generating condition (see for example §V.38 in Rogers & Williams ()), $(X^{1},X^{2},X^{3})$ has a density in $\mathbb{R}^{3}$ for any choice of an initial condition, if $(\sigma^{1}(x),0,0),(0,\sigma^{2}(x),0)$ span $\mathbb{R}^{2}$ for all $x\in\mathbb{R}^{3}$ and $[(\sigma^{1}(x),0,0),\mu]=\sigma^{1}\partial_{1}\mu-\mu^{i}\partial_{i}\sigma^
{1}$ or $[(\sigma^{2}(x),0,0),\mu]$ span the third direction. Arrows going from the first three nodes to the fourth will be necessary for Markovianity as long as the coefficients of $X^{4}$ depend on the respective variables on the support of such density. The point is that even though $(X^{1},X^{2},X^{3})$ is only driven by a two-dimensional Brownian motion, the coefficients of the SDE (at the initial condition) impart a ‘twist’ to the solution, making causal minimality a trivial condition on $(\sigma^{4},\mu^{4})$ which is verified as soon as there is any dependence.

###### Definition A.5 ((Global) Faithfulness of ⟂ ⟂ s , h + superscript subscript perpendicular-to absent perpendicular-to 𝑠 ℎ \operatorname{\perp\perp}_{s,h}^{+} start_OPFUNCTION ⟂ ⟂ end_OPFUNCTION start_POSTSUBSCRIPT italic_s , italic_h end_POSTSUBSCRIPT start_POSTSUPERSCRIPT + end_POSTSUPERSCRIPT in 𝒢 ~ ~ 𝒢 \tilde{\mathcal{G}} over~ start_ARG caligraphic_G end_ARG ) .

###### Example A.6 (Causal minimality) .

### A.8 Global Markov Property for the Symmetric Independence Criterion

In this section, we provide a proof for the global Markov property of the symmetric independence criterion in order for the PC-algorithm to be applicable. It is based on the notion of independence models:

###### Definition A.7 .

Let $V\cong[d]$ be a set. An *independence model* $\mathcal{J}(V)$ over $V$ is a ternary relation over disjoint subsets of $V$,

$$ $\displaystyle\mathcal{J}(V)\subset\{(A,B,C)\mid A,B,C\subset V\>\text{disjoint}\}$ $$

When $(A,B,C)$ is a triple in $\mathcal{J}(V)$, $(A,B,C)\in\mathcal{J}(V)$, we also use $\langle A,B\mid C\rangle$ to denote $(A,B,C)$. This notation highlights the fact that $C$ is a conditioning set. An independence model $\mathcal{J}(V)$ is called a *semigraphoid* if 1.-4. hold for all disjoint $A,B,C\subset V$.

- 1.
(symmetry) $\langle A,B\mid C\rangle\in\mathcal{J}(V)$ $\Rightarrow$ $\langle B,A\mid C\rangle\in\mathcal{J}(V)$
- 2.
(decomposition) $\langle A,B\mid C\rangle\in\mathcal{J}(V)$ and $D\subset B$ $\Rightarrow$ $\langle A,D\mid C\rangle\in\mathcal{J}(V)$
- 3.
(weak union) $\langle A,B\cup D\mid C\rangle\in\mathcal{J}(V)$ $\implies$ $\langle A,B\mid C\cup D\rangle\in\mathcal{J}(V)$
- 4.
(contraction) $\langle A,B\mid C\rangle\in\mathcal{J}(V)$ and $\langle A,D\mid B\cup C\rangle\in\mathcal{J}(V)$ $\implies$ $\langle A,B\cup D\mid C\rangle\in\mathcal{J}(V)$

A semigraphoid $\mathcal{J}(V)$ is called *graphoid* if 5. holds.

- 5.
(intersection) $\langle A,B\mid C\cup D\rangle\in\mathcal{J}(V)$ and $\langle A,D\mid B\cup C\rangle\in\mathcal{J}(V)\implies\langle A,B\cup D\mid C
\rangle\in\mathcal{J}(V)$

An independence model $\mathcal{J}(V)$ called *compositional* if 6. holds.

- 6.
(composition) $\langle A,B\mid C\rangle\in\mathcal{J}(V)$ and $\langle A,D\mid C\rangle\in\mathcal{J}(V)\implies\langle A,B\cup D\mid C
\rangle\in\mathcal{J}(V)$

The graphical criterion of d-separation defines a independence model on a DAG $\mathcal{G}=(V,\mathcal{E})$ by

$$ $\displaystyle\mathcal{J}_{d-sep}(\mathcal{G})=\left\{(A,B,C)\>\mid\>A,B,C \subset V\>\text{disjoint},A\operatorname{\perp\perp}_{d-sep}B\mid C\right\}$ (9) $$

###### Theorem A.8 ( Lauritzen & Sadeghi ( 2018 ) ) .

For any DAG $\mathcal{G}=(V,\mathcal{E})$, the independence model $\mathcal{J}_{d-sep}(\mathcal{G})$ is a compositional graphoid.

Let $V=[d]$. Given random variables $X^{i}:(\Omega,\mathcal{A},P)\rightarrow(\mathcal{X}^{i},\mathcal{A}_{i}),\ i
\in V,$ and $A,B,C\subset V$ disjoint, conditional independence

$$ $\displaystyle X^{A}\operatorname{\perp\perp}X^{B}\mid X^{C}:\Leftrightarrow \sigma(X^{A})\operatorname{\perp\perp}\sigma(X^{B})\mid\sigma(X^{C})$ $$

defines the *probabilistic independence model*

$$ $\displaystyle\mathcal{J}(P)=\left\{(A,B,C)\>\mid\>A,B,C\subset V\>\text{ disjoint},X_{A}\operatorname{\perp\perp}X_{B}\mid X_{C}\right\}.$ (10) $$

We will use this as a *symmetric conditional independence* relation. Note that $X^{i}$ may take values in a path-space in which case it is a stochastic process. The independence model $\mathcal{J}(P)$ is a semigraphoid which is an immediate consequence of sub-$\sigma$-algebra properties.

For $v\in V$, we let $\mathrm{nd}_{v}^{\mathcal{G}}$ denote the set of *nondescendants* of $v$, i.e., the set of nodes $i$ such that there is no directed path from $v$ to $i$.

###### Definition A.9 (Directed global and local Markov properties) .

Let $\mathcal{G}=(V,\mathcal{E})$ be a DAG and $\mathcal{J}(V)$ be an independence model over $V$. The independence model $\mathcal{J}(V)$ satisfies the *global Markov property w.r.t. $\mathcal{G}$* :$\Leftrightarrow$

$$ $\displaystyle\langle A,B\mid C\rangle\in\mathcal{J}_{d-sep}(\mathcal{G})\quad \implies\quad\langle A,B\mid C\rangle\in\mathcal{J}(V)\qquad\forall A,B,C \subset V\>\text{disjoint}$ (11) $$

The independence model $\mathcal{J}(V)$ satisfies the *directed local Markov property w.r.t. $\mathcal{G}$* :$\Leftrightarrow$

$$ $\displaystyle\langle\{v\},\left(\mathrm{nd}^{\mathcal{G}}_{v}\backslash\mathrm {pa}^{\mathcal{G}}_{v}\right)\mid\mathrm{pa}^{\mathcal{G}}_{v}\rangle\in \mathcal{J}(V)\text{ for all }v\in V$ (12) $$

###### Theorem A.10 ( Lauritzen et al. ( 1990 ) ) .

Let $\mathcal{G}=(V,\mathcal{E})$ be a DAG, let and $\mathcal{J}(V)$ a semigraphoid over $V$. The directed global and local Markov properties are equivalent, that is,

$$ $\displaystyle\mathcal{J}(V)\quad\text{satisfies \lx@cref{creftype~refnum}{glob al_Markov} w.r.t.\ }\mathcal{G}\>\Leftrightarrow\>\mathcal{J}(V)\quad\text{ satisfies \lx@cref{creftype~refnum}{local_Markov} w.r.t.\ }\mathcal{G}$ (13) $$

We therefore only have to establish the local Markov property for the symmetric independence model $\operatorname{\perp\perp}_{\mathrm{sym}}$:

###### Proposition A.11 .

Let $X$ be a set of variables induces by the model in eq. 7 with a constant and diagonal diffusion matrix such that it induces the DAG $\mathcal{G}=(V,\mathcal{E})$. Then the system satisfies the *directed local Markov property* w.r.t. $\mathcal{G}$,

$$ $\displaystyle X^{i}\operatorname{\perp\perp}_{\mathrm{sym}}\mathrm{nd}^{ \mathcal{G}}_{X^{i}}\setminus\mathrm{pa}^{\mathcal{G}}_{X^{i}}\mid\mathrm{pa}^ {\mathcal{G}}_{X^{i}}\qquad\forall i\in[d]$ (14) $$

###### Proof.

The idea of the proof is that all information about node $X^{i}$ is contained in its parents $pa_{i}:=\mathrm{pa}^{\mathcal{G}}_{X^{i}}$, the Brownian motion $W^{i}$, and the initial condition $X_{0}^{i}$. A similar idea can be used to show the result in the case of a SCM Pearl ().
We let $\tilde{\mathcal{J}}(P)$ denote the semigraphoid induced by $\operatorname{\perp\perp}_{\mathrm{sym}}$ on the variable set $\tilde{V}=\{X^{1},\ldots,X^{n}\}\cup\{W^{1},\ldots,W^{n}\}\cup\{X^{1}_{0},
\ldots,X^{n}_{0}\}$, and we let ${\mathcal{J}}(P)$ denote the semigraphoid induced by $\operatorname{\perp\perp}_{\mathrm{sym}}$ on the variable set ${V}=\{X^{1},\ldots,X^{n}\}$. We define the DAG

$$ $\displaystyle\tilde{\mathcal{G}}=\left(\tilde{V}=\{X^{1},\ldots,X^{n}\}\cup\{W ^{1},\ldots,W^{n}\}\cup\{X^{1}_{0},\ldots,X^{n}_{0}\},\>\tilde{E}=E\cup\{X^{j} \leftarrow W^{j}\}\cup\{X^{j}\leftarrow X^{j}_{0}\}\right).$ (15) $$

By construction, $\mathcal{F}_{t}(X^{i})\subset\mathcal{F}_{t}(W^{i})\vee\mathcal{F}_{t}(X^{
\mathrm{pa}^{\mathcal{G}}_{i}})\vee\sigma(X^{i}_{0})$, and therefore

$$ $\displaystyle\mathcal{F}_{t}(X^{i})\operatorname{\perp\perp}\mathcal{F}_{t}(X^ {\mathrm{nd}^{\mathcal{G}}_{i}\setminus\mathrm{pa}_{i}}\vee W^{-i}\vee X_{0}^{ -i})\mid\mathcal{F}_{t}(X^{\mathrm{pa}_{i}})\vee\mathcal{F}_{t}(W^{i})\vee \sigma(X^{i}_{0})$ $$

where the parent and nondescendant sets are to be read in the graph $\mathcal{G}$, $W^{-i}=\{W_{1},\ldots,W_{n}\}\setminus\{W_{i}\}$, and $X_{0}^{-i}=\{X_{0}^{1},\ldots,X_{0}^{n}\}\setminus\{X_{0}^{i}\}$. The above statement corresponds to eq. 12 when $v\in\{X_{1},\ldots,X_{n}\}$. When $v=W_{i}$ or $v=X_{0}^{i}$, $v$ has no parents in $\tilde{\mathcal{G}}$ and eq. 12 also holds. Therefore, the independence model $\tilde{\mathcal{J}}(P)$ satisfies the local directed Markov property w.r.t. to $\tilde{\mathcal{G}}$.

The independence model $\tilde{\mathcal{J}}(P)$ is a semigraphoid, and by Theorem A.10 it satisfies eq. 11 w.r.t. $\tilde{\mathcal{G}}$. We have $X^{i}\operatorname{\perp\perp}_{d-sep}X^{\mathrm{nd}_{i}\setminus\mathrm{pa}_{
i}}\mid X^{\mathrm{pa}_{i}}$ in $\tilde{\mathcal{G}}$, and one therefore has

$$ $\displaystyle\mathcal{F}_{t}(X^{i})\operatorname{\perp\perp}\mathcal{F}_{t}(X^ {\mathrm{nd}_{i}\setminus\mathrm{pa}_{i}})\mid\mathcal{F}_{t}(X^{\mathrm{pa}^{ i}})$ $$

establishing eq. 12 for $\mathcal{J}(P)$ w.r.t. $\mathcal{G}$.
∎

###### Definition A.7 .

###### Theorem A.8 ( Lauritzen & Sadeghi ( 2018 ) ) .

###### Definition A.9 (Directed global and local Markov properties) .

###### Theorem A.10 ( Lauritzen et al. ( 1990 ) ) .

###### Proposition A.11 .

###### Proof.

### A.9 Proof of Corollary 3.3 and Algorithm for Robust Causal Discovery

###### Proof.

This is clear from the data generating mechanism in eq. 1 where all Brownian motions $\mathrm{d}W^{k}_{t}$ and initial conditions $X_{0}^{k}$ are jointly independent.
Since $X^{j}_{[0,T]}$ is fully determined by $\{X_{0}^{\mathrm{an}_{j}},\mathrm{d}W^{\mathrm{an}_{j}}_{t}\}$ where $\mathrm{an}_{j}$ are the ancestors of $j$ in $\mathcal{G}$, it follows that $i\in\mathrm{an}_{j}$ if and only if $X^{j}_{[0,T]}\perp\not\hskip 1.00006pt\perp X_{0}^{i}$. Since we already know from the CPDAG that $i$ and $j$ are adjacent and the graph is acyclic, $i\in\mathrm{an}_{j}$ implies $i\to j$.
∎

The full written out algorithm then looks like Algorithm 2 and is sound and complete by Corollary 3.3.

Figure: Algorithm 2 Robust Causal discovery for acyclic SDE models leveraging the initial values

###### Proof.

### A.10 Post-Processing Without Jointly Independent Initial Values

In real-world settings where there is no natural ‘starting point’ of the processes, we can typically not assume joint independence of initial values $X^{i}_{0}$ as ‘the process also happened before’.
In this subsection we extend Algorithm 2 to a version that is applicable in scenarios where we cannot assume the initial values $(X_{0}^{i})_{i\in[d]}$ to be jointly independent.

Figure: Algorithm 3 Robust causal discovery for acyclic SDE models.

The algorithm works similar as to Algorithm 2, except for the fact that for each pair $i,j\in V$ with an unoriented edge, we have to condition on the parents of $j$ when testing for the existence of edge $i\to j$ to prevent information flowing from $i$ to $j$. As it can happen that $j$ is connected to $k\neq i$ with another undirected edge, the set $\mathcal{K}_{i\to j}$ is required:

$$ $\displaystyle\mathcal{K}_{i\to j}:=\mathrm{pa}^{\mathcal{G}_{post}}_{j,known} \times\mathcal{P}(V_{j,un}),$ $$

where $\mathrm{pa}^{\mathcal{G}_{post}}_{j,known}:=\{k\in V\backslash\{i,j\}:(k_{0},j
_{1})\in\tilde{E}_{post},(j_{0},k_{1})\not\in\tilde{E}_{post}\}$, $V_{j,un}:=\{k\in V\backslash\{i\}:(k_{0},j_{1})\in\tilde{E}_{post},(j_{0},k_{1
})\in\tilde{E}_{post}\}$. $\mathcal{K}_{i\to j}$ is required to implement the procedure of conditioning on all known parents of j ($\mathrm{pa}^{\mathcal{G}_{post}}_{j,known}$) and in case there are unspecified other neighbors $V_{j,un}$ of $j$, loops through the options of either condition on them or don’t condition on them. In case $j$ has no other undirected adjacencies, this set only contains one element, the parents of $j$ except $i$.

###### Proposition A.12 .

Algorithm 3 is sound and complete assuming faithfulness, meaning if we denote the output graph of Algorithm 3 by $\mathcal{G}_{ext}=(V,E_{ext})$, the output edges $E_{ext}$ coincide with the edges of the true DAG $\mathcal{G}=(V,E)$.

###### Proof.

- ”$\supset$”:
We have to show $(i,j)\in E\implies(i,j)\in E_{ext}$. This direction is clear but will still be stated here. Since the original PC algorithm (see e.g. Spirtes et al. ()) is sound, $(i,j)\in\tilde{E}_{post}$
Since $(i,j)\in E$, $i_{0}$ and $j_{1}$ cannot be separated in the lifted graph $\forall K\subset V\backslash\{i,j\}$. Thus by the faithfulness assumption $X^{i}\perp\not\hskip 1.00006pt\perp_{s,h}^{+}X^{j}\mid X^{K}$ $\forall K\subset V\backslash\{i,j\}$ which holds especially for $K\in\mathcal{K}_{i\to j}$, hence $(i,j)\in E_{ext}$.
- ”$\subset$”:
Proof by contraposition: $(i,j)\in E_{ext}\implies(i,j)\in E$ is equivalent to $(i,j)\not\in E\implies(i,j)\not\in E_{ext}$.
Since the PC result $\mathcal{G}_{pre}=(V,E_{pre})$ has the correct skeleton and v-structures, the only option would be that wlog $i\neq j$ s.t. $j\to i\in\mathcal{G}$, $(i_{0},j_{1}),(j_{0},i_{1})\in\tilde{E}_{post}$ meaning the original PC was to detect an adjacency but unable to discard $(i,j)$. But then $j_{0}$ is separated from $i_{1}$ by the parents $\mathrm{pa}_{i}^{\mathcal{G}}$ of $i$, which do not contain $j$. Therefore, they are contained in $\mathcal{K}_{j\to i}$ as a single element and by the global Markov property, we obtain $X^{j}\operatorname{\perp\perp}_{s,h}^{+}X^{i}\mid X^{\mathrm{pa}_{i}^{\mathcal
{G}}}$, thus $(i,j)\not\in E_{ext}$.

∎

###### Proposition A.12 .

###### Proof.

### A.11 Consistency of the Conditional Independence Test

In this section, we argue, that under certain assumptions on our system, Theorem 1 of Doran et al. () still holds even without assuming the existence of a density.
At first, it is to be noted that the signature kernel eq. 3 is bounded on compact sets $C\subset\operatorname{\mathrm{BV}}(C(I,\mathbb{R}^{d_{x}}))$ of BV paths, an assumption that can be made as we evaluate our SDEs over finite, discrete time steps. Satisfying this boundedness, it defines an RKHS over these compact sets. In addition, the signature kernel is characteristic on these bounded sets (see Proposition 3.3 Cass et al. ()).
We use the notation $X:\left(\Omega,\mathcal{A},P\right)\rightarrow\left(\mathcal{X},\mathcal{A}_{
\mathcal{X}}\right)$, $Y:\left(\Omega,\mathcal{A},P\right)\rightarrow\left(\mathcal{Y},\mathcal{A}_{
\mathcal{Y}}\right)$, $Z:\left(\Omega,\mathcal{A},P\right)\rightarrow\left(\mathcal{Z},\mathcal{A}_{
\mathcal{Z}}\right)$ with $\mathcal{X},\mathcal{Y},\mathcal{Z}$ being compact subsets in $\operatorname{\mathrm{BV}}(C(I,\mathbb{R}^{n_{x}})),BV(C(I,\mathbb{R}^{n_{y}})
),BV(C(I,\mathbb{R}^{n_{z}}))$ and denote the RKHS’s defined by the signature kernel as $(\mathcal{H}_{\mathcal{X}},k_{\mathcal{X}}),(\mathcal{H}_{\mathcal{Y}},k_{
\mathcal{Y}}),(\mathcal{H}_{\mathcal{Z}},k_{\mathcal{Z}})$.

Our measure-theoretical argumentation now follows the definitions in Park & Muandet () and the original proof in Doran et al. ().
Let $H$ be a Banach space, $\mathcal{A}^{\prime}_{H}\subset\mathcal{A}_{H}$ a sub-$\sigma$-algebra of the Borel-$\sigma$-algebra.

###### Definition A.13 .

Let $F:(\Omega,\mathcal{A},P)\rightarrow(H,\mathcal{A}_{H})$ a Bochner $P$-integrable random variable. The *conditional expectation of $H$* w.r.t. $E$ is a Bochner $P$-integrable RV $F:(\Omega,\mathcal{A},P)\rightarrow(H,\mathcal{A}^{\prime}_{H})$ s.t.

$$ $\displaystyle\forall A\in\mathcal{A}^{\prime}_{H}\quad\int_{A}F\mathrm{d}P= \int_{A}F^{\prime}\mathrm{d}P\,.$ $$

Due to the existence and (almost sure) uniqueness, such an $F^{\prime}$ is denoted $\mathbb{E}[F\mid\mathcal{A}^{\prime}_{H}]$ with the notation $\mathbb{E}[F\mid Z]:=\mathbb{E}[F\mid\sigma(Z)]$ for a general $Z:(\Omega,\mathcal{A},P)\rightarrow(\mathcal{Z},\mathcal{A}_{\mathcal{Z}})$ which in turn can be used to define the *conditional expectation* $P(A\mid\mathcal{A}^{\prime}_{H}):=\mathbb{E}[\mathbbm{1}_{A}\mid\mathcal{A}^{
\prime}_{H}]$ and can be used to define the *conditional kernel mean embedding*

$$ $\displaystyle\mu_{X\mid Z}:=\mu_{P_{X\mid Z}}:=\mathbb{E}_{X\mid Z}[k_{ \mathcal{X}}(X,\,\cdot\,)\mid Z]$ $$

As proven in Theorem 5.8 of Park & Muandet () under the assumption, that $P_{X\mid Z}:=\mathbb{E}[\,\cdot,\mid Z]$ admits a regular version (meaning its can be written as a transition kernel, see Park & Muandet ()) and $k_{\mathcal{X}}\otimes k_{\mathcal{Y}}$
is characteristic, it holds that

$$ $\displaystyle X\operatorname{\perp\perp}Y\mid Z\>\Leftrightarrow\>\lVert\mu_{ XY\mid Z}-\mu_{X\mid Z}\otimes\mu_{Y\mid Z}\rVert_{\mathcal{H}_{\mathcal{X}} \otimes\mathcal{H}_{\mathcal{Y}}}=0\>\Leftrightarrow\>\lVert\mu_{XY\mid Z} \otimes\mu_{Z}-\mu_{X\mid Z}\otimes\mu_{Y\mid Z}\otimes\mu_{Z}\rVert_{\mathcal {H}_{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{\mathcal{ Z}}}=0$ $$

where the second equivalence is trivial. Let now $\sigma\in\mathcal{S}_{n}$ be a permutation such that its permutation matrix $M_{\sigma}=(m_{ij})$

$$ $m_{ij}=\begin{cases}1&\>\text{if }j=\sigma(i)\>,\\ 0&\>\text{if }j\neq\sigma(i)\>.\end{cases}$ $$

satisfies $\mathrm{Tr}(M_{\sigma})=0$. The empirical estimates for the given sample $\{(x_{i},y_{i},z_{i})\}_{i\in[n]}$ and its permuted version under $\sigma$

$$ $\displaystyle\hat{\mu}_{P_{X,Y,Z}}$ $\displaystyle\coloneqq\frac{1}{n}\sum_{i=1}^{n}k_{\mathcal{X}}(x_{i},\cdot) \otimes k_{\mathcal{Y}}(y_{i},\cdot)\otimes k_{\mathcal{Z}}(z_{i},\cdot),$ (16) $\displaystyle\hat{\mu}_{P^{\prime}_{X,Y,Z}}$ $\displaystyle\coloneqq\frac{1}{n}\sum_{i=1}^{n}k_{\mathcal{X}}(x_{i},\cdot) \otimes k_{\mathcal{Y}}(y_{\sigma(i)},\cdot)\otimes k_{\mathcal{Z}}(z_{i}, \cdot).$ (17) $$

Similar to Laumann et al. (), can now write the distance between test statistic and the approximated test statistic in terms of the sum of the distance between the empirical estimates and its respective counterparts by using the inverse and the regular triangle inequality

$$ $\displaystyle\left\lvert\lVert\mu_{XY\mid Z}\otimes\mu_{Z}-\mu_{X\mid Z} \otimes\mu_{Y\mid Z}\otimes\mu_{Z}\rVert_{\mathcal{H}_{\mathcal{X}}\otimes \mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{\mathcal{Z}}}-\lVert\hat{\mu}_{P_ {X,Y,Z}}-\hat{\mu}_{P^{\prime}_{X,Y,Z}}\rVert_{\mathcal{H}_{\mathcal{X}} \otimes\mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{\mathcal{Z}}}\right\rvert$ $\displaystyle\leq{}$ $\displaystyle\left\lvert\lVert\mu_{XY\mid Z}\otimes\mu_{Z}-\mu_{X\mid Z} \otimes\mu_{Y\mid Z}\otimes\mu_{Z}-\hat{\mu}_{P_{X,Y,Z}}+\hat{\mu}_{P^{\prime} _{X,Y,Z}}\rVert_{\mathcal{H}_{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}} \otimes\mathcal{H}_{\mathcal{Z}}}\right\rvert$ $\displaystyle\leq{}$ $\displaystyle\lVert\mu_{XY\mid Z}\otimes\mu_{Z}-\hat{\mu}_{P_{X,Y,Z}}\rVert_{ \mathcal{H}_{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{ \mathcal{Z}}}+\lVert\mu_{X\mid Z}\otimes\mu_{Y\mid Z}\otimes\mu_{Z}-\hat{\mu}_ {P^{\prime}_{X,Y,Z}}\rVert_{\mathcal{H}_{\mathcal{X}}\otimes\mathcal{H}_{ \mathcal{Y}}\otimes\mathcal{H}_{\mathcal{Z}}}.$ $$

By Park & Muandet (), the first term tends to zero for $n\rightarrow\infty$. With a similar line of argumentation as in Doran et al. () on can establish a majorant for the estimated HSCIC, namely

$$ $\displaystyle\lVert\hat{\mu}_{P_{X,Y,Z}}-\hat{\mu}_{P^{\prime}_{X,Y,Z}}\rVert_ {\mathcal{H}_{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{ \mathcal{Z}}}$ $\displaystyle={}$ $\displaystyle\left\lVert\frac{1}{n}\sum_{i=1}^{n}k_{\mathcal{X}}(x_{i},\cdot) \otimes k_{\mathcal{Y}}(y_{i},\cdot)\otimes k_{\mathcal{Z}}(z_{i},\cdot)-\frac {1}{n}\sum_{i=1}^{n}k_{\mathcal{X}}(x_{i},\cdot)\otimes k_{\mathcal{Y}}(y_{ \sigma(i)},\cdot)\otimes k_{\mathcal{Z}}(z_{i},\cdot)\right\rVert_{\mathcal{H} _{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{\mathcal{Z}}}$ $\displaystyle={}$ $\displaystyle\left\lVert\frac{1}{n}\sum_{i=1}^{n}k_{\mathcal{X}}(x_{i},\cdot) \otimes k_{\mathcal{Y}}(y_{i},\cdot)\otimes k_{\mathcal{Z}}(z_{i},\cdot)-\frac {1}{n}\sum_{i=\sigma(j),j=1}^{n}k_{\mathcal{X}}(x_{\sigma^{-1}(i)},\cdot) \otimes k_{\mathcal{Y}}(y_{i},\cdot)\otimes k_{\mathcal{Z}}(z_{\sigma^{-1}(i)} ,\cdot)\right\rVert_{\mathcal{H}_{\mathcal{X}}\otimes\mathcal{H}_{\mathcal{Y}} \otimes\mathcal{H}_{\mathcal{Z}}}$ $\displaystyle\leq{}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\left\lVert\left(k_{\mathcal{X}}(x_{i}, \cdot)-k_{\mathcal{X}}(x_{\sigma^{-1}(i)},\cdot)\right)\otimes k_{\mathcal{Y}} (y_{i},\cdot)\otimes\left(k_{\mathcal{Z}}(z_{i},\cdot)-k_{\mathcal{Z}}(z_{ \sigma^{-1}(i)},\cdot)\right)\right\rVert_{\mathcal{H}_{\mathcal{X}}\otimes \mathcal{H}_{\mathcal{Y}}\otimes\mathcal{H}_{\mathcal{Z}}}$ $\displaystyle\leq{}$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\underbrace{\left\lVert k_{\mathcal{X}}( x_{i},\cdot)-k_{\mathcal{X}}(x_{\sigma^{-1}(i)},\cdot)\right\rVert_{\mathcal{H }_{\mathcal{X}}}}_{\leq 2M_{x}}\underbrace{\left\lVert k_{\mathcal{Y}}(y_{i}, \cdot)\right\rVert_{\mathcal{H}_{\mathcal{Y}}}}_{\leq M_{y}}\left\lVert\left(k _{\mathcal{Z}}(z_{i},\cdot)-k_{\mathcal{Z}}(z_{\sigma^{-1}(i)},\cdot)\right) \right\rVert_{\mathcal{H}_{\mathcal{Z}}}$ $\displaystyle\leq{}$ $\displaystyle 2M_{x}M_{y}\left(\frac{1}{n}\operatorname{\mathrm{Tr}}(PD^{RKHS} )\right)$ $$

where the assumptions can be made due to the boundedness of the kernels.
Hence, under the assumption that

$$ $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\delta_{(x_{i},y_{\sigma(i)},z_{i})} \rightarrow P_{X\mid Z}\otimes P_{Y\mid Z}\otimes P_{Z}$ (18) $$

a statement that holds true in the case of densities (see supplementary material of Janzing et al. ()), since $P_{\mid Z}$ is regular, $\hat{\mu}_{P^{\prime}_{X,Y,Z}}\rightarrow\mu_{X\mid Z}\otimes\mu_{Y\mid Z}
\otimes\mu_{Z}$.

Hence under the assumptions that $P_{\cdot\mid Z}$ admits a regular version (an assumption also made by Laumann et al. ()), the considered random variables operate on compact sets of BV-paths and eq. 18 holds, the test is consistent if and only if $\frac{1}{n}\operatorname{\mathrm{Tr}}(PD^{RKHS})\rightarrow 0$.
This is the exact same statement as in Theorem 1 of Doran et al. ().

###### Definition A.13 .

## Appendix B Experiments and Evaluation

This section provides details about the experimental framework, elaborating on the data generation process, the parameters and architectures employed, and presents additional results.

**Table 4: Performance of *future-extended h-locally conditionally independence criterion* in the bivariate setting for the two heuristic $\sigma_{1}$ and $\sigma_{2}$ for different sample sizes. The SHD is measured for different values of $s$ in $\operatorname{\perp\perp}_{s,t}^{+}$ with $t=1$, defining different lengths of intervals in the past and future. Overall the performance is best of shorter values of the past intervals compared to the future ($s<\frac{t}{2}$).**
|  | $\sigma_{1}$ | $\sigma_{2}$ |  |  |
| --- | --- | --- | --- | --- |
| $(\times 10^{2})$ | $n=400$ | $n=600$ | $n=400$ | $n=600$ |
| $s=0.1$ | $\mathbf{76.2\pm 2.5}$ | $\mathbf{58.4\pm 2.7}$ | $14.9\pm 2.0$ | $15.8\pm 2.0$ |
| $s=0.2$ | $77.3\pm 2.7$ | $72.2\pm 2.8$ | $\mathbf{12.9\pm 1.7}$ | $\mathbf{10.9\pm 1.7}$ |
| $s=0.3$ | $80.2\pm 2.3$ | $76.2\pm 2.6$ | $21.8\pm 2.2$ | $17.8\pm 2.0$ |
| $s=0.4$ | $89.1\pm 2.0$ | $92.0\pm 2.1$ | $31.7\pm 2.4$ | $25.7\pm 2.2$ |
| $s=0.5$ | $89.1\pm 2.0$ | $98.1\pm 2.0$ | $38.6\pm 2.5$ | $27.7\pm 2.2$ |
| $s=0.6$ | $82.4\pm 1.9$ | $96.0\pm 2.1$ | $60.0\pm 2.7$ | $48.5\pm 2.9$ |
| $s=0.7$ | $102.0\pm 1.2$ | $104.0\pm 2.0$ | $80.2\pm 2.2$ | $61.4\pm 2.7$ |
| $s=0.8$ | $104.0\pm 1.4$ | $95.0\pm 1.5$ | $81.1\pm 2.0$ | $87.1\pm 2.4$ |
| $s=0.9$ | $95.0\pm 1.6$ | $104.4\pm 1.8$ | $102.0\pm 1.2$ | $101.0\pm 2.1$ |

### B.1 Kernel Heuristic and Interval Choice

In kernel methods, besides selecting the right kernel function, the correct choice of its parameters is vital.
A reasonable choice for radially symmetric kernels $k\left(x_{i},x_{j}\right)=f(\lVert x_{i}-x_{j}\rVert/\sigma)$, $f:\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}$ is to choose $\sigma$ to lie in the same order of magnitude as the pairwise distances $\left(\lVert x_{i}-x_{j}\rVert\right)_{i,j}$.
Choosing the bandwidth for the RBF kernel in our setting with samples $\left(x_{i}\right)_{i\in[n]}$, $x_{i}\in\mathcal{C}^{0}(I,\mathbb{R}^{d})$ is not straight forward, as the RBF kernel ultimately acts on signatures and there is no immediate way in obtaining a ‘typical scale’ of pairwise distances of signatures. To still set the bandwidth purely based on observed data, there are two natural ways to implement a median heuristic based on two different collections of pairwise distances. The first option is to choose

$$ $\sigma_{1}=\mathrm{Median}\bigl{\{}\lVert x_{i}-x_{j}\rVert_{L^{2}}\,\big{|}\, i,j\in[n]\bigr{\}}\>,$ $$

where $\|\cdot\|_{L^{2}}$ is the is the $L^{2}$ norm of entire paths, i.e., integrated over the time interval $I$.
Given that in practice we observe the paths $x_{i}$ still at fixed time points $t_{1}<\ldots<t_{m}$ within the interval $I$, we can also consider the heuristic

$$ $\sigma_{2}=\mathrm{Median}\bigl{\{}\lVert x_{i,t_{l}}-x_{j,t_{k}}\rVert_{2}\, \big{|}\,i,j\in[n],l,k\in[m]\bigr{\}}\>,$ $$

where $\|\cdot\|_{2}$ is the Euclidean norm in $\mathbb{R}^{d}$ and we consider pairwise distances of all observations across all paths and time points.
In early benchmarks, we have empirically found $\sigma_{2}$ to yield a better heuristic by comparing the performance of both heuristics by evaluation of the SHD in a causal discovery task in the linear bivariate setting $X^{1}\rightarrow X^{2}$ without diffusion interaction with $a_{21}\sim\mathcal{U}\left((-2.5,-1]\cup[1,2.5)\right)$, $a_{11},a_{22}\sim\mathcal{U}\left((-0.5,0.5)\right)$ and $\sigma_{i}\sim\mathcal{U}\left([0,0.5)\right)$. The experiment is performed for different choices of intervals $[0,s],[s,s+h]$ in the future-extended $h$-local conditional independence. As can be seen from Table 4, $\sigma_{2}$ and the intervals $[0,0.1],[0.1,1]$ seem to perform best on average and are therefore used in the rest of the experiments. In all the experiments of the paper, in the implementation of the signature kernel we use a depth parameter of $4$, the RBF kernel and we add time as an extra dimension to the signature kernel.

Figure: Figure 6: Performance in terms of type I and type II errors of different unconditional (HSICb and HSICγ) and conditional tests (KCIPT, KCITb, KCITγ, SDCIT). Among the conditional tests, SDCIT outperforms KCIPT in terms of type I error, while performig in a similar range in terms of type II error. Among the unconditional tests, HSICb outperforms HSICγ
Refer to caption: x5.png

### B.2 Empirical Comparison of Different Tests

Different kernel-based conditional and unconditional tests are compared. HSIC with bootstrap (HSICb) and with $\gamma$-distribution approximation (HSICγ) (Gretton et al., ) were implemented for the unconditional setting, while KCIPT (Doran et al., ), KCIT with bootstrap (KCITb), KCIT with $\gamma$ approximation (KCITγ) (Zhang et al., ) and SDCIT for the conditional setting (Lee & Honavar, ). Figure 6 shows that SDCIT and HSICb are overall the best performing ones in terms of type I and type II error. In this experiment, we sampled 100 different linear SDEs and we sample 400 different trajectories within time interval [0,1]. The number of bootstrap samples over permutation test is set to 100, the number of permutations for a single permutation test to 20000, the number of null samples via Monte Carlo from all values in permutation test to 20000 and 1000 for HSIC, the number of null samples for KCIT-bootstrap to 20000 and the number of null samples for SDCIT is set to 1000. These are the tests parameters also used in the rest of the experiments.

### B.3 Additional Experimental Results for the Bivariate Case

Figure 7 shows that in the simple unconditional bivariate setting of the power analysis in Section 4 our test maintains high power up to substantial levels of missingness and only starts degrading once more than 80% of observations are dropped at random starting from $64$ time steps over the interval $[0,1]$ for $1000$ different SDEs.

Figure: Figure 7: Test power over the missingness level (percentage of observations that is being dropped at random). While both our method as well as Laumann et al. () maintain their power up to high levels of missingness, our test still consistently outperforms the baseline. Lines show means over 1000 settings in the power analysis in Section 4 and shaded regions show standard errors.
Refer to caption: x6.png

### B.4 Building Blocks for Causal Discovery

The performance of the conditional and unconditional independence tests was evaluated using widely recognized causal structures that are at the foundation of causal discovery. Specifically, we examined scenarios involving two variables, both connected and unconnected, and three-variable configurations comprising chain, fork, and collider structures. The analysis of type I and type II errors, as shown in Table 5, indicates robust performance for both conditional SDCIT and unconditional HSIC tests. These results refer to the symmetric conditional independence $\operatorname{\perp\perp}_{sym}$, while future-extended h-local conditional independence are presented in Table 7. The tables show that the test is able to reconstruct the presence and the directionality of the edges overall all building blocks of causal structure learning.
We also perform the same experiments for the setting in which the causal dependence only comes from the diffusion. Referring to the linear SDE in eq. 5, $a_{ij}=a_{ji}=0$, $a_{11},a_{22}\sim\mathcal{U}[-0.5,0.5]$, $\sigma_{12},\sigma_{21}\sim\mathcal{U}[(-2.5,-1.0)\cup(1.0,2.5)]$ and $\sigma_{11},\sigma_{22}\sim\mathcal{U}[-0.5,0.5]$. $b_{i}$ and $d_{i}$ are sampled from $\mathcal{U}[-0.1,0.1]$ and $\mathcal{U}[-0.2,0.2]$. Table 6 includes the results for all the different causal structures.

**Table 5: Test performance on the key building blocks for causal discovery. Here, $\operatorname{\perp\perp}$ refers to $\operatorname{\perp\perp}_{\mathrm{sym}}$ except for $X^{1}_{p}\operatorname{\perp\perp}X^{1}_{f}$ which means $X^{1}_{[0,t/2]}\operatorname{\perp\perp}X^{1}_{[t/2,t]}$. We use SDCIT for CI tests and bootstrapped HSIC unconditional independence. We simulate $40$ samples of $512$ trajectories for $10$ different SDEs (400 tests per column) with $25\%$ of the 128 original observations dropped uniformly at random. The error represent type I or type II error when the null hypothesis $H_{0}$ should be rejected (✗) or accepted (✓), respectively.**
| graph | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$ | $X^{1}$$X^{2}$ | $X^{1}$$X^{2}$ |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $H_{0}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}_{p}\!\operatorname{\perp\perp}\!X^{1}_{f}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}$ |
| should reject | ✓ | ✗ | ✓ | ✓ | ✗ | ✓ | ✗ | ✓ |
| error | $0.0000$ | $0.0075$ | $0.0000$ | $0.0000$ | $0.9600$ | $0.0000$ | $0.0450$ | $0.0000$ |

**Table 6: Test performance on the key building blocks for causal discovery. Here, $\operatorname{\perp\perp}$ refers to $\operatorname{\perp\perp}_{\mathrm{sym}}$ except for $X^{1}_{p}\operatorname{\perp\perp}X^{1}_{f}$ which means $X^{1}_{[0,t/2]}\operatorname{\perp\perp}X^{1}_{[t/2,t]}$ in a setting where the dependence only comes from the diffusion. We use SDCIT for CI tests and bootstrapped HSIC unconditional independence. We simulate $40$ samples of $512$ trajectories for $10$ different SDEs (400 tests per column) with $25\%$ of the 128 original observations dropped uniformly at random. The error represent type I or type II error when the null hypothesis $H_{0}$ should be rejected (✗) or accepted (✓), respectively.**
| graph | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$$X^{2}$$X^{3}$ | $X^{1}$ | $X^{1}$$X^{2}$ | $X^{1}$$X^{2}$ |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $H_{0}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{1}_{p}\!\operatorname{\perp\perp}\!X^{1}_{f}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}$ |
| should reject | ✓ | ✗ | ✓ | ✓ | ✗ | ✓ | ✗ | ✓ |
| error | $0.0000$ | $0.0072$ | $0.0000$ | $0.0000$ | $0.9600$ | $0.00$ | $0.07$ | $0.01$ |

**Table 7: Test performance on the chain, fork and collider causal structures for the *future-extended h-locally conditionally independence criterion*. Here, $\operatorname{\perp\perp}$ refers to $\operatorname{\perp\perp}_{s,h}^{+}$ and we use SDCIT. We simulate $10$ samples of $512$ trajectories for $10$ different SDEs ($100$ tests per column). The error represent type I or type II error when the null hypothesis $H_{0}$ should be rejected (✗) or accepted (✓), respectively.**
| graph | $X^{1}$$X^{2}$$X^{3}$ |  |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| $H_{0}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{1}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{1}|X^{2}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{3}|X^{1}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{2}|X^{1}$ |
| should reject | ✗ | ✓ | ✗ | ✗ | ✗ | ✓ |
| error | 0.07 | 0.33 | 0.04 | 0.04 | 0.03 | 0.00 |
| graph | $X^{1}$$X^{2}$$X^{3}$ |  |  |  |  |  |
| $H_{0}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{1}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{1}|X^{2}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{3}|X^{1}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{2}|X^{1}$ |
| should reject | ✓ | ✗ | ✗ | ✗ | ✗ | ✓ |
| error | 0.12 | 0.02 | 0.04 | 0.02 | 0.04 | 0.06 |
| graph | $X^{1}$$X^{2}$$X^{3}$ |  |  |  |  |  |
| $H_{0}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{2}|X^{3}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{1}|X^{3}$ | $X^{1}\!\operatorname{\perp\perp}\!X^{3}|X^{2}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{1}|X^{2}$ | $X^{2}\!\operatorname{\perp\perp}\!X^{3}|X^{1}$ | $X^{3}\!\operatorname{\perp\perp}\!X^{2}|X^{1}$ |
| should reject | ✗ | ✓ | ✓ | ✓ | ✓ | ✗ |
| error | 0.08 | 0.4 | 0.56 | 0.46 | 0.36 | 0.16 |

### B.5 Causal Discovery

We tested SCOTCH using various sparsity parameters and epochs to identify the optimal configuration. Table 8 confirms that the configuration with $\lambda=200$ and $n_{e}=2000$ outperforms others in the tested settings, and hence is chosen for the comparison with our method in Table 2. The table also reveals that SCOTCH’s performance is not robust to changes in $\lambda$, showing substantial variations under the same number of epochs, optimizer, and learning rate with different $\lambda$ values. These results demonstrate that the SCOTCH’s performance can significantly deteriorate with changes in the parameter $\lambda$, while keeping the number of epochs constant. For example, when $\lambda$ is adjusted from 200 to 100, the Structural Hamming Distance (SHD, scaled by $10^{2}$) increases from an average of 538 to 10,275 for $n_{e}=2000$.

**Table 8: SHD ($\times 10^{2}$) of discovered graphs, for $200$ samples and $d$ nodes in the SDE model. The mean and standard error of SCOTCH for different values of $\lambda$ and $n_{e}$ for 20 different SDEs is presented. The best performing setting for higher dimensional graphs is $\lambda=200$ and $n_{e}=2k$.**
| $\times 10^{2}$ | $(\lambda,n_{e})$ | $d=3$ | $d=5$ | $d=10$ | $d=20$ | $d=50$ |
| --- | --- | --- | --- | --- | --- | --- |
| SCOTCH | 100, 2k | $188\pm 28$ | $417\pm 86$ | $250\pm 61$ | $1525\pm 1160$ | $10275\pm 6176$ |
| SCOTCH | 200, 2k | $110\pm 21$ | $270\pm 48$ | $530\pm 223$ | $\mathbf{370\pm 174}$ | $\mathbf{538\pm 70}$ |
| SCOTCH | 1, 1k | $\mathbf{20\pm 11}$ | $\mathbf{85\pm 20}$ | $\mathbf{375\pm 28}$ | $1455\pm 103$ | $6095\pm 153$ |
| SCOTCH | 50, 1k | $330\pm 27$ | $1255\pm 54$ | $5895\pm 107$ | $1340\pm 243$ | $2994\pm 941$ |
| SCOTCH | 200, 1k | $400\pm 17$ | $1370\pm 29$ | $6425\pm 80$ | $705\pm 77$ | $7863\pm 2670$ |

### B.6 Conditional Independence Testing Beyond the SDE Model

In this section, we aim to demonstrate that our proposed non-parametric conditional independence (CI) test extends beyond the setting of SDEs and even beyond semi-martingales, rendering it more broadly applicable than most existing methods, which usually assume independent (additive) noise.

#### Fractional Brownian Motion.

To illustrate this, we start by applying our test to a stochastic process driven by fractional Brownian motions (fBMs).
Fractional Brownian motion is a generalization of standard Brownian motion, a continuous time Gaussian process $B_{H}(t)$ with zero mean and $\mathbb{E}[B_{H}(t)B_{H}(s)]=\frac{1}{2}\left(\lvert t\rvert^{2H}+\lvert s
\rvert^{2H}-\lvert t-s\rvert^{2H}\right)$ that incorporates a parameter $H\in(0,1)$, the so called *Hurst parameter*, which governs long-range dependencies. For

- •
$H=0.5$ it reduces to the regular Brownian Motion.
- •
$H>0.5$ increments are positively correlated.
- •
$H<0.5$ increments are negatively correlated.

Specifically, we consider the following system of stochastic differential equations driven by fractional Brownian motions

$$ $\displaystyle\begin{pmatrix}dX_{t}\\ dY_{t}\end{pmatrix}=\begin{pmatrix}0&0\\ a_{21}&0\end{pmatrix}\begin{pmatrix}X_{t}\\ Y_{t}\end{pmatrix}dt+\begin{pmatrix}d_{1}&0\\ 0&d_{2}\end{pmatrix}\begin{pmatrix}dW^{1}_{t}\\ dW^{2}_{t}\end{pmatrix}\>,$ (19) $$

with $W^{1}_{t}$, $W^{2}_{t}$ independent fBMs with Hurst parameter $H$ and $a_{21},d_{1},d_{2}\sim\mathcal{U}([-2,2])$.
To quantify the effectiveness of our test in this context, we measure the test power for different strengths of interaction for different Hurst parameters.
Figure 8 shows this power analysis for $500$ random settings of eq. 19.
The power of our test steadily approaches 1 as the interaction strength increases to 1 regardless of $H$.
We note that since the above example does not have independent increments, Algorithm 1 is not applicable anymore for causal discovery in such settings.
Instead, the fractional Brownian motion example merely highlights the usefulness of our CI test independently from its application in causal discovery in the SDE model.

#### Functional data.

For further corroboration of our claims, we turn to the following example from functional data analysis (FDA), taken from Laumann et al. (), to test causal discovery in a non SDE setting.
After drawing a DAG $\mathcal{G}=(V,E)$, processes for source vertices, i.e., vertices $i=v_{i}$ without parents in $\mathcal{G}$, are generated according to

$$ $\displaystyle X^{i}(t)=\sum_{m=1}^{M}c^{i}_{m}\phi_{m}(t)+\epsilon_{t}^{i}$ $$

with Fourier basis functions $\phi_{1}(t)=1$, $\phi_{2}(t)=\sqrt{2}\sin\left(2\pi t\right)$, $\phi_{3}(t)=\sqrt{2}\cos\left(2\pi t\right)$ and so on, the weights $c^{i}_{m}\sim\mathcal{N}(0,1)$, and the additive noise $\epsilon_{t}^{i}\sim\mathcal{N}(0,1)$. For all vertices $i$ with at least one parent, we set

$$ $\displaystyle X^{i}(t)=a\sum_{k\in\mathrm{pa}_{i}^{\mathcal{G}}}\int_{0}^{t}X^ {k}(s)\beta^{k}(s,t)ds+\epsilon_{t}^{i}$ $$

with $\beta^{k}(s,t)=8(s-c^{k}_{1})^{2}-8(t-c^{k}_{2}))^{2}$, $c^{k}_{1},c^{k}_{2}\sim\mathcal{U}_{[0,1]}$ and $a=1$.

We draw 40 DAGs with 40 distinct FD generating mechanisms for each dimension $d\in\{3,5,10,20,50\}$.
Table 9 shows that SCOTCH cannot capture the underlying dependencies, as it is tailored specifically to the SDE model. It is thus misspecified for such a function generating mechanism.
Under such (arguably mild) misspecification, SCOTCH is even outperformed by PCMCI, whereas our method still performs much better than both baselines.

Figure: Figure 8: Test power over for different fBM-driven processes eq. 19 for different H. The CI test picks up dependence for as little as $100$ samples given a certain amount of interaction-strength. Lines show means over $500$ settings and shaded regions show standard errors.
Refer to caption: extracted/5924253/figures/fbm_final.png

**Table 9: SHD ($\times 10^{2}$) comparison of SigKer to the baselines PCMCI and SCOTCH (with $\lambda=200$ and $n_{e}=2000$) in the functional data setting. We report mean and standard errors over 40 instances.**
| $\times 10^{2}$ | $d=3$ | $d=5$ | $d=10$ | $d=20$ | $d=50$ |
| --- | --- | --- | --- | --- | --- |
| $\operatorname{\perp\perp}_{\mathrm{sym}}$ | $40\pm 8$ | $22\pm 10$ | $509\pm 25$ | $1542\pm 60$ | $6786\pm 223$ |
| $\operatorname{\perp\perp}_{\mathrm{sym}}$+p.p. | $42\pm 8$ | $141\pm 13$ | $489\pm 28$ | $1498\pm 56$ | $6529\pm 250$ |
| PCMCI | $110\pm 30$ | $320\pm 42$ | $1150\pm 121$ | $6518\pm 293$ | $7799\pm 250$ |
| SCOTCH | $558\pm 14$ | $1861\pm 43$ | $8567\pm 63$ | $35467\pm 833$ | $105165\pm 10934$ |

### B.7 Real-World Pairs Trading Example

In the pairs trading experiment, stock price data is downloaded from Yahoo Finance for a predefined list of stocks over a specific period, divided into training (1st January 2010 to 31st December 2011) and trading intervals (1st January 2012 to 31st December 2012). The chose stocks are Trinity Industries (TRN), Brandywine Realty Trust (BDN), Commercial Metals Company (CMC), The New York Times Company (NYT), New York Community Bancorp (NYCB), The Wendy’s Company (WEN), CNX Resources Corporation (CNX). Logarithms of the stock prices are computed to stabilize variance and normalize the prices. During the training period, pairs of stocks are selected based on various statistical tests: Cointegration Test (Engle-Granger) to determine if a pair of stocks is cointegrated, suggesting a long-term equilibrium relationship; Augmented Dickey-Fuller (ADF) Test to assess the stationarity of the spread (ratio) of a pair’s prices; and Granger Causality Test to check if the price of one stock in the pair can predict the price of the other. Pairs are selected based on the p-values from these tests, with a threshold of $0.05$ determining significance. In the trading phase, a rolling window is used to continuously recalculate the mean and standard deviation of the spread between the selected pairs. Trades are initiated based on the z-score of the spread, where a high z-score indicates a short position and a low z-score indicates a long position. Positions are managed and closed when the spread returns to a predefined z-score threshold. The strategy’s performance is evaluated based on total return, annual percentage rate (APR), Sharpe ratio, maximum drawdown (maxDD), and maximum drawdown duration (maxDDD).

### B.8 Evaluation Metrics

Following common conventions in the causal discovery literature, we evaluate our algorithms using the following metrics.

#### Structural Hamming Distance (SHD).

Given two directed graphs $\mathcal{G}_{1}=\left(V,E_{1}\right)$, $\mathcal{G}_{2}=\left(V,E_{2}\right)$ over the same nodes with different sets of edges and let $A_{1},A_{2}\in\{0,1\}^{n\times n}$ be their adjacency matrices. Then we define the structural hamming distance (SHD) as

$$ $\displaystyle\operatorname{SHD}=\sum_{i,j\in[n]}\lvert(A_{1})_{ij}-(A_{2})_{ij}\rvert$ $$

#### Normalized SHD.

In order to compare the recovery performance of our algorithms for different types of graph-sizes, the normalized SHD is defined as

$$ $\displaystyle\operatorname{SHD_{norm}}=\frac{\operatorname{SHD}}{d(d-1)}$ $$

## Appendix C Baselines

In Section 4, we benchmark our method against other baselines, including Granger causality, CCM, PCMCI, and Laumann (Laumann et al., ). For the latter, we used their implementation provided in the causal-fda package.
For the Granger-implementation for two variables ($d=2$), we used Seabold & Perktold (), for CCM we used Javier (), and for PCMCI we used the tigramite package (Runge et al., ).
In PCMCI, tests for edges are conducted by applying distance correlation-based independence tests (Székely et al., ) between the variables’ residuals after regressing out other nodes using Gaussian processes. For SCOTCH implementation (Wang et al., ), we use the package causica. For SCOTCH, we always use a learning rate of $0.001$ and keep the same default parameters for the learning algorithm.

## Appendix D Details on the Partially Observed Setting

Since in the partially observed setting there could in principle be infinitely many unobserved variables, special tools are required.
Usually, the graphical framework of *maximal ancestral graphs* (MAGs) (loosely speaking DAGs with also bidirected edges) is used, which
encodes ancestral relations in the DAG $\mathcal{G}$ and aims at graphically representing conditional independencies implied by $\mathcal{G}$ involving only the (marginal of the) observed variables.
The unique MAG $M^{\mathcal{G}}$ corresponding to the partially observed DAG $\mathcal{G}$ can be constructed via Zhang () the following rules.

- •
For $v_{1},v_{2}\in V_{obs}$, $v_{1}\operatorname{\text{---}}v_{2}$ in $M^{\mathcal{G}}$ if and only if there exists an inducing path relative to $V_{L}$ between $v_{1},v_{2}$ in $\mathcal{G}$.
- •
For each pair of adjacent vertices $v_{1}\operatorname{\text{---}}v_{2}$ in $M^{\mathcal{G}}$ we orient the edge between them as follows:
–
$v_{1}\rightarrow v_{2}$ if $v_{1}\in\mathrm{an}^{\mathcal{G}}_{v_{2}}$ and $v_{2}\not\in\mathrm{an}^{\mathcal{G}}_{v_{1}}$
–
$v_{1}\leftarrow v_{2}$ if $v_{2}\in\mathrm{an}^{\mathcal{G}}_{v_{1}}$ and $v_{1}\not\in\mathrm{an}^{\mathcal{G}}_{v_{2}}$
–
$v_{1}\leftrightarrow v_{2}$ if $v_{2}\not\in\mathrm{an}^{\mathcal{G}}_{v_{1}}$ and $v_{1}\not\in\mathrm{an}^{\mathcal{G}}_{v_{2}}$

An *inducing path* is a path with all colliders on the path being in $\mathrm{an}^{\mathcal{G}}_{\{v_{1},v_{2}\}}\cap V_{obs}$ and all other nodes on the path in $V_{L}$.
The obtained MAG $M^{\mathcal{G}}$ therefore preserves causal ancestral relations with respect to $\mathcal{G}$.
The *Fast Causal Inference* (FCI) algorithm (Zhang, ) run on the observed marginal outputs an equivalence class of MAGs, known as a *Partial Ancestral Graph* (PAG), which captures the same adjacencies but leaves some endpoints unoriented as there are multiple MAGs with the same conditional independence relations.
For example, in Figure 4 FCI cannot rule out an additional latent confounder between $A$ and $C$ such that the induced MAG has the same m-separation properties.
Unlike FCI, by leveraging the direction of time, we can recover the full underlying MAG $M^{\mathcal{G}}$, which is substantially more informative than FCI’s PAG.
First, we establish adjacencies analogous to the first part of Algorithm 2 (or simply by running FCI with our symmetric criterion) and then follow the steps in the construction of the MAG from a DAG to orient these adjacent edges $v_{i}\operatorname{\text{---}}v_{j}$, $v_{i}\neq v_{j}\in V_{obs}$:

- •
$v_{i}\rightarrow v_{j}$ if $X_{0}^{i}\perp\not\hskip 1.00006pt\perp_{\mathrm{sym}}X_{[0,T]}^{j}$ and $X_{0}^{j}\operatorname{\perp\perp}_{\mathrm{sym}}X_{[0,T]}^{i}$
- •
$v_{i}\leftarrow v_{j}$ if $X_{0}^{i}\operatorname{\perp\perp}_{\mathrm{sym}}X_{[0,T]}^{j}$ and $X_{0}^{j}\perp\not\hskip 1.00006pt\perp_{\mathrm{sym}}X_{[0,T]}^{i}$
- •
$v_{i}\leftrightarrow v_{j}$ if $X_{0}^{i}\perp\not\hskip 1.00006pt\perp_{\mathrm{sym}}X_{[0,T]}^{j}$ and $X_{0}^{j}\perp\not\hskip 1.00006pt\perp_{\mathrm{sym}}X_{[0,T]}^{i}$

We note that these criteria precisely encode the (non-)ancestral relationships in the MAG construction rules.
Uniqueness of the obtained result follows from Richardson & Spirtes (, Corollary 3.10) as two ancestral graphs $\mathcal{G}_{1}$, $\mathcal{G}_{2}$ are equal, if they share the same adjacencies and ancestral relations.

We have thus constructed a sound and complete algorithm (assuming access to a CI oracle) to recover the unique MAG $M^{\mathcal{G}}$ for the ground truth DAG $\mathcal{G}$ in the partially observed SDE model.
To illustrate our causal discovery algorithm in practice, we conduct experiments on an SDE model with a graph $\mathcal{G}$ as depicted in Figure 4 where $V_{obs}=\{A,B,C,D\}$, $V_{L}=\{U\}$, $E=\{(A,C),(B,C),(U,C),(U,D)\}$.
While FCI run on data coming from a Markov model with this DAG is able to detect the adjacency structure, it is unable to decide whether $A$ is truly an ancestor of $C$.
It will correctly infer the arrow-head into $C$, but cannot rule out an arrow-head at $A$, see Figure 4.
Neural network approaches like SCOTCH that directly infer the functional relationships by model fitting are not able to deal with such partially observed settings and are typically expected to predict edges that do not exists in the ground truth graph, since they can not model potentially infinitely many unobserved processes.
To demonstrate this concretely, we ran SCOTCH as well as our FCI-inspired algorithm on $100$ SDEs with the adjacency from our example graph $\mathcal{G}$ and compare how often $A$ and $D$ are (falsely) adjacent in the output. SCOTCH predicted an edge between $A$ and $D$ in 88% and ours algorithm in only 8% of all settings, see also Figure 4.

## Peer Reviews

confidence: 4
contribution: 4
initial_rating: 8
presentation: 4
questions: - The  Oracle used for this algorithm, how likely it is to be wrong? How much it would affect the performance?
- What happens if the underlying causal graph has a cycle? How does the algorithm handle this situation?
- How does the algorithm handles the case where the observation in time is limited (subsampled) ?
soundness: 4
strengths: - Provides a solid framework for causal inference in continuous-time SDEs, going beyond the limitations of traditional discrete-time models.
- Introduces a practical CI test using signature kernels, suited for handling path-dependent random variables.
- Shows strong performance with incomplete data and irregularly sampled time series.
- Empirical tests and real-world examples, like pairs trading, demonstrate the algorithm’s practical effectiveness and accuracy.
- Method is not heavily dependent on hyper parameters.
summary: This paper introduces an innovative approach to uncovering causal relationships in stochastic processes, using conditional independence (CI) tests based on signature kernels to detect causal links within stochastic differential equation (SDE) models. It presents a comprehensive algorithm to reconstruct causal graphs, even with partially observed data and irregular sampling patterns. Benchmark tests show this method outperforms existing causal discovery techniques (in small graphs) in continuous-time settings, showing particular strength in cases with incomplete data and path-dependence (without need for hyper parameter tuning)
weaknesses: - The paper assumes stationarity and acyclicity, which may restrict its use in scenarios where causal relationships change over time.
- Limited discussion on performance of the algorithm if the IC Oracle is wrong.
- While I appreciate the rigor in the paper, the length of the paper, considering the appendix is more than $2 /3$ of the page limit. Might be more suitable for a journal setting.

---

confidence: 3
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. Could you briefly explain how to incorporate or partially incorporate PC and FCI into the proposed algorithm, given that both assume IID samples? Is any adjustment needed for non-IID samples?

2. Do KCIPT and SDCIT require IID samples as well? If so, a similar question arises as in item 1.

3. Are the simulated datasets used in the experiments 'continuous-time'? If so, how do you choose the discretization interval for PCMCI and other baselines that assume 'discrete-time'? Does using different discretization intervals influence their performance? How do you compare the output of PCMCI and the proposed algorithm, given that PCMCI includes time lags and may cause 'cycles' according to the definitions in this context?

4. Are you considering using more baselines designed for time series, particularly non-stationary time series? The number of baselines included is limited. Is there a specific reason for not using metrics such as F1 score, precision, and recall but just SHD?

5. Is there a power analysis for the conditional test, and how does it perform with different conditioning sets?

6. Could you explain how to obtain bi-directed edges as shown in Figure 4?
soundness: 2
strengths: 1. The effect working on the 'continuous-time' path sequence is impressive.

2. The signature kernel conditional independence test is a novel CI test for path sequences.

3. The paper is well written, and the framework of the paper is straightforward.

4. The proposed algorithm has been applied to a series of simulations and one real-world dataset.
summary: The authors propose a sound and complete causal discovery algorithm for stochastic processes, incorporating a consistent signature kernel conditional independence test. The stochastic process is modeled using a stochastic differential equation (SDE), and the entire process is segmented into intervals that exhibit the Markov property with respect to an acyclic dependence graph.
weaknesses: 1. There seems to be no assumption sections. The acyclic assumption seems like a common assumption used in many causal discovery methods; however, it is more restricted in this paper as a cycle could be easily created if the causal relations exist for both $X^i_{0,s}$ to $X^j_{s,s+h}$ and $X^j_{s+h,s+2h}$ to $X^i_{s+2h,s+3h}$. Such pairs of causal relations are allowed in many causal discovery methods for time series, such as PCMCI. Therefore, though the authors claim that the proposed algorithm will not rely on the 'discrete-time' assumption, they did not discuss the impact of not assuming a time lag based on the 'discrete-time' assumption and the additional limitations from the acyclic assumption in this paper, compared to many previous causal discovery algorithms. Hence, the limitation discussion and comparison is not comprehensive. The acyclic assumption in this context is not very practical.

2. Does the proposed CI test only work in the specific setting of this paper? Could it be utilized outside of this setting, for more general time series?

3. In the limitations and requirements section, it is stated that the proposed algorithm far exceeds other existing causal discovery methods for time series data, which may not be the best way to frame the discussion. For instance, in part (b), it is mentioned that the proposed algorithm can handle confounders; however, there are already many algorithms that allow for confounders in time series, such as LPCMCI [1] and tsFCI [2]. Additionally, there are algorithms for non-stationary time series, such as CD-NOD [3], and those with special periodic patterns, such as PCMCI$_{\Omega}$ [4]. Therefore, it is difficult to conclude that the proposed algorithm far exceeds other related work, given different settings and assumptions.

4. The number of baselines in the experiment results is limited for both causal discovery and the CI test. Please refer to the questions section for more details.

5. Please correct me if I am mistaken, but it seems there is no computational complexity analysis or running time results provided.

6. I may have misunderstood, but does the first selected interval of $[0, s], [s, h]$ have to start at the beginning? Based on line 233, there are two copies of vertex $V$. Does this imply $h = T$? If not, how many intervals are possible? If multiple intervals are allowed, can the acyclic assumption be relaxed since different time-ordered intervals resemble the concept of "time lag" in a 'discrete-time' setting?

7. By intuition, is having $V_0$ and $V_1$ essentially a sub-sampling technique, where samples in $[0,s]$ and samples in $[s,h]$ are considered? The full causal graph that the algorithm aims to discover is restricted to this $V_0$ and $V_1$, and the estimated causal graph will be influenced by the intervals selected. Again, the full causal graph discussed here is different from the one in many causal discovery algorithms, referred to as the full time causal graph, which does not require sub-sampling and is more comprehensive. Therefore, this discussion needs to be handled with care.

[1]Gerhardus, Andreas, and Jakob Runge. "High-recall causal discovery for autocorrelated time series with latent confounders." Advances in Neural Information Processing Systems 33 (2020): 12615-12625.

[2]Entner, Doris, and Patrik O. Hoyer. "On causal discovery from time series data using FCI." Probabilistic graphical models 16 (2010).

[3]Huang, Biwei, et al. "Causal discovery from heterogeneous/nonstationary data." Journal of Machine Learning Research 21.89 (2020): 1-53.

[4]Gao, Shanyun, et al. "Causal discovery in semi-stationary time series." Advances in Neural Information Processing Systems 36 (2024).

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: Some minor typos and unclear sentences:

- Line 421: “even in the settings it was tailored to” -> what was tailored to these settings, SigKer or the state of the art?
- For figure 2, maybe explicitly say which graph is the lifted graph (right) and which is G (left)
- Line 157: “is inapplicable” -> “it is inapplicable”?
- Some pedantic styling comments:
- Line 105: you have double parenthesis with the citation to Laumaann et al.
- Lines 170,200,420:  you appear to be using hyphens (-) instead of em dashes (---) for interjections. See the JMLR formatting guide (under dashes): https://www.jmlr.org/format/format.html
soundness: 4
strengths: The paper is well-written and a pleasure to read. All theoretical claims are backed by careful proofs, and care is taken to ensure reproducibility of their experimental results. The authors openly address the limitations of their methodology,  and additional interesting results and discussions are provided in a well-organized and comprehensive appendix.

I believe the paper earns its place in the suite of tools for causal inference with time-series data from dynamical systems. In my opinion, this is a high-quality paper that deserves acceptance.
summary: The paper makes contributions to causal discovery in stochastic dynamical systems. In particular, it introduces a novel framework of conditional independence constraints in SDE processes, which are then leveraged to provide a causal discovery (CD) algorithm. The authors also provide a conditional independence (CI) test to evaluate these constraints from data. The authors evaluate the CD algorithm on synthetic data and compare it to other state-of-the-art baselines. They also evaluate their CI test on synthetic data and in a small case study on stock trading pairs.
weaknesses: The only major weakness is the lack of real-world experiments for the causal discovery algorithm, which I understand to be the main contribution of the paper. Only the CI test has a real-world data experiment on a downstream task, which I found creative and is well-documented in the appendix.

Naturally, because we are in causal inference, real-world data with a ground truth can be difficult to find. However, I would like to point the authors to two recent papers that provide real-world data with a causal ground truth and whose settings appear to be a good fit for the method in this paper:

**“[Causal discovery in a complex industrial system: A time series benchmark](https://arxiv.org/abs/2310.18654)” by Mogensen et al. (2023).**

The paper presents a real-world dataset from a dynamical system with partially observed data. The authors provide a ground-truth causal graph (section 2.4). The paper comes with a website (https://soerenwengel.github.io/essdata) with links to the dataset and preprocessing code.

**“[The Causal Chambers: Real Physical Systems as a Testbed for AI Methodology](https://arxiv.org/abs/2404.11341)” by Gamella et al. (2024)**

The authors build two physical devices, one of which (wind tunnel) produces real-world, time-series data from a dynamical system. There is a causal ground-truth graph for this system (Figure 3), which the authors use to benchmark the PCMCI+ algorithm (Figure 6a). I believe this is an extension of one of your baselines (PCMCI), which makes me suspect your method is also applicable. The authors provide a well-documented [notebook ](https://github.com/juangamella/causal-chamber-paper/blob/main/case_studies/causal_discovery_time.ipynb) to download the dataset and reproduce the PCMCI+ experiment. Using your method may be plug-and-play in this case.

There may be other suitable real-world datasets, but I found none after this search. A real-world experiment for the main contribution of the paper (the causal discovery algorithm) would further elevate the value of the paper, and I would be happy to raise my score as a result. This is only a suggestion, and my decision to accept is independent of whether the authors do this or not.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: The experimental setup involves setting small parameters for self-loops (e.g.,  $a_{ii} \sim  U([-0.5, 0.5])$), while the parameters influencing causal effects between different variables are set higher (e.g., $a_{ij} \sim U([-2,-1] \cup [1, 2])$). While I understand this may help to emphasize inter-variable causal relationships, practical applications may not always exhibit this distinction so clearly. Additionally, could the authors clarify whether errors related to self-loops are reported in the Structural Hamming Distance (SHD) results?
soundness: 3
strengths: 1. The paper is well-motivated, presenting a clear need for improved causal discovery methods within the context of stochastic processes.
2. Extensive experimental evaluations support the effectiveness of the proposed approach. The method demonstrates consistent superiority over several baseline models, underscoring its potential contributions to the field.
summary: This paper introduces a novel approach for causal discovery in stochastic processes modeled by stochastic differential equations (SDEs). The authors propose a conditional independence (CI) test based on signature kernels, which enables the identification of causal relationships in continuous-time dynamical systems, even under irregular or partially observed data. Key contributions include a causal discovery algorithm that leverages the directionality of time, as well as an efficient signature kernel-based CI test.
weaknesses: 1.	The paper’s claims regarding diffusion-dependence cases may be somewhat overstated. In particular, for SDE models involving "driving noise" (i.e., cases where the diffusion coefficient depends on $X_t$), certain **causal graphs may not be identifiable from observational data**. See example 5.5 in [1]. In such cases, identifying the generator of the SDE model—so as to identify the post-interventional distribution—may be a more reasonable goal, as explored in [2] for linear SDEs. This may also provide some insight into the **diffusion-dependence** results presented in Table 1, where the proposed method shows comparatively lower performance.

[1] Hansen, Niels, and Alexander Sokol. "Causal interpretation of stochastic differential equations." (2014): 1-24.

[2] Wang, Yuanyuan, et al. "Generator identification for linear SDEs with additive and multiplicative noise." Advances in Neural Information Processing Systems 36 (2024).

2.	The section on signature kernels in Section 2 is mathematically dense and may be challenging for readers unfamiliar with the topic. Including more intuitive explanations would enhance accessibility.
3.	Given the extensive use of mathematical notations, a summary table listing the key notations could be helpful for readers in following the manuscript’s developments.

Meta-review: This paper tackles the challenging and important problem in causal discovery of extracting causal relationships from temporal data governed by a system of stochastic differential equations (SDEs). The authors show how to extract a causal graph from a prototypical system, and propose methods for learning this graph. Reviewers were unanimously in favour of acceptance, with 3/4 reviewers scoring the paper >=8. I recommend the paper be accepted as an oral presentation given the novelty and significance of its contributions.

Decision: Accept (Spotlight)