# Visualization for Dichotomous Variables,
the Independence and Markov chains

ZHANG Yan
<EMAIL_ADDRESS>
##  1 Introduction

In probability theory, the independence is a very fundamental concept, but
with a little mystery. People can always easily manipulate it logistically but
not geometrically, especially when it comes to the independence relationships
among more that two variables, which may also involve conditional
independence. Here I am particularly interested in visualizing Markov chains
which have the well known memoryless property. I am not talking about drawing
the transition graph, instead, I will draw all events of the Markov process in
a single plot. Here, to simplify the question, this work will only consider
dichotomous variables, but all the methods actually can be generalized to
arbitrary set of discrete variables. Let me end the introduction with one
final work below:

This is a non-time-homogeneous Markov chain with 10 dichotomous variables, and
the conditional probabilities are randomly generated. We can actually observe
some sort of fractal pattern in this graph, which is just the reflection of
the Markov property.

##  2 Visualization for Dichotomous Variables

Assume probability space $(\Omega,\mathcal{F},P)$, a dichotomous variable is
just a function $X$: $\Omega\mapsto\\{0,1\\}$, which means a dichotomous
variable can split the sample space $\Omega$ into two parts, $X^{-1}(0)$ and
$X^{-1}(1)$. So, if I assume the points in a plane geometry like a circle, a
triangle or a rectangle as $\Omega$, then I can use anther plane geometry or
simply a line to represent a dichotomous variable if it split the $\Omega$
into two parts. Like what we can see below:

And actually, we don’t set the constraint that the two parts should be two
intact subspaces, which means the follow dichotomous variables are also valid,
but in these cases, we need to set a rule to specify which subspaces are parts
of $X^{-1}(0)$ and which are parts of $X^{-1}(1)$.

For the above example, we can assume the outer ones are parts of $X^{-1}(0)$,
the inner one is $X^{-1}(1)$. I will avoid this kind of phenomenon if
possible, but in the visualization of long Markov chain, it’s not avoidable.

##  3 Visualization for the Independence

With the assistance of squares, circles and lines, we can already visualize
all kinds of ternary independence-dependence relationships quite elegantly.

Firstly, I visualize the common effect relationship:

In the above graph, we can see that $A\nVbar C$, $B\nVbar C$, $A\Vbar B$ but
$A\nVbar B|C$, which is the selection bias phenomenon in causal inference.

And then, the common cause relationship:

This time, $A\nVbar B$, $A\nVbar C$, $B\Vbar C|A$ but $B\nVbar C$, which is
the confounding phenomenon in causal inference. And furthermore, this graph
can also represent a Markov chain of length three: $B\to A\to C$, because
$B\nVbar A$, $A\nVbar C$, $B\nVbar C$ but $B\Vbar C|A$. So, for this three
variables, we actually can’t tell the difference between common cause
relationship and causal chain relationship. And by symmetry, the above graph
can also be $C\to A\to B$, which is a well-known phenomenon (Markov
equivalence graphs [1]) that the following causal diagrams are equivalent to
each other in a sense.

##  4 Visualization for Markov chains

To draw the Markov chains, we need to adopt some more complicated rules.
Firstly, for simplicity, I set our aim to be generating a Markov chain
$X_{0}\to X_{1}\to\dots\to X_{n}$ with symmetric time-homogeneous conditional
probabilities: $P(X_{i+1}=0|X_{i}=0)=1/3$, $P(X_{i+1}=0|X_{i}=1)=2/3$, and the
initial probability $P(X_{0}=1)=1/2$. Let me show you the first three graphs:

(a) $X_{0}$

(b) $X_{0}\to X_{1}$

(c) $X_{0}\to X_{1}\to X_{2}$

In (a), one single vertical line was used to represent $X_{0}$, splitting the
whole space into two equal parts, where the left one is $X_{0}^{-1}(0)$ and
the right one is $X_{0}^{-1}(1)$.

In (b), based on $P(X_{1}=0|X_{0}=0)=1/3$, I used a horizontal line to split
subspace $X_{0}^{-1}(0)$ into two parts, the lower one $X_{1}^{-1}(0)\cap
X_{0}^{-1}(0)$ and the upper one $X_{1}^{-1}(1)\cap X_{0}^{-1}(0)$. Similarly,
by $P(X_{1}=0|X_{0}=1)=2/3$, two more parts, $X_{1}^{-1}(0)\cap X_{0}^{-1}(1)$
and $X_{1}^{-1}(1)\cap X_{0}^{-1}(1)$ were obtained. And $X_{1}^{-1}(0)$ is
just the union of $X_{1}^{-1}(0)\cap X_{0}^{-1}(0)$ and $X_{1}^{-1}(0)\cap
X_{0}^{-1}(1)$, and $X_{1}^{-1}(1)=(X_{1}^{-1}(1)\cap
X_{0}^{-1}(0))\cup(X_{1}^{-1}(1)\cap X_{0}^{-1}(1))$. We can see in this plot,
$X_{0}\nVbar X_{1}$.

In (c), by $P(X_{2}=0|X_{1}=0)=1/3$ and $P(X_{2}=0|X_{1}=1)=2/3$, there are
four subspaces for $X_{2}^{-1}(0)$ and $X_{2}^{-1}(1)$ formed respectively by
four new vertical lines. Easy to see that $X_{1}\nVbar X_{2}$, and because
$P(X_{2}=0|X_{0}=0)=5/9$ which is not equal to $P(X_{2}=0|X_{0}=1)=4/9$,
$X_{0}\nVbar X_{2}$ can also be deduced. But this time, we actually have
$X_{0}\Vbar X_{2}|X_{1}$, which makes the graph a Markov chain.

In next page, I drew the plots from 2-chain to 10-chain, each color
representing one single dichotomous variable.

Finally, if I change the conditional probabilities and the initial
probability, I can actually draw all kinds of Markov chains made by
dichotomous variables as I like, such as the graph I showed in the
introduction. And my current method is not the only way to visualize such
complicated independence relationships, one can easily find another one if he
understand the fundamental, and sometime it’s just a matter of aesthetics.

## References

  * [1] Christopher Meek.  Causal inference and causal explanation with background knowledge.  arXiv preprint arXiv:1302.4972, 2013.


# Castor: Causal Temporal Regime Structure Learning

Antiquus S. Hippocampus, Natalia Cerebro & Amelie P. Amygdale
Department of Computer Science
Cranberry-Lemon University
Pittsburgh, PA 15213, USA
<EMAIL_ADDRESS>
&Ji Q. Ren & Yevgeny LeNet
Department of Computational Neuroscience
University of the Witwatersrand
Joburg, South Africa
<EMAIL_ADDRESS>
Coauthor
Affiliation
Address
email Use footnote for providing further information about author (webpage,
alternative address)— _not_ for acknowledging funding agencies. Funding
acknowledgements go at the end of the paper.

###### Abstract

The task of uncovering causal relationships among multivariate time series
data stands as an essential and challenging objective that cuts across a broad
array of disciplines ranging from climate science to healthcare. Such data
entails linear or non-linear relationships, and usually follow multiple a
priori unknown regimes. Existing causal discovery methods can infer summary
causal graphs from heterogeneous data with known regimes, but they fall short
in comprehensively learning both regimes and the corresponding causal graph.
In this paper, we introduce CASTOR, a novel framework designed to learn causal
relationships in heterogeneous time series data composed of various regimes,
each governed by a distinct causal graph. Through the maximization of a score
function via the EM algorithm, CASTOR infers the number of regimes and learns
linear or non-linear causal relationships in each regime. We demonstrate the
robust convergence properties of CASTOR, specifically highlighting its
proficiency in accurately identifying unique regimes. Empirical evidence,
garnered from exhaustive synthetic experiments and two real-world benchmarks,
confirm CASTOR’s superior performance in causal discovery compared to baseline
methods. By learning a full temporal causal graph for each regime, CASTOR
establishes itself as a distinctly interpretable method for causal discovery
in heterogeneous time series.

##  1 Introduction

††Contact: {abdellah.rahmani<EMAIL_ADDRESS>
Multivariate Time Series (MTS) is a very common type of data in a wide variety
of fields. Uncovering the causal relationships among MTS variables and
understanding how they evolve over time is crucial in numerous fields, such as
climate science and health care. Although randomized controlled trials are
widely recognized as the definitive method for determining causal
relationships (Hariton & Locascio, 2018; McCoy, 2017), they often present
challenges in terms of cost, ethics, or feasibility. For example: learning
gene regulatory networks via gene knockout experiments would be prohibitively
expensive on a large scale. Consequently, a multitude of causal discovery
approaches now focus on extracting causality from observational data sources
(Löwe et al., 2022; Bussmann et al., 2021; Pamfil et al., 2020; Moraffah et
al., 2021; Runge, 2018; Wu et al., 2020).
Yet, the task of temporal causal discovery can be challenging. Firstly, the
causal relationships between variables in the real-world context can manifest
in either linear or non-linear forms. Secondly, certain causal effects may be
instantaneous, making the differentiation between cause and effect
particularly challenging. Lastly, the presence of history-dependent noise,
where noise levels depend on past observations, adds yet another layer of
complexity. Various research works, such as those introducing models like
Dynotears (Pamfil et al., 2020) and Rhino (Gong et al., 2022), have aimed to
address these challenges. However, a common assumption in most existing
methods is that time series observations are homogeneous under a singular
”domain/regime” described by one causal graph. This assumption often falls
short in real-world scenarios, like Electroencephalographic (EEG) time series
in epilepsy setting (Tang et al., 2021; Rahmani et al., 2023), or a genomic
cancer data (Saeed et al., 2020), for example. In such scenarios, the samples
originate from a mixture of different causal models on the same set of
variables, and the causal effects of the mixture distribution can in general
not be faithfully represented by a single Directed Acyclic Graph (DAG). The
authors in Thiesson et al. (2013), address causal discovery for heterogeneous
data in the context of independent identically distributed (iid) variables.
They develop a method grounded in the Expectation-Maximization (EM) algorithm
(Dempster et al., 1977) and designed to unearth and learn from a mixture of
DAGs. The above work however learn one causal graph per regime and does not
incorporate the notion of time during causal graph learning, potentially
leading to missed or inaccurately predicted causal relations. Recent work, CD-
NOD, developed by Huang et al. (2020), tackled the issue of heterogeneity
differently, and addresses heterogeneous time series by modulating causal
relationships through a regime index. While this provides a summary graph
highlighting behavioral changes across regimes, it cannot infer individual
causal graphs and necessitates prior knowledge about the number of regimes.
Overall, most of the existing methods make certain assumptions that may limit
their applicability, such as a prior knowledge of the number of regimes, prior
information about regime indices, which are often not readily available.
We introduce CASTOR, a framework for regime-dependent causal structure
learning. CASTOR is designed to learn both linear and nonlinear temporal
causal relationships from heterogeneous multivariate time series (MTS),
without necessitating prior knowledge of regime indexes or the total number of
regimes involved. Using a predetermined window for index initialization,
CASTOR maximizes a score function via EM algorithm to infer the number of
regimes and learns linear or non-linear causal relationships inherent in each
regime. To the best of our awareness, our framework represents an
unprecedented initiative in the time series landscape for the simultaneous
learning of causal graphs, regime numbers, and regime indices. Our key
contributions include:

  * •

In this paper, we introduce CASTOR, a score based method for temporal causal
structure learning tailored for heterogeneous time series.

  * •

CASTOR employs the EM algorithm in conjunction with the Notears (Zheng et al.,
2020) penalization method to deduce regime partitions and infer causal graphs.

  * •

Building upon the intuition of LIN (Liu & Kuang, 2023), we demonstrate that
regime indices can indeed be recovered, subject to Kullback–Leibler divergence
constraints. Furthermore, we establish that within this framework, the
underlying causal structure is identifiable, thus confirming the robustness
and efficacy of our approach.

  * •

We conduct extensive synthetic experiments with ablation studies to
demonstrate the advantages of CASTOR and its robustness across different
settings. Additionally, we demonstrate its performance compared to baselines
in real-world discovery benchmarks.

##  2 Framework

Consider the multrivariate time series (MTS)
$\bm{(X_{t})}_{t\in\mathcal{T}}=(X_{t}^{i})_{i\in\mathbf{V},t\in\mathcal{T}}$
where $|\mathbf{V}|=d$ number of components of MTS
$\bm{(X_{t})}_{t\in\mathcal{T}}$ and $\mathcal{T}$ is time index set. The MTS
$\bm{(X_{t})}_{t\in\mathcal{T}}$ is aligned with a temporal causal graph
elucidating its generative process. Each component of the MTS at time t is
represented by a singular node within the graph $\mathcal{G}$, leading us to
define $\mathcal{G}$ as follows:

###### Definition 1

(Temporal Causal Graph, Assaad et al. (2022)) $\mathcal{G}=(V,E)$ is the
temporal causal graph associated with the MTS $\bm{(X_{t})}_{t\in\mathcal{T}}$
is a DAG and the set of vertices in that graph consists of the set of
components $X^{1},\ldots,X^{d}$ at each time $t\in\mathbb{N}$. The edges $E$
of the graph are defined as follows: variables $X_{t-\tau}^{i}$ and
$X_{t}^{j}$ are connected by a lag-specific directed link
$X_{t-\tau}^{i}\rightarrow X_{t}^{j}$ in $\mathcal{G}$ pointing forward in
time if and only if $X^{i}$ at time $t-\tau$ causes $X^{j}$ at time $t$ with a
time lag of $\tau>0$ for $i=j$ and with a time lag of $\tau\geq 0$ for $i\neq
j$.

In the rest of our paper, we refer to the parents of node $i$ in $\mathcal{G}$
at previous time (lagged parents) by $\mathbf{Pa}_{\mathcal{G}}^{i}(<t)$ and
to the parents at the current time $t$ (instantaneous parents) by
$\mathbf{Pa}_{\mathcal{G}}^{i}(t)$. To clarify our notation, if $X^{i}$ at
time $t-\tau$ causes $X^{j}$ at time $t$ with a time lag of $\tau>0$, we have
an edge $X_{t-\tau}^{i}\rightarrow X_{t}^{j}$ and
$X^{i}\in\mathbf{Pa}_{\mathcal{G}}^{j}(<t)$ and if $\tau=0,\text{ we have an
edge }X_{t}^{i}\rightarrow X_{t}^{j}\text{ and
}X^{i}\in\mathbf{Pa}_{\mathcal{G}}^{j}(t)$.
Structural Equation Models (SEMs). Schölkopf et al. (2021) For a MTS
$\bm{(X_{t})}_{t\in\mathcal{T}}$, the associated SEM for a given temporal
causal graph, denoted as $\mathcal{G}$ for additive noise models (ANMs):

$X_{t}^{i}=g_{i}\left(\mathbf{Pa}_{\mathcal{G}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}}^{i}(t)\right)+\epsilon_{t}^{i},$
(1)

where $g_{i}$ is a deterministic function and $\epsilon$ are mutually
independent noise variables. SEM describes the causal relationships between
MTS components given a causal graph $\mathcal{G}$. For the same MTS and by
assuming causal Markov property (Definition 4 in the Appendix B), we can
further define a

Figure 1: Overview of CASTOR. The dashed edges refer to the time lagged links,
the normal arrows represent the instantaneous links.

causal graphical model (CGM) by a distribution $p$ as follows:

$p\left(\bm{X}_{t}\mid\mathcal{G}\right)=\prod_{i=1}^{d}p\left(X_{t}^{i}\mid\mathbf{Pa}_{\mathcal{G}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}}^{i}(t)\right),$
(2)

with $p\left(\bm{X}_{t}\mid\mathcal{G}\right)$, called observational
distribution, capturing the joint probability between $X^{i}$ and
$p\left(X_{t}^{i}\mid\mathbf{Pa}_{\mathcal{G}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}}^{i}(t)\right)$,
refering to the conditional distribution of $X_{t}^{i}$ given parents of
$X_{t}^{i}$ in the temporal causal graph $\mathcal{G}$. A MTS may have one or
multiple regimes. In our context, we characterize an MTS with multiple regimes
as follows:

###### Definition 2

(MTS with multiple regimes) We say that MTS
$(\bm{X}_{t})_{t\in\mathcal{T}}=(X_{t}^{i})_{i\in\mathbf{V},t\leq\mathcal{T}}$
is composed of $K$ disjoint regimes if it exists a unique time partition
$\mathcal{E}=(\mathcal{E}_{u})_{u\in\\{1,...,K\\}},\text{ such that:
}\cap_{u=1}^{K}\mathcal{E}_{u}=\emptyset$,
$\mathcal{T}=\cup_{u=1}^{K}\mathcal{E}_{u}$ and $\forall u\in\\{1,...,K\\},$
the MTS $(\bm{X}_{t})_{t\in\mathcal{E}_{u}}$ is stationary and associated with
a unique temporal causal graph $\mathcal{G}^{u}$.

It is crucial to note that $\mathcal{G}^{u}$ are regime-dependent, meaning
that graphs vary across different regimes i.e,
$\mathcal{G}^{u}\neq\mathcal{G}^{v}$. In the rest of our paper, we note
$\mathcal{G}=\\{\mathcal{G}^{u}\\}_{u\in\\{1,...,K\\}}$ the set of temporal
causal graph associated with a MTS $(\bm{X}_{t})_{t\in\mathcal{T}}$ composed
of $K$ multiple regimes. Within each regime $u$, the temporal causal graph
$\mathcal{G}^{u}$, can be represented by a collection of adjacency matrices,
collectively denoted as $\bm{G}^{u}_{\\{0:L\\}}$ with maximum lag $L$. Here,
$\bm{G}^{u}_{\tau\in\\{1:L\\}}$ refers to the adjacency matrix of the temporal
causal graph $\mathcal{G}^{u}$ at a specific time lag $\tau\in\\{1:L\\}$. In
other words, if $X_{t-\tau}^{i}\rightarrow X_{t}^{j}$, then the coefficient
$[G^{u}_{\tau}]_{ij}\neq 0$, and 0 otherwise. When $\tau=0$, we have
instantaneous links within the temporal causal graph $\mathcal{G}^{u}$.
We now propose a novel functional form for SEM that incorporates linear/non-
linear relations, instantaneous links and multiple regimes:

$\forall u\in\\{1,...,K\\},\forall
t\in\mathcal{E}_{u}:X_{t}^{i}=g^{u}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right)+\epsilon_{t}^{i},$
(3)

where $g^{u}_{i}$ is a general differentiable linear/non-linear function and
$\epsilon_{t}^{i}\sim\mathcal{N}(0,1)$, follows to a normal distribution. By
assuming causal markov property, we can define the associated CGM:

$\forall u\in\\{1,...,K\\},\forall
t\in\mathcal{E}_{u}:p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right)=\prod_{i=1}^{d}p^{(u)}\left(X_{t}^{i}\mid\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right),$
(4)

As we explained above, the temporal causal dynamics of an MTS with $K$
multiple regimes typically cannot be captured by a single DAG; instead, a
mixture of DAGs is more appropriate. We define this mixture distribution as
follows:

$p\left(\bm{X}_{t}\right)=\sum_{u=1}^{K}\pi_{t,u}(\alpha)p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right),$
(5)

where $p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right)$ is defined in Eq (4)
and $\pi_{t,u}(\alpha)$ represent mixing coefficients. These coefficients can
be learned freely; however, in our work, we define $\pi_{t,u}(\alpha)$ based
on the work by Samé et al. (2011) as the following:
$\pi_{t,u}(\alpha)=\frac{\exp\left(\bm{\alpha}_{u,1}t+\bm{\alpha}_{u,0}\right)}{\sum_{u=1}^{K}\exp\left(\bm{\alpha}_{u,1}t+\bm{\alpha}_{u,0}\right)}$.
The soft-max function defined in this way ensures a smooth transition between
the different regimes. Previous work assumes prior knowledge of time partition
$\mathcal{E}$ or report a summary causal graph, falling short of elucidating
the full temporal causal graph. In the next section we will present CASTOR a
score based method for temporal causal structure learning tailored for
heterogeneous time series.

##  3 CASTOR: Causal temporal regime structure learning

In this section, we present a novel approach designed to learn the set of
temporal causal graphs $\mathcal{G}=\\{\mathcal{G}^{u}\\}_{u\in\\{1,...,K\\}}$
associated with MTS $(\bm{X}_{t})_{t\in\mathcal{T}}$ composed of multiple
regimes, the total number of regimes $K$, and the associated time partition
$\mathcal{E}=(\mathcal{E}_{u})_{u\in\\{1,...,K\\}}$. Sections 3.1 and 3.2
present solutions tailored to both linear and nonlinear contexts,
respectively.

###  3.1 Castor for linear causal relationships

We first consider the case where the causal relationships between the
components are linear. Based on Eq (3), the SEM can then be articulated as
follows: $\forall u\in\\{1,...,K\\},\forall
t\in\mathcal{E}_{u}:\bm{X}_{t}=\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}+\bm{\epsilon}_{t}.$
From the causal Markov assumption (4) in the Appendix B and using Eq (5) we
have :

$\log
p\left(\bm{X}_{0:|\mathcal{T}|}\right)=\sum_{t=0}^{|\mathcal{T}|}\log\sum_{u=1}^{K}\pi_{t,u}(\alpha)p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right),$
(6)

Our objective is to learn simultaneously $K$,
$\mathcal{E}=(\mathcal{E}_{u})_{u\in\\{1,...,K\\}}$ and
$\mathcal{G}^{u}=\bm{G}^{u}_{\\{0:L\\}},\forall u\in\\{1,...,K\\}$ that
maximise the $\log p\left(\bm{X}_{0:|\mathcal{T}|}\right)$ defined in Eq (6).
First learning the temporal causal graphs, represented as $\forall
u\in{1,...,K}:\mathcal{G}^{u}=\bm{G}^{u}_{\\{0:L\\}}$, concurrently entails
the estimation of the regime distribution
$p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right)$. In line with the
definition of CGM in Equation (2), we model CASTOR’s joint density of the uth
regime by:

$f^{u}\left(\bm{X}_{t}\right):=\prod_{i=1}^{d}f^{u}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right),$
(7)

where $f^{u}$ is a distribution family. It is important to highlight that
while $f^{u}$ can, in theory, be any distribution, in this particular study,
we focus on Gaussian noise, used by many works and for which they showed the
identifiability of causal graphs for one regime (Huang et al., 2020; Peters &
Bühlmann, 2014). As a result, our distribution $f^{u}$ will be a normal
distribution. Estimating the aforementioned parameters ($K,\mathcal{E}\text{
and }\mathcal{G}^{u},\forall u\in\\{1,...,K\\}$) is challenging as the regime
to which each time-step belongs remains unknown. The presence of the sum
within the log terms in Eq (6) renders the log likelihood intractable. The EM
algorithm (Dempster et al., 1977) is well-suited to address such problems; it
introduces a new variable to model regime participation $\gamma_{t,u}$,
thereby transforming the log likelihood to the following form: $\log
p\left(\bm{X}_{0:|\mathcal{T}|}\right)=\sum_{t=0}^{|\mathcal{T}|}\sum_{u=1}^{K}\gamma_{t,u}\log\pi_{t,u}(\alpha)p^{(u)}\left(\bm{X}_{t}\mid\mathcal{G}^{u}\right).$
To learn the regime indexes
$\mathcal{E}=(\mathcal{E}_{u})_{u\in\\{1,...,K\\}}$ and the number of regimes
$K$, CASTOR initially divides the MTS into $N_{w}>K$ equal time windows in the
initialisation step, where each window represents one regime. CASTOR uses this
initialisation and the formulation (7) to estimate the temporal causal graphs
by maximising the following equation:

$\displaystyle\sup_{\theta,\alpha}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=0}^{|\mathcal{T}|}\gamma_{t,u}\log\pi_{t,u}(\alpha)f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{
s.t., }\bm{G}^{u}_{0}\text{ is a DAG,}$ (8)
$\displaystyle\iff\left\\{\begin{array}[]{l}\max_{\alpha}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\ln\left(\pi_{t,u}(\alpha)\right)\\\
\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{ s.t.,
}\bm{G}^{u}_{0}\text{ is a DAG,}\end{array}\right.$

Figure 2: CASTOR’s regime learning procedure

where $\theta$ stands for $\forall u\in\\{1,...,K\\}:\bm{G}^{u}_{\\{0:L\\}}$
and $|\mathcal{G}^{u}|$ is the number of edges in the temporal causal graph of
regime $u$. The first term is the averaged log-likelihood over data, the
second is a penalty term with positive small coefficient $\lambda$ that
controls the sparsity constraint. We impose an acyclicity constraint solely on
the adjacency matrix $\bm{G}^{u}_{0}$ of instantaneous links, as the other
adjacency matrices $\bm{G}^{u}_{\\{1:L\\}}$ are inherently acyclic by
definition, because these matrices establish links between variables at time
$t$ and their time-lagged parents at time $t-\tau$. From the Eq (8), We can
note $\mathcal{S}$, the score function of CASTOR as follows:

$\mathcal{S}(\mathcal{G},\mathcal{E}):=\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{ s.t
}\bm{G}^{u}_{0}\text{ is a DAG,}$ (9)

To address the optimization challenge that incorporates the acyclicity
constraints, we employ an augmented Lagrangian method Zheng et al. (2018);
Pamfil et al. (2020); Brouillard et al. (2020); Liu & Kuang (2023). Our Eq (9)
can be succinctly written as:

$\min_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\left\|\bm{X}_{t}-\left(\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}\right)\right\|_{F}^{2}+\lambda|\mathcal{G}^{u}|+\frac{\rho}{2}h\left(\bm{G}^{u}_{0}\right)^{2}+\alpha
h\left(\bm{G}^{u}_{0}\right),$ (10)

where $\alpha,\rho$ characterize the strength of the DAG penalty. The function
$h(\bm{G})=\operatorname{tr}\left(e^{\bm{G}\odot\bm{G}}\right)-d$ corresponds
to the acyclicity constraints proposed in Zheng et al. (2018) ($\odot$ is the
Hadamard product). After this initial step of learning the graphs with $N_{w}$
equal windows, our method alternates between updating regime indexes
$\mathcal{E}=(\mathcal{E}_{u})_{u\in\\{1,...,K\\}}$, during the expectation
phase (E-step), and inferring the temporal causal graphs
$\bm{G}^{u}_{\\{0:L\\}}$ during the maximisation phase (M-step). In the
expectation step, CASTOR updates the probability $\gamma_{t,u}$, that refers
to the probability of $\bm{X}_{t}$ belonging to regime $u$, using the
following equation:

$\gamma_{t,u}=\frac{\pi_{t,u}(\alpha)f^{u}\left(\bm{X}_{t}\right)}{\sum_{j=1}^{N_{w}}\pi_{t,j}(\alpha)f^{j}\left(\bm{X}_{t}\right)}$
(11)

Following each E-step, CASTOR filters out the regimes characterized by an
insufficiency of samples. In other words, CASTOR eliminates all regimes with
fewer samples than $\zeta$ that is the minimum regime duration, denoted as
$\zeta$ and defined as a hyper-parameter. Subsequently, samples from the
discarded regimes are allocated to the closest regime in terms of probability
(if the discarded regime is $u$, the sample $\bm{X}_{t}$ will be allocated in
the next iteration to regime $v$ with the highest $\gamma_{t,v}$). Figure 2
illustrates the regime learning process.

###  3.2 CASTOR for non-linear causal relation

We now describe our approach to discerning non-linear causal relationships
from MTS, starting with the case $K=1$ for ease of description. Subsequently,
we will outline how this methodology extends to address MTS with multiple
regimes. In this setting, the SEM is the same as the one expressed in Eq (1):
$X_{t}^{i}=g_{i}\left(\mathbf{Pa}_{\mathcal{G}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}}^{i}(t)\right)+\epsilon_{t}^{i},$
where $g_{i}$ represents a differentiable function that carries the non
linearity property for the causal relationships and
$\epsilon_{t}^{i}\sim\mathcal{N}(0,1)$. The associated CGM is defined
similarly to Eq (2).
Our objective is to recover the temporal causal graph
$\mathcal{G}=\bm{G}_{\\{0:L\\}}$, which is equivalent to estimate the true
distribution defined in Eq (2):
$p\left(\bm{X}_{t}\mid\mathcal{G}\right)=\prod_{i=1}^{d}p\left(X_{t}^{i}\mid\mathbf{Pa}_{\mathcal{G}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}}^{i}(t)\right)$.
We employ Neural Networks (NN) to accommodate the non-linearity introduced in
our problem formulation, similarly to what is done in Brouillard et al.
(2020); Liu & Kuang (2023), while utilizing NOTEARS (Zheng et al., 2018)
constraints to enforce acyclicity for the graph modeling instantaneous links.
The NN is used to estimate the parameters of our distribution, defined
as111Our work builds upon the foundational ideas presented in Zheng’s work
(Zheng et al. (2020)) on causal discovery in non-linear iid data.:

$f\left(\bm{X}_{t}\right):=\prod_{i=1}^{d}f_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right),$
(12)

with $f$ a normal distribution family. The parameters of the distribution, are
determined by a neural networks NN based on variables’ parents’ value. For
each time step $t$, we aggregate all time-lagged variables to form a single
time lag vector, denoted as
$\bm{X}_{t}^{\text{lag}}=\left[\bm{X}_{t-1}|\cdots|\bm{X}_{t-L}\right]$, which
encompasses the lagged data. We employ $\bm{X}_{t}$ and
$\bm{X}_{t}^{\text{lag}}$ as inputs for different neural networks $\text{
NN}_{i},\forall i\in\\{1,...,d\\},$ with the objective of estimating the
parameters of the distributions $f_{i}$. Mathematically, our Neural Networks
are formulated as follows: $\forall i\in\\{1,...,d\\}:\text{
NN}_{i}(\bm{X}_{t},\bm{X}_{t}^{\text{lag}})=\psi_{i}\left(\phi_{i}(\bm{X}_{t}),\phi^{\text{lag}}_{i}(\bm{X}_{t}^{\text{lag}}))\right),$
where $\psi_{i}$ are neural networks composed of locally connected layers
introduced by Zheng et al. (2020) and activation functions
$\phi_{i},\phi^{\text{lag}}_{i}$ $(i\in\\{1,...,d\\})$ composed of linear
layers and sigmoid activation functions. The locally connected layers help to
encode the variable dependencies in the first layer. Let $\Theta_{i}$ and
$\Theta_{i}^{\text{lag}}$ represent the parameters of the first layer for
$\phi_{i},\phi^{\text{lag}}_{i}$ respectively. For a given node $i$, the
instantaneous and time-lagged interactions with another node $j$ can be
succinctly captured by examining the norms of the corresponding columns $j$ in
the weight matrices of the initial layers such that:
$[\bm{G}_{0}]_{ij}=\left\|\Theta_{i}(\text{column
j})\right\|_{2},[\bm{G}_{\tau}]_{ij}=\left\|\Theta_{i}^{\text{lag}}(\text{column
j})\right\|_{2}$. As we know the matrix
$\Theta_{i}^{\text{lag}}\in\mathbb{R}^{d\times dL}$, recovering lag matrices
by the above equation requires a reshape formulation, we keep the notation
above for simplicity. The objective is thus to learn $\Theta_{i}$ and
$\Theta_{i}^{\text{lag}}$ that will encode the causal relationships by
maximizing the following score function:

$\mathcal{S}(\mathcal{G}):=\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{t=1}^{|\mathcal{T}|}f(\bm{X}_{t})-\lambda|\mathcal{G}|,\text{
s.t }\bm{G}_{0}\text{ is a DAG,}$ (13)

Given a loss function $\mathcal{L}$ such as least squares, maximising Eq (13)
is equivalent to:

$\displaystyle\min_{\theta,\mathcal{G}}\frac{1}{|\mathcal{T}|}\sum_{t=1}^{|\mathcal{T}|}\sum_{i=1}^{d}\mathcal{L}(X^{i}_{t},\psi_{i}\left(\phi_{i}(\bm{X}_{t}),\phi^{\text{lag}}_{i}(\bm{X}_{t}^{\text{lag}})\right))+\lambda|\mathcal{G}|+\frac{\rho}{2}h\left(\bm{G}_{0}\right)^{2}+\alpha
h\left(\bm{G}_{0}\right)$ (14)

where $\theta$ are the neural networks parameters and $\mathcal{G}$ is a
notation of $\bm{G}_{\\{0:L\\}}$. For the acyclicity constraint, we add the
augmented Lagrangian term $\frac{\rho}{2}h\left(\bm{G}_{0}\right)^{2}+\alpha
h\left(\bm{G}_{0}\right)$ to our optimization problem.
Similarly to the linear case, CASTOR for non linear case uses the EM algorithm
that alternates between distribution parameters learning (these parameters
include the causal graphs) and regime learning. The sole change between the
linear and non-linear scenarios resides in the methodology employed for
estimating the distribution parameters. In the linear instance, we have
demonstrated that the parameters are directly linked to the temporal causal
graphs through the adjacency matrices. Conversely, in the non-linear scenario,
we leverage Neural Networks (NNs) to estimate the distribution parameters.
Hence, the Maximisation step of the EM algorithm, specifically tailored for
the non-linear case, is as follows:

$\displaystyle\min_{\theta,\mathcal{G}}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\sum_{i=1}^{d}\gamma_{t,u}\mathcal{L}(X^{i}_{t},\psi^{u}_{i}\left(\phi^{u}_{i}(\bm{X}_{t}),\phi^{u,\text{lag}}_{i}(\bm{X}_{t}^{\text{lag}})\right))+\lambda|\mathcal{G}^{u}|+\frac{\rho}{2}h\left(\bm{G}^{u}_{0}\right)^{2}+\alpha
h\left(\bm{G}^{u}_{0}\right)$ (15)

where $\theta$ summarises all the network parameters and
$\mathcal{G}=\\{\mathcal{G}^{u}\\}_{u\in\\{0,...,K\\}}$. The formulations of
the E-step and regime adjustment remains the same as in the linear case.
Algorithm (1) outlines our CASTOR model for both linear and nonlinear causal
relationships. It details the process for updating parameters at each
iteration and employs a minimum regime duration $\zeta$ to eliminate
unnecessary regimes, thereby determining the optimal number of regimes.

Algorithm 1 CASTOR algorithm

procedure CASTOR($\bm{X},W,L,N_{\text{iter}},\zeta$)$\triangleright$
$\bm{X}\in\mathbb{R}^{T\times d}$,$W$: window size

while $\text{iter}\leq N_{\text{iter}}$ do

$\gamma_{t,u}\leftarrow\frac{\pi_{t,u}(\alpha)f^{u}\left(\bm{X}_{t}\right)}{\sum_{j=1}^{N_{w}}\pi_{t,j}(\alpha)f^{j}\left(\bm{X}_{t}\right)}$
$\triangleright$ E step 11

$\alpha\leftarrow\operatorname*{arg\,min}_{\alpha}\sum_{u=1}^{K}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\ln\left(\pi_{t,u}(\alpha)\right)$
$\triangleright$ M step, 8

$\mathcal{G}\leftarrow\text{argmin of Eq \ref{1313} or
\ref{16}}$$\triangleright$ M step

if $\sum_{t}^{|\mathcal{T}|}\gamma_{t,u}\leq\zeta$ then

$\forall t:\gamma_{t,u}\leftarrow 0$ $\triangleright$ Cancel the regime $u$

return $\gamma,\mathcal{G}$

##  4 Theoretical guarantees

In this section, we present the theoretical underpinnings that ensure the
robustness of CASTOR as a causal discovery method for time series data
encompassing multiple regimes. Specifically, we provide guarantees on two
fronts: the unambiguous identification of regime indices, and the accurate
inference of causal relationships. Together, these guarantees confirm the
ability of CASTOR to uniquely discern causal relations across time series with
varying regimes. CASTOR learns the temporal causal graphs by optimizing a
score function (Eq (9):
$\mathcal{S}(\mathcal{G},\mathcal{E})=\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{\mathrm{u}}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|$) during the M-step.
However, the score function also depends on a time index partition
$\mathcal{E}$, which CASTOR learns in the expectation step. Consequently, the
convergence of CASTOR and its capability to accurately discern the true regime
partition, as well as construct appropriate temporal causal graphs, are
contingent upon effectively optimizing the score function until it attains its
optimal value. Let $\mathcal{G}^{*}$ be the set of ground truth causal graphs
i.e
$\mathcal{G}^{*}=\\{\mathcal{G}^{*,u}\\}_{u\in\\{1,...,K\\}}=\\{\bm{G}^{*,u}_{\\{0:L\\}}\\}_{u}$
and $\mathcal{E}^{*}$ the correct partition based on regime index. The
following theorem summarizes our main theoretical contribution.

###### Theorem 1

We assume that each regime has enough data and the penalty coefficients in
Equation (15) are sufficiently small, and all the assumption 1, 2,3,4,5,6,7
hold (Appendix B for precise statements)), we have for any estimation
$(\hat{\mathcal{G}},\hat{\mathcal{E}}):\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)>\mathcal{S}(\hat{\mathcal{G}},\hat{\mathcal{E}}),$
if $\exists u\in\\{1,...,K\\}\text{ s.t }\hat{\mathcal{G}^{u}}\text{ disagrees
with }$ $\mathcal{G}^{*,u}$ on instantaneous or/and time lagged link, or any
regime $\hat{\mathcal{E}}_{u}\in\hat{\mathcal{E}}$ is close to none of the
true regimes in the sense of Kullback–Leibler divergence.

Proof.  Details of the proof in the Appendix B.1.
The theorem indicates that inaccuracies in identifying causal structures or
real regimes yield a sub-optimal estimation score. By optimizing according to
Equation (9) we will asymptotically identify the actual regimes and recover
true causal graphs, leading to the subsequent corollary.

###### Corollary 1.1

Given the conditions stated in Theorem 1, if the score proposed in Equation
(9) is optimized, then samples in each regime would approach one of the true
regimes asymptotically in the sense of Kullback-Leibler divergence and the
causal structure of each regime is identifiable.

##  5 Experiments

###  5.1 Synthetic data

$K=2$ $K=3$ $K=4$ $d$ Method Type SHD F1 SHD F1 SHD F1 10 VARLINGAM Inst Lag
26 16 $18.2_{\pm 3.8}$ $10.4_{\pm 7.8}$ 39 28 $11.0_{\pm 6.4}$ $5.01_{\pm
1.4}$ 59 33 $9.70_{\pm 1.8}$ $5.10_{\pm 3.1}$ RPCMCI Inst Lag - 38 -
$42.3_{\pm 11.1}$ - 41 - $18.8_{\pm 2.5}$ - - - - CASTOR Inst Lag 0 0 100±0
100±0 1 0 97.3±2.5 100±0 1 2 99.3±0.9 98.0±1.4 40 VARLINGAM Inst. Lag 134 102
$8.40_{\pm 1.2}$ $1.2_{\pm 0.1}$ 202 137 $9.83_{\pm 0.9}$ $1.13_{\pm 0.8}$ 248
225 $10.9_{\pm 3.1}$ $1.43_{\pm 0.6}$ RPCMCI Inst. Lag - 155 - $42.1_{\pm
3.5}$ - 321 - $18.4_{\pm 14}$ - - - - CASTOR Inst. Lag 2 0 98.3±1.7 100±0.0 9
1 98.2±1.2 99.8±0.2 7 2 98.3±0.4 98.9±0.9

Table 1: F1 and SHD Scores by Models and Settings: $d$ indicates number of
nodes, $K$ refers to the number of regimes, and ’Type’ refers to the causal
links as either instantaneous or time-lagged.

We perform extensive experiments to evaluate the performance of our method,
CASTOR, in both linear and non-linear causal relationships. For groundtruth
graph generation, we employ the Barabási-Albert (Barabási & Albert, 1999)
model with a degree of 4 to establish instantaneous links, while we utilize
the Erdős–Rényi (Newman, 2018) random graph model with degrees ranging from 1
to 2 for time-lagged relationships. We focus on scenarios with a single time
lag $L=1$, although additional experiments involving multiple lags are
detailed in the Appendix. The duration of each regime is chosen randomly from
the set {300, 400, 500, 600}. Each experiment (combination choice of $K$ and
$d$) was repeated three times under multiple settings, all combinations yield
to more than 60 different datasets.
Linear Relationships.  We examine varying numbers of nodes, specifically
$\\{5,10,20,40\\}$, and generated time series with different regime counts
$\\{2,3,4,5\\}$. Our model’s performance is benchmarked against multiple
baselines, namely RPCMCI (Saggioro et al., 2020) and VARLINGAM (Hyvärinen et
al., 2010) and the results are presented in Table 1. RPCMCI represents the
sole baseline tailored to address a similar setting. RPCMCI necessitates
specific parameters, including the number of regimes and the maximum number of
transitions, and with this input, it only infers time-lagged relations. Even
with this detailed information, RPCMCI struggles to achieve convergence,
particularly in settings with more than 3 different regimes. In contrast,
CASTOR does not only surpass RPCMCI in performance but also converges
consistently, correctly identifying the number of regimes and recovering both
the regime indices and the underlying causal graphs of each regime. We can
notice that CASTOR successively infers the regime indexes and learns as well
the instantaneous links (more than 95% F1 score in different settings) as well
as time lagged relations (more than 97% in almost all the settings). When we
compare CASTOR with VARLINGAM (which performs causal discovery method for MTS
data which can model both lagged and instantaneous links), the former also
demonstrates markedly superior performance. To provide context, we manually
partition our generated data into discrete regimes to facilitate the
evaluation of VARLINGAM. We then executed VARLINGAM on each segmented regime,
synthesized the graphs, and compared these composite structures against their
true counterparts. Even when executing VARLINGAM separately on each regime,
CASTOR still surpasses VARLINGAM, all without access to any prior information,
such as the number of regimes or the indices of the regimes. Additional
results that confirm the above results are in the Appendix C.
Non-linear Relationships.  For non-linear relationships, the functions
$g^{u}_{i}$ defined in Eq (3) include random weights generated from a uniform
distribution over the interval ]0,2], coupled with activation functions
selected randomly from the set {Tanh, LeakyReLU, ReLU}. In this case, we
compare our model against various baseline models, namely DYNOTEARS (Pamfil et
al., 2020) and VARLINGAM (Hyvärinen et al., 2010) and conducted experiments
with different numbers of regimes {2, 3, 4} and nodes{10, 20}. We can see from
Figure (3) that both CASTOR and DYNOTEARS exhibit superior performance to
VARLINGAM. It is important to outline that DYNOTEARS and VARLINGAM are each
applied to individual regimes separately; neither is designed to learn or
infer the number or indices of regimes. CASTOR demonstrates comparable
performance to DYNOTEARS in modeling time-lagged relations for non-linear
scenarios while learning also the regime indexes. It also succeeds in
identifying instantaneous links. DYNOTEARS achieves inferior results in
identifying instantaneous links due to his formulation that takes into
consideration only linear relationships. This is much clearer in high
dimensions due to the increasing complexity of the problem. Additional
experiments are available in Appendix C, and show the same trends as explained
above.

Figure 3: F1 by Models and Settings, Inst means instantaneous links and lag
time-lagged ones.

###  5.2 Web activity dataset

Model | F1 Reg1 | F1 Reg 2 | Reg Acc
---|---|---|---
PWGC | 53.8 | 20.0 | -
VARLINGAM | 66.7 | 0 | -
CASTOR | 18.2 | 28.5 | 100
Table 2: F1 Scores across Two IT Datasets: ’Reg’ indicates regime and ’Reg
Acc’ refers to regime identification accuracy.

Model | F1 FMRI | Reg Acc
---|---|---
PWGC | 66.7 | -
VARLINGAM | 47.6 | -
CASTOR | 66.5 | 80
Table 3: F1 Scores for FMRI data, ’Reg’ indicates regime and ’Reg Acc’ refers
to regime identification accuracy.

We now evaluate222Our evaluation methodology aligns with the approach outlined
in Assaad et al. (2022) CASTOR on two stacked IT monitoring time series, each
comprising 1106 timestamps and 7 nodes, sourced from EasyVista
333https://github.com/ckassaad/causal_discovery_for_time_series. These series
capture activity metrics from a web server. We compared our method with PWGC
(Granger Causality (Granger, 1969)) and VARLINGAM (Hyvärinen et al., 2010). As
it is evident from Table 3, CASTOR proficiently identifies the exact number of
regimes and regime indices. It is pertinent to note that PWGC and VARLINGAM
are not inherently designed to infer regime indices, we evaluate both models
separately on each regime. On regime 2, CASTOR outperforms PWGC and VARLINGAM
in learning causal relationships. While VARLINGAM exhibits superior results
compared to CASTOR in regime 1, it is not designed to learn the causal graph
and the indices. Also, given that our data pertains to IT monitoring, the
likelihood of the presence of instantaneous links is relatively low, which
could account for the good performance of VARLINGAM and PWGC that models only
time-lagged links.

###  5.3 Nestim brain connectivity

Finally we test the efficacy of CASTOR on FMRI imaging dataset Nestim (Smith
et al., 2011), a resource commonly utilized as a benchmark in the field of
temporal causal discovery (Löwe et al. (2022); Khanna & Tan (2019); Assaad et
al. (2022)). Each time series in the dataset simulates the neural signals for
an individual human subject, encompassing $d=5$ distinct nodes. Notably, the
majority of time series in the Nestim dataset have the same summary causal
graph. In order to create MTS with multiple regimes, we generated additional
synthetic data using the generative process outlined in the previous section
5.1. Our choice of baseline models aligns with those utilized in the IT
monitoring data evaluation, as detailed in Section 6.2. Table 3 shows that
CASTOR successfully infers exact number of regimes and it succeeds to detect
80% of the real FMRI data. Even if the regime accuracy is 80%, CASTOR shows
robustness and versatility in identifying complex temporal relationships
within this type of data by achieving similar results to the baselines that
have access to the true regime indices.

##  6 Related works

Causal structure learning from time series.  Assaad et al. (2022) offer an
extensive survey of methods for learning temporal causal relationships. Most
notably, Granger causality is the primary approach used for causal discovery
from time series (Amornbunchornvej et al., 2019; Wu et al., 2020; Löwe et al.,
2022; Bussmann et al., 2021; Xu et al., 2019). However, it is unable to
accommodate instantaneous effects. DYNOTEARS (Pamfil et al., 2020), on the
other hand, leverages the acyclicity constraint established by Zheng et al.
(2018) to continuously relax the DAG and differentiably learn instantaneous
and time lagged structures. However, DYNOTEARS is still limited to linear
functional forms. TiMINo (Peters et al., 2013) provides a general theoretical
framework for temporal causal discovery with functional causal models and also
a practical algorithm that learns casual relationships that can be non linear.
However, the aforementioned methods assume that MTS are composed of a single
regime.
Causal structure learning from heterogeneous data. Several studies have sought
to tackle the challenge of causal discovery in heterogeneous data (Huang et
al., 2020; Saeed et al., 2020; Zhou et al., 2022; Günther et al., 2023;
Saggioro et al., 2020). Remarkably, Huang et al. (2020) address heterogeneous
time series by modulating causal relationships through a regime index. While
it provides a summary graph highlighting behavioral changes across regimes,
they cannot infer individual causal graphs and necessitates prior knowledge of
number of regime. Meanwhile, LIN (Liu & Kuang, 2023) investigates the problem
of causal structure learning from MTS with interventional data but in the
absence of domain (observational or interventional) data. LIN cannot learn
different graphs: it learns the indices of different domains and one causal
graph (represents the instantaneous links) per MTS. Finally, Saggioro et al.
(2020) assume knowledge of the number of regimes and propose the inference of
only time-lagged links. Furthermore, they evaluate their algorithm on graphs
with a limited number of nodes.

##  7 Conclusion

The task of inferring temporal causal graphs from observational time series is
essential in numerous fields. It involves modeling linear and non-linear
relationships and identifying multiple regimes, often without prior knowledge
of regime indices. We introduce CASTOR, a new score-based with proven
convergence properties, preventing the need for prior knowledge of regimes.
Our method demonstrates superior performance in handling both linear and non-
linear relationships across multiple regimes in synthetic and real datasets.
By learning a full temporal causal graph for each regime, CASTOR distinguishes
itself as a notably interpretable method for causal discovery in heterogeneous
time series. CASTOR holds significant promise, and we envision future
applications in real-world scenarios such as epilepsy detection, alongside
exploring its adaptability to non-Gaussian noise.

## References

  * Amornbunchornvej et al. (2019) Chainarong Amornbunchornvej, Elena Zheleva, and Tanya Y Berger-Wolf.  Variable-lag granger causality for time series analysis.  In _IEEE International Conference on Data Science and Advanced Analytics (DSAA)_ , 2019.
  * Assaad et al. (2022) Charles K Assaad, Emilie Devijver, and Eric Gaussier.  Survey and evaluation of causal discovery methods for time series.  _Journal of Artificial Intelligence Research_ , 73:767–819, 2022.
  * Barabási & Albert (1999) Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  _Science_ , 286(5439):509–512, 1999.
  * Brouillard et al. (2020) Philippe Brouillard, Sébastien Lachapelle, Alexandre Lacoste, Simon Lacoste-Julien, and Alexandre Drouin.  Differentiable causal discovery from interventional data.  _In Advances in Neural Information Processing Systems_ , 33:21865–21877, 2020.
  * Bussmann et al. (2021) Bart Bussmann, Jannes Nys, and Steven Latré.  Neural additive vector autoregression models for causal discovery in time series.  In _Discovery Science_ , 2021.
  * Dempster et al. (1977) Arthur P Dempster, Nan M Laird, and Donald B Rubin.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the royal statistical society: series B (methodological)_ , 39(1):1–22, 1977.
  * Gong et al. (2022) Wenbo Gong, Joel Jennings, Cheng Zhang, and Nick Pawlowski.  Rhino: Deep causal temporal relationship learning with history-dependent noise.  _arXiv preprint arXiv:2210.14706_ , 2022.
  * Granger (1969) Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  _Econometrica: journal of the Econometric Society_ , pp. 424–438, 1969.
  * Günther et al. (2023) Wiebke Günther, Urmi Ninad, and Jakob Runge.  Causal discovery for time series from multiple datasets with latent contexts.  _arXiv preprint arXiv:2306.12896_ , 2023.
  * Hariton & Locascio (2018) Eduardo Hariton and Joseph J Locascio.  Randomised controlled trials—the gold standard for effectiveness research.  _BJOG: an international journal of obstetrics and gynaecology_ , 125(13):1716, 2018.
  * Huang et al. (2020) Biwei Huang, Kun Zhang, Jiji Zhang, Joseph Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf.  Causal discovery from heterogeneous/nonstationary data.  _Journal of Machine Learning Research_ , 21(1), 2020.
  * Hyvärinen et al. (2010) Aapo Hyvärinen, Kun Zhang, Shohei Shimizu, and Patrik O Hoyer.  Estimation of a structural vector autoregression model using non-gaussianity.  _Journal of Machine Learning Research_ , 11(5), 2010.
  * Khanna & Tan (2019) Saurabh Khanna and Vincent YF Tan.  Economy statistical recurrent units for inferring nonlinear granger causality.  _arXiv preprint arXiv:1911.09879_ , 2019.
  * Liu & Kuang (2023) Chenxi Liu and Kun Kuang.  Causal structure learning for latent intervened non-stationary data.  In _International Conference on Machine Learning_ , 2023.
  * Löwe et al. (2022) Sindy Löwe, David Madras, Richard Zemel, and Max Welling.  Amortized causal discovery: Learning to infer causal graphs from time-series data.  In _Conference on Causal Learning and Reasoning_ , 2022.
  * McCoy (2017) C Eric McCoy.  Understanding the intention-to-treat principle in randomized controlled trials.  _Western Journal of Emergency Medicine_ , 18(6):1075, 2017.
  * Moraffah et al. (2021) Raha Moraffah, Paras Sheth, Mansooreh Karami, Anchit Bhattacharya, Qianru Wang, Anique Tahir, Adrienne Raglin, and Huan Liu.  Causal inference for time series analysis: Problems, methods and evaluation.  _Knowledge and Information Systems_ , 63:3041–3085, 2021.
  * Newman (2018) Mark Newman.  _Networks_.  Oxford university press, 2018.
  * Pamfil et al. (2020) Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam.  Dynotears: Structure learning from time-series data.  In _International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Peters & Bühlmann (2014) Jonas Peters and Peter Bühlmann.  Identifiability of gaussian structural equation models with equal error variances.  _Biometrika_ , 101(1):219–228, 2014.
  * Peters et al. (2012) Jonas Peters, Joris Mooij, Dominik Janzing, and Bernhard Schölkopf.  Identifiability of causal graphs using functional models.  _arXiv preprint arXiv:1202.3757_ , 2012.
  * Peters et al. (2013) Jonas Peters, Dominik Janzing, and Bernhard Schölkopf.  Causal inference on time series using restricted structural equation models.  _In Advances in neural information processing systems_ , 2013.
  * Peters et al. (2017) Jonas Peters, Dominik Janzing, and Bernhard Schölkopf.  _Elements of causal inference: foundations and learning algorithms_.  The MIT Press, 2017.
  * Rahmani et al. (2023) Abdellah Rahmani, Arun Venkitaraman, and Pascal Frossard.  A meta-gnn approach to personalized seizure detection and classification.  In _IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2023.
  * Runge (2018) Jakob Runge.  Causal network reconstruction from time series: From theoretical assumptions to practical estimation.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 28(7), 2018.
  * Saeed et al. (2020) Basil Saeed, Snigdha Panigrahi, and Caroline Uhler.  Causal structure discovery from distributions arising from mixtures of dags.  In _In International Conference on Machine Learning_ , 2020.
  * Saggioro et al. (2020) Elena Saggioro, Jana de Wiljes, Marlene Kretschmer, and Jakob Runge.  Reconstructing regime-dependent causal relationships from observational time series.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 30(11), 2020.
  * Samé et al. (2011) Allou Samé, Faicel Chamroukhi, Gérard Govaert, and Patrice Aknin.  Model-based clustering and segmentation of time series with changes in regime.  _In Advances in Data Analysis and Classification_ , 5:301–321, 2011.
  * Schölkopf et al. (2021) Bernhard Schölkopf, Francesco Locatello, Stefan Bauer, Nan Rosemary Ke, Nal Kalchbrenner, Anirudh Goyal, and Yoshua Bengio.  Towards causal representation learning 2021.  _arXiv preprint arXiv:2102.11107_ , 2021.
  * Smith et al. (2011) Stephen M Smith, Karla L Miller, Gholamreza Salimi-Khorshidi, Matthew Webster, Christian F Beckmann, Thomas E Nichols, Joseph D Ramsey, and Mark W Woolrich.  Network modelling methods for fmri.  _Neuroimage_ , 54(2):875–891, 2011.
  * Tang et al. (2021) Siyi Tang, Jared A Dunnmon, Khaled Saab, Xuan Zhang, Qianying Huang, Florian Dubost, Daniel L Rubin, and Christopher Lee-Messer.  Self-supervised graph neural networks for improved electroencephalographic seizure analysis.  _arXiv preprint arXiv:2104.08336_ , 2021.
  * Thiesson et al. (2013) Bo Thiesson, Christopher Meek, David Maxwell Chickering, and David Heckerman.  Learning mixtures of dag models.  _arXiv preprint arXiv:1301.7415_ , 2013.
  * Wu et al. (2020) Tailin Wu, Thomas Breuel, Michael Skuhersky, and Jan Kautz.  Discovering nonlinear relations with minimum predictive information regularization.  _arXiv preprint arXiv:2001.01885_ , 2020.
  * Xu et al. (2019) Chenxiao Xu, Hao Huang, and Shinjae Yoo.  Scalable causal graph learning through a deep neural network.  In _ACM international conference on information and knowledge management_ , 2019.
  * Zheng et al. (2018) Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  _In Advances in neural information processing systems_ , 2018.
  * Zheng et al. (2020) Xun Zheng, Chen Dan, Bryon Aragam, Pradeep Ravikumar, and Eric Xing.  Learning sparse nonparametric dags.  In _International Conference on Artificial Intelligence and Statistics_ , 2020.
  * Zhou et al. (2022) Fangting Zhou, Kejun He, and Yang Ni.  Causal discovery with heterogeneous observational data.  In _Uncertainty in Artificial Intelligence_ , 2022.
  * Zhu et al. (1997) Ciyou Zhu, Richard H Byrd, Peihuang Lu, and Jorge Nocedal.  Algorithm 778: L-bfgs-b: Fortran subroutines for large-scale bound-constrained optimization.  _ACM Transactions on mathematical software (TOMS)_ , 1997.

##  Appendix A Expectation-Maximization derivation

In this section, we shall elucidate the computational details surrounding the
resolution of our optimization problem. Specifically, we will provide clarity
on the various equations introduced in Section 3, namely, Eq (8, 10, 11, 15).
E-step.  We model regime participation through a binary latent variable
$z_{t}\in\mathbf{R}^{N_{w}}$; $\bm{X}_{t}$ belongs to regime $u\Rightarrow
z_{t,u}=1$.

$\displaystyle\gamma_{t,u}$
$\displaystyle=p\left(z_{t,u}=1\mid\bm{X}_{t},\bm{G}^{u}_{\\{0:L\\}}\right)$
(16) $\displaystyle=\frac{p\left(z_{t,u}=1\right)p\left(\bm{X}_{t}\mid
z_{t,u}=1,\bm{G}^{u}_{\\{0:L\\}}\right)}{\sum_{j=1}^{N_{w}}p\left(z_{t,j}=1\right)p\left(\bm{X}_{t}\mid
z_{t,j}=1,\bm{G}^{j}_{\\{0:L\\}}\right)}$
$\displaystyle=\frac{\pi_{t,u}(\alpha)f^{u}\left(\bm{X}_{t}\right)}{\sum_{j=1}^{N_{w}}\pi_{t,j}(\alpha)f^{j}\left(\bm{X}_{t}\right)}$

M-step. Having estimated probabilities $\gamma_{t,u}$ in the E-step, we can
now maximise the expected posterior distribution given the MTS
$(\bm{X}_{t})_{t\in\mathcal{T}}$ and we have:

$\displaystyle\sup_{\theta,\alpha}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=0}^{|\mathcal{T}|}\gamma_{t,u}\log\pi_{t,u}(\alpha)f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{
s.t }\bm{G}^{u}_{0}\text{ is a DAG,}$ (17)
$\displaystyle\iff\left\\{\begin{array}[]{l}\max_{\alpha}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\ln\left(\pi_{t,u}(\alpha)\right)\\\
\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{ s.t
}\bm{G}^{u}_{0}\text{ is a DAG,}\end{array}\right.$

We know
$f^{u}\left(\bm{X}_{t}\right)=\mathcal{N}\left(\bm{X}_{t}\bm{G}_{0}^{u}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}_{\tau}^{u},I\right)$,
hence:

$\displaystyle\iff\sup_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,\text{ s.t
}\bm{G}^{u}_{0}\text{ is a DAG,}$ (18)
$\displaystyle\iff\min_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t\in\mathcal{E}_{u}}\left\|\bm{X}_{t}-\left(\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}\right)\right\|_{F}^{2}+\lambda|\mathcal{G}^{u}|,\text{
s.t }\bm{G}^{u}_{0}\text{ is a DAG}$
$\displaystyle\iff\min_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\left\|\bm{X}_{t}-\left(\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}\right)\right\|_{F}^{2}+\lambda|\mathcal{G}^{u}|,\text{
s.t }\bm{G}^{u}_{0}\text{ is a DAG}$
$\displaystyle\iff\min_{\theta}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\gamma_{t,u}\left\|\bm{X}_{t}-\left(\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}\right)\right\|_{F}^{2}+\lambda|\mathcal{G}^{u}|+\frac{\rho}{2}h\left(\bm{G}^{u}_{0}\right)^{2}+\alpha
h\left(\bm{G}^{u}_{0}\right),$

The only difference between the linear and the non linear cases is how we
estimate the mean of the normal distribution $f^{u}$ for every regime $u$. As
we mentioned in section 3.2, we estimate these means using NNs and we have
$f^{u}_{i}\left(\bm{X}_{t}\right)=\mathcal{N}\left(\psi^{u}_{i}\left(\phi^{u}_{i}(\bm{X}_{t}),\phi^{u,\text{lag}}_{i}(\bm{X}_{t}^{\text{lag}})\right),1\right)$,
Hence, our M-step for non-linear CASTOR:

$\displaystyle\min_{\theta,\mathcal{G}}\frac{1}{|\mathcal{T}|}\sum_{u=1}^{N_{w}}\sum_{t=1}^{|\mathcal{T}|}\sum_{i=1}^{d}\gamma_{t,u}\mathcal{L}(X^{i}_{t},\psi^{u}_{i}\left(\phi^{u}_{i}(\bm{X}_{t}),\phi^{u,\text{lag}}_{i}(\bm{X}_{t}^{\text{lag}})\right))+\lambda|\mathcal{G}^{u}|+\frac{\rho}{2}h\left(\bm{G}^{u}_{0}\right)^{2}+\alpha
h\left(\bm{G}^{u}_{0}\right)$ (19)

##  Appendix B Regime and causal graphs identifiability

In this section, we concentrate on establishing the identifiability of regimes
and causal graphs within the CASTOR framework. Before diving into the details,
let us set and clarify the required assumptions.

###### Definition 3

(Causal Stationarity Runge (2018)). The time series (that has one regime)
process $(\bm{X}_{t})_{t\in\mathcal{T}}$ with a graph $\mathcal{G}$ is called
causally stationary over a time index set $\mathcal{T}$ if and only if for all
links $X_{t-\tau}^{i}\rightarrow X_{t}^{j}$ in the graph

$X_{t-\tau}^{i}\not\\!\perp\\!\\!\\!\perp
X_{t}^{j}\mid\bm{X}_{<t}\backslash\left\\{X_{t-\tau}^{i}\right\\}\text{ holds
for all }t\in\mathcal{T}$

This elucidates the inherent characteristics of the time-series data
generation mechanism, thereby validating the choice of the auto-regressive
model. In our setting, we generalize Causal Stationarity as follows:

###### Assumption 1

(Causal Stationarity for time series with multiple regimes). The time series
process $(\bm{X}_{t})_{t\in\mathcal{T}}$ comprise multiple regimes $K$, where
$K$ is the number of regime, we note $\mathcal{E}_{u}=\\{t|\gamma_{t,u}=1\\}$
the set of time indexes where the regime $u$ is active, and
$\mathcal{T}=\cup_{u}\mathcal{E}_{u}$. $(\bm{X}_{t})_{t\in\mathcal{T}}$ with a
graph $\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,K\\}}$ is called causally stationary
over a time index set $\mathcal{T}$ if and only if for all $u\in\\{1,..,K\\}$,
$(\bm{X}_{t})_{t\in\mathcal{E}_{u}}$ is causal stationary with graph
$\mathcal{G}^{u}$ for time index set $\mathcal{E}_{u}$.

###### Definition 4

(Causal Markov Property, Peters et al. (2017)). Given a DAG $\mathcal{G}$ and
a joint distribution $p$, this distribution is said to satisfy causal Markov
property w.r.t. the DAG $\mathcal{G}$ if each variable is independent of its
non-descendants given its parents.

This is a common assumptions for the distribution induced by an SEM. With this
assumption, one can deduce conditional independence between variables from the
graph.

###### Assumption 2

(Causal Markov Property for multiple regimes). Given a set of DAGs
$\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,K\\}}$ and a set of joint distribution
$\\{p^{(u)}\\}_{u\in\\{1,..,K\\}}$, we say that this set of distributions
satisfies causal Markov property w.r.t. the set of DAGs
$\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,K\\}}$ if for every $u$: $p^{(u)}$ satisfy
causal Markov property w.r.t the DAG $\mathcal{G}^{u}$.

###### Definition 5

(Causal Minimality, Gong et al. (2022)). Consider a distribution $p$ and a DAG
$\mathcal{G}$, we say this distribution satisfies causal minimality w.r.t.
$\mathcal{G}$ if it is Markovian w.r.t. $\mathcal{G}$ but not to any proper
subgraph of $\mathcal{G}$.

###### Assumption 3

(Causal Minimality for multiple regimes). Given a set of DAGs
$\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,K\\}}$ and a set of joint distribution
$\\{p^{(u)}\\}_{u\in\\{1,..,K\\}}$, we say that this set of distributions
satisfies causal minimality w.r.t. the set of DAGs
$\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,K\\}}$ if for every $u$: $p^{(u)}$ satisfy
causal minimality w.r.t the DAG $\mathcal{G}^{u}$.

###### Assumption 4

(Causal Sufficiency). A set of observed variables $V$ is causally sufficient
for a process $\bm{X}_{t}$ if and only if in the process every common cause of
any two or more variables in $\bm{V}$ is in $\bm{V}$ or has the same value for
all units in the population.

This assumption implies there are no latent confounders present in the time-
series data.

###### Assumption 5

(Well-defined Density). We assume the joint likelihood induced by the CASTOR
SEM (Eq. (3)) is absolutely continuous w.r.t. a Lebesgue or counting measure
and $\forall u:\left|\log
p^{(u)}\left(\bm{X}_{t\in\mathcal{E}_{u}};\mathcal{G}^{u}\right)\right|<\infty$
for all possible $\mathcal{G}^{u}$.

This presumption ensures that the resulting distribution possesses a well-
defined probability density function. It is also necessary for Markov
factorization property.

###### Assumption 6

(CASTOR in DAG space). We assume the CASTOR framework can only return the
solutions from DAG space.

###### Definition 6

(Source model Gong et al. (2022)) For any MTS defined as follows:

$X_{t}^{i}=f_{i,t}\left(\mathbf{Pa}_{G}^{i}(<t),\mathbf{Pa}_{G}^{i}(t)\right)+\epsilon_{t}^{i}$
(20)

A source model characterizes the initial conditions:

$X_{s}^{i}=f_{i,s}\left(\mathbf{Pa}_{G}^{i}\right)+\epsilon_{s}^{i}$ (21)

for $s\in[0,\mathcal{S}]$, where $\mathcal{S}$ is the length for the initial
conditions and $\mathbf{Pa}_{G}^{i}$ contains the parents for node $i$. We
define $p_{s}\left(\bm{X}_{0:\mathcal{S}}\right)$ as the induced joint
distribution for the initial conditions.

###### Assumption 7

(CASTOR initial condition). We assume that the initial conditions are known
and source model is identifiable.

###  B.1 Proof of theorem 1

Assuming the aforementioned assumptions we want to prove the theorem 1.
We consider $\mathcal{G}=\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,N_{w}\\}}$,
$\mathcal{E}=\cup_{u=1}^{N_{w}}\mathcal{E}_{u}$ where $N_{w}$ is the number of
window, $\mathcal{G}^{*}=\\{\mathcal{G}^{*,u}\\}_{u\in\\{1,..,K\\}}$, $K$ is
the exact number of true regimes and
$\mathcal{E}^{*}=\cup_{u=1}^{K}\mathcal{E}^{*}_{u}$. We denote
$\mathcal{E}_{c}\mathcal{E}_{\ell}^{*}$ the set of time indexes that is shared
between regime $c$ of our model estimation and the true regime $\ell$ and
$\pi_{\ell}^{*}:=\frac{\left|\mathcal{E}_{\ell}^{*}\right|}{T},\pi_{c}:=\frac{\left|\mathcal{E}_{c}\right|}{T}$,
$\pi_{c\ell}:=\frac{\left|\mathcal{E}_{c}\mathcal{E}_{\ell}^{*}\right|}{T}$.
Our objective is to prove that for any estimation
$(\hat{\mathcal{G}},\hat{\mathcal{E}})$ : if $\exists u\in\\{1,...,K\\}\text{
s.t. }\hat{\mathcal{G}^{u}}\text{ disagree with}$ $\mathcal{G}^{*,u}$ on
instantaneous or/and time lagged link, or any regime
$\hat{\mathcal{E}}_{u}\in\hat{\mathcal{E}}$ is close to none of the true
regimes in the sense of Kullback–Leibler divergence:
$\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)>\mathcal{S}(\hat{\mathcal{G}},\hat{\mathcal{E}}).$
We have by Eq (9)

$\mathcal{S}(\mathcal{G},\mathcal{E}):=\sup_{\theta,\mathcal{G}}\frac{1}{T}\sum_{u=1}^{N_{\mathrm{u}}}\sum_{t\in\mathcal{E}_{u}}\log
f^{u}\left(\bm{X}_{t}\right)-\lambda|\mathcal{G}^{u}|,$

where $\lambda$ is the sparsity penalty coefficient and
$f^{u}\left(\bm{X}_{t}\right):=\prod_{j=1}^{d}f^{u}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right)$
with
$f^{u}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right)$
the function used to describe the distribution family in Eq (7).
We will structure the proof as follows:

  * •

Prove that if the score is optimized, then all the estimated regimes will be
pure (have only elements of the same true regime).

  * •

Prove that, when the regimes are pure and $N_{w}=K$, we have
$\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)>\mathcal{S}(\hat{\mathcal{G}},\hat{\mathcal{E}}).$
for any estimation $\hat{\mathcal{G}}$ where $\exists u\in\\{1,...,K\\}\text{
s.t }\hat{\mathcal{G}^{u}}\text{ disagrees with}$ $\mathcal{G}^{*,u}$ on
instantaneous or/and time lagged link.

####  B.1.1 Optimizing the score will lead to pure regimes

Ignoring penalty terms, we have:

$\displaystyle-\mathcal{S}(\mathcal{G},\mathcal{E})$
$\displaystyle=-\sup_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\pi_{c\ell}\frac{1}{\left|\mathcal{E}_{e}\mathcal{E}_{\ell}^{*}\right|}\sum_{t\in\mathcal{E}_{e}\mathcal{E}_{\ell}^{*}}\left[\log
f^{c}(\bm{X}_{t})\right]$ (22)
$\displaystyle\rightarrow-\sup_{\phi}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\log f^{c}\right]$
$\displaystyle=-\sup_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\sum_{i=1}^{d}\log
f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right]$
$\displaystyle=-\sup_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}$
$\displaystyle\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[-\log\frac{p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)}{f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)}+\log
p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right]$
$\displaystyle=-\sup_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}$
$\displaystyle\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\Big{[}-\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)$
$\displaystyle-\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\Big{]}$
$\displaystyle=-\sup_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[-\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
$\displaystyle-\sup_{\theta}\sum_{c=1}^{N_{\mathrm{c}}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[-\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\right]$
$\displaystyle=\inf_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{j=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
$\displaystyle+\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{j=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\right]$
$\displaystyle=\inf_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
$\displaystyle+\sum_{\ell=1}^{K}\sum_{i=1}^{d}\left(\sum_{c=1}^{N_{\mathrm{c}}}\pi_{c\ell}\right)\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\right]$
$\displaystyle=\inf_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
$\displaystyle+\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{\ell}^{*}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\right]$

Note that $\theta$ could be the parameters of the neural networks used in Eq
(15) for non linear causal relationship or
$\theta=\\{\mathcal{G}^{u}\\}_{u\in\\{1,..,N_{w}\\}}$ for linear case Eq 10.

For the score of ground truth (ignoring penalty terms):

$\displaystyle-\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)$
$\displaystyle\rightarrow
0+\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{\ell}^{*}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{H}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\right)\right](\text{
by Assumption 5 \ref{ass5}) }$ (23)

Combining Equation (22) and Equation (23) , we have (considering penalty
terms):

$\displaystyle\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)-\mathcal{S}(\mathcal{G},\mathcal{E})=$
$\displaystyle\inf_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}$
(24) $\displaystyle\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
$\displaystyle+\lambda\left(\sum_{c=1}^{N_{w}}|\mathcal{G}^{c}|-\sum_{\ell=1}^{K}\left|\mathcal{G}^{*,\ell}\right|\right)$

The first term in Equation (24) is the score term, others are penalty term.
In the following lines, our goal is to demonstrate that optimizing the score
term ensures that all identified regimes will accurately match the real
regimes. In other words, each estimated regime will be a true representation
of an actual one. Additionally, by shifting samples from less significant
regimes (regimes with few samples) to the most similar significant regimes,
our variable $N_{w}$ will eventually stabilize at the value of K. To do this,
we will proceed by contradiction:
Suppose the score term in Eq (24) is optimized and there exists a regime $e$
that is not pure, i.e., there exist $a,b\in[K]$ with $a\neq b$ but
$\pi_{ea}>0$ and $\pi_{eb}>0$. Since they are different distributions for two
different regimes with two different causal graphs, there exists
$i\in\\{1,...,d\\}$ such that
$p^{(a)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,a}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,a}}^{i}(t)\right)\right)\neq
p^{(b)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,b}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,b}}^{i}(t)\right)\right)$.
Then the score term in Equation (24) has the following lower bound:

$\displaystyle\inf_{\theta}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{e\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)\right]$
(25) $\displaystyle\geq$
$\displaystyle\inf_{\theta}\Big{\\{}\pi_{ea}\underset{\bm{X}_{t}\sim
p^{(a)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(a)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)\right]$
$\displaystyle+\pi_{eb}\underset{\bm{X}_{t}\sim
p^{(b)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(b)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)\right]\Big{\\}}$

As we assumed that the score term in Eq (24) is optimized, it means that:

$\displaystyle
0=\inf_{\theta}\sum_{c=1}^{N_{w}}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{c\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{c}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{c}}^{i}(t)\right)\right)\right]$
(26) $\displaystyle\Rightarrow
0=\inf_{\theta}\sum_{\ell=1}^{K}\sum_{i=1}^{d}\pi_{e\ell}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\left[\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)\right]$
$\displaystyle\Rightarrow\left\\{\begin{array}[]{l}\mathrm{D}_{\mathrm{KL}}\left(p^{(a)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)=0\\\
\mathrm{D}_{\mathrm{KL}}\left(p^{(b)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)\|f^{e}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{e}}^{i}(t)\right)\right)=0\par\end{array}\right.$
$\displaystyle\Rightarrow
p^{(a)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)=p^{(b)}\left(X_{t}^{i}\mid\left(\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{*,\ell}}^{i}(t)\right)\right)$

and the last line, Eq (26), is a contradiction because the two distributions
represent two different regimes with two different graphs. Hence, if the score
term of Eq (24) is optimized all the estimated regimes will be pure.
First case: if we matched the samples of less significant regimes to the wrong
regimes, the regime is not pure and then the score term is not optimized
(contradiction).
Second case: If we eliminate a lot of regimes such that $N_{w}\leq K-1$, at
least one of our estimated regimes will not be pure and this contradicts the
assumption of optimized score term (same reasoning).
Based on this reasoning, optimizing the score term of Equation 24 will ensure
convergence to the true number of regimes and also every regime will be pure.

####  B.1.2 In case of edge disagreement
$\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)>\mathcal{S}(\hat{\mathcal{G}},\hat{\mathcal{E}})$

Now we will show that Eq (24 )is positive, if $\exists u\in\\{1,...,K\\}\text{
s.t }\hat{\mathcal{G}^{u}}\text{ disagrees with}$ $\mathcal{G}^{*,u}$ on
instantaneous or/and time lagged link.
We assume that each estimated regime
$\hat{\mathcal{E}}_{c}\left(c\in\left\\{1,\ldots,N_{w}\right\\},N_{w}\geq
K\right)$ contains samples from same true regime. Then Equation (24) has lower
bound:

$\begin{array}[]{r}\inf_{\theta}\sum_{\ell=1}^{K}\pi_{\ell}^{*}\underset{\bm{X}_{t}\sim
p^{(\ell)}}{\mathbb{E}}\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\|f^{\ell}\right)\\\
\geq\left(\min_{\ell}\pi_{\ell}^{*}\right)\inf_{\theta}\sum_{\ell=1}^{K}\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\|f^{\ell}\right)\end{array}$
(27)

Equation (27) is positive if and only if $\eta(\mathcal{G})$ is positive.

$\eta(\mathcal{G}):=\inf_{\theta}\sum_{\ell=1}^{K}\mathrm{D}_{\mathrm{KL}}\left(p^{(\ell)}\|f^{\ell}\right)$
(28)

Let assume that $\exists r\in\in\\{1,...,K\\}\text{ s.t
}\hat{\mathcal{G}^{r}}\text{ disagrees with}$ $\mathcal{G}^{*,r}$ on
instantaneous or/and time lagged link. We follow the same intuition as Gong et
al. (2022); Peters et al. (2013; 2017):

$\displaystyle\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)$ (29)
$\displaystyle=$
$\displaystyle\mathrm{D}_{\mathrm{KL}}\left[p^{(r)}_{s}\|f^{r}_{s}\right]+\sum_{t=\mathcal{S}+1}^{\left|\mathcal{E}_{r}^{*}\right|}\mathbb{E}_{p^{(r)}\left(\bm{X}_{0:t-1}\right)}\left[\mathrm{D}_{\mathrm{KL}}\left[p^{(r)}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\|f^{r}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\right]\right]$

Based on Assumption 7, we know that the source model is known, which leads to
the following result:

$\displaystyle\text{Eq
(\ref{27})}\iff\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)=\sum_{t=\mathcal{S}+1}^{\left|\mathcal{E}_{r}^{*}\right|}\mathbb{E}_{p^{(r)}\left(\bm{X}_{0:t-1}\right)}\left[\mathrm{D}_{\mathrm{KL}}\left[p^{(r)}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\|f^{r}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\right]\right]$
(30)

We will show that $\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)$ is
positive in two cases:

  * •

Disagreement on lagged parents only.  This means that for all
$t\in[\mathcal{S}+1,T]$, the instantaneous connections at $t$ for
$\hat{\mathcal{G}^{r}}$ and $\mathcal{G}^{*,r}$ are the same, and $\exists
t\in[\mathcal{S}+1,T]\text{ and }i\in\\{1,...,d\\}$ such that
$\mathbf{Pa}_{\mathcal{G}^{*,r}}^{X_{t}^{i}}(<t)\neq\overline{\mathbf{Pa}}_{\hat{\mathcal{G}^{r}}}^{X_{t}^{i}}(<t)$.
We can use a similar argument as the theorem 1 in Peters et al. (2013).
Without loss of generality, we assume under $\hat{\mathcal{G}^{r}}$, we have
$X^{j}_{t-\tau}\rightarrow X^{i}_{t}$ and there is no connections between them
under $\mathcal{G}^{*,r}$. Thus, from Markov conditions, we have

$X^{i}_{t}\perp\\!\\!\\!\perp
X^{j}_{t-\tau}\mid\bm{X}_{0:t-1}\cup\mathrm{ND}_{t}^{X^{i}}\backslash\left\\{X^{i}_{t},X^{j}_{t-\tau}\right\\}$

under $\mathcal{G}^{*,r}$, where $\mathrm{ND}_{t}^{X^{i}}$ are the non-
descendants of node $X^{i}_{t}$ at some time $t$. However, from the causal
minimality and Proposition 6.16 in Peters et al. (2017), we have

$X^{i}_{t}\not\\!\perp\\!\\!\\!\perp
X^{j}_{t-\tau}\mid\bm{X}_{0:t-1}\cup\mathrm{ND}_{t}^{X^{i}}\backslash\left\\{X^{i}_{t},X^{j}_{t-\tau}\right\\}$

under $\hat{\mathcal{G}^{r}}$, and we have
$\mathrm{D}_{\mathrm{KL}}\left[p^{(r)}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\|f^{r}\left(\bm{X}_{t}\mid\bm{X}_{0:t-1}\right)\right]\neq
0$. Hence, $\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)\neq 0$

  * •

Disagreement on instantaneous parents.  In this Section we will use two
different results one for the linear and the other one for the non linear
case.

    * –

Linear case.  For this case, we will use Theorem 1 in Peters & Bühlmann
(2014). In this theorem, the author confirms that the graph is identifiable
for linear models with Gaussian additive noise, if for each
$j\in\\{1,\ldots,d\\}$, the weights of the causal relations $\beta_{jk}\neq 0$
for all $k\in\mathbf{PA}_{j}^{\mathcal{G}_{0}}$. For our instantaneous links,
we have all the weights of the parents are non null. Hence, the instantaneous
links are identifiable. Otherwise if
$\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)\neq 0$

    * –

Non linear case.  Using Theorem 2 from Peters et al. (2012), we can notice
that our instantaneous links are Identifiable Functional Model Class,
$(\mathcal{B},\mathcal{F})$-IFMOC, they belongs exactly to the 3rd class of
Lemma 3 in Peters et al. (2012): nonlinear ANMs:
$\mathcal{F}_{3}=\\{f(X,n)=\phi(X)+n\\}$ $\mathcal{B}_{3}=\\{(\phi,X,N)$ not
lin., Gauss, Gauss $\\}\backslash\tilde{B}_{3}$. Hence our instantaneous links
are identifiable, otherwise,
$\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)\neq 0$.

Based on the above reasoning, we can show that if $\exists
r\in\in\\{1,...,K\\}\text{ s.t., }\hat{\mathcal{G}^{r}}\text{ disagree with}$
$\mathcal{G}^{*,r}$ on instantaneous or/and time lagged links,
$\mathrm{D}_{\mathrm{KL}}\left(p^{(r)}\|f^{r}\right)\neq 0$.
Thus,$\eta(\mathcal{G})>0$. Then as we assume in Theorem 1 that $\lambda$ is
sufficient small would implies Equation (26) is positive.
If $|\hat{\mathcal{G}^{r}}|\geq\left|\mathcal{G}^{*,r}\right|$ then clearly Eq
26 is positive. Let
$\mathbb{G}^{+}:=\left\\{\hat{\mathcal{G}^{r}}\in\mathbb{G}||\hat{\mathcal{G}^{r}}|<|\mathcal{G}^{*,r}\mid\right\\}$.
To make sure that we have
$\mathcal{S}\left(\mathcal{G}^{*},\mathcal{E}^{*}\right)-\mathcal{S}(\mathcal{G},\mathcal{E})>0$
for all $\mathcal{G}\in\mathbb{G}^{+}$, we need to pick $\lambda$ sufficiently
small. Choosing
$0<\lambda<\min_{\mathcal{G}\in\mathbb{G}^{+}}\frac{\eta(\mathcal{G})}{\left(\sum_{c=1}^{N_{w}}|\mathcal{G}^{c}|-\sum_{\ell=1}^{K}\left|\mathcal{G}^{*,\ell}\right|\right)}$
is sufficient.

##  Appendix C Further experimental results

###  C.1 Synthetic data

We employ the Erdos–Rényi (ER) (Newman, 2018) model with mean degrees of 1 or
2 to generate lagged graphs, and the Barabasi–Albert (BA) (Barabási & Albert,
1999) model with mean degrees 4 for instantaneous graphs. The maximum number
of lags, $L$, is set at 1 or 2. We experiment with varying numbers of nodes
$\\{5,10,20,40\\}$ and different numbers of regimes $\\{2,3,4,5\\}$, each
representing diverse causal graphs. The length of each regime is randomly
sampled from the set $\\{300,400,500,600\\}$.

  * •

Linear case.  Data is generated as follows:

$\forall u\in\\{1,...,K\\},\forall
t\in\mathcal{E}_{u}:\bm{X}_{t}=\bm{X}_{t}\bm{G}^{u}_{0}+\sum_{\tau=1}^{L}\bm{X}_{t-\tau}\bm{G}^{u}_{\tau}+\bm{\epsilon}_{t},$

with $\bm{G}^{u}_{0}$ is adjacency matrix of the generated graph by BA model,
$\forall\tau\in\\{1,..,L\\}:\bm{G}^{u}_{\tau}$ are the adjacency of the time
lagged graphs generated by ER and $\epsilon_{t}\sim\mathcal{N}(0,I)$, follows
to a normal distribution.

  * •

Non-linear case.  The formulation used to generated the data is:

$\forall u\in\\{1,...,K\\},\forall
t\in\mathcal{E}_{u}:X_{t}^{i}=g^{u}_{i}\left(\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(<t),\mathbf{Pa}_{\mathcal{G}^{u}}^{i}(t)\right)+\epsilon_{t}^{i},$

where $g^{u}_{i}$ is a general differentiable linear/non-linear function and
$\epsilon_{t}^{i}\sim\mathcal{N}(0,1)$, follows a normal distribution. The
function $g^{u}_{i}$ is a random combination between a linear transformation
and a randomly chosen function from the set: $\\{\texttt{Tanh, LeakyReLU,
ReLU}\\}$.

###  C.2 Baselines

All used benchmarks for the synthetic experiments are run by using publicly
available libraries: VARLINGAM (Hyvärinen et al., 2010) is implemenented in
the lingam444https://lingam.readthedocs.io/en/latest/ python package. RPCMCI
(Saggioro et al., 2020) is implemented in
Tigramite555https://jakobrunge.github.io/tigramite/ and DYNOTEARS (Pamfil et
al., 2020) on causalnex666https://causalnex.readthedocs.io/en/latest/ package.
We fine tuned the parameters to achieve the optimal graph for each model.
For CASTOR, an edge threshold of 0.4 is selected. In the linear scenario, we
establish $\zeta=100$ as the minimum regime duration, while in the non-linear
context, $\zeta$ is set at 200. To demonstrate the model’s robustness to the
choice of the window size, we train CASTOR using diverse window sizes,
specifically $w=200$ or $w=300$. For the sparsity coefficient, we use
$\lambda=0.05$. In order to optimise our M-step, we use L-BFGS-B algorithm
(Zhu et al., 1997).

###  C.3 Further experiments: Linear case

|  |  | $K=2$ | $K=3$ | $K=4$
---|---|---|---|---|---
$d$ | Method | Type | SHD | F1 | SHD | F1 | SHD | F1
5 | VARLINGAM |  | Inst.
---
Lag
| 7
---
9
| $22.9_{\pm 5.6}$
---
$8.01_{\pm 12}$
| 13
---
14
| $18.4_{\pm 3.6}$
---
$15.9_{\pm 4.2}$
| 14
---
14
| $24.1_{\pm 6.5}$
---
$8.20_{\pm 6.5}$
RPCMCI |  | Inst.
---
Lag
| -
---
-
| -
---
-
| -
---
-
| -
---
-
| -
---
-
| -
---
-
CASTOR |  | Inst.
---
Lag
| 0
---
0
|  100±0.0
---
100±0.0
| 0
---
0
|  100±0.0
---
100±0.0
| 2
---
1
|  97.3±1.9
---
97.2±2.0
20 | VARLINGAM |  | Inst
---
Lag
| 62
---
48
| $8.83\pm 1.7$
---
$2.96\pm 0.2$
| 90
---
61
| $10.3\pm 2.4$
---
$1.66\pm 1.24$
| 121
---
93
| $14.0\pm 2.4$
---
$2.30\pm 0.7$
RPCMCI |  | Inst
---
Lag
| -
---
67
| -
---
$46.6\pm 9.14$
| -
---
178
| -
---
14.1±2.3
| -
---
-
| -
---
-
CASTOR |  | Inst
---
Lag
| 3
---
0
|  97.0±2.7
---
100±0
| 28
---
9
|  88.1±6.2
---
89.9${{}_{\pm}5.1}$
| 4
---
1
|  98.3±1.4
---
99.7±0.4
Table 4: F1 and SHD Scores by Models and Settings: $d$ indicates number of nodes, $K$ refers to the number of regimes, and ’Type’ refers to the causal links as either instantaneous or time-lagged. |  |  | $K=5$
---|---|---|---
$d$ | Method | Type | SHD | F1
10 | VARLINGAM |  | Inst.
---
Lag
| 59
---
51
| $11.7_{\pm 2.3}$
---
$4.69_{\pm 2.4}$
RPCMCI |  | Inst.
---
Lag
| -
---
-
| -
---
-
CASTOR |  | Inst.
---
Lag
| 4
---
1
|  98.4±1.4
---
97.1±1.4
Table 5: F1 and SHD Scores by Models and Settings: $d$ indicates number of
nodes, $K$ refers to the number of regimes, and ’Type’ refers to the causal
links as either instantaneous or time-lagged.

In this section, we present additional results using linear synthetic data
with varying numbers of nodes, specifically $\\{5,20\\}$, and diverse numbers
of regimes $\\{2,3,4,5\\}$. Our model’s performance is compared against RPCMCI
(Saggioro et al., 2020) and VARLINGAM (Hyvärinen et al., 2010), and the
results are displayed in Table 4, 5. CASTOR not only outperforms RPCMCI but
consistently converges, accurately identifying the number of regimes and
recovering both the regime indices and their respective underlying causal
graphs. In scenarios with 5 different regimes, where RPCMCI fails to converge,
CASTOR infers the true number of regimes, their partitions, and efficiently
learns the causal graphs, achieving more than a 98% F1 score. Furthermore, in
comparison to VARLINGAM (a method capable of modeling both lagged and
instantaneous links in MTS data), CASTOR demonstrates markedly superior
performance.
Figure 4 displays a comparison between the graphs estimated by CASTOR and the
ground truth graphs. CASTOR effectively learns the causal graphs and the
distinct regime partitions. In figure 5 we test CASTOR on data with $L=2$
lags, hence CASTOR needs to estimate 3 adjacency matrices in every regimes
(the matrix for instantaneous links, first time lagged links and second time
lagged links). We can notice that CASTOR estimates well the graphs and infers
the regime partitions.

Figure 4:  The estimated temporal causal graphs for five regimes (Linear
case), with one matrix representing instantaneous links and another
delineating time-lagged relations. The second column showcases the actual
causal graphs, while the final column highlights the discrepancies between the
estimated and true graphs.

Figure 5:  The estimated temporal causal graphs For MTS with 2 lags and 10
nodes and 2 regimes for the linear case

###  C.4 Further experiments: nonlinear case

In this section, we present additional results using non-linear synthetic data
with varying numbers of nodes, specifically $\\{10,20\\}$, and diverse numbers
of regimes $\\{3,5\\}$. In this case, we compare our model against various
baseline models, namely DYNOTEARS (Pamfil et al., 2020) and VARLINGAM
(Hyvärinen et al., 2010). CASTOR showcases performance similar to DYNOTEARS in
modeling non-linear time-lagged relations while simultaneously learning the
regime indexes. Additionally, it effectively identifies instantaneous links.
As depicted in figures 7 and 8, the performance of CASTOR diminishes with an
increase in the number of regimes in non-linear scenarios. This outcome is
understandable since non-linearity poses a more challenging causal discovery
problem, and an increase in the number of regimes augments the number of
parameters, consequently affecting our model’s performance.

Figure 6:  F1 by Models for K=3 Setting, Inst means instantaneous links and
lag time-lagged ones. Figure 7:  The estimated temporal causal graphs for 3
regimes (non Linear case), with one matrix representing instantaneous links
and another delineating time-lagged relations. The second column showcases the
actual causal graphs, while the final column highlights the discrepancies
between the estimated and true graphs.

Figure 8:  The estimated temporal causal graphs for 5 regimes (non Linear
case), with one matrix representing instantaneous links and another
delineating time-lagged relations. The second column showcases the actual
causal graphs, while the final column highlights the discrepancies between the
estimated and true graphs.