Flow based approach for Dynamic Temporal Causal models with non-Gaussian or Heteroscedastic Noises
Authors: Abdellah Rahmani
Pascal Frossard
Venue: NeurIPS 2025 poster ()
Keywords: Bayesian structure learning
Bayesian expectation maximisation
non stationarity
time series
heteroscedasticity
Primary area: probabilistic_methods

Abstract: Understanding causal relationships in multivariate time series is crucial in many scenarios, such as those dealing with financial or neurological data. Many such time series exhibit multiple regimes, i.e., consecutive temporal segments with a priori unknown boundaries, with each regime having its own causal structure. Inferring causal dependencies and regime shifts is critical for analyzing the underlying processes. However, causal structure learning in this setting is challenging due to (1) non-stationarity, i.e., each regime can have its own causal graph and mixing function, and (2) complex noise distributions, which may be non-Gaussian or heteroscedastic. Existing causal discovery approaches cannot address these challenges, since generally assume stationarity or Gaussian noise with constant variance. Hence, we introduce FANTOM, a unified framework for causal discovery that handles non-stationary processes along with non-Gaussian and heteroscedastic noises. FANTOM simultaneously infers the number of regimes and their corresponding indices and learns each regime’s Directed Acyclic Graph. It uses a Bayesian Expectation Maximization algorithm that maximizes the evidence lower bound of the data log-likelihood. On the theoretical side, we prove, under mild assumptions, that temporal heteroscedastic causal models, introduced in FANTOM's formulation, are identifiable in both stationary and non-stationary settings. In addition, extensive experiments on synthetic and real data show that FANTOM outperforms existing methods.

Flow based approach for Dynamic Temporal Causal models with non-Gaussian or Heteroscedastic Noises
Abdellah Rahmani
Pascal Frossard LTS4, EPFL
LTS4, EPFL
Lausanne, Switzerland abdellah. rahmani wew.nieeellf..o.. ch
Lausanne, Switzerland pascal www.ssrddeprolll.cooooooo ch
Abstract
Understanding causal relationships in multivariate time series crucial many scenarios, such as those dealing with financial neurological data. Many such time series exhibit multiple regimes, i.e., consecutive temporal segments with priori unknown boundaries, with each regime having its own causal structure. Inferring causal dependencies and regime shifts is critical for analyzing the under- lying processes. However, causal structure learning in this setting is challenging due to (1) non-stationarity. i.e., each regime can have its own causal graph and mixing function, and (2) complex noise distributions, which may be non-Gaussian or heteroscedastic. Existing causal discovery approaches cannot address these challenges, since generally assume stationarity or Gaussian noise with constant variance. Hence, we introduce FANTOM, a unified framework for causal discovery that handles non-stationary processes along with non-Gaussian and heteroscedastic noises. FANTOM simultaneously infers the number of regimes and their cor- responding indices and learns each regime's 's Directed Acyelic Graph. It uses a Bayesian Expectation Maximization algorithm that maximizes the evidence lower bound of the data log- likelihoooo. On the theoretical side, we prove, under mild assumptions, that temporal heteroscedastic causal models, introduced in FAN- TOM's formulation, are identifiable in both stationary and non-stationary settings. In addition, extensive experiments synthetic and real data show that FANTOM outperforms existing methods.
Introduction
Causal structure learning from multivariate time series (MTS) is fundamental problem with diverse applications in traffic modeling [9], biology [44], climate science [43], or healthcare [47]. However, identifying causal relationships MTS poses several challenges. First, real-world time series are often non-stationary, exhibiting multiple unknown regimes, each potentially governed by different causal relationships. Examples include changing dependencies across climate conditions [27], financial markets [21]. and epileptic seizure stages [56]. Second, many MTS display complex noises, e.g., non-Gaussian or even heteroscedastic noise. whose variance depends on both instantaneous and lagged causes. This occurs in fMRI data [46], EEG measurements [22], or financial data [18].
Recent causal discovery methods capture linear and nonlinear interactions with instantaneous and lagged effects 35, 51]. More recently, Gong et al. [15] and Wang et al. (55) explored structural equation models for single stationary regime governed by one causal graph with historically dependent noise, where noise variance depends solely on time-lagged variables, neglecting heteroscedasticity. Existing multi-regime methods include RPCMCI [45], which identifies only linear, time-lagged interactions and requires prior knowledge of regime numbers and transitions: and CD-NOD [21], which handles causal discovery from non-stationary MTS, but is limited to homoscedastic noise, cannot infer individual causal graph for each regime, and is incapable of
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

FANTOM
Reginn
Regime
Cuical guich ind
Cunsal guphi g
E-step: Regime learning
M-stop: Craph learning
G
Posterior prohal probabilitis sla to r/s,, Rcll
Normalizing Flow
Nu
for ss n
Temporal GNN
Figure I: Illustration of FANTOM processing a MTS with two ground truth regimes (K = 2). The algorithm recovers the regime indices Σī and Σद and learns temporal causal graph for each regime (dashed edges represent time-lagged links; solid arrows indicate instantaneous links). In the E-step, posterior probabilities p The M-step then infers causal graphs within each regime. Here, Now is the number of regimes that converges K- 2.
Ft F=()) are estimated, where 2 = means at, belongs regime
identifying recurring regimes. CASTOR [39], infers both regime indices and separate causal graphs regime, accommodating instantaneous and lagged causal relationships, yet it is still restricted to normal noise. Consequently, the previous methods cannot jointly infer the number of regimes, their indices, their causal graphs, nor effectively manage either Gaussan or heteroscedastic noises.
To address these limitations, we propose FANTOM. new framework for Structural Equation Models (SEMs) in multi-regime MTS under either non-Gaussian or heteroscedastic noises. FANTOM is, to the best of our knowledge, the first method to handle heteroscedasticity in both stationary and non-stationary MTS, as well as non-Gaussianity in non-stationary settings. Given a MTS with multiple regimes, FANTOM simultaneously learns each regime's causal graph and mixing function, determines the number of regimes, and infers their indices (Figure 1). It uses Bayesian Expectation Maximization (BEM) [12] procedure to optimize the evidence lower bound (ELBO), alternatively assigning regime indices (Expectation step) and inferring causal relationships in each regime (Maximization step). Unlike Bausssan-bassed approaches, FANTOM employs conditional normalizing flows [12] to handle complex distributions and compute regime membership probabilities. It uses Bayesian structure learning that averages across all plausible graphs and naturally filters out spurious edges. Under mild assumptions, prove that temporal heteroscedastic causal models identifiable for both stationary and multi-regime MTS. Across extensive comparisons with existing multi-regime causal discovery methods, show that FANTOM consistently achieves superior performance in structure learning and regime detection. Moreover, it outperforms stationary models, even when they are provided with ground-truth regime partitions, synthetic and real-world datasets. The main contributions of this paper can be summarized as follows:
We introduce FANTOM, unified framework for causal discovery in multi-regime MTS that handles both homoscedastic non- Gaussian and heteroscedastic noises while simultaneously discovering the number of regimes, their indices, and their corresponding causal graphs.
Under mild assumptions in causal discovery, we prove identifiability of the temporal het- eroscedastic causal models in the stationary case and show that the number of regimes, their indices, and their graphs are identifiable (up to permutation) in the non-stationary setting.
We demonstrate, via extensive experiments, that FANTOM outperforms state-of- state-of-the-ar methods both synthetic and real-world datasets.
Related work. Many works tackle causal structure learning from stationary MTS, Granger causality is the primary approach used for this purpose [32, 8]. However, it unable to accommodate

---

instantaneous effects. DYNOTEARS [35], learns instantaneous and time lagged structures and leverages the acyclicity constraint, introduced by Zheng et al. in [60], to turn the DAG learning problem to a purely continuous optimization problem. However, DYNOTEARS is limited to linear SEMs. Runge et al. [42] proposed a two-stage algorithm PCMCI+ that can scale to large time series, PCMCI+ is able also to handle non linear relationships. Neverthless, DYNOTEARS and PCMCI+ are restricted to homoscedastic noises where variance is a constant over time. For this reason, Rhino [15] and SCOTCH [55] introduced models that tackle stationary MTS with historical noise, where the noise variance is function of solely time lagged parents. However Rhino, DYNOTEARS and PCMCI+, cannot handle heteroscedastic noise and are limited to stationary MTS.
Several studies have sought to tackle the challenge of causal discovery in non-stationary MTS [21. 17, 45, 33]. Remarkably, Huang et al. [21] address the setting of time series composed of different regimes by modulating causal relationships through a regime index. CD-NOD detects change points and outputs single summary causal graph, but it overlooks recurring regimes and provides neither regime-specific graphs nor their count. RPCMCI [45] provides graph for each regime, yet it assumes prior knowledge of the number of regimes, restricts edges to time-lagged links, and offers no identifiability guarantees. Balsells-Rodas et al. [5] establish identifiability for first-order, regime-dependent causal discovery in multi-regime MTS with Gaussian noise and offer practical algorithm, but their framework allows only a single time-lagged edge and excludes instantaneous links. Finally, CASTOR [39] jointly infers regime labels and their causal graphs, capturing both instantancous and lagged links, under an equivariant Gaussian noise assumption. However, none of these models can simultancously learn the number of regimes, their indices, and their structures under non-Gaussian noise, nor do they offer identifiability guarantees. FANTOM fills this gap and even generalises richer heteroscedastic noise settings while providing identifiability results. More detailed related work is provided in Appendix A).
Problem formulation
In this Section, we introduce our notation, define a temporal causal graph, and then we present a new Structural Equation Models (SEMs) for multi-regime MTS with aa.GGausiiinn/eteoooddssssiiiii noise.
Notation. Matrices, vectors, and scalars are denoted by uppercase bold G,, lowercase bold x, and lowercase letters zi_,, respectively. Ground-truth variables are indicated with an asterisk, such as G We assume that all distributions have densities p (xx)) w.r.t. the Lebesgue measure. The notation [|0 L|] represents the set of integers {0, Y0,.... L} and |- denotes set cardinality. The notation (ET) = (x))sev represents a MTS of V| = d components and length [T], T is the time index set and x<< refers to (æI-L, .... 1.
Definition 2.1 (Temporal Causal Graph [39]). The temporal causal graph, associated with the MTS (#:))eT, is defined by a Directed Acyclic Graph (DAG) g = (V, E). represented by a collection of adjacency matrices G] {Go, ... Gr), and fixed maximum lag L. Its vertices V consist of the set of components z....... x for each t' ||t are defined as follows: Vr E [[1 : L]] the variables ri_, and x are connected by lag-specific directed link z + x in G pointing forward in time if and only if x' at time t - causes x' at time t. Then the coefficient [Gy associated with the adjacency matrix G, E (I will be non-zero and z² (a t); the lagged parents of node i in G. For instantaneous links (r = 0), we cannot have self loops i.e. + j. If r O, we have an edge mt -> x! and xf E Pag (f); the instantancous parents of j at the current time t, if and only if x² at time causes x at time t.
: t|]. The edges E of the graph
In many real world scenarios, a non-stationary MTS (æ:)IET may exhibit K distinct, non-overlapping
regimes, where non-stationarity is not modeled by continuous changes but rather by sequences of piecewise-constant regimes, as in climate science [27], finance [18], or epileptic recordings [56]. Every regime r is stationary MTS block, has its own temporal causal graph G" (Detinition 2.1). At each time = 1, 2, . [T] there is a discrete latent state 2, E {1, K} partition, i.c., 2 set of all time indices at which regime appears. We gather these sets into I = (I))=e(1:KT yielding a unique time partition of the MTS (ax))ee JeeT composed of K regimes. In addition, the observation Te follows, novel and general SEM that takes into account non stationarity, and handles both
} that models the regime means that x; belongs to regime r and we denote I, [f| = Vt T the
3

---

non-Gaussian case and heteroscedastic setting, that we introduce as follows, Vr E [|1: KI], Vt E I,:
-" (Pag- (< (<)), Page (1) + : (Pag- (), Pag (1) "6"
(I)
where for and gh are general differentiable functions, with gf strietly positive, and 6." following an arbitrary probability density. We assume E(()") 0 and (() = 1 without loss of generality. We denote the set of these temporal causal graphs as 9 = (g" >re]]l In the case of non-stationary MTS, regimes appear sequentially with at least minimum duration Ç. and a subsequent regime v (where u = r +++ 1) begins only after at least ( time units from the start of regime r. Additionally, if regime r reoccurs, its duration in the second appearance is also no less than Ç samples. We refer to the phenomenon, where each regime persists for at least ( consecutive time steps, as the regime persistent dynamics assumption and define it as follows:
Assumption 2.2. We consider MTS with K multiple regimes (x2)) where the SEM is defined in (q) Given variable i, we assume that regime r E [|| K]] is (-persistent if the parents (Pag. (< (), Pag- (1) and functional dependencies (f"r, ghh :6 .e)") are stationary for consecutive time steps I.
The persistence assumption permits to capture different regime dynamics, whether arising from changes in the causal graph across regimes, commonly observed in climate science [27] and epileptic recordings [56], or from shifts in functional dependencies, which correspond soft interventions in causal discovery. Our newly introduced SEM in (1.(1) generalizes several existing approaches in three novel aspects: (1) when K = 1, if (Pagi (< 1), Pagi (4)) = 1 for all i E [[1 : d we recover the classical additive noise models in causal discovery. Thus, allowing g²1 to be a strictly positive and differentiable function, not only extends Rhino's SEM [15] but also yields a new, general SEM for stationary multivariate time series with heteroscedastic noise, (2) When K > 1, if gB "Pag (< t), Pagr (t)) = 1. and f.r ~ (0, σ,) for all i E [|| d[], we recover the setting introduced in [39, 4]. Then, allowing e." to follow an arbitrary probability density then yields the first SEM for non-stationary MTS composed of multiple regimes with non-Gaussian noise. Finally, (3) permitting gن be a strictly positive, differentiable function leads, to the best of our knowledge, to the first general SEM for non-stationary MTS with heteroscedastic noise.
3 Flow based approach for Dynamic Temporal Causal models
3.1 ELBO formulation
Many real-world time series e.g., EEG data [22, 40]. climate data [27, 16], and financial data [18] are non-stationary and exhibit complex noise distributions. Existing causal discovery methods cannot jointly recover the number of regimes, their indices, and their structures under either non- non-Gaussian or heteroscedastic noises. With FANTOM, our objective to close this gap. simultaneously learning the number of regimes K, their indices I (re)r 1:K]], and DAGS g (9") [[]]k]]] in both homoscedastic non-Gaussian and general heteroscedastic settings. Because integrating over latent regimes makes the data log-likelihood intractable, we instead maximise its Evidence Lower BOund (ELBO). Proposition 3. proved in Appendix G, formalises this ELBO for Nw > K provisional regimes. Section 3.2 outlines the initialization trick that instantiates these No, regimes, and the E-step that progressively merges them until N₂ settles at K.
Proposition 3.1. Let (#t)ee be a MTS composed of multiple regimes and following the SEM described in (1). The data likelihood admits the following evidence lower bound (ELBO):
log pe (HIET) 2 Σ (g) (p((ele [log pas (xt æ2t, G")) + log (2+)
No.
H (plz) De, xxd))l +
Σ
E Egge (0r) [log (gg) + H (gor (g")) =E)OO
where 4I E TT:T:: 2, E [[1 Nwl] are the discrete latent variables and N., is the number of regimes.

---

Here, è are all the learnable parameters of our model (detailed in Section 3.2), log po=: (at Iccs G represents the observational log-likelihood of 2t belonging to regime 2, [[1, No, qo (g") represents the variational distribution that approximates the intractable posterior pa- (9" FIEI, and p(ze [2t, 2+ represents the posterior distribution of the latent variables 2₁. The distribution p(g") is the graph prior and (2+) represents our prior belief about the membership of samples to the causal models; typically we model it as a time varying function.
3.2 Model parametrization
In this Section, we will define and motivate FANTOM design choices. FANTOM maximises the ELBO (2). with Bayesian Expectation Maximization (BEM) scheme. Because BEM normally needs the number of regimes a priori, we first describe an initialisation trick that removes this require- ment. Next, we motivate our prior over the latent regime indicator 2, and show how Temporal Graph Neural Network, combined with CNFs, models the
likelihood pa*. (æt (
Initialization trick. FANTOM initially divides the MTS into N > K equal time windows in the initialisation step (the length of the initialized windows is greater than ( minimum regime duration), where each window represents one initial regime estimate. Our initialization scheme builds some initial pure regimes (regimes composed of samples from the same ground truth regime) and other impure ones (regimes composed of samples from two neighboring ground truth regimes). After initialization, FANTOM alternates between two phases. The E-step (subsection 3.4) updates the regime indices I = (Ir)r E []1 Nw and removes regimes with too few samples, updating N The M-step (subsection 3.3.1) learns the graphs (Gr)r E [|I Nw/] and models heteroscedastic noise with Bayesian structure-leaming method that uses conditional normalizing flows (CNFs). These phases repeat until the algorithm reaches the maximum number of iterations.
BEM choice motivation. We argue that inferring regimes and learning their associated DAGs are interdependent tasks, making the BEM algorithm particularly well-suited for this problem. A two-step alternative, that detects change points with KCP [1]. then runs causal discovery on each segment breaks down: First, change point detection methods like KCP [1] fail detect regime shifts driven by changes in causal mechanisms, because those involve shifts in conditional distributions (See Appendix F.3.3). It also treats recurring regimes as distinct, forcing redundant model fits and raising computation costs. Second, heteroscedastic noise further degrades existing causal methods (Table 2). FANTOM addresses all three issues by coupling CNFs, which capture heteroscedasticity, with Bayesian structure learning that prunes spurious edges.
Time varying weight modeling. We use time-varying weights, initially proposed for financial data modeling [59, 57], as priors for the discrete latent variables zt. To support smooth regime transitions consistent with our persistence assumption (Assumption 2.2), we adopt a flexible functional form based on the softmax transformation of learnable parameter w" R2 and time index t:
P (2t = Te (ew"") Σ>=, exp(w]- " This formulation encourages that, if x; belongs to regime in the current iteration, it is only allowed to remain the same regime r or smoothly transition to neighboring regimes (r 1, + 1) in the next iteration. See Section 3.4 and Appendix B for details.
exp(=] (t+)))
Bayesian structure learning. Following [14, 15, 55J, FANTOM employs Bayesian structure learning. We approximate the intractable posterior (G (IET) using the variational distribution (G) I (g") 8(0"), where denotes the Dirac delta function. Following [14, 15, 55]. we model I (9") as a product of independent Bernoulli variables and compute its expectation with a single Monte Carlo sample using the Gumbel-Softmax trick [26]. Additional details are in Appendix C.
Likelihood of SEM. Using the functional form Eq.(1), we have
yt for = 9" (r (Pag- (< 1), Pa (1) of" x - f (Pag- (< 1), Pay- (1)),
then we can write the observational likelihood:
(3) where Phetern refers to the density function of the heteroscedastic conditions. To estimate fir we build upon the model formulation of [14], which uses neural networks to describe the functional
pa- (x)) Pag- (< t), Pa (t)) Phetero
(vi" Pag- (< t). Pag
5

---

relationship for R + R. Specifically, we propose flexible functional designs for fo, which must respect the relations encapsulated in g". Namely, if x & Pa (< t) Pag (t), then Of:r 10xl 0. We design
for (Pag- f), Pa) - 60"
ΣΣ r=
G
O'
e
c
(4)
where 1/" and O" are neural network blocks illustrated in Figure 2 with all the other colored blocks. Instead of using a neural network block per node, we adopt weight-sharing mechanism by using trainable em- beddings 6r,, for E [|0 : Li] and i E (1,... d). For the heterosçedastic term, we introduce an invert- ible mapping fg
Neural Network blok
Layer
IR + R such that:
Concnt
Pot (9"" (Pag- 1), Pag- (1)) -4" -+""==""" =
Lineur Roile
and
Layer Norm
Zey
where ni (), 1). The design of th."" needs properly balance the flexibility and tractability of the transformed noise density. We choose conditional normalizing flows, for heteroscedastic noise, called conditional spline flow [10], that transforms our het- eroscedastic noise distribution to a fixed normal noise n,"" for regime r. The spline parameters are predicted using a hyper-network with similar form (q. to incorporate heteroscedasticity. Due the invert- ibility of pi.". parents is:
Linenr Leaky ReLa
Neural Network block
Nenral Network block
Conceul
Inputs
Figure 2: Temporal graph neural network
(TGNN) used by FANTOM.
"rr the noise likelihood conditioned on all
Phetero
(0) ir Pag- Pag-((())) = Pater
de
(5)
Be
where (-) the standard normal density. In the non-Gaussian, non-heteroscedaasie case, FAN- TOM learns base distribution based on composite affine-spline transformation of standard Gaussian. Finally, the system parameters comprise the learnable parameters of the time varying weights w", of the variational inference o" and of the neural networks 6". We use θ, introduced in Eq.(2), to note all the learnable parameters of FANTOM, and we have the set of parameters is: ((w", or 0r) No.
(-)
3.3 BEM algorithm
3.3.1
M step: graph learning
FANTOM applies BEM to maximize the ELBO described in Proposition 2. It begins by initializing the regime partitions 57 = p(z r|æ,, xx(e) using equally sized windows, selected via hyperparameter. Note that these binary regime indices 3% are updated during the E-step. Then, in the M-step, FANTOM incorporates 8r the ELBO (2) to estimate the DAGS for each regime and learn the parameters w" that align ne (w=" with B3. We have:
No. 3r log pa, (we, wee,
N, E gn (9r) [logp (g") + H (go. (G"))
arg max 2
t=1
IT No.
of log ((w)))
(6)
The maximization the above equation can be decomposed into two distinct and separate maximiza- tion problems. The first problem, regime alignment, focuses aligning Tr(w") with 3r, While the

---

second one, graph learning, involves estimating DAGs for every regime.
For the graph prior (9 for all r []1 N have to combine two components: DAG constraint and graph sparseness prior. Inspired by [60, 15, 14], we propose the following un- normalised prior p (g") x expl (->a G. phh (G)) - ar (Gi) Using E3.(3), we have
log (0t, | C<,, 9-) = Σ
(vt" Pag- (< t), Pa)-(() where Phetero and 1/2" are defined in Eq.(5). The parameters H", o² are learned by maximizing the graph learning maximiza- tion problem (6), where the Gumbel-softmax gradient estimator is used [26]. We also leverage augmented Lagrangian training similar to [35, 39, 14], to anneal a, p.
i-l log Phetero
3.4 E step: Regime learning
In the E-step, FANTOM updates the posterior probability B = p(2z = 2e. xee) (see Eq.(7), with derivation provided in Appendix B):
B (2+ = r) Ict,2 - g*)
D = p)p (x: Tct, J. Go)
X (w" ) (ot &et, 9"'
where p(æ;, Ezts r,) denotes the observational likelihood of being generated by the SEM from Eq (1) for regime т. This probability is computed using the CNFs trained during the M-step, following the same reasoning as in (3).
The probability of x belonging to regime is influ- enced by two main factors: the observation's 's position within its current regime and whether that regime is designated as pure or impure. For example, if æ; is in a pure regime but near the boundary in the current iteration, ((("") and [1("+] are nearly
1.0
0.8
ne(w")
equal (e.g.. #se[1100, 1500 (wo*) vs. Tre[1100,1500 (600)(w> (w²)
nitw?)
in Figure 3). Nonetheless, since regime was learned from pure data, p(æ; Vett &t = r, 9t) stays high, keeping 8 its maximum value and maintaining æ, in regime r for the next iteration. In the other hand, if æ, is in an impure regime + 1 near the boundary during the current iteration, Te(w") and Ty (w"+1) also close in value (e.g., (te[1501,1800) (w ) T1e01150.18000(44 (w²) in Figure 3). However, because the causal graph for regime r is more reliable (hav- ing been derived from pure data). p(æ;, &F 2t
0.4
m((w))
0.2
0.0
1000
Time
Figure 3:
Illustration of n. (w") after Fan-
tom's 's first iteration with equal windows of 1500 samples for an MTS of 4500 samples with two ground-truth regimes: Ií [|0 1999 and Iž = [| |000 4500|].
(") > p(æ; æ<t; 24 + 1. G )). 2, moves from regime +
). As result, in the next iteration.
For simplicity reasons, we explicit these cases from one border but the same thing happens in the other border which accelerates convergence. More details about other cases and Figures that illustrate the idea could be found in Appendix B.
After updating P. for each sample æ,, FANTOM assigns value of 1 the most probable regime (with the highest (,), and 0 to others. Additionally, FANTOM filters out regimes with insufficient samples (fewer than Ç, the minimum regime duration, defined as hyper-parameter). Discarded regime samples are then reassigned to the nearest regime in terms of probability 0 in the subsequent iteration which in general neighboring regime ensured by the way we set up the probability 37 x Te (w("")) (p'")((( (xx) T<t, 2t r").
4 Identifiability results
Identifiability an important statistical property to ensure that the causal discovery problem is meaningful. In causal analysis, the whole point is to find out which variable causes others; if the model is not identifiable, the analysis is not possible at all. This section proves that causal discovery from multi-regime MTS is identifiable in the FANTOM framework namely, when (i) the noise is
7

---

non-Gaussian or heteroscedastic and (ii) the latent variable has time- varying-parem parent prior. We first formalize identifiability for this setting, then state three theorems covering both stationary and multi-regime MTS under the two noise assumptions.
Definition 4.1. The conditional distribution of regli- MTS with a time varying prior is: p(a,(æx<)) Σ (w"")"" (2, 9*). We say this model is identifiable up to permutation and translation, if:
.rr E [[I K|] the causal model (0", ("""")) is identifiable.
For any two models with parameters (w", 00" "r)))== and (ō", ("))"") such that for any E T:
p(a, (et))) (t)), we have K = K and it exists a permutation a and translation function 0 2 such that fr = (r())) and w" = e((0(r))
Following [15, 39, 7]. we use the common assumptions of causal discovery settings (Causal Markov property H.2, stationarity H.2, minimality H.2, sufficiency H.2), see Appendix H.1 for precise statements. We present our first theoretical results, Theorem 4.2, states that for any stationary MTS, composed of K = 1 regime and following (1.(1) with ~ NN( 1) the ground truth solution G is uniquely identifiable, the detailed proof can be found in Appendix H.
Theorem 4.2 (Identifiability Temporal Heteroscedastic Gaussian noise model (THGNM)). Assume Causal Markov property, stationarity, minimality, sufficiency and let (æI)IET be MTS following THGNM, E TT:
t -' (Pab(<<)) Pa)(()) + (Pag(< 1), Pa((()) G.
(8)
where f and g¹ are differentiable functions, with 9i strictly positive and et ~ (0, 1) are mutually independent normal noises. The THGNM is identifiable if 2 is not a polynomial of degree two.
Identifiability of Temporal Restricted Heteroscedastic noise model (TRHNM), We present and prove in Appendix H.3 our second identifiability results of a TRHNM (Theorem H. where ef follow any arbitrary density distribution. We states the results for bivariate time series, in which we show that, a backward model exists, a differential equation will always hold. Then inspired
from Peters et al. [38], Immer al. [25], and Strobl et al. [50]. we define TRHNM and show its identifiability.
Our last theoretical result states that the mixture of identifiable temporal causal models with either non-Gaussian or heteroscedastic noises is also identifiable as defined in definition 4.1.
Theorem 4.3 (Identifiability of the mixture of identifiable temporal causal models). Let F be family of K identifiable temporal causal models, F = (per (-I-, (1)))))) that are linearly independent and let Mr be the family of all -finite mixtures of elements from F, i.e.,
EK (aol:w"'o-(eeee (x,
E EF :e(w"))))) ) 0 and EKI ne(w") = =1}
Then the family MK is identifiable as defined in definition 4.1.
Identifiability is important in causal discovery as shown in several papers [7, 15, 35, 5]. We extend these guarantces to FANTOM's s settings by proving identifiability for stationary temporal models with heteroscedastie noise, showing that weaker assumptions suffice when the noise ef is simplified, recovering identifiability under explicit restrictions for arbitrary noise, and covering multi-regime MTS scenarios in both cases. Although convergence rates or finite finite-sample bounds remain elusive because the BEM objective is non-convex, our experiments demonstrate that FANTOM still converges in non-Gaussian and heteroscedastic settings. Further convergence rates or finite-data bounds are however extremely challenging due to the non-convexity of the acyclicity constraint in BEM procedure. Yet, we empirically demonstrate, in the experiments section, that FANTOM converges in both non- Gansaia and heteroscedastic cases.
8

---

Table I: Average SHD, Fl scores, NHD and Ratio for different models with d 10 nodes and K regimes. Split denotes whether regime separation is automatic (v) or manual (x), and Type classifies the graph as window (W) or summary (S). Inst. refers to instantaneous links, and Lag to time-lagged
edges.
Hoummeddasii nou- Caamoiun pose
Hecroscolaaiic muiic
Model
Split Type|
Regime
SHDI
Fif
NIIDI
SHHD.
Rati
SHLD
NEED
Ratio:
[SID]
F17
NHID
Ranoo
ROMCH+
17.
74.9
01.12
14.3
88.3
0:01
0.11
46.1
11.1
0005
0.35
66.0
0.00
0.80
DYNOTEARS
W
4200
54.4
0.06
82.1
0.00
0.12
88.5
31.5
0.14
0.68
118.0
37.5
0.17
0,6
CASTOR
W
22.0
66.2
0.03
0.35
17.0
84.4
0.01
0.15
104.0
25.4
0.19
0.76
133.5
34.8
0.24
0.64
RICMET
W
CASTOR
W
47.0
34.8
0.05
0.65
90.5
9944
0.07
0.40
51A
BANTOMI
0.33
99.5
89.1
0.01
0.10
5,67
93.3
0,006
0.06
12.3
90.9
0.012
0.08
SHD]
F11
NHD]
Rato.
Asc:
SHD)
FT7
NHD
Rizao )
Acc.
COD
42.5
31
0.42
1.67
42
6.15
061
0.93
PANTOM
4.5
95.6
0.04
0.84
96.6
4.5
96.5
0084
0.03
96.8
5 Experiments
5.1 Synthetic data
Data generation. We conduct extensive experiments to evaluate FANTOM's performance on synthetic datasets. For ground truth graph generation, use the Barabási- Alberr model (degree 4) for instantancous links and the Erdos Erdos-Rényi model (degree 1-2) for time-lagged relationships. For data generation process, J Df are chosen to be randomly initialized MLPs with one hidden layer and activation functions randomly chosen from (Tanh, Exp}. We evaluate the different models on multiple complex noise distributions; non-stationary MTS with either (1) heteroscedastic noise or (2) non non-Gaussian homoscedastic noise, details in Appendix 1. We consider L= 1, while additional experiments with multiple lags are provided in the Appendix F.1. Regime durations are randomly selected from {1000, 1500, 2000, 2500}. We test different numbers of nodes {5, 10, 20, 40) and varying regime counts (K E 3)) Each combination of K and d nodes repeated three times, resulting in over 24 distinct datasets!
Benchmarks. We benchmark our model against several baselines, including causal discovery methods for MTS with multiple regimes, such as CASTOR [39], CD-NOD [21] and RPCMCI [45]. Since CD-NOD returns summary graph (see Appendix F.1), compute a comparable summary graph from FANTOM's output for fair evaluation. FANTOM is also compared with models for single-regime MTS, including Rhino [15], PCMCI+ [42], DYNOTEARS [35]. Given that these models cannol deal with multiple regimes, we put them in more favorable position than ours and provide these models with the true regime partition information. This is done by training the aforementioned models on each pure regime separately (regime governed by the same causal model).
Evaluation Metrics. We assess the performance of our proposed method for learning the DAGs using four key metrics: 1) Fl score, representing the harmonic mean of precision and recall; 2) Structural Hamming Distance (SHD), which counts discrepancies (e.g.. reversed, missing, redundant edges) between two DAGs; 3) Normalized Hamming Distance (NHD) measures how many edges differ normalized by the total number of possible edges; 4) Ratio NHD computes the ratio between the NHD and the baseline NHD of an output with the same number of edges but with all of them incorrect. For the regime detection task, we use Accuracy (Reg Acc) metric.
Results and discussion. Table shows results on MTS with multiple regimes under heteroscedastic noise (right part of the table) and homoscedastic non-Gaussian noise (left part of the table). In the homoscedastic, non-Gaussian scenario, baselines generally perform better in the graph learning task, yet FANTOM still surpasses them regime detection, with 99.4% accuracy, and instantancous links, 99.5% F1. For the regime detection task, CASTOR succeeds to detect the regimes but with low accuracy (51.6%) compared to FANTOM (99.4%) and this due the fact that CASTOR assumes equivariance normal noise with, however RPCMCI does not converge in this case too. For time- lagged connections, Rhino (95.2% FI). that has access the ground-truth regime labels by training it on pure regime separately, slightly outperforms FANTOM (89. 19.1% 1% F1) that learns simultaneously the number of regimes, their indices and structures. In the heteroscedastic setting, FANTOM consistently outperforms both multi-regime baselines (CASTOR, RPCMCI, CD-NOD) and stationary approaches. It achieves the top scores on all metries: for instantaneous links, an Fl of 93.3%, 60%
hntpss//eithbb..ommaaaaaaa 1/fantom.ii
9

---

improvement over the second-best, and ratio of 0.06, 0.64 lower than the next-best DYNOTEARS. For time-lagged links, an FI of 90.9%, 25% higher than Rhino, and a ratio of 0.08. FANTOM also detects the correct number of regimes and their indices with over 96% accuracy. By contrast, RPCMCI struggles to converge due to its homoscedastic assumption and time-lag-only dependencies, CASTOR relies on Gaussian noise and cannot detect regime labels in the absence of ground truth, DYNOTEARS, PCMCI+, CD-NOD, and Rhino likewise fail in this heteroscedastic scenario, even when regime labels are given priori. Although Rhino models history-dependent noise, it does not handle general heteroscedasticity. Overall, the table shows that FANTOM is the only method that remains robust across both noise scenarios, matching or surpassing specialised baselines while simultaneously discovering regimes and graphs.
Appendices F.2.1 and F.2.2 report additional results with standard deviations, confirming that FAN- TOM sustains its performance when scaled to graphs of 20 and 40 nodes. Ablation study in the Appendix F.1 further shows FANTOM's 's robustness towards the choice of initialized window and ç.
5.2 Real world data
Wind Tunnel. We use the wind tunnel datasets from Gamella et al. [13], featuring two controllable fans pushing air through a chamber, barometers measuring air pres- sure at various locations, and hatch con- trolling an external opening. The dataset comprises 16 variables across two regimes of 10,000 samples each: the first is obser- vational, while the second involves soft in- terventions on five variables (see Appendix E.4.1). We compare FANTOM to the afore- mentioned baselines, with results in Table 2. FANTOM is the only model that detects the regime with 99.9% accuracy and outperforms all baselines the causal graph learning task, achieving 38.5% on Fl score. Notably, FANTOM surpasses all baselines even when they are given the ground-truth regime partitions.
Table 2: Performance on Wind Tunnel data evaluated summary causal graph.
Split
SIDD
FI1
Ratio
Reg Acc.
PCMCI+ DYNOTEARS CASTOR
X
37
22.9
0.77
X
X
34
0
x
X
104
17.2
0.82
X
CD-NOD
20.0
0.80
X
Rhino
X
47
32.0
0.68
x
CASTOR FANTOM
V
120
19.5
0.80
49.9
38.5
0,61
99.9
Epilepsy detection. Huizenga et al. [22] show that scalp potential fields are contaminated by heteroscedastic noise in EEG measurements. We evaluate FANTOM's performance in detecting epileptic regimes using EEG signals from 10 different patients the Temple University Hospital EEG Seizure Corpus (TUSZ) dataset [40, 52]. We treat this an unsupervised regime detection problem, analyzing roughly 100 seconds of recordings at a 250 Hz sampling rate for each patient, capturing both normal and seizure states (see Appendix E.4.2). The recordings consist of 19 electrodes, each considered a causal variable. FANTOM detects the correct regime partitions with an average 82.7% accuracy across all patients. The seizure regime's learned graph is denser and more connected than that the normal state, which aligns with the generalized seizures affecting multiple brain regions. Full details and illustrations are provided in Appendix E.4.2.
6 Conclusion
We introduced FANTOM, a unified framework for multi-regime MTS that jointly infers (i) the number of regimes, (ii) their boundaries, and (iii) their corresponding causal DAG, under either non-Gaussian or heteroscedastic noises. Under mild assumptions in causal discovery, we prove identifiability of the temporal heteroscedastic causal models in the stationary case and show that the number of regimes, their indices, and their graphs are identifiable (up to permutation) in the non-stationary setting. Extensive experiments on synthetic and real-world data show consistent gains over strong baselines. FANTOM offers principled means to uncover regime-specific causal dynamics, enhancing regime detection, and causal discovery with potential applications in various domains such as finance, climate science, and neuroscience.
10

---

Acknowledgment
We thank Alessandro Favero, Nikolaos Dimitriadis, and Ortal Senouf for helpful feedbacks and comments. This work was supported by the SNSF Sinergia project 'PEDESITE: Personalized Detection of Epileptic Seizure the Internet Things (IoT) Era'
References
[1] Sylvain Arlot, Alain Celisse, and Zaid Harchaoui. A kernel multiple change- algorithm via model selection. Journal of machine learning research, 2019.
[2] Charles K Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. Journal of Artificial Intelligence Research, 2022.
[3] Karim Assaad, Emilie Devijver, Eric Gaussier, and Ali Ait-Bachir. mixed noise and constraint- based approach to causal inference in time series. In ECML PKDD, 2021.
[4] Carles Malsells-Rodas, Yixin Wang, and Yingzhen Li. On the identifiability of switching dynamical systems. arXiv preprint (xXix:2305.15922 2023.
[5] Carles Baisells-Rodas, Yixin Wang, and Yingzhen Li. On the identifiability of switching dynamical systems. In Forty-first International Conference on Machine Learning, 2024.
[6] Albert-László Barabási and Réka Albert. Emergence of scaling random networks. Science,
1999.
[7] 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, 2020.
[8] Bart Bussmann, Jannes Nys, and Steven Latré. Neural additive vector autoregression models for causal discovery in time series. In Discovery Science, 2021.
[9] Yuxiao Cheng, Ziqian Wang, Tingxiong Xiao, Qin Zhong, Jinli Suo, and Kunlun He. Causaltime: Realistically generated time-series for benchmarking of causal discovery. arXiv preprint arXiv:2310.01753, 2023.
[10] Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Advances neural information processing systems, 2019.
[11) Doris Entner and Patrik O Hoyer. On causal discovery from time series data using fci. Proba- bilistic graphical models, 2010.
[12] Nir Friedman. The bayesian structural em algorithm. arXiv preprint 1rXXiv1301.73:11 2013.
[13] Juan I. Gamella, Jonas Peters, and Peter Bühlmann. Causal chambers as a real-world physical testbed for ai methodology. Nature Machine Intelligence 7(1):107-118; 18, 2025.
[14] Tomas Geffner, Javier Antoran, Adam Foster, Wenbo Gong, Chao Ma, Emre Kiciman, Amit Sharma, Angus Lamb, Martin Kukla, Nick Pawlowski, et al. Deep end-to-end causal inference. arXiv preprint arXiv:2202.00 02/50 2022.
[15] Wenbo Gong. Joel Jennings, Cheng Zhang, and Nick Pawlowski. Rhino: Deep causal temporal relationship learning with history history-dependent noise. Preprint arXiv:2210.14706, 2022.
[16] Wiebke Günther, Urmi Ninad, Jonas Wahl, and Jakob Runge. Conditional independence testing with heteroskedastic data and applications to causal discovery. Advances in Neural Information Processing Systems, 2022.
[17] Wiebke Günther, Oana-luliana Popescu, Martin Rabel, Urmi Ninad, Andreas Gerhardus, and Jakob Runge Causal discovery with endogenous context variables. Advances in Neural Information Processing Systems, pages 36243-36284, 2024.

---

[18] James D Hamilton and Raul Susmel. Autoregressive conditional heteroskedasticity and changes in regime. Journal of econometrics, 1994.
[19] Uzma Hasan, Emam Hossain, and Md Osman Gani. A survey on causal discovery methods for iid and time series data. arXiv preprint arXiv:2303. 2023.
[20] Stefan Haufe, Klaus-Robert Müller, Guido Nolte, and Nicole Krämer. Sparse causal discovery in multivariate time series. In causality: objectives and assessment. PMLR, 2010.
[21] Biwei Huang, Kun Zhang, Jiji Zhang, Joseph Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from heterngenneuus/oosstatinnnn data. In Journal of Machine Learning Research, 2020.
[22] Hilde M Huizenga and Peter CM Molenaar. Equivalent source estimation of scalp potential fields contaminated by heteroscedastic and correlated noise. Brain topography, 1995.
[23] Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo. Constraint- causal discovery: Conflict resolution with answer set programming. In UAI, 2014.
[24] 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, 2010.
[25] Alexander Immer, Christoph Schultheiss, Julia E Vogt, Bernhard Schölkopf, Peter Bühlmann, and Alexander Marx. On the identifiability and estimation of causal location-scale noise models. In International Conference on Machine Learning. PMLR, 2023.
[26] Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterizztion with gumbel-softmat arXiv preprint arXiv: 16..111111011 01 144, 2016.
[27] Soufiane Karmouche, Evgenia Galytska, Jakob Runge, Gerald A Meehl, Adam S Phillips, Katja Weigel, and Veronika Eyring. Regime-orienter causal model evaluation of Illantcc-pccife teleconnections cmip6. Earth System Dynamics, 2023.
[28] Nan Rosemary Ke, Olexa Bilaniuk, Anirudh Goyal, Stefan Bauer, Hugo Larochelle, Bernhard Schölkopf. Michael C Mozer, Chris Pal, and Yoshua Bengio. Learning neural causal models 2019.
from unknown interventions. Preprint
[29] Ilyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International conference on artificial intelligence and statistics. PMLR, 2021.
[30] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[31] Lars Lorch, Jonas Rothfuss, Bernhard Schölkopf, and Andreas Krause. Dibs: Differentiable bayesian structure learning. In Advances in Neural Information Processing Systems, 2021.
[32] 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.
[33] Sarah Mameche, Lénaig Cornanguer, Urmi Ninad, and Jilles Vreeken. Spacetime: Causal discovery from non-stationary time series. arXiv preprint 1rXiv:2501.102355 2025.
[34] Mark Newman, Networks. Oxford university press, 2018.
[35] 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.
[36] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference on time series using restricted structural equation models. Advances neural information processing systems, 2013.
[37] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Flements of causal inference: founda- tions and learning algorithms. The MIT Press, 2017.
12

---

[38] Jonas Peters, Joris M Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery with continuous additive noise models. The Journal of Machine Learning Research, 2014.
[39] Abdellah Rahmani and Pascal Frossard. Causal temporal regime structure learning. In The 28th International Conference on Artificial Intelligence and Statistics, 2025.
[40] Abdellah Rahmani, Arun Venkitaraman, and Pascal Frossard. A meta-gnn approach to per- sonalized seizure detection and classification. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023.
[41] Jakob Runge. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2018.
[42] Jakob Runge.
Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Conference on Uncertainty in Artificial Intelligence. PMLR, 2020.
[43] Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science advances, 2019.
[44] Karen Sachs, Omar Perez, Dana Pe'er, Douglas A Lauffenburger, and Garry P Nolan. Causal protein-signaling networks derived from multiparameter single-cell data. Science, 2005.
[45] 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, 2020.
[46] Christof Seiler and Susan Holmes. Multivariate heteroscedasticity models for functional brain connectivity. Frontiers in neuroscience, 2017.
[47] Xinpeng Shen, Sisi Ma, Prashanthi Vemuri, and Gyorgy Simon. Challenges and opportunities with causal discovery algorithms: application to alzheimer's pathophysiology. Scientific reports, 2020.
[48] Le Song, Mladen Kolar, and Eric Xing. Time-varying dynamic bayesian networks. In Advances in Neural Information Processing Systems, 2009.
[49] Peter Spirtes, Clark Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
[50] Eric V Strobl and Thomas A Lasko. Identifying patient-specific root causes with the het- croscedastic noise model. Journal of Computational Science, 2023.
[51] Xiangyu Sun, Oliver Schulte, Guiliang Liu, and Pascal Poupart. Nts-notears: Learning nonpara- metric dbns with prior knowledge. arXiv preprint (tti::20090022880 2021.
[52] 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 electroencephaloggaahic seizure analysis. arXiv preprint arXiv:2104.08336, 2021.
[53] Sofia Triantafillou and Ioannis Tsamardinos. Constraint-based causal discovery from multiple interventions over overlapping variable sets. The Journal of Machine Learning Research, 2015.
[54] Sumanth Varambally, Yi-An Ma, and Rose Yu. Discovering mixtures of structural causal models from time series data. arXiv preprint 120.032300 2023.
[55] Benjie Wang, Joel Jennings, and Wenbo Gong. Neural structure learning with stochastic differential equations. arXiv preprint arXiv:2311.03309 2023.
[56] Xiaojia Wang, Yanchao Liu, and Chunfeng Yang. Ictal-onset localization through effective connectivity analysis based on rnn-ge with intracranial eeg signals in patients with epilepsy. Brain Informatics, 2024.
[57] Chun S Wong and Wai K Li. On a logistic mixture autoregressive model. Biometrika, 2001.
13

---

[58] Sidney Yakowitz and John D Spragins. On the identifiability of finite mixtures. The Annals of Mathematical Statistics, 1968.
[59] Shuguang Zhang, Minjing Tao, Xu-Feng Niu, and Fred Huffer. Time-varying gaussian-cauchy mixture models for financial risk management. arXiv preprint arXiv:2002.06102, 2020.
[60] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances in neural information processing systems, 31, 2018
[61] Shengyu Zhu, Ignavier Ng. and Zhitang Chen. Causal discovery with reinforcement learning. Preprint arXiv:1906.0444/ 2019.
14

---

Table of content
A Detailed related works
16
Expectation step: Derivation, intuition and illustration
17
C Maximization step:
18
C.I Mathematical Derivation: Equation (2) to (6)
18
Variational Inference Details
19
D Limitations & Risk of spurious causality
19
Data generation and Baselines
20
E.1
Synthetic data
20
E.2
Baselines
20
E.3
Optimization parameters
E,4
Real world data
E.4.1
Causal Chambers data
21
E,4.2
Epilepsy data
22
Additional Experiments
24
F.1
Ablation studies
24
F.1.1
Impact of Heteroscedastic Modeling and Flow Architecture
24
E.1.2
Robustness to the pruning step, minimum regime duration, and regime ini- tialization
25
R11.3
Robustness data standardization
26
F22
Additional results on synthetic data
27
F.2.1
Heteroscedastic noise with different number of nodes and regimes
27
F.2.2
Non-Gaussian noise with different number of nodes and regimes
29
F3
Additional experiments for L=2
E33.1
Illustrations of learned graphs
32
F.3.2
Time complexity analysis
35
F3.33
Regime detection experiments
36
G Proof of proposition
37
Proofs of our theoretical contributions
38
H.I
Assumptions
38
H.2 Proof of theorem 4.2
39
H.3 Identifiability results in the case of Temporal General Heteroscedastic Noise Models
40
H.4 Proof of theorem 4.3
42
15

---

θ
2t
x <t
Et
t t=
L,
IT
Figure 4: Graphical model of FANTOM. Observed variables (æx) are in gray, while latent variables (2²) and parameters (9) are in white. Blue edges represent parameter variiabee interactions.
A Detailed related works
Causal structure learning from IID data. Causal structure learning has become an active area of research. hasan et al. [19] recently presented a comprchensive review of causal discovery methods for IID data and time series. For IID data, several approaches rely on conditional independence to infer causal relationships from observational data, such as the classical PC algorithm [49]. Additionally, some methods extend beyond observational data, incorporating interventional data to enhance causal inference, including COmbINE [53] and HEJ [23]. These approaches utilize data collected from controlled interventions to uncover causal relationships. novel research direction introduced in [60] addresses the combinatorial challenges of structure learning by formulating it as continuous constrained optimization problem, thus avoiding computationally costly combinatorial searches. Similarly, Zhu et al. [61] utilize the acyclicity constraint but employ reinforcement learning techniques for estimating directed acyclic graphs (DAGs). In contrast, Ke et al. [28] propose an approach that learns DAGs from interventional data by optimizing unconstrained objective function. [7] provide a comprehensive analysis of continuous-constrainee applicable interventional data scenarios. Another significant method, DiBS [31], estimates the full posterior distribution over Bayesian networks from limited observations, enabling quantification of uncertainty and assessment of confidence in causal discovery.
methods, offering a generalized framework
Causal structure learning from stationary MTS. The previously mentioned state-of-the-ar methods primarily target independent observations rather than temporal dependencies. Assaad et al. [2] provide comprehensive review of approaches specifically designed for causal discovery from MTS. To model causal relationships involving time dependencies, researchers frequeetly employ Dynamic Bayesian Networks (DBNs), which effectively capture discrete-time temporal dynamics within directed graphical frameworks. Some methods neglect contemporancous (instantaneous) dependencies and focus exclusively on recovering time-lagged causal links [20, 48), and tsFCI [11], the latter adapting the Fast Causal Inference algorithm [49] for time series data. Runge et al. [43] introduced PCMCI, a scalable two-stage algorithm for time series, initially focusing only on time-lagged relationships. They subsequently extended it PCMCI+ [42], enabling the identification of contemporaneous causal connections. Additionally, models addressing non-Gaussian instantaneous effects have been developed, such VARLINGAM [24], which integrates non- Gaussian instantaneous models with autoregressive components. Another significant method Time- series Models with Independent Noise (TiMINo) [36], which studies nonlinear and instantancous effects using constrained SEMs. Pamfil et al. [35] recently proposed DYNOTEARS, leveraging an algebraic characterization of graph acyclicity from [60] to estimate both instantaneous and time- lagged relationships from time series data. DYNOTEARS utilizes a score- based DBN learning approach optimized via an augmented Lagrangian framework, enabling causal graph inference without assumptions on the underlying topology. In contrast, methods like NBCB [3], a noise- based//onstraint-basee approach, aim to learn summary causal graph directly from observational time series data, going beyond Markov equivalence constraints even in the presence of instantaneous relationships. Rhino [15] introduced the first model that tackle stationary MTS with historical noise, where they assume that the noise variance changes over time a function of solely time lagged
16

---

1.0
00
m(w+))
0.6
Zz
m(wr))
0.4
n(w"))
0.2
0.0
60
1000 2000
5002
4200
Iz Figure 6: Example of initialization with N = 3 windows, I1 and Σ3 are pure regimes while I2 is impure (composed of sam- ples from two ground-truth regimes Zí and I =)).
I I2
Time
Illustration of Te (w") after Fan- tom's first iteration with equal windows of 1500 samples for an MTS of 4500 samples with two ground-truth regimes: Ij = 110 1999[] and In [|0000: 4500
Figure 5:
parents. All the aforementioned methods does not handle heteroscedastie setting and also fail short in the case of MTS with multiple regimes.
Causal structure learning from MTS with multiple regimes. Some research have aimed to address this challenge by developing methods for causal discovery heterogeneous data. An example of such method is CD-NOD [21]. tackles time series with various regimes. By using the time stamp IDss as a surrogate variable, CD-NOD output one summary causal graph where the parents of each variable are identified as the union of all its parents in graphs from different regimes. Then it detects the change points by using a non-stationary driving force that estimates the variability of the conditional distribution p(z, (union parents of x;;) over the time index surrogate. While CD-NOD provides summary graph capturing behavioral changes across regimes, falls short in inferring individual causal graphs, also CD-NOD cannot handle either non Gaussiann heteroscedastic noise. The overall summary graph does not effectively highlight changes between regimes. Additionally, CD-NOD detects the change points but fails to determine the regime indices, rendering it incapable of inferring the precise number of regimes. In scenarios involving recurring regimes, CD-NOD is unable to detect this crucial information. Another relevant work dealing with MTS composed of multiple regimes is RPCMCI [45]. In this paper, the model [45] learns temporal graph for each regime. However, they focus initially on inferring only time-lagged relationships and require prior knowledge of the number of regimes and transitions between them. [5] addresses first-order regime-dependent causal discovery from MTS with multiple regimes. They proved that first-order Markov switching models with non-linear Gaussian transitions are identifiable up to permutations. Their work offers also a causal discovery in time series data. However, its primary
practical algorithms for limitation is the assumption of solely time-lagged relationships, with the theory being restricted to a single time lag. CASTOR [39] learns regime indices and their corresponding causal graphs, including instantaneous and time-lagged relationships, under the assumption of normally distributed noise with equivariance. However, like other causal discovery methods for non- stationatr MTS, does not address non-Gaussian or heteroscedastic noise. [54] tackle different setting in which they aim to discover a mixtures of Structural Causal Models from a datasets of MTS. They assume that they have different stationary MTS in the same dataset, one regime per MTS, and each one could be explained by one causal model in the mixture. In their case, they assume that every MTS in the dataset is stationary but the whole dataset is mixture. In our case, we assume that we have only one non-stationary MTS and it is composed of different regimes where we do not know when the regime starts and ends, and our goal is to identify the regimes and the corresponding causal graphs.
B Expectation step: Derivation, intuition and illustration
B = 2p = 1 Zt, T2L, G((oL)
p (z+ = r) p (x: Fets 21 = =. )r)
Nu ΣNHi((e = 1) (x:
Iet, = j. gi)
X (T(W")) (x) 3c = r,)",,
17

---

where p(z; Ict,
G denotes the likelihood of Is being generated by the SEM from Equation (1) for regime r. This probability is computed using the normalizing flows trained during the M-step, following the same reasoning as in Eq(3).
The probability of at belonging to regime r is influenced by two main factors: the observation's position within its current regime and whether that regime is designated as pure or impure. In order to elarify the intuition behind pure and impure regimes, Figure 6 shows an example of such case. presents an initialization of three equal windows while the MTS is composed of two ground truth regimes presented by green color Ií and red one Zz. such case, the regimes I1 and Σ3 are pure, because they are composed of samples coming from the same ground truth regime (I₁ for the regime
I and Ig for the regime Iz), while Zg is an impure regime and has samples from the two ground truth regimes.
We highlight all the different cases for sample at either near the border or not of regime and also either the causal graph learned for r is on pure or impure data:
Case 1: If 2; is in pure regime and is far from the boundary in the current iteration, Te(w") takes a high value (for example, 10000 (w²) ) in Figure 3). Because regime r was trained on pure data, its causal graph is more accurate, leading to a high likelihood p(z; Ø<e 2, = r, 9*). Consequently, J5 x TH (w") p(ax; ci Ft r, 9r) remains dominant, causing to stay in regime at the next iteration.
Case 2: If Ze is in a pure regime but is near the boundary in the current iteration, Ty (ww")) and (a) are nearly equal (e.g., [1e[1100,1000] (w²) vs. T (w>²) in Figure 3)). Nonetheless, since regime r was learned from pure data, p(æ; Ict, maximum value and maintaining 2; in regime r for the next iteration.
g") stays high, keeping 3 at its
Case 3: If is in an impure regime + 1 near the boundary during the current iteration, Te Te(w") and Ty (1) are also close in value (e.g., #se[1501, 10000 (w= vs. [re[1501,1800 (w²) in Figure 3). However, because the causal graph for regime is more reliable (having been derived from pure data), p(æ; Tgt 2t = r, g") > p( E<tt: &t = r + 1, G)+ As a result, D, moves from regime + 1 to in the next iteration.
Case Dy belongs to impure regime r + 1 and is far from the border (e.g., [1801, 2500] in Figure 3). In this case, it's uncertain whether Is will switch regimes in the next iteration. However, as the pure regime expands with each iteration, H, will eventually be near the border of regime r + 1, bringing us back to Case 3.
C Maximization step:
C.1 Mathematical Derivation: Equation (2) to (6)
We perform the maximization of ELBO presented in proposition 3.1, using a BEM procedure. where we alternate between E-step (updating posterior probabilities while fixing all the parameters = {0', o², w") No and M step (Updating the parameters while using the posteriors learned in the E-step), we can summarize the process as follows:
In the Estep: we learn 3 = (2e 21, Quta) ) where Ouss is O of the previous iteration.
In the M-step: we fix the learned posterior probabilities the values (r and we update O.
By fixing the posterior probabilities the M-step, the entropy of these probabilities H (p (Zt Tee Ools is constant, which allow us to discard it.
18

---

The detail derivation from Eq.(2) to (q) is the following:
E) (9) [Bp(z)lenxee)) Plog (Ie Get + logp (2+)) H (p (Fe oe +<)))
Nw
B
(Or) [logp(9")] + H (qp, (9"))
IT
Σ (S)
Σ Ep(er (ze 2z, Tec Oold (logpe,, (a) I2t + log (ze))
N Bo (0+) [log (G") 11 (go, (9"))
+ H (p (Zt I tee Oold )]
Cte
(9)
We replace p (zz Fe I<t, Oold ) by 8 (zz) by our prior choice ne (w"), and we discard the entropy of the posterior because it is constant in these steps. Hence, we got the (q)
Nw Se (log (zt I<t, g") + log Ty (w"))
ELBO(0) = = E (C))
(10)
N
E
(0)) (ogp((")) H (go. (9")).
After the maximization of (Q) FANTOM update the posteriors (z+ Wi,Tech Ooia) following Eq.(7).
C.2 Variational Inference Details
We provide the detailed formulations for gør (9"). In order to model the temporal adjacency matrices
G where E [| K]], we use two learnable matrices U Q2 Rdxd such that:
exp (uR.i5)
Pk.y
(11)
exp (ux,is) + exp (4k.a))
For instantaneous graphs G we used the same trick in [15, 54], in which we employ three lower triangular learnable matrices Uo, Qus Eo E d to characterise three scenarios: (1) (2)
j + i; (3) no edge between them. For node i j:
exp (Way)
p(i +)) exp (ux) + exp (gis) + exp (eig)
exp (giz)
p(j + exp (u:)) + exp (gij) ) + exp (eig)
(12)
exp (eiy)
p( no edge exp (uy) + exp (g:s) ) + exp (eiz)
With this formulation, the instantaneous adjacency matrix free of self-loops, eliminating any length-1 cycles.
D Limitations & Risk of spurious causality
FANTOM's performance deteriorates when regime contains only handful of samples or is recorded an extremely low sampling rate. This shortcoming is not surprising, estimating separate causal graph for each regime is intrinsically difficult in the presence of multiple regimes and heteroscedastic noise. Yet this ability to pinpoint which edges vanish emerge from one regime to the next what makes FANTOM valuable in domains such as healthcare and climate science, where regime shifts
19

---

earry substantive meaning. Importantly, in realistic settings where each regime offers sufficient data, for example, epileptic seizures that last several minutes 250 Hz, FANTOM delivers strong results and provides insights unattainable with existing methods.
FANTOM requires suitable initial segmentation effectively learn the regime indices. Our ablation studies demonstrate that selecting a reasonable initial window size establishes a basis for accurate regime detection, which is crucial for achieving high performance. However, we highlight a critical limitation regarding the subsequent pruning step: if the initial window is too high (over-pruning), it can drastically reduce the number of initial regimes below the ground truth. This initialization leads to loss of essential information and significant deterioration of FANTOM's final performance.
Spurious causality is practical risk, so advocate using the learned graphs decision support rather than as ground truth. Before acting on them, especially in medical or financial settings, practitioners should adopt safeguards: (i) expert vetting of proposed edges and mechanisms; (ii) robustness checks (e.g., stability under resampling or perturbations, alternative specifications); and (iii) interventional validation when feasible.
E
Data generation and Baselines
E.1
Synthetic data
We employ the Rényi (ER) [34] model with mean degrees of 1 2 to generate lagged graphs and the Barabasi-Albert (BA) [6] model with mean degrees for instantancous graphs The maximum number of lags, L. is set at 1. We experiment with varying numbers of nodes {10, 20, 40} and different numbers of regimes (2,3) 3), each representing diverse causal graphs or mixing functions. The length of each regime is randomly sampled from the set (1000, 1500, 2000, 2500, 3000).
Heteroscedastic case. heteroscedastic settings, noise variance shifts across both variables and observations, making the underlying DAG much harder to recover from data. Given random set if directed acyclic graphs g (g") r[[])][]]]]]- we generate observations from the SEMs in Eq as follows:
Vr (1,..., K),
for (Pa (< t), Pag. (t)) erplgi" (Pai (< t), Pal (t)) 4",
where fhrr, gf are chosen to be randomly initialized ML.s with one hidden layer of size number of nodes and tanh activation functions. é: follows either normal distribution (0,, 1) or a more complex one obtained by transforming samples from a standard Gaussian with an MLP with random weights and sin activation function.
Homoscedastic non-Gaussian case. The formulation used to generated the data is:
ir Vr (1,..., E = f (Pag- (< Pag. + Et'
where fir is general differentiable non-linear function. The function far is random combination between a linear transformation and a randomly chosen function from the set: (Tanh, Exp]. e. follows either Triangular distribution or more complex one obtained by transforming samples from standard Gaussian with an MLP with random weights and sin activation function.
E.2 Baselines
DYNOTEARS [35]. DYNOTEARS formulates causal discovery for multivariate time series through linear vector autoregressive (VAR) model that simultancously captures lagged and instantaneous causal effects. Its key innovation is the DAGness penalty a smooth, continuously differentiable relaxation of the acyclicity constraint optimized via an augmented Lagrangian scheme alongside mean squared error loss. DYNOTEARS emerges as the special case of FANTOM obtained by setting
K = 1, using linear component functions fd, fixing the noise scaling to g = I and d ~ (0, in Eq(I). For comparing with this model, we use publicly available package causalnex*.
PCMCI+ [42]. a scalable two-stage algorithm for time series, enabling the identification of contem- poraneous causal connection. As DYNOTEARS, PCMCI+ is special case of FANTOM, obtained
hattps:/ccasaalne: crr.aaaeeseeeerrrrrttooooooo.os
20

---

by setting K = 1, using linear or non linear component functions ft., fixing the noise scaling to t.1
ghr = 1 and allowing ç to follow any distribution. For the comparison, we use publicly available package Tigramite"
Rhino [15]. Gong et al. propose the first structural equation models with historically dependent noise, where noise variance depends solely on time-lagged variables, the Rhino's SEM is as follow:
क = (Pac(< (), (a) + g (Pac(< (),
(13)
where of (0, Rhino neglects heteroscedasticity, and assumes single stationary regime governed by one causal graph. By our SEM proposed in Eq(l), we can recover the Rhino's SEM by setting K and making g' a function of only time lagged parents. In Rhino, they took a non linear transformation of normal noise which equivalent in our case to allowing cf² follow any distribution. To compare with Rhino, we used the open package causica".
RPCMCI [45]. RPCMCI learns regime indices and time lagged causal relationships from multi- regime MTS. FANTOM's s SEM in Eq(1) generalize it. We can recover RPCMCI settings by making frr depends only on time lagged relations and gf 1. For the comparison, available package Tigramite'.
publicly
CASTOR [39]. CASTOR learns number of regimes their indices and also their corresponding DAGs including instantaneous and time lagged causal relationships from multi-regime MTS. But they assume that they only have gaussian noise with equivariance. FANTOM's SEM in Eq(1) generalize it. We recover CASTOR settings by making g²¹ 1 and ct" (0, 1). For the comparison, use publicly available code CASTOR
F.3 Optimization parameters
Heteroscedastic settings. Unless noted (i.e., in the synthetic-data study), we set the model lag to the true value of I and allow FANTOM to capture instantaneous effects. The variational posterior 9a (9") is initialized to prefer sparse graphs (edge probability < 0.5). Heteroscedastic noise is modeled with conditional normalizing flows (CNFs). Every neural block is a two-layer MLP with 32 hidden units, residual connections, and layer normalization.
Gradients for discrete edges are estimated via the Gumbel-Softmax trick, using hard forward pass and a soft backward pass with temperature 0.25, All spline flows employ 128 bins, and each transformation uses an embedding dimension equal to the number of nodes.
The sparsity penalty is fixed at X, = 50. For graphs with 10 or 20 nodes we use p p= and a = 0. whereas for 40 nodes set p 0.001. Models are optimized with Adam [30] at learning rate of 0.005. We establish Ç = 900 as the minimum regime duration, and we use 1000 as initial window size.
Homoscedastic non-Gaussian settings. Unless noted (i.e., in the synthetic-data study), we set the model lag to the true value of and allow FANTOM to capture instantaneous effects. We use the same parameters as the heteroscedastic settings. The main difference is the use composite of affine-spline transformation. All spline flows employ 16 bins, and each transformation uses an embedding dimension equal to the number of nodes.
The sparsity penalty is fixed at X = 5. For graphs with 10 or 20 nodes we use 1 and a = 0. whereas for 40 nodes we set p 0.001. Models are optimized with Adam [30] a learning rate of 0.005.
E.4
Real world data
E.4.1 Causal Chambers data
We the wind tunnel datasets from Gamella al. [13], featuring two controllable fans pushing air through a chamber, barometers measuring air pressure at various locations, and a hatch controlling an
/tw. http:///ww.....eeeeee///aaaaaaaa btatpa://jakobrrngegggne htt.ge/ggggtabb.o/iiggmaat.....
md
21

---

Puch. Outch berometer
baromeder
Spenker
"Epoaker amplifiatiio cirout
Signal F
Ampliter
source
Polentomelers As
Mossuurd signaks Sr, Or, Re.
Microphane
Exhaust
Inner barometers Pa.. Pap
Intake Pi.. Cu. Th. Ru. On
Mr Ryg Out
LustoRunt Cosa Tea. Rust Oce
Oue.0
Figure 7: Figure taken from Gamella et al. [13]. Diagrams of the wind tunel causal chamber and its main components, including the amplification circuit that drives the speaker of the wind tunnel. The variables measured by the chamber are displayed in black math print. Sensor measurements denoted by a tilde. Manipulable variables, that is, actuators and sensor parameters, are shown in bold symbols.
external opening. The tunnel is chamber with two controllable fans that push air through it and barometers that measure air pressure different locations. A hatch precisely controls the area of an additional opening to the outside (see Figure 7). The dataset comprises 16 variables: controllable load of the two fans Lins Lout. their measurable speed (w/ (ow ). the current draw by the fans
Öi : ( (Cin
Cau ), the resulting air pressure inside the chamber (Paw, Pup or at its intake (
Pint ), and the hatch H. In the circuit that drives the speaker, we can manipulate the potentiometers ( (A1, that control the amplification, monitoring the resulting signal different points of the
circuit ($,, S2 and through the microphone output (M Pamb is the ambient pressure measure by the outer barometer.
We evaluate all the models on two regimes of 10,000 samples each: the first is observational, while the second involves soft interventions on five variables:
Tout the resolution of the tachometer timer that measures the elapsed time between suc- cessive revolutions of the fan. Choosing microseconds yields higher resolution in the fan-speed measurement. Hence intervention on Tou yields change Woat
Oup the oversampling rates when taking measurements of the current (in Cou ), amplifier
(S,,3, ) and microphone signals (M), and of air pressure at the different barometers (Pup: Pamb Pint).
Rm Rout R2 the reference voltages, in volts, of the sensors used to measure the current ( Cin Com ), and amplifier ( ($2 ), respectively.
For the training procedure. we start by a window of 6000 samples, which gives us three different initial regimes then FANTOM converges smoothly to the exact number of regimes K = 2. We set our lag time lagged and we allow the presence of instantaneous parents. Regarding the parameters, we use sparsity coefficient equal to 50, spline of 8 bins and MLPs of size 32.
We compare FANTOM to the baselines mentioned in the main text, with results in Table 2. FANTOM is the only model that detects the regime with 99.9% accuracy and outperforms all baselines on the graph learning task, achieving 38.5% Fl score. Notably, FANTOM surpasses all the models tailored for stationary MTS, even when they are given the ground-truth regime partitions.
E.4.2 Epilepsy data
Huizenga et al. [22] show that scalp potential fields are contaminated by heteroscedastic noise in EEG measurements. We evaluate FANTOM's performance in detecting epileptic regimes using EEG
22

---

EBG Signals
1p1
02
Amplnues
102
Figure 8: EEG signals for patient of id 7170, session the TUSZ data.
signals from 10 different patients in the Temple University Hospital EEG Seizure Corpus (TUSZ) dataset [52]. The dataset encompasses multiple seizure types; in this study, foeus generalized seizures, which engage the entire brain. Each patient's 's record contains scalp EEG signals from 19 channels (Figure 8), each considered causal variable.
State-of-the-art methods [52] use graph neural networks (GNNs) for seizure detection, typically building fixed distance graph and feeding Fast Fourier Transform (FFT) coefficients from 10-12 seconds EEG windows as node features. These approaches (i) train a single model for all patients, offering personalization; (ii) cannot operate in zero-shot unsupervised settings; and (ii) reuse an identical graph for both seizure and normal periods. FANTOM addresses these limitations: detects seizures in personalized manner without any training data and learns distinct temporal causal graphs for normal and seizure states.
Preprocessing, We apply FANTOM to 10 different patients from TUSZ dataset, we treat this as unsupervised regime detection problem, analyzing roughly 100 seconds recordings at 250 Hz sampling rate for each patient, capturing both normal and seizure states. Before running FANTOM. we filter out the EEG signals by a band pass filter of order 6, the lower frequency is 0.5Hz while the highest frequency is 50Hz.
Table 3: Seizure detection accuracy using FANTOM for 10 different patients
Patient id
Regime Acc
0002
84.8
0021
80.8
0302
76.6
0492
86,6
6440
86.1
6520
80.2
7128
94.1
7170
81.1
7936
82.0
8303
75.1
Avg
82.7 82.7+5.2
23

---

Seizure
Background
Seizure
Background
Patient id: 7128
Patient id: 6440
Patient id: 0492
Patient id: 8303
Figure 9: Figure illustrates the summary causal graph per regime learned by FANTOM for different patients
We segment each EEG recording into fixed 12 windows (3000 samples). FANTOM initialized with eight initial regimes and reliably converges to the two ground-truth states; quantitative results appear Table 3. We employ temporal lag of eight samples and allow instantaneous parental links. Averaged over all patients, FANTOM achieves 82.7% -asiimmassinmmnt accuracy. The graph learned for the seizure state is substantially denser and more interconnected than that of the normal state (Figure 9), consistent with the widespread neural involvement of generalized seizures.
F Additional Experiments
F.1
Ablation studies
F.111 Impact of Heteroscedastic Modeling and Flow Architecture
To assess the contribution of heteroscedastic modelling and the flow architecture, we ran an ablation on MTS comparing three FANTOM variants: (i) Gaussian output (no flow, FANTOM (ANTOM-GGuss), (ii) a simple coupling-spline flow (FANTOM-Spline), and (iii) the full model with a conditional normalizing flow (CNF, FANTOM). the single-regime setting, removing or simplifying the CNF consistently degrades performance. In the multi-regime setting, only the CNF variant remains stable and converges to the two ground truth regimes; the alternatives collapse to one regime, Results are reported in the table below:
Gauss),,,
24

---

Table Importance of heteroscedastic modelling and flow architecture for 10 nodes. means the method collapses to regime.
Heteroscedastic noise
Inst
Lag
Regime
K
model
SHD
F1
SHD
F1
Accuracy
FANTOM-
39
0.0
44
42.5
FANTOM-Spline
17
0.0
18
0
FANTOM
0
100
0
100
FANTOM-GAA Gauss
2
FANTOM-Spie Spllne
FANTOM
1 97.8
1
98.9
98.6
F.1.2 Robustness to the pruning step, minimum regime duration, and regime initialization
We evaluated robustness the pruning step, minimum regime duration ((), and regime initialization. Across these variations, FANTOM's predictive performance remains stable. We find that changing (window size. O primarily affects runtime, especially long horizons, without materially impacting graph and regime learning performances. To show efficiency under long-horizon conditions, we used regimes with lengths [2000, 3000, 2000, 2500].
Table Performance of FANTOM with varying hyperparameters (window size and () on 10 node graphs with 2, 3, and 4 regimes.
Heteroscedastic noise
Inst.
Lag
Regime
(window, ()
K
Numb
Running
SHD
F1
SHD
Fl
Acc.
(1000,900)
17'8s
91.2
13
85.7
96.5
(1500,
10' 12s
1
97.9
9
89.8
95.1
(1000,900)
30'33s
5
94.2
8
94.4
96.4
(1500, 1200)
16'2s
96.7
10
91.9
91.8
(1000,900)
40'2s
5
95.8
12
92.7
92.9
(1500, 1200)
25'38s
97.5
16
90.5
91.5
Table Performance of FANTOM with varying hyperparameters (window size and () on node graphs with 2, 3, and 4 regimes.
Heteroscedastic noise
Inst.
Lag
Regime
(window, ()
Numb
Running
SHD
F1
SHD
Fl
Acc.
(1000,900)
22' 15s
33
81.9
27
91.9
100
(1500,
10'20s
35
83.8
25
91.7
100
(1000,900)
3
4
34'25s
49
85.6
58
88.6
99.8
(1500,1200)
15'52s
43
87.3
72
87.3
99.9
(1000,900)
66'30s
34
86.5
53
91.9
99.7
(1500,1200)
21'38s
32
87.1
54
91.8
99.6
From the table, FANTOM scales to 40 nodes and converges to four regimes under different ini- tializations. With a window initialization of 1500 and 1200, it starts with six initial regimes and converges smoothly to four with 99.6% accuracy in 2 m With window initialization of
25

---

1000 and Ç = 900, it starts with nine initial regimes and converges to four with 99.7% accuracy in 66min30s. All rebuttal experiments were run on a Tesla V100-SXM2 (26GB). These results indicate that FANTOM scales without requiring large GPU memory.
To test. scalability, (lengths 2000, 3000, 2000, 2500). As expected, additional regimes and longer horizons increase runtime. With 40 nodes and two regimes, training finishes in 10 min 20 s, whereas the four regime setting requires 21 min 38 s.
FANTOM
40 node dataset with four regimes and long horizons
We showed in the table above (Table 5 and 6) that FANTOM is robust to the initialization parameter (window size that impacts directly the initial regime count and the pruning threshold), e.g. in the case of 40 nodes with 4 regimes, with window initialization of 1500 and 1200, it starts with six initial regimes and converges smoothly to four in 21 min38s. With window initialization of 1000 and Ç 900, it starts with nine initial regimes and converges to four with 99.6% accuracy in 66 min 30s.
E.1.3 Robustness to data standardization
We investigated FANTOM's robustness to data standardization, common challenge for most optimization-base causal discovery models. The results in Table clearly demonstrate FANTOM's resilience. This robustness stems from its of Bayesian structure learning estimate distribution over plausible graphs and conditional normalizing flows to effectively model complex data distri- butions. In contrast, models such as CASTOR are sensitive standardization and even fail when provided with ground-trute regime labels.
Table 7: Robustness of FANTOM to data standardization (10 nodes).
Heteroscedastic noise
Inst.
Lag
Regime
model
data
SHD
F1
SHI)
F1
Acc.
CASTOR
Raw
90
26.9
28
29.7
x
standardized
28
0.0
33
0.0
X
Raw
4
93.3
7
90.9
97.1
FANTOM
standardized
5
91.5
7
90.9
91.1
CASTOR
Raw
135
24.1
135
29.6
x
standardized
37
0.0
48
0.0
Raw
3
96.1
10
91.5
92.2
FANTOM
standardized
2
97.3
9
92.3
92.2
From the Table 7, FANTOM achieves the same results on raw and standardized data while CASTOR fails in both cases even when give it access to the ground truth regime labels. In the case of standardized data, CASTOR predicts adjacency matrices full of zeros due its incapability of handling such scenarios.
26

---

F.2 Additional results on synthetic data
F.2.1
Heteroscedastic noise with different number of nodes and regimes
Table 8: Average SHD, Fl scores, NHD and Ratio for different models with d = 10 nodes and K regimes, Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
leteroscedastic noise, K 2 and d
10
Model
Split
Inst.
Lag
Regime
SHD
FIf
NHD4
Ratio!
SHD
FIf
NHD]
Ratio,
Acc.
PCMCI+
X
34.3+3.
14.2+2s
0.09
0.85+0.1
26.0-43
25.6448
0.07
0.74+0.
X
Rhino
x
28.025.1
8.56149
0.09
0.924am
38.5
58.4
0.13-00
0.43
x
DYNOTEARS
58.5
31.2a28
0.22.0.0
0.68
79.5-aa
33,9+44
0.27=6.
0.66
CASTOR
90.0 X
28.1
0.38
0.74
90.5_a2
35.6134
0.42
0.64
x
RPCMCT
CASTOR
FANTOM
4.67 21.
91.7 132
0.005
10.0
0.04
13.
11.6 88.2
0.02 0.0 0.08 10.0
96.6.11
Table 9: Average SHD, Fl scores, NHD and Ratio for different models with d = 10 nodes and K = regimes. Split denotes whether regime separation is automatic (v) manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Heteroscodastic noise, K and d
Model
Split
Inst.
Lag
Regime
SHI
F11
NHD]
Ratio
SHD
Fif
NHD.
Ratio
Acc.
PCMCI+
46.
11.1
0.05
0.88
46.0
19.0
0.05
0.80
X
Rhino
X
44.5
5.11
0.06..00
0.94-00
53.5-1.
64.7146
0.07
0.35
X
DYNOTEARS
89.5+15
31.5+04
0.14
0.68-oo
118.0=4m
37.5+12
0.17+aa
0.61+00
CASTOR
104.
23.4
0.19
0.76
133.5
34.8
0.24
0.64
x
RPCMCI
CASTOR
FANTOM
5.67 93.3+12
0.006
0.06
12.3
90.9+1.
0.012+aa
0.08
97.1
Under heteroscedastic conditions, d = 10 node graphs with either two or three regimes, FANTOM outperforms every baseline. Among methods explicitly designed for multi-regime MTS, it is the only one that converges to the true regime partitions, attaining 97.1 1% (Table 9) in regime detection accuracy. CASTOR and RPCMCI break down in heteroscedastic settings.
To give stationary MTS baselines the best possible chance, we supply them with the ground truth regime partitions. This is done by training the aforementioned models on each pure regime separately (regime governed by the same causal model). Yet, FANTOM still dominates the structure learning task, while learning the number of regime and their indices well, achieving an NHD of 0.006 (Table 9). Because NHD penalizes every missing, extra, or mis-oriented edge, a value of 0.006 implies that FANTOM not only recovers the graph skeleton but orients edges with high precision.
Fl of 93.3 % and
Table 10: Average SHD, Fl scores, NHD and Ratio for different models with d 20 nodes and K regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Heteroscedastic noise, and
Inst.
Lag
Regime
Model
Split
SHD|
F11
NHDJ
Ratiol.
SHDI
FIf
NHD.
Ratiol,
Acc.
PCMCI+
84.0
39.7
0.03
0.6
42.5
76.9
0.1410a 0.22
X
Rhino
x
70.0+60
1.26+12
0.04
0.98-no
188.0.3.
26.3+as
0.14+3a
0.73+ao
X
DYNOTEARS
221.5
26.9-12
0.22
0.72
45.0gmm
61.52a9
0.02
0.38
CASTOR
377.511s
14.9101
0.84
379.511s
19.1
0.41
0,80_ao
x
RPCMCI
CASTOR
FANTOM
9,00
89.1 137
0.006
0.10.eo
26.0
85.4
0.01
0.14..o
97.8.02
27

---

Table 11: Average SHD, Fl scores, NHD and Ratio for different models with d 20 nodes and K 3 regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Heteroscedastic noise, K and
20
Model
Split
Inst.
Lag
Regime
SHD4
FIf
NHD
Ratio]
SHD
F1f
NHD,
Ratio.
Acc.
PCMCI+
X
83.5
39.9+16
0.03+00
0,59
38.0,2a
77.8.09
0.01 0.22
X
Rhino
X
108.5=11
1.6+12
0.02+0.
0,98+a0
292.0+a2
25.8+8.
0.09400
0.74+0.0
X
DYNOTEARS
325 -4.
27.6+10
0.14+0.0
0.72+ao
55.5+0.5
69.8+7.
+.01
0.29+0.0
X
CASTOR
412.3
13.8
0.20
0.86
570.
17.5
0.28
0.81
X
RPCMCI
CASTOR
FANTOM
13.5 6.5
88.2+62
0.05+0.0
0.11 +00 46.5+a5
84.6.25
0.01+00
0.14+0.0 97.8+1.
Under heteroscedastic conditions, d = 20 node graphs with either two three regimes, FANTOM outperforms all the baselines. Among methods explicitly designed for multi-regime MTS, it is the only one that converges to the true regime partitions, attaining 97.8% in regime detection accuracy (Table 10). CASTOR and RPCMCI break down in heteroscedastic settings, also in the case of 20 nodes.
To give stationary MTS baselines the best possible chance, we supply them with the ground truth regime partitions. This is done by training the aforementioned models on each pure regime separately (regime governed by the same causal model). Yet, FANTOM still dominates the structure learning task, while learning the number of regime and their indices as well, achieving Fl of 89.1 % and an NHD of 0.006 for instantaneous links and an Fl of 95.4% and and an NHD 0.01% % for time lagged (Table 10).
Table 12: Average SHD, Fl scores, NHD and Ratio for different models with d 40 nodes and K = 2 regimes. Split denotes whether regime separation is automatic (V) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Heteroscedastic noise, K 2 and d 40
Model
Split
Inst.
Lag
Regime
SHD4
F1f
NHD.
Ratiol SHD!
Fif
NHD.
Ratio.
Acc.
PCMCI+
X 146.5=5.0
25.8+a.1
0.02=m.
0.74
109.0+16.
38.0+4.
0.01=0.0
0.42-o.
x
Rhino
X
137.0
0.0
0.02
1.00
700.0sst.
28.215.7
0.12
0.71
x
DYNOTEARS
X
137.5192
17.1
0.0201aa
0.82
129.0.1s.
36.6149
0.01
0.63
x
CASTOR
x 151.0+1.
0.0+a.
0.03
1.00=oc
333.0+93
19.6+56
0.050+aa
0.80
x
RPCMCI
CASTOR
V
FANTOM
27.00.6.0
85.6137
0.005100
0.14..0 .00 26.5
91.9
0.004100
0.08 0.0 99.9...1
Table 13: Average SHD, Fl scores, NHD and Ratio for different models with d = 40 nodes and K regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to tame-lagged edges.
Heteroscedastic noise, 3 and 40
Model
Split
Inst.
Lag
Regime
SHD
FI1
NHD
Ratio
SHD
F1f
NHD]
Ratio
Ave.
PCMCI+
x
222.0+14
25,2+04
0.01+a.
0.75+no
142.0+17.
62.4+10
0.01+6o
0.37+a.
x
Rhino
x
210.5102
0.0
0.01
1.00
1005.
29.3154
0,08
0.7.0.
x
DYNOTEARS
203.0++4
15.2+22
0.01=aa
0.85+ao
161.0=14
51.3+15
0.01
400
0.49+a.1
x
CASTOR
X
224.0.102
0.00
0.01
1.00
501.0.a.
23.1
0.03
0.76
x
RPCMCI
CASTOR
FANTOM
52.0 82.2
0.004.0.0
65.0
17.
87.9
0.004...
0.11 inn 99.8
Under heteroscedastic conditions, d = 40 node graphs with either two or three regimes, FANTOM outperforms every baselineand shows that it can scale to large graphs even in this complex setting. Among methods explicitly designed for multi-regime MTS, it is the only one that converges to the true regime partitions, attaining 99.9% in regime detection accuracy. CASTOR and RPCMCI break
28

---

down in heteroscedastic settings (Table 12). Large scale graphs helps FANTOM to differentiate between the different regimes and increases the regime detection accuracy by 2% compared to 20 node graphs settings.
To give stationary MTS baselines the best possible chance, we supply them with the ground truth regime partitions. This is done by training the aforementioned models on each pure regime separately (regime governed by the same causal model). Yet, FANTOM still dominates the structure learning task, while learning the number of regime and their indices as well, achieving an FI of 85.6 % and an NHD of 0.14 (Table 12).
F.2.2 Non-Gaussian noise with different number of nodes and regimes
Table 14: Average SHD, Fl scores, NHD and Ratio for different models with d = 10 nodes and K regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Homoscedastic non- Gaussians noise, K d 10
Model
Split
Type
Inst.
Lag
Regime
SHD
F17
NHD.
Ratio,
SHDJ
FIf
NHD]
Ratiol
Acc.
PCMCI+
X
6.00+a6
88.1-05
0.01+a0
0.12-00
8.00+5.
90.7+6.1
0.01
0.09+00
X
Rhino
W
2.50
96.0.a.
0.004
0.04
4.00
96.0
0.006.0.0
0.04.00
X
DYNOTEARS
X
W
42.0.n. 51.8279
0.12
0.48
.n0
8.00
86.8
0.02
0.12
X
CASTOR
W
17.0ga0 62.6285
0.055
0.37
11.0g1
85.2.o.
0.02ao0
0.15
RPCMCI
CASTOR
W
34.0
20.
42.5
0.09 a00
0.57_na
45.0
47.0
0.13
0.52
77.0
FANTOM
W
1.00.
98.2....
0.002
.0.0
0.01 0.0
8,00,40
91.0..aa
0.02.co
0.08 +0:36
98.6.00
SHD
Fif
NHDJ
Ratio
Acc.
CD.NOD
33.0gam
47.0as.
0.41 son
0.52
x
FANTOM
S
7.5
93.1+22
0.07
0,06
99.6=0.0
Table 15: Average SHD, Fl scores, NHD and Ratio for different models with d 10 nodes and K regimes. Split denotes whether regime separation is automatic (V) or manual (x). Inst. refers to instantancous links, and Lag time-lagged edges.
Homoscedastic non Cou-saas noise, K 3 and d
10
Model
Split
Type
Inst.
Lag
Regime
SHD
FI1
NHD,
Ratio
SHD
F1f
NHD
Ratio
Acc.
PCMCI+
x
W
17.5
74.9 5.
0.02
0.24
14.5
88.3116
0.01
0.11,00
Rhino
W
2.50
96.8_13
0.002ano
0.03
6.00
95.2
0.006
0.04
DYNOTEARS
X
W
42.0.s..
54.4=m1
0.06
0.45
21.0=1.4
82.
0.02
0.17
CASTOR
22.0+28 66.2-25
0.030+no
0.33+ae
17.0-12
84.4+18
0.01-mo
0.15+on
RPCMCI
W
CASTOR
W
47.0
34.8.m.
0.05
0.65
59.5
39.4.11.
0.07
0.60
51.6289
FANTOM
W
0.33+a4
99.5
0.00
0.00
12.5
89.1
0.01 0.01=n.0
0.10+au
99.4.01
SHD..
F1
NHD.
Ratio 4
Acc.
CD-NOD
X
S
42.5
31.8-111
0.42-no
0.67 0.04:09
x
FANTOM
4.5+20
95.6 1.5
0.04.
96.6.0.2
Under homoscedastic, non non-Gaussian noise with d = 10 node graphs and either two or three regimes, FANTOM outperforms every bascline instantaneous-link inference, reaching Fl of 98.2 % and an NHD of 0.002. Among methods explicitly designed for multi-regime MTS, both FANTOM and CASTOR recover the exact number of regimes, whereas RPCMCI fails to converge to the true partitions. FANTOM further surpasses CASTOR in regime detection (98.6 % vs. 77.0 %) and in DAG learning (F1 98.2 % vs. 42.5 %)..
To give the stationary-Mo
baselines their best chance, we supply them with the ground-truth regime labels and train each model on the corresponding pure regime. Even in this favorable setting, only Rhino exceeds FANTOM on time-lagged links, achieving an F1 of 96.0 % compared with FANTOM's 91.0% %.
29

---

Table 16: Average SHD, Fl scores, NHD and Ratio for different models with d 20 nodes and K 2 regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Homoscedastic non- Gaussian noise. K 2 and 20
Model
Split Type
Inst.
Lag
Regime
SHD
FIf
NHDJ
Ratio
SHDI
Flf
NHDL
Ratio
Acc.
PCMCI+
W
46.0.ns
54,9-u6
0,03+00
0,45+00
17.0-4a
88.7+a5
0.009+0.0
0.11+0.0
X
Rhino
W
17.5eu.
82.5=16
0.007
0.17
29.5=as
83.5-21
0.01
0.16
DYNOTEARS
W
44,5-a5
44,9-31
0.03+00
0.55+0.0
46,5-12
55,8-59
0,025
0.44+00
CASTOR
W
139.5 41.8
0.10
0.58
186.0
40.1
0.13
0.60
RPCMCI
W
CASTOR
W
68.0gas
37.2-12
0.04
0.62
86.0.n9
37.4
0.06
0.64
46.6
FANTOM
W
3.50
97.2
0.002
0.02aco
11.5
93.1=16
0.006
0.06
100..0.
SHD
FI7
NHD
Ratio
Acc.
CD-NOD
S
106.
26.0=a1
0.31
0.73-na
FANTOM
S
4.0-0.
98.3
0.01
0.01
100.
Under homoscedastic, non-Gaussian noise with d = 20 node graphs and two regimes, FANTOM outperforms every baseline instantaneous-ling link and time lagged link inference, reaching an F1 of 97.2 % and an NHD of 0.002 for instantaneous links and an F1 of 93.1 % on time lagged relationships. Among methods explicitly designed for multi-regime MTS, both FANTOM and CASTOR recover the exact number of regimes, whereas RPCMCI fails to converge to the true partitions. FANTOM further surpasses CASTOR in regime detection (100. % vs. 46.6 %) and in DAG learning (F1 97.2 % Vs. 37.2 %).
To give the stationary- -TTS basclines their best chance, we supply them with the ground regime labels and train each model on the corresponding pure regime. Even in this favorable setting, FANTOM out performs all the baselines, achieving an F1 of 97.2 %.
Table 17: Average SHD, Fl scores, NHD and Ratio for different models with d 20 nodes and K = 3 regimes. Split denotes whether regime separation is automatic (V) or manual (x). Inst. refers instantaneous links, and Lag time-lagged edges.
Homoscedastic non Gaussian noise, K 3 and d
20
Model
Split Type
Inst.
Lag
Regime
SHD4
Flf
NHD4
Ratio
SHDI
FIf
NHDI
Ratio_
Acc.
PCMCI+
X
W
66.5-21
55.1
0.02
0.45
24.5
88.9
0.007
0.11
x
Rhino
W
27.5
82.1
0.007
0.17
50.52as
81.8
0.01
100
0.18
DYNOTEARS
X
W
69.0aaa
43.6
0.02aae
0.56aco
54.5
66.3a24
0.01
0.34
x
CASTOR
X
W
167.0_99
38.3246
0.05
0.61aco
264.51s4
41.2215
0.08
0.58
RPCMCI
W
CASTOR
W
121.5
34.3-
0.03
0.65
+0.1
181.5
33.7
0.05
0.66
79.8
FANTOM
W
7.00
96.1
0.001
0.03
45.51250
85.2
0.01
0.14
100..o.m
SHDJ
F1
NHDJ
Ratio
Acc.
CD-NOD
S
133.
31.5
0.43
0.61_ap
FANTOM
S
6.0=1.
98.2.m.2
0.01
0.01
100.
Under homoscedastic, non-Gaussian noise with = 20 node graphs and three regimes, FANTOM outperforms every baseline on instantancous- link inference, reaching an Fl of 96.1 % and an NHD of 0.001. Among methods explicitly designed for multi-regime MTS, both FANTOM and CASTOR recover the exact number of regimes, whereas RPCMCI fails to converge to the true partitions. FANTOM further surpasses CASTOR in regime detection (100. % vs. 79.8 %) and DAG learning (F1 96.1 % vs. 34.3 %).
To give the stationary- baselines their best chance, we supply them with the ground-truth regime labels and train each model on the corresponding pure regime. For this scenario of 20 nodes and 3 regimes and benefiting from regime labels, PCMCI+ exceeds FANTOM and Rhino on time-lagged links, achieving an Fl of 88.9 % compared with FANTOM's 85.2%.
30

---

Table 18: Average SHD, Fl scores, NHD and Ratio for different models with d nodes and K regimes. Split denotes whether regime separation is automatic (v) or manual (x). Inst. refers instantaneous links, and Lag to time-lagged edges.
Homoscedastic non Gaussian noise. K 2 and 40
Model
Split Type
Inst.
Lag
Regime
SHD1
FIf
NHDJ
Ratio
SHDI
FIf
NHDL
Ratio
Acc.
PCMCI+
W
91.0.as
56.7-a1
0.01+00
0,4340.
25.0-14
91.9-os
0.003+0.
0.08
X
Rhino
24.5.2.
90.7.0s
a.004ac.
0.09 10.0
34.5.50
89.8.as
0.005 1000
2.10
X
DYNOTEARS
W
100.0-57
40.8-37
0.015+6a
0.59+6.
58.5-17
76,4-54
0,009+00
0.21
CASTOR
94.0 67.3
0.010
0.33
100.0
78.6.54
3.01
0.21 30.0
RPCMCI
W
CASTOR
W
FANTOM
W
29.5
1
0.004
0.1
32.3
90.8
0.005
0.08
100..o.m
SHD
FL
NHD
Ratio
Acc.
CD-NOD
X
260.a50
17.8
0.18
0.82
x
FANTOM
S
39.5
92,4
0.02
0.07
100..oo
Under homoscedastic, non- Gaussian noise with 40-node graphs and two regimes. FANTOM is the only method that recovers the correct number of regimes, whereas CASTOR and RPCMCI fail to converge. FANTOM achieves 100% regime-detection accuracy and an 88.1% F1 score DAG recovery.
For fair comparison with stationary-MTS baselines, we provide these models with the ground-truth regime labels and train them separately on each pure regime. In this advantaged setting, Rhino surpasses FANTOM instantaneous links (F1 90.7% vs. 88.1%), while PCMCI+ leads time-lagged links (F1 = 91.9% vs. 90.8%).
Table 19: Average SHD, F1 scores, NHD and Ratio for different models with d 40 nodes and K regimes. Split denotes whether regime separation is automatic (v) manual (x). Inst. refers to instantaneous links, and Lag to time-lagged edges.
Homoscedastic non-Gaussian noise, K and 40
Model
Split Type
Inst.
Lag
Regime
SHDI
FIf
NHD.
Ratio
SHDI
F11
NHD.
Ratio]
Acc.
PCMCI+
W
141.5
56.5
0.008
0.43
37.5192
91.6128
0.003
0.08
Rhino
X
W
53.5
85.5
0.0042no
0.14
52.01s.7
89.2122
0.004
0.1 Laa
DYNOTEARS
W
152.0+12
40.8+a7
0.01-00
0.59-am
91.0+17.
75.8
0.006-co
0.24
CASTOR
X
W
93.5+22 65.9+76
0.009
0.34
94.0.a5
78.6+15
0.009
0.21+an
X
RPCMCI
W
CASTOR
W
FANTOM
W
35.3+
.0.
0.002 000
0.09
40.0
46.3+0.
90.7+13
0.003
-0.0
0.09
100..0.0
SHD
F1
NHD
Ratio |
Acc.
CD-NOD
S
342.
17.1
0.24
0.82_m.
FANTOM
S
57.0+26
92.8-n7
0.03
0.07
100..o.
F3 Additional experiments for L-=2
Table 20: Performance of FANTOM on L 2 time series under varying node counts d {10,20, 403 and regime settings
and 3
Heteroscedastic noise
Inst.
Lag
Regime
k
SHD
F1
SHD
Fl
Acc.
10
97.7
0
100.0
94.3
2
20
14
87.0
4
98.7
99.0
40
15
94.5
3
99.5
98.9
10
1
98.5
1
99.5
91.5
3
5
97.5
11
97.4
97.9
40
43
88.6
3
99.0
99.0
31

---

FANTOM maintains similar performance with a larger lag (L = 2). It achieves an Fl score above 87% on instantaneous links for 10, 20, and 40 nodes MTS with or regimes. For time-lagged links, we evaluate using the summary graph of lagged relations; the table shows that FANTOM detects them effectively.
F.3.1 Illustrations of learned graphs
Estimased de Regime ]
True Ce. Regine
Diffrrnce 16 Gsl Regine
Estrmatedd di hegmm 2
True Cr. Rerime
Difermnce 16 Gil Regme 3
Extimated fir Regime? 2
True Go Regime
Diferencce (6 (6-66) Reg me
Estimated 6, Regime 2
True Cy Rugime
Difference 16 -G,) Regime
Figure 10: The estimated temporal causal graphs for two regimes, in Heteroscedastic case, consist of one matrix of 10 rows and 10 columns representing instantaneous links and another of 10 rows and 10 columns delineating time-lagged relations (with a maximum lag T. in this case). Dark red indicates value of one (presence of an edge), while gray symbolizes value of 0 (absence of an edge). The second column displays the ground-truth causal graphs, and the final column highlights the difference between the estimated and true graphs.
32

---

Regime
Regime
(6 Regime
Estimated 6, Reaime
Difference i
Estimated G, Regime
Difference 16
Estimated G Regime 2
True Regine
Difference i G.I. Regine
Extimated Go Regime
True Gr. Regim
Difference (6 Gol. Regime
Estimated G Regime 3
Trae G., Regine
Differance id G.l. Rogine
Figure 11: The estimated temporal causal graphs for three regimes, in Heteroscedastic case, consist of one matrix of 10 rows and 10 columns representing instantaneous links and another of 10 rows and 10 columns delineating time-lagged relations (with maximum lag I. = in this case). Dark red indicates value of one (presence of an edge), while gray symbolizes value of (absence of an edge). The second column displays the ground-truth causal graphs, and the final column highlights difference between the estimated and true graphs.
33

---

Estimated
Diferrnce 16,
Regime
Regine
J6, Gel Rogin
Negime
C.. Regim
Lifferecce JC, -a. Nogime
The estimated temporal causal graphs for three regimes, in Heteroscedastic case, consist of one matrix of 40 rows and 40 columns representing instantaneous links and another of 40 rows and 40 columns delineating time- time-lagged relations (with maximum lag L = in this case). Dark red indicates value of one (presence of an edge), while gray symbolizes value of (absence of an edge). The second column displays the ground-truth causal graphs, and the final column highlights the difference between the estimated and true graphs.
Figure 12:
34

---

Regime
Regim1
(6 Regime
2
Estimated G, Regime
Difference C
Extimated G. Regime
Difference GG
Estimated G. Regime 2
Trus Regine
Difference C G.I. Regine
Estimated Co Regime
True Regie
Differncce 16 Gol. Regim
Estinated G Regime 3
Tras G.,, Regine
Differance C G.I. Rogine
Figure 13: The estimated temporal causal graphs for three regimes, in Non- Gaussian case, consist of one matrix of 40 rows and 40 columns representing instantancous links and another of 40 rows and 40 columns delineating time-lagged relations (with maximum lag I. = 1 in this case). Dark red indicates value of one (presence of an edge), while gray symbolizes a value of (absence of an edge). The second column displays the ground-truth causal graphs, and the final column highlights the difference between the estimated and true graphs.
F.3.2 Time complexity analysis
We start first by computing the time complexity of our Temporal Graph Neural network illustrated in Figure 2. We note the input size, e the embedding size used for eij and h hidden layer size. After the first NN block,
TNN-1 O(dd + e)h hd) O(d) ((h(++)))
35

---

then after the matrix multiplication block and the second NN block, we have:
Forwardd
O (Lde 2h(d e)).
where L. is the maximum lag.
Using the same architecture for Conditional normalizing flow has time complexity of:
TeNp = O (Ld + h + dK)).
The complexity of FANTOM per iteration is O (N (T (2Ld + h (e + dk)), where K is the number of bins, N. is the number of regimes, and [7] is the number of samples.
F.3.3 Regime detection experiments
We compare FANTOM CASTOR [39] and KCP [1] in the task of regime detection. KCP is multiple change-point detection method designed to handle univariate, multivariate, or complex data. Being non-parametric, KCP does necessitate knowing the true number of change points in advance. detects abrupt changes in the complete distribution the data by employing a characteristic kernel.
CASTOR is causal discovery model specifically designed for multi-regime MTS. CASTOR is learns number of regime and their indices and their corresponding causal graphs without any prior knowledge. But it is limited to normal noise with equivariance. FANTOM learns the regime indices, while handling heteroscedastic noises.
We opted to perform regime detection with 10 nodes and three different regimes. For a fair comparison, we chose three regimes without re-occurrence, as KCP only detect change points and cannot identify the re-occurrence of specific regime.
Regarding the models employed, by the authors?. For KCP, we employ the Rupture package".
the open-source code of CASTOR implemented in Python
Regime Accuracy
80
40
30
KCP
CASTOR Methods
FANTOM
Figure 14: Comparison between FANTOM, CASTOR and KCP on regime detection for a MTS with heteroscedastic noise and composed of regimes using accuracy metric. The number of nodes is ==11 10.
From Figure 14, evident that FANTOM outperform CASTOR and the change-point detection method KCP. This outcome can be attributed to the limitation of KCP in detecting changing points within causal mechanisms that are represented by conditional distributions. FANTOM outperforms
1/CASTRR
"attps://centze-bbeelll git.ub.tthbb.......... htw.huuuooo.ouuuuuuudoooo.ss
36

---

CASTOR in detecting regime indices. This result can be explained by the fact that CASTOR fails in handling heteroscedastic noises and fails to learn meaningful graphs which also lead to poor performance in regime detection.
From this analysis and the other experiments shown in the different tables, in scenarios involving MTS with multiple regimes with non-Gaussian or Heteroscedastic noises, FANTOM offers robust solution. Additionally, employing other methods split the regimes and learn the causal graph through traditional causal discovery methods may not be an optimal solution:
conclude that
We demonstrate that regime indices are not well recoverable by other state-of of-the-art change point detection method KCP. Therefore, employing KCP to learn the regimes and subse- quently using methods like DYNOTEARS, PCMCI+, or Rhino to learn the graph may not constitute optimal solution.
In cases of regime recurrence, the aforementioned methods are unable to accurately detect the exact number of regimes. Therefore, if employs KCP and subsequently the regime partitions revealed by KCP as an input to a causal discovery method (such as PCMCI+, DYNOTEARS, Rhino, etc.), the running time will be significantly high.
We show throughout all our experiments, that FANTOM outperforms all causal discovery method in DAG learning task, even when these models are more favorable scenarios, by having access the regime labels beforehand.
G Proof of proposition
In this section, we are going to provide the entire proof of our proposition 3.1. As we state before, we note G = (9") [e[]iell:
Proof. we have:
log (IIET)
log pe (x), x22))
po (2e ce ((9 (9)
log Σ
(g)
g
E (y (logre (a. 9) + PG) log g0 (9))
Let's focus on the first term of our last inequality logpe (x, >< 9), we have:
past (x G, 24 p(z)) (( Ør, x a.))
logpe Pe (z 9) = loog Σ
p Øt, x xx1)
> Epte)ler
plog (xa g.t) ) + log logp(z;) 2x) - 1x ac))
2 (p(z)m))ee)) plog po= (xx Eet. gz1) + log (2)) + (p(2)/zr,
Including this result in the previous equation gives the following:
log pe (FIET)
B9)( (ppe [log pa (x,, (eces (za) H (p(zalue, xce))]
Eggr (9r) [logp (g") + H (gor (g"))
FLBOO(O),
note that the priors p(g") and the variational estimations gor (g") are independents.
37

---

H Proofs of our theoretical contributions
In this section, concentrate establishing the identifiability of regimes and causal graphs within the FANTOM framework. Before diving into the details, let us set and clarify the required assumptions.
H.1 Assumptions
Definition H.1 (Causal Stationarity, [41]. A stationary time series process (#t))TT with graph is called causally stationary over a time index set T if and only if for all links oi-, + of in the graph
x/-r MC of
Ce {
This elucidates the inherent characteristics of the time-series data generation mechanism, thereby validating the choice of the auto-regressive model.
Assumption H.2 (Causal Stationarity for MTS with multiple regime). A MTS (#t)te with K regimes, graph set (g") [[]: and regime partition & = (Em)ue[ll:K]) causally stationary over the time index set if, for each regime u E [[] K]], the sub-series (It)tes, causally stationary with graph 9 as defined in definition H.1.
Definition H.3. (Causal Markov Property, [37]). Given DAG & and joint distribution p, this distribution is said to satisfy causal Markov property w.r.t. the DAG G if each variable is independent of its non- 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 H.4 (Causal Markov Property (CMP)). A set of joint distributions (p((l9)))
satisfies the CMP with respect the DAGs (g") ([[[][]]]] if, for each E [|] K]]. the distribution
p(-19") satisfies the CMP relative the DAG G". Specifically, in every regime r, each variable is independent of its non-descendants given its parents.
joint distribu- tion (p(-I9")) (rr))))))((:)))) we say that this of distributions satisfies causal minimality w.r.t. the set of DAGs (g"), [E[]IK]]] if for every r: p(-1g") is Markovian w.r.t the DAG G" but not any proper subgraph of g".
Assumption H.5 (Causal Minimality). Given a set of DAGs (9") (r"))e[[: and
Assumption H.6 (Causal Sufficiency). process if and only if in the process every common cause of any two or more variables in V is in V has the same value for all units the population.
of observed variables V is causally sufficient for
This assumption implics there are no latent confounders present in the time-series data.
Causal graph
Causal Markov
Causal sufficiency
Faithfulness
Heteroscedastic noise
Stationarity per regime
/ Minimality
DYNOTEARS
W
X
X
PCMCI+
W
RPCMCI
W
F M
X
Rhino
V
X
X
CD-NOD
S
V
F
X
CASTOR
X
FANTOM
W
M
Table 21: Summary of the main assumptions of algorithms considered in the paper. For causal graphs, means that the algorithm provides summary causal graph and W means that the algorithm provides a window causal graph; corresponds to faithfulness and M to minimality. An empty cell mean that the information given in the corresponding column was not discussed by the authors the corresponding algorithm.
The table 21 illustrates that most assumptions (causal sufficiency, causal Markov, faithfulness/mini- mality) commonly shared among various state-of-the-arr models in causal discovery.
38

---

However, FANTOM, CASTOR, RPCMCI, and CD-NOD relax the assumption of stationarity and instead assume that the MTS (Multivariate Time Series) are composed of different regimes. While CD-NOD predicts only summary causal graph, FANTOM, CASTOR and RPCMCI predict window causal graph, which can subsequently be used to reconstruct summary graph. FANTOM the only model that can handle heteroscedastic noise.
H.2 Proof of theorem 4.2
We start first by proving theorem 4.2. so, we will prove identifiability the case of bivariate time series, Lemma H.7. Then we will prove identifiability in the case of MTS.
Lemma H.7 Assume Causal Markov property, minimality, stationarity, sufficiency and (st x KET
be bivariate time series such that (st, "EET following Eq( I) where K and e N NN007 1). We have o) ot follow Gaussian distribution, if f" f2 are non linear and 1. gt g are not a polyno- mial of degree two then the bivariate Temporal heteroscedastic Gaussian noise (THGNM) model is identifiable.
Proof of the Lemma. Let's assume we have two temporal causal graph G and G' for the bivariate TGHNM.
Disagreement in time lagged relationships. Assume that and g' do not differ in the instantaneous effects. Vi {1,2): (t) Pag:(t). Hence and Wlog, there is some k 0 and edge
=/=x 2 in 9 but not in G'. From g' and the Causal Markov property, we have that s/_k IL 2? S, where (= S S L,iE (1, U UND)), {=!_ and ND, are all X; that are non- descendants (wrt instantancous effects) 2². Applied to G, causal minimality leads to a contradiction because r/-a 1d xt S in g, and the above reasoning shows that it exists subgraph G' of g that is Markovian to the joint distribution of the data.
Disagreement on instantaneous parents. Now, let's assume we have forward model, Vt E T:
2 = f (Pa) (), 2), + g (Pab(< f), d.
We will prove by contradiction that a backward model
x? = 2 (Pag 2)) + g² (Pa) f) x|=)) e2
can not exists.
We note h. Pat(< t) U Pağ(<< t). We know that of I1 (zt Pab(<< t), Pag(<< t)) and
ef IL (x., Pap(< f), ( < t)), using Lemma 36 Peters al. [38], we have:
x = f (pa a8lm)) +
2 = f² (pa2 (pa2.: =bln.) + g² (pa2 1n, E7
where wi] h, is xt conditioned on ht This last result contradicts the theorem states by Khemakhem et al. [29]. Hence, our bivariate TGHNM identifiable.
Let's prove this results in the case of MTS (Theorem 4.2). In the case of Disagreement time lagged relationships, we can use the same proof for the bivariate case.
Disagreement instantaneous parents. Let's assume we have two temporal causal graph G and G' such that G + g'. According to the Propostion 28 in Peters et al [38], for 9 and g' be two different DAGS over a of variables V. such that w, is generated by our HNM and satisfies the Markov condition and causal minimality with respect to G and G'. Then there are variables o7, x7 E V such
that for the set Q = (t){ Y == Pag (1) {x!) } and S = Q U Y, we have: 1) x? > x! in & and / - 2 in g'. 2) S C ND variables of a in graph G. ND the set of non-descnndant (wrt instantaneous effects) of x) in graph g.
#f Maf} and S C ND d' #?
Mæ!) (M+)) Pa) is the set of parent
We consider = with p(a) > 0. Denote a.! : =ff S = and ? == x? S= 8. Lemma 37 in Peters et al. [38] states that if p(æ,) is generated according to the SEM models as follows:
39

---

= Ji (Pag(<< 1), (a)(U), e) i ,d) EV
with corresponding DAG G, then for variable zi E V. S S CND then ej 1L S. Our TGHNM can be viewed one specific class of the SEM in the aforementioned equation. Hence. Lemma 37 holds under our TGHNM and renders e IL (2t. S) and ef IL (at 5): using Lemma 36 in Peters et al. [38], we have:
46 a/ h. = f1 (q.paker (n.) 18 + gt (q, 2f c
2 a =
P-(( (P. aa + (y, pales mn)) E
where rile,, z/ conditioned on ht, This results contradict our previous proved Lemma, then THGNM is identifiable model under the conditions stated in the theorem.
Theorem H.8 (Identifiability of Temporal Non Gaussian noise model (TNGNM)). Assume Causal Markov property, stationarity, minimality, sufficiency and let (II)ET be MTS following TNGNM, Vt S
= (Pag(< 1), Pa) +
(14)
where f is differentiable function, and ei are mutually independent noises and follow non Gaussian distribution. The TNGNM identifiable
Proof. The proof of this theorem could be concluded from theorem 1 in Rhino [15]. Eq(14) is special case of Rhino SEMs.
H.3 Identifiability results in the case of Temporal General Heteroscedastic Noise Models
In this section, we will present our identifiability results for the case of Temporal General Het- eroscedastic Noise, where MTS has the following SEM :i E T.
f (Pat(<< t), Pa)(()) + (Pag(< t), (t)(())) c4.
(15)
where f and g' are differentiable functions, with g, strictly positive and ei are mutually independent normal noises and can have any arbitrary density distribution. We assume E(c)) = 0 and ((({)))= =1 without loss of generality.
We will start first by showing that if backward model, respects to instantaneous links, exists in the bivariate case then, the data generating mechanism must fulfill the Partial Differential Equation (PDE). Then, following Peters al. [38] and Strobl et al. [50] for defining Restricted SEM on iid, we will define a Temporal Restricted Heteroscedastic Noise model and show its identifiability
Lemma H.9. Assume Causal Markov property, minimality, stationarity, sufficiency and (xt x DEEETT be bivariate time series. Then we have time lagged parents are identifiable, and backward model with respect to instantaneous links can be fit i.e. Vt E T:
x = = (5 + g (i c
(16)
. - ff (
g² ( .
We note h t< t), and let v. (-) and vz(.) be the twice differentiable log densities of X7 and of respectively. For compact notation, define
(x? x.) log
( (2) 2!))
J - ( /B)
log
P2
g2(x!)
i - --(
-1880)))) and
V/
(
G(Z
+ - ((()))
40

---

Assume that g','', fulfill the following PDE for all ([,,)) with G -.) + O.
and g² are twice differentiable, Then, the data generating mechanism must . f,
(
0
0 ((
G
(17)
(0 G 3H)
0²
( 0
(o)
G कt) oā! X (0) x/).
We drop the time-lagged parent in Eq (16) to simplify the notation, since conditioning the history makes it redundant.
Proof of the Lemma. Let's assume we have two temporal causal graph 9 and G' for the bivariate temporal heteroscedastic causal models where the noise distribution could follow any arbitrary distribution.
Disagreement in time lagged relationships. Assume that and G' do not differ in the instantaneous effects. Vi E {1, Pag (t) Pagy((), (t). Hence and Wlog, there is some k > 0 and edge "/ + 2 in G but not in G.. From G' and the Causal Markov property, we have that r/ IL zf S, and ND, are all X that are non-
where S ((=)-,, S 1 S L, {1. ( UDD)) descendants (wrt instantaneous effects) of z²- Applied to G. causal minimality leads to contradiction because _/_x 1L. x; S in g, and the above reasoning shows that it exists subgraph G' of G that is Markovian to the joint distribution of the data.
{=]_
Disagreement in instantaneous parents. Now, let's assume we have a forward model, after conditioning on h, = Pag(< P U 1). 44 E T:
= l,. f² (5)) + -
We want to prove that if backward model
27 = 2 In. f (āt) 9² (5!) -
exists then the PDE in Eq(17) fulfilled.
Our conditioning trick on time lagged parents makes the use of Immer et al. [25] theorem 1 feasible
in our case. We employ the change of variables from { + to [ and the proof will be the same as Immer et al.. Hence, we leverage Theorem of Immer el al. and we conclude that if backward model exists the PDE E17) is verified.
Theorem H.10 (Identifiabilityy of Temporal Restricted Heteroscedastic noise model (TRHNM)).
Assume Causal Markov property, minimality, sufficiency and let (II) be MTS following
Eq(1) where K 1. The graph G is uniquely identified if Vi E [|] dj, x² E Pa (t) and S such that (Pat(t))al) S (Nd (z:) vel).. there exists S = where p(s) > 0
h; Paig(< t) U Pa t) and p (ai, 1 |s, he) do not satisfy PDE of Equation 17, and we call the model that verify this condition, the Temporal Restricted Heteroscedastic noise model.
Proof of the theorem. We will follow the same steps in the proof of theorem 4.2. Let's assume have two temporal causal graph and G' for the multivariate TRHNM. We assume also that
Vi E [|1 ]. Vj x E Pag (1) and S such that (P)(())) CS S (Nd (ai) ver).. there exists
S where p(s) > 0 h, Pa t) U Pa
(< t) and p (x), s, h.)) do not satisfy PDE of
Equation 17.
We will start by showing that time lagged parents are identifiable. Same reasoning in the bivariate case.
Disagreement in time lagged relationships. Assume that G and g' do not differ in the instantaneous effects. Vi E 21; (t) = (l). (1). Hence and Wlog, there is some k > 0 and an edge s/=x + 2 in 9 but not in g'. From g' and the Causal Markov property, we have that x:-. II. x S, where S = (f=)_,, SIS L, {1, 27) U ND.)) 2?}.. and ND, are all X that are non- descendants (wrt instantancous effects) of 2f. Applied to G. causal minimality Ieads to contradiction
41

---

because T/_x +t z; S in G, and the above reasoning shows that it exists subgraph G' of G that is Markovian the joint distribution of the data.
Disagreement on instantaneous parents. Let's now assume we have two temporal causal graph g and 9" such that 9 + G'. According the Propostion 29 in Peters et al [38], for G and g' be two different DAGs over a set of variables V, such that æ; is generated by our TRHNM and satisfies the Markov condition and causal minimality with respect
and g'. Then there are variables Y == Pag (t){x; and S QUY, in G and o + x in g'. 2) S C CNDS
27 V such that for the set Q (ag(t))
if {aī} and S ND Mat) (a)(() is the set of parent variables of a in graph G. ND is the set of non-descendant (wrt instantaneous effects) of z in graph g.
have: 1) 27 + xt
We consider S = s with p(s) > O. Lemma 37 in Peters et al. [38] states that if p(æ:) is generated according to the SEM models as follows:
o = Se (Pag(< (), (a)()))))))) i "d) EV
with corresponding DAG G, then for variable af E V. V. ND then of 11. S. Our TRHNM at can be viewed one specific class of the SEM in the aforementioned equation. Hence, Lemma 37 holds under TRHNM and applying it to 2 renders ef Il (x S) and x es 11 (x)
Using now Lemma 36 in Peters al. [38], and we denote x.!" at Ss and x,"" = xt 8- We have:
4 1.* I, = = (2]7"(n)) +g ("7 (1/1mm (m)
(18)
1. (mc) is at conditioned on h,, S and h, Pag(< 1)
2 f² (xt m.) + g²
62
Pağ(<< 1) in this case, which also equal to h, Pag. 1) U (a t), because we proved identifiability of time lagged parents. Eq(18 8) raise contradiction because, having these forward and backward models imply the verification of PDE 17. But we chose & such that this PDE is not verified hence contradiction. Then, the identifiability of our Temporal Restricted Heteroscedastic noise.
where "I Ib.
H.4 Proof of theorem 4.3
In this section, we want to prove the identifiability of mixture of Temporal causal models either in the case of Temporal Heteroscedastic Gaussian noise, Temporal Restricted Heteroscedastic noise and Homoscedastic NonGaussian noise.
Proof. Let F be a family of identifiable temporal causal models either from TGHNM or TNGNM, F = (po- (-1:, (")))))) that are linearly independent and let Mr be the family of all K-finite mixtures of elements from F, i.e.
Mx (pl2)/=2)) ER (p(w"))o xree (eee.""), P' (1-(-1))))) E FF mp(w") > 0 and )=w")
First, we introduce a result from Yakowitz & Spragins [58] that established a necessary and sufficient condition for the identifiability of finite mixtures of multivariate distributions.
Theorem H.11 (Identifiability of finite mixtures of distributions, Yakowitz. & Spragins [58]). Let F (F(x; a),aa E T x E Re be a finite mixture of distributions. Then F identifiable if and only if F is linearly independent set over the field of real numbers.
We will further assume that it exists two distribution such that V( Fee) covering the space value of random variables (X,, X xet):
K
p(xx(e))) =
Σ
Te(w") (0x (e
(19)
k
p(xx/xc))) =
Σ
me (5" (0 (2ee
42

---

Our objective is to show first that K = K, it exists permutation o and translation function 2: IR2 - 2 (0", 9r) = ( and w" = ((()
Using that J (po- (-1,, 9T) );= are linearly independent and fixing
R
(w") por (B g")
Te
(w0"
Pg. (20a 9*)
K
k
O,per (2 x<1 9")
Σ
bg
and this true V(a,,, xxee) covering the space value of the random variable = (,,, <t) By
using theorem H.11, we can conclude that: K K and exists a permutation a such that: (0",9") (") = (0- (r) (r(r) and A4 E TT: mT(W"") ((("))
To proof our identifiability as defined in definition 4.1, we still need to prove that w" ((a)()).. We have V/ E T vT(w")= (e(""), we take two indices r, E [[1: K]]
Te(ww)) = m((30(r) TT(w²)= (e)
(20)
To handle time varying weights identifiability, we will consider ratios of mixture weights:
explw] (+wo)) 2i exp(11 (w) ++ad))
TT TT(W")
Tel (w=)
expl (w] (twr)
E exp(ur ( )+4a))
expl ت ,t, utug (1))
Te (())
Σx exp (wارسس) (w 16
(wo(s))
exp (up (ن (+/
Ex expl (ما 1
By Equation 20:
Te (eww))
Te (r))
2e (w+)
Te (x)
<> exp [['] - wf) 1 (w] wū)] = exp [(w)(t] wg((s) 1 (w)(r) - w =
* T (Wi wf) t + (un wi)) (wg ((0) t ++ (ugit) w/ w6(())
As consequence of the last equation, we have for all the indices:
wi w = w/ (r)) (s))
w (3)
w/ w w/
wh
w w/ (() w wi w((() A
wi w (r) w wi (+) Ao
Hence it exists a translation function @ _ + R2a such that Vr []] K]]:
w = (())
Hence our mixture of temporal causal models is identifiable as defined in definition 4.1.
43

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
Yes the main claims are clear in the abstract, introduction and in the whole paper
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Yes, we have section in appendix D in which we talk about the limitation of our work.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (c.g.. independence assumptions, noiseless settings, model sell-speciiiiatinn asymptotic approximations only holding locally). ). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex, system might not be used reliably to provide closed captions for online lectures because it fails handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and faimesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
44

---

Answer: [Yes]
Justification: Our work offers three principal theoretical contributions. All assumptions and complete proofs appear in Appendix H. Due space constraint, we state two theorems in the main text and the third in the appendix, with each theorem explicitly citing its underlying assumptions.
Guidelines:
The answer NA means that the paper does not include theoretical results. All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theorems.
The proof's can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Yes, we provide in the main text figures that illustrate the exact architecture used for our model. In the appendices, we provide the exact typer-parameter needed to reproduce our results in all the presented experiments. Our code is also provided in the supplementary material.
Justification:
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g.. in the case of large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make clear how to reproduce that algorithm. (b) the contribution is primarily a new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (c.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibility.
45

---

In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: Our code, with all the instruction reproduce our results, is provided in the supplementary materials
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed run reproduce the results. See the NeurIPS code and data submission guidelines (https:
//niss. cctpublcc/aaaeecccccccccccccccct /ece//giiee /wwedeeee/oeebmmeeeeeee cebbbbisssionii ssisnnolicyyyycocoo for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) is recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: All the detailed needed presented in Appendix E.3.
code and reproduce results
Guidelines:
The answer NA means that the paper does include experiments.
The experimental setting should be presented in the core the paper to a level detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, appendix, as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: In the main text, we presented only average scores for the different metric, But, in the appendix we provide erro bars for our model and all the other baselines.
Guidelines:
The answer NA means that the paper does not include experiments.
46

---

The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables figures symmetric error bars that would yield results that are out range (e.g. negative rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We provide time complexity and running time in appendices F.3.2 and F.1
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs as well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https: httpss/////u..../.ee/ cc Spppblic/tthissuiddeeepoppppppee
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special cireumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is a special consid- eration due to laws or regulations their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: In our conclusion, talk about paper's broader impact.
Guidelines:
47

---

The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g.. gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, possible, include a URL
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
48

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. cww.coo////oo/// datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. If this information is available online, the authors encouraged to reach out the asset's s creators.
13. New assets
Question: Are assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification: We presented all the details about model: parameters, architecture, data. Also the detailed proof' are provided
Guidelines:
The answer NA means that the paper does not release assets.
Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, curation, other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the paper does involve crowdsourcing nor research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
49

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology. scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification:
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our LLM policy (https: :ottpa:/on /tt////eissp.///ee cc/Conferencees22225 /LLM) for what should or should be described.
50


Reviews:
[{"review_id": "4u9A0rBQkY", "reviewer": "Reviewer_yhgx", "summary": "The paper presents FANTOM, a method for causal discovery in multivariate time series that may switch among an unknown number of regimes and exhibit non-Gaussian or heteroscedastic noise. FANTOM combines a Bayesian EM procedure with conditional normalising flows. The authors prove identifiability for both single-regime and multi-regime settings. Experiments on synthetic data and two real-world benchmarks show that FANTOM recovers regime boundaries and regime-specific causal graphs more accurately than several existing methods.", "questions": "Can the authors provide additional information mentioned in weakness?", "limitations": "- Scalability / Efficiency – report wall-clock training time, GPU hours, and memory for the 40-variable experiment, and discuss how cost grows with more variables, longer horizons, or >3 regimes.\n- Hyper-parameter sensitivity – show how results change with the initial regime count, pruning threshold, and normalising-flow depth; document failure cases (e.g., over-pruning).\n- Risk of spurious causality – recommend practical safeguards (expert vetting, interventional checks) before acting on the learned graphs, especially in medical or financial domains.\n- Privacy & misuse – clarify whether the method can train on anonymised / differentially-private data and outline precautions when applied to sensitive EEG or industrial logs.", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "## Strengths \n- Rigorous methodology. The Bayesian EM routine cleanly separates regime assignment (E-step) and DAG/noise learning (M-step), with conditional normalising flows providing flexible likelihoods.\n- Provable guarantees. The paper supplies identifiability theorems for both single- and multi-regime settings, lending theoretical credibility to the approach.\n- Provable guarantees. The paper supplies identifiability theorems for both single- and multi-regime settings, lending theoretical credibility to the approach.\n\n## Weaknesses\n- Ablation gaps. There is no study isolating the importance of heteroscedastic modelling, flow architecture, or the pruning step, so robustness remains uncertain.\n- Compute cost unreported. Training time and GPU usage for the flow-based M-step are omitted, making efficiency hard to judge.", "was_revised": true, "final_justification": "During the rebuttal, most of the concerns were resolved. Although I still have concerns about the risk of spurious causality, especially in real-world scenarios, I will keep my positive comment on this work.", "rebuttal": "We thank the reviewer for their valuable feedback and suggestions. We are pleased they find our paper rigorous and supported by theoretical results for the single and multi regime settings. Below, we address their specific comments and suggestions. We hope that our responses are satisfactory so that the rating can yet be increased.\n\n>Ablation gaps. There is no study isolating the importance of heteroscedastic modelling, flow architecture, or the pruning step, so robustness remains uncertain.\n\nWe thank the reviewer for the helpful suggestions. To assess the contribution of heteroscedastic modelling and the flow architecture, we ran an ablation on MTS comparing three FANTOM variants: (i) Gaussian output (no flow, FANTOM-Gauss), (ii) a simple coupling‑spline flow (FANTOM-Spline), and (iii) the full model with a conditional normalizing flow (CNF, FANTOM). In the single‑regime setting, removing or simplifying the CNF consistently degrades performance. In the multi‑regime setting, only the CNF variant remains stable and converges to the two ground truth regimes; the alternatives collapse to one regime. Results are reported in the table below: \n\n**Table 1: Importance of heteroscedastic modelling and flow architecture. - means the method collapses to 1 regime.**\n\n|   |               |      | Heteroscedastic |     |      |          |\n|:-:|---------------|------|:---------------:|:---:|:----:|----------|\n|   |               | Inst |                 | Lag |      | Regime   |\n| K | model         | SHD  | F1              | SHD | F1   | Accuracy |\n| 1 | FANTOM-Gauss  | 39   | 0.0             | 44  | 42.5 | only 1 regime        |\n|   | FANTOM-Spline | 17   | 0.0             | 18  | 0    | only 1 regime         |\n|   | FANTOM        | 0    | 100             | 0   | 100  | only 1 regime        |\n| 2 | FANTOM-Gauss  | -    | -               | -   | -    | -        |\n|   | FANTOM-Spline | -    | -               | -   | -    | -        |\n|   | FANTOM        | 1    | 97.8            | 1   | 98.9 | 98.6     |\n\nWe evaluated robustness to the pruning step, minimum regime duration (ζ), and regime initialization. Across these variations, FANTOM’s predictive performance remains stable. We find that changing (window size, ζ) primarily affects runtime, especially at long horizons, without materially impacting graph and regime learning performances. To show efficiency under long-horizon conditions, we used regimes with lengths [2000, 3000, 2000, 2500].\n\n**Table 2: Performance of FANTOM with varying hyperparameters (window size and ζ) on 10‑node graphs with 2, 3, and 4 regimes.**\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Number | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | of iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   10  | 4    | 17'8s   |           5           | 91.2 |  13 | 85.7 |  96.5  |\n|    (1500,1200)    | 2 |   10  | 4    | 10'12s  |           1           | 97.9 |  9  | 89.8 |  95.1  |\n|     (1000,900)    | 3 |   10  | 4    | 30'33s  |           5           | 94.2 |  8  | 94.4 |  96.4  |\n|    (1500,1200)    | 3 |   10  | 4    | 16'2s   |           3           | 96.7 |  10 | 91.9 |  91.8  |\n|     (1000,900)    | 4 |   10  | 4    | 40'2s   |           5           | 95.8 |  12 | 92.7 |  92.9  |\n|    (1500,1200)    | 4 |   10  | 4    | 25'38s  |           3           | 97.5 |  16 | 90.5 |  91.5  |\n\n**Table 3: Performance of FANTOM with varying hyperparameters (window size and ζ) on 40‑node graphs with 2, 3, and 4 regimes.**\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Number | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | of iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   40  | 4    | 22'15s  |           33          | 81.9 |  27 | 91.9 |   100  |\n|    (1500,1200)    | 2 |   40  | 4    | 10'20s  |           35          | 83.8 |  25 | 91.7 |   100  |\n|     (1000,900)    | 3 |   40  | 4    | 34'25s  |           49          | 85.6 |  58 | 88.6 |  99.8  |\n|    (1500,1200)    | 3 |   40  | 4    | 15'52s  |           43          | 87.3 |  72 | 87.3 |  99.9  |\n|     (1000,900)    | 4 |   40  | 5    | 66'30s  |           34          | 86.5 |  53 | 91.9 |  99.7  |\n|    (1500,1200)    | 4 |   40  | 4    | 21'38s  |           32          | 87.1 |  54 | 91.8 |  99.6  |\n\nFrom the table, FANTOM scales to 40 nodes and converges to four regimes under different initializations. With initialization = 1500 and ζ =1200, it starts with **six initial regimes** and converges smoothly to **four with 99.6% accuracy** in 21 min 38 s. With initialization = 1000 and ζ =900, it starts with **nine initial regimes** and converges to **four with 99.7% accuracy** in 66 min 30 s. All rebuttal experiments were run on a Tesla V100‑SXM2 (26 GB), whereas the results in the paper were obtained on an NVIDIA GTX TITAN (12 GB). These results indicate that FANTOM scales without requiring large GPU memory.\n\n>Scalability / Efficiency – report wall-clock training time, GPU hours, and memory for the 40-variable experiment, and discuss how cost grows with more variables, longer horizons, or >3 regimes.\n\n**Table 4: Runtime for different number of nodes and under different number of regimes.**\n\n|       |      |         | (window, $\\zeta$)=(1500,1200) |         |      |         |\n|:-----:|------|---------|:-----------:|:-------:|:----:|:-------:|\n|       | K=2  |         | K=3         |         | K=4  |         |\n| Nodes | Number | Running | Number | Running | Number | Running |\n|       | of iter | Time    | of iter        | Time    | of iter | Time    |\n| 10    | 4    | 10'12s  | 4           | 16'2s   | 5    | 25'38s  |\n| 20    | 4    | 9'50s   | 4           | 15'3s   | 5    | 21'10s  |\n|   40  | 4    | 10'20s  | 4           | 15'52   | 5    | 21'38s  |\n\nWe thank the reviewer for these suggestions. To test scalability, we ran FANTOM on a **40-node dataset with four regimes and long horizons** (lengths 2000, 3000, 2000, 2500). As expected, additional regimes and longer horizons increase runtime. With **40 nodes and two regimes**, training finishes in **10 min 20 s**, whereas the **four-regime setting requires 21 min 38 s**.\n\n>Hyper-parameter sensitivity – show how results change with the initial regime count, pruning threshold, and normalising-flow depth; document failure cases (e.g., over-pruning).\n\nWe thank the reviewer for suggesting enriching experiments. We showed in the table above (Table 2 and 3) that FANTOM is robust to the initialization parameter (**window size that impacts directly the initial regime count and $\\zeta$ the pruning threshold**), e.g. in the case of **40 nodes with 4 regimes**, with window initialization = 1500 and ζ =1200, it starts with **six initial regimes** and converges smoothly to four in 21 min 38 s. With initialization = 1000 and ζ =900, it starts with **nine initial regimes** and converges to four with 99.6% accuracy in 66 min 30 s. \n\nWe also conducted new experiments to explore the sensitivity of FANTOM to normalizing flow depth. We pick K=2 regimes and d=10 nodes.\n\n**Table 4: FANTOM performance with different normalizing flow depth.**\n\n| K=2, d=10               |      | Heteroscedastic |     |      |          |\n|-------------------------|------|:---------------:|:---:|:----:|----------|\n|                         | Inst |                 | Lag |      | Regime   |\n| Normalizing flows depth | SHD  | F1              | SHD | F1   | Accuracy |\n| 16                      | 0    | 100             | 1   | 98.9 | 97.1     |\n| 32                      | 0    | 100             | 0   | 100  | 96.6     |\n| 64                      | 0    | 100             | 1   | 98.9 | 94.8     |\n| 128                     | 0    | 100             | 0   | 100  | 97.2     |\n\nFANTOM shows robustness to the choice of the depth of NF, regardless of the number of bins, it achieves the same range of performances.\n\nRegarding over pruning, we want to thank you for this suggestion, we will add a discussion in our manuscript, where we will clarify that over pruning will affect the performance of FANTOM.\n\n>Risk of spurious causality – recommend practical safeguards (expert vetting, interventional checks) before acting on the learned graphs, especially in medical or financial domains.\n\nWe thank the reviewer for this suggestion. We will add the following paragraph to our updated version:\nSpurious causality is a practical risk, so we advocate using the learned graphs as decision support rather than as ground truth. Before acting on them, especially in medical or financial settings, practitioners should adopt safeguards: (i) expert vetting of proposed edges and mechanisms; (ii) robustness checks (e.g., stability under resampling or perturbations, alternative specifications); and (iii) interventional validation when feasible.\n\n>Privacy & misuse – clarify whether the method can train on anonymised / differentially-private data and outline precautions when applied to sensitive EEG or industrial logs.\n\nWe will add a paragraph to our updated paper to clarify that our method trains on anonymised data and clarifies the precautions when applied to EEG data."}, {"review_id": "BPwqkxpx64", "reviewer": "Reviewer_4BQU", "summary": "The paper introduces FANTOM, a causal discovery framework that handles multi-regime multivariate time series (MTS), homoscedastic non-Gaussian and heteroscedastic noise. FANTOM is capable of learning the regime indices even when its unknown. The authors also provide identifiability guarantees for both heterscedastic and homoscedastic noise settings. Given the MTS multi-regime data, FANTOM learns the regimes and the causal graph via Bayesian expectation maximization (BEM) optimizing evidence lower bound (ELBO).", "questions": "1. In assumption 2.2, the MTS system is assumed to be stationary for $\\zeta$ consecutive steps. Is the $\\zeta$ assumed to known apriori? Is the initial regime assignments based on the knowledge of $\\zeta$? \n2. How would the initialization of $N_w$ guaranteed to be greater than $K$ when $K$ unknown? \n3. Given that FANTOM runs until a maximum number of iterations are reached (mentioned in line 176), how is the maximum number of iterations chosen? \n4. How does FANTOM compare to the baseline (specifically CASTOR) in terms of computational cost? How much does modeling heteroscedastic and non-Gaussian noise affect the computational cost when compared to CASTOR? \n\nIf the authors can satisfactorily address the questions above, I would be inclined to raise my overall recommendation.\n\n(Minor comment) It would help to define the function $H$ after introducing it in equation (2).", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Quality: The paper is technically sound, the main claims of the paper (causal discovery under MTS multi regime setting) is supported by the theoretical results established in section 4, and experimentally in section 5 where the proposed method shows improved performance over the baselines in certain settings. \n\nClarity: The presentation is sound. However, the paper is dense with technical details making it a bit challenging to read. \n\nSignificance: The paper addresses the difficult task of causal discovery from time series data multiple regime changes and heteroscedastic noise. The experimental results shows that the the proposed method, FANTOM, outperforms the baselines in the heteroscedastic setting, and in the real world data experiments. However, the proposed work assumes that the minimum regime width is assumed to known beforehand, which may not always be practical.  \n\nOriginality: The proposed work extends the results of Rahmani and Frossard (2025) [1] (CASTOR), both follow a similar framework with FANTOM allowing for heteroscedastic noise as well as non-Gaussian noise compared to CASTOR. \n\n[1] Rahmani, Abdellah, and Pascal Frossard. \"Causal Temporal Regime Structure Learning\".", "was_revised": true, "final_justification": "All the points I raised were addressed during the rebuttal.", "rebuttal": "We thank the reviewer for their valuable feedback and interesting questions. We are pleased that they find our paper technically sound, supported by theoretical results, deem the task we address difficult and our performance superior to existing baselines. Below, we address their specific comments and questions. We hope that our responses are satisfactory so that the rating can yet be increased.\n\n>Originality: The proposed work extends the results of Rahmani and Frossard (2025) [1] (CASTOR)\n\nWe thank the reviewer for their feedback, It s true that FANTOM and CASTOR uses the same time varying prior for the latent factors  as well as the initialization trick, however FANTOM first employs Conditional Normalizing flows (CNF) to handle complex noise distributions (heteroscedastic noises and Non-Gaussian noises). In such cases CASTOR fails to detect the regimes as well as learning the graphs (**see Table 1 in our paper**). Moreover FANTOM uses Bayesian structure learning instead of structure learning (used by CASTOR) and this with CNF helps solving **the problem of data standardization** faced by CASTOR, Dynotears and other models. We conducted new experiments to show this phenomena, here are the results:\n\n**Table 1: Robustness of FANTOM to data standardization**\n| data         | K | model  | Inst SHD | Inst F1 | Lag SHD | Lag F1 | Reg Acc |\n|--------------|---|--------|----------|---------|---------|--------|---------|\n| Raw          | 2 | CASTOR with regime labels| 90       | 26.9    | 28      | 29.7   | x       |\n| standardized | 2 | CASTOR with regime labels| 28       | 0.0     | 33      | 0.0    | x    |\n| Raw          | 2 | FANTOM | 4        | 93.3    | 7       | 90.9   | 97.1    |\n| standardized | 2 | FANTOM | 5        | 91.5    | 7       | 90.9   | 91.1    |\n| Raw          | 3 | CASTOR with regime labels| 135      | 24.1    | 135     | 29.6   | x       |\n| standardized | 3 | CASTOR with regime labels| 37       | 0.0     | 48      | 0.0    | x       |\n| Raw          | 3 | FANTOM | 3        | 96.1    | 10      | 91.5   | 92.2    |\n| standardized | 3 | FANTOM | 2        | 97.3    | 9       | 92.3   | 92.2    |\n\nFrom the Table, FANTOM achieves the same results on raw and standardized data while CASTOR fails in both cases even when we give it access to the ground truth regime labels. In the case of standardized data, CASTOR predicts adjacency matrices full of zeros due to its incapability of handling such scenarios. We also want to highlight that in real world scenarios where the noise distribution is complex CASTOR fails to learn meaningful graphs (In both our presented real world data sets CASTOR collapses to one regime).\n\n> In assumption 2.2, the MTS system is assumed to be stationary for ζ consecutive steps. Is the ζ assumed to known apriori? Is the initial regime assignments based on the knowledge of ζ?\n\n* **Minimum regime duration (ζ).** As in CASTOR, we treat ζ as a hyperparameter that sets the minimum length of a regime. When domain knowledge is available, it should guide this choice. For example, epileptic seizures typically last minutes to hours; setting ζ≈3 or 4 minutes will allow FANTOM to detect the meaningful regimes. In climate applications, seasons span months; setting ζ≈1 month is appropriate.\n\n* **Initialization and pruning.** FANTOM uses two hyperparameters: (i) the window initialization that gives rise to initial regimes, and (ii) ζ which cancels short regimes during the BEM training loop. We choose the window slightly larger than ζ to avoid early collapse to a single regime (e.g., window = 1000, ζ = 900 in our experiments).\n\n* **No domain knowledge available.** If no prior knowledge is available, taking a reasonable small (window, ζ) will help learning the meaningful regime smoothly. This choice will not affect the performance of FANTOM but it will affect running time. \n\nHere we conducted some ablation studies that show that FANTOM is robust to the (window size, ζ ) choices, the regime length are in this order [2000,3000,2000,2500], we tried to choose long horizons in order to show that FANTOM is scalable to long horizons regardless of the choice of the hyper parameters:\n\n**Table 2: Performance of FANTOM with varying hyperparameters (window size and ζ) on 10‑node graphs with 2, 3, and 4 regimes.**\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Number | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | of iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   10  | 4    | 17'8s   |           5           | 91.2 |  13 | 85.7 |  96.5  |\n|    (1500,1200)    | 2 |   10  | 4    | 10'12s  |           1           | 97.9 |  9  | 89.8 |  95.1  |\n|     (1000,900)    | 3 |   10  | 4    | 30'33s  |           5           | 94.2 |  8  | 94.4 |  96.4  |\n|    (1500,1200)    | 3 |   10  | 4    | 16'2s   |           3           | 96.7 |  10 | 91.9 |  91.8  |\n|     (1000,900)    | 4 |   10  | 4    | 40'2s   |           5           | 95.8 |  12 | 92.7 |  92.9  |\n|    (1500,1200)    | 4 |   10  | 4    | 25'38s  |           3           | 97.5 |  16 | 90.5 |  91.5  |\n\n**Table 3: Performance of FANTOM with varying hyperparameters (window size and ζ) on 40‑node graphs with 2, 3, and 4 regimes.**\n\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Number | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | of iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   40  | 4    | 22'15s  |           33          | 81.9 |  27 | 91.9 |   100  |\n|    (1500,1200)    | 2 |   40  | 4    | 10'20s  |           35          | 83.8 |  25 | 91.7 |   100  |\n|     (1000,900)    | 3 |   40  | 4    | 34'25s  |           49          | 85.6 |  58 | 88.6 |  99.8  |\n|    (1500,1200)    | 3 |   40  | 4    | 15'52s  |           43          | 87.3 |  72 | 87.3 |  99.9  |\n|     (1000,900)    | 4 |   40  | 5    | 66'30s  |           34          | 86.5 |  53 | 91.9 |  99.7  |\n|    (1500,1200)    | 4 |   40  | 4    | 21'38s  |           32          | 87.1 |  54 | 91.8 |  99.6  |\n\nIn the **four‑regime setting** with ζ=900 and window = 1000, FANTOM initializes with **nine regimes: [1000,1000,1000,1000,1000,1000,1000,1000,1500]** and converges smoothly to the **four ground‑truth regimes: [2000,3000,2000,2500], with 99.7%**, this configuration requires 66 min 30 s. With initialization = 1500 and ζ =1200, it starts with **six initial regimes: [1500, 1500, 1500, 1500, 1500, 2000]** and converges smoothly to **four: [2000,3000,2000,2500], with 99.6%** accuracy in 21 min 38 s. Overall, Regime and Graph learning performances has not been affected by hyperparameters choices.\n\n>How would the initialization of Nw guaranteed to be greater than K when K unknown?\n\nWe thank the reviewer for the question. As discussed earlier, our method has two hyperparameters: the window size and ζ. By choosing them carefully, you can ensure $N_w \\geq K$. In practice, set (window,ζ) as small as possible while keeping enough data per window will both satisfy $N_w \\geq K$ and help FANTOM converge reliably to the true regimes. In our earlier ablations, any (window,ζ) smaller than (2000,1900) and larger than ~(800,700) will work, since the Temporal GNN and CNF need few data to learn the graph effectively. We also showed that initialization does not affect FANTOM’s performance; it only changes runtime (see the previous table 2 and 3).\n\n>Given that FANTOM runs until a maximum number of iterations are reached (mentioned in line 176), how is the maximum number of iterations chosen?\n\nThe number of iterations is a hyperparameter. In general, increasing it does not degrade FANTOM’s performance: once FANTOM reaches the optimum, it stabilizes, and additional iterations leave the solution unchanged.\n\n>How does FANTOM compare to the baseline (specifically CASTOR) in terms of computational cost? How much does modeling heteroscedastic and non-Gaussian noise affect the computational cost when compared to CASTOR?\n\nWe thank the reviewer for the suggestion. We ran additional experiments comparing the computational cost of CASTOR and FANTOM in the case of heteroscedastic noise. We set the regime lengths to [2000, 3000, 2000], with ζ = 1200, window = 1500, and d = 10. The results are as follows:\n\n| K=3, d=10 |         |\n|-----------|---------|\n| model     | runtime |\n| CASTOR    | 14'59   |\n| FANTOM    | 16'2s   |\n\nCASTOR runs ~2 minutes faster than FANTOM, but it collapses to a single regime at the first iteration which makes it model the MTS as only one regime and reduces the computational cost in the next iterations. Moreover, CASTOR cannot learn the ground-truth regimes and graphs in heteroscedastic settings, where FANTOM performs significantly better.\n\nIn terms of memory: All results in the paper were obtained on an NVIDIA GTX TITAN with 12 GB RAM, indicating that FANTOM has modest memory requirements.\n\n>Reviewer comment: (Minor comment) It would help to define the function H after introducing it in equation (2).\n\nWe thank the reviewer for their comment. We will introduce the entropy H with its mathematical formulation in our updated version."}, {"review_id": "fYsiI2ROi7", "reviewer": "Reviewer_KqQF", "summary": "The authors present a novel method for causal discovery in multivariate time series. It can cope with heteroscedastic models (for arbitrary noise distributions), lagged and instantaneous effects as well as stationarity and non-stationarity. The method is based on a Bayesian Expectation Maximization algorithm which learns both the different regimes of the time series and the underlying causal structure(s). Moreover, identifiability results, a detailed empirical comparison with further causal discovery methods, and real data examples are provided.", "questions": "-\tDid you investigate whether the method is robust to scaling? E.g., standardizing the data, which is often a problem in likelihood-based methods (Reisach et al. *Beware of the simulated dag! causal discovery benchmarks may be easy to game*, 2021).\n- In Eq. (9) (Appendix C) should small $k$ be $\\tau$?\n-\tCan you elaborate a bit more on how you carried out the maximization of the ELBO in the M-step?\n-\tWhat about $L > 1$? Line 323 references Appendix F.1 for such cases, but I could not find results in that section.", "limitations": "Yes (except scaling; see questions).", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "Strengths:\n\n-\tThe authors fill a gap in the literature developing a causal structure learning algorithm for time series that can handle (non-) stationary heteroscedastic and non-Gaussian settings.\n-\tThe authors back up their method with convincing empirical (experiments) and theoretical results (identifiability, building upon proofs for the stationary setting).\n-\tA detailed explanation of the methodology is provided, including sufficient intuition, and supported by clear and illustrative figures. \n\nWeaknesses (minor):\n\n-\tMore details on how exactly the “joint” maximization of Eq. (6) in 3.3 is carried out could be incorporated.\n-\tWithin the main text, a more thorough explanation of Theorem 4.3 would be appreciated.\n-\tPlease explain why the entropy of $p(z_t|x_t,x_{<t})$ does not show up in Eq (6).\n-\tThe writing can be improved. For instance, the authors alternate between “non-Gaussian”, “non Gaussian”, “log likelihood”, “log-likelihood” and “non-stationary” and “non stationary”, … There should be consistency in this regard. Also usage of “can not” and “definition H.1” instead of “Definition”, …\n-\tThe placement of Table 2 should be changed to satisfy the formatting requirements.", "was_revised": true, "final_justification": "Since the authors could answer my questions in a satisfactory manner, I increased my score.", "rebuttal": "We thank the reviewer for the valuable feedback. The reviewer observes that the proposed method fills an important gap in the literature, is rigorously benchmarked, and delivers substantial empirical and theoretical contributions. They also acknowledge the clarity of the exposition, including intuitive explanations and illustrative figures. Below, we address their specific comments and questions. We hope that our responses are satisfactory so that the rating can yet be increased.\n\n>More details on how exactly the “joint” maximization of Eq. (6) in 3.3 is carried out could be incorporated.\n\nWe thank the reviewer for the suggestion and the question about computing Eq. (6). We will add the following discussion to the revised manuscript: Eq. (6) splits into two **independent maximization problems**: (i) graph learning and (ii) regime alignment. The colors in Eq. (6) marks the terms belonging to each subproblem. **We solve these two maximizations separately.** We would like to mention that the regime‑alignment step is computationally light and runs in only a few seconds.\n\n>Can you elaborate a bit more on how you carried out the maximization of the ELBO in the M-step (Eq 6)?\n\nAfter estimating the posterior probabilities $\\beta^r_t$ in the E‑step, we apply hard thresholding, which assigns each timestamp to a single regime and makes $\\\\beta^r_t$ binary ($\\\\forall t \\in \\\\mathcal{T}: \\\\sum\\_{r=1}^{N_w}\\beta^r_t =1$ and $\\\\beta^r_t \\\\in \\\\{0,1\\\\}$). We then optimize the two objectives separately. \n* For regime alignment, we use SGD to learn the parameter $\\omega^r$ in a softmax of an affine function of the time index. \n* For graph learning, we use an augmented Lagrangian to enforce the acyclicity prior and optimize with Adam, following Geffner et al. (2022) [1] and Gong et al. [2]. The procedure alternates between: (i) optimizing the graph objective with fixed penalty parameters $\\rho$ and $\\alpha$ (defined in our graph prior line 240 and 247) for about 2000 steps, and (ii) updating $\\rho$ and $\\alpha$. We repeat (i) (ii) until convergence or until hitting the maximum number of outer updates (typically 3).\n\n[1] Geffner, Tomas, et al. \"Deep End-to-end Causal Inference.\" Transactions on Machine Learning Research.\n\n[2] Gong, Wenbo, et al. \"Rhino: Deep Causal Temporal Relationship Learning with History-dependent Noise.\" The Eleventh International Conference on Learning Representations.\n\n>Within the main text, a more thorough explanation of Theorem 4.3 would be appreciated.\n\nWe thank the reviewer for this constructive suggestion, which improves the quality of our paper. We will add the following explanation to the revised manuscript:\n\n*The theorem states that for identifiable temporal causal models, if we choose K linearly independent component models, then their mixture with time‑varying weights is identifiable, up to a permutation of the mixture components and up to a translation (additive shift) of the time‑varying weights. This establishes identifiability for mixtures with time‑varying weights and generalizes the standard case with static weights.*\n\n>Please explain why the entropy of $p(z_t|x_t,x_{<t})$ does not show up in Eq (6).\n\nWe thank the reviewer for pointing out this question, we will add a section in our appendices in which we will clarify the mathematical derivation and the intuition behind it.\n\nWe perform the maximization of ELBO presented in proposition 3.1, using a BEM procedure, where we alternate between E-step (updating posterior probabilities while fixing all the parameters $\\Theta = \\\\{\\theta^r,\\phi^r, \\omega^r\\\\}\\_{r=1}\\^{N_w}$) and M step (Updating the parameters while using the posteriors learned in the E-step), we can summarize the process as follows:\n\n* In the Estep: we learn $\\beta\\_t^r = p(z_t|x_t,x_{<t}, \\Theta\\_{old})$, where $\\Theta\\_{old}$ is $\\Theta$ of the previous iteration.\n* In the M-step: we fix the learned posterior probabilities to the values $\\beta\\_t^r $ and we update $\\Theta$.\n* By fixing the posterior probabilities in the M-step, the entropy of these probabilities is a constant, which allow us to discard it.\n\nThe detail derivation from Eq.(2) to Eq.(6) is the following:\n\n$\\\\mathrm{ELBO}(\\Theta) = \\\\sum_{t=1}^{|T|}\\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})} \\Big[\\\\mathbb{E}\\_{p(z_t \\mid x_t, x_{<t})}\\big(\\log p_{\\theta_{z_t}}(x_t \\mid x_{<t}, \\mathcal{G}^{z_t})+ \\log p(z_t) \\big)+ H\\big(p(z_t \\mid x_t, x_{<t})\\big)\\Big]+ \\sum\\_{r=1}^{N_w}\n\\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q_{\\phi_r}(\\mathcal{G}^r)\\right)$\n \n$= \\sum_{t=1}^{|T|} \\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})}\\Big[\\sum_{z_t} p(z_t \\mid x_t, x_{<t}, \\Theta_{\\mathrm{old}})\\big(\\log p_{\\theta_{z_t}}(x_t \\mid x_{<t}, \\mathcal{G}^{z_t}) + \\log p(z_t)\\big) \\Big] + \\\\underbrace{H\\left(p(z_t \\mid x_t, x_{<t},\\Theta_{\\mathrm{old}})\\right)}_{\\text{Cte}} + \\\\sum\\_{r=1}^{N_w} \\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q\\_{\\phi_r}(\\mathcal{G}^r)\\right)$\n\nWe replace $p(z_t \\mid x_t, x_{<t}, \\Theta_{\\mathrm{old}})$ by $\\beta_t^{r}$, $p(z_t)$ by our prior choice $\\pi_t(\\omega^r)$, and we discard the entropy of the posterior because it is a constant in these steps. Hence, we got the Eq (6):\n\n$=\\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})}\\Big[\\\\sum\\_{t=1}^{|T|}\\sum_{r=1}^{N_w}\\beta_t^r \\big(\\log p_{\\theta_r}(x_t \\mid x_{<t},\\mathcal{G}^r) + \\log \\pi_t(\\omega^r)\\big)\\Big]  + \\sum_{r=1}^{N_w} \\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q_{\\phi_r}(\\mathcal{G}^r)\\right).$ Eq(6)\n\nAfter the maximization of Eq.(6), FANTOM update the posteriors $p(z_t \\mid x_t, x_{<t}, \\theta_{\\mathrm{old}})$ following Eq.(7).\n\n>The writing can be improved. Format\n\nWe thank the reviewer for their careful reading, their remarks and suggestions that enhance the quality of our paper. We will correct all these inconsistencies in our updated version. Regarding table 2 we will correct the format.\n\n>Did you investigate whether the method is robust to scaling? E.g., standardizing the data, which is often a problem in likelihood-based methods (Reisach et al. Beware of the simulated dag! causal discovery benchmarks may be easy to game, 2021).\n\nThank you for raising this important and challenging question about causal discovery with optimization‑based models. We investigated it. **FANTOM is robust to data standardisation** because it uses **Bayesian structure learning** to estimate a distribution over plausible graphs, and **conditional normalizing flows** to model complex distributions. In contrast, other models like CASTOR are sensitive to standardization and fails even when we provide ground truth regime labels. The table below summarizes the results.\n\n**Table 1: Robustness of FANTOM to data standardization**\n| data         | K | model  | Inst SHD | Inst F1 | Lag SHD | Lag F1 | Reg Acc |\n|--------------|---|--------|----------|---------|---------|--------|---------|\n| Raw          | 2 | CASTOR with regime labels| 90       | 26.9    | 28      | 29.7   | x       |\n| standardized | 2 | CASTOR with regime labels| 28       | 0.0     | 33      | 0.0    | x       |\n| Raw          | 2 | FANTOM | 4        | 93.3    | 7       | 90.9   | 97.1    |\n| standardized | 2 | FANTOM | 5        | 91.5    | 7       | 90.9   | 91.1    |\n| Raw          | 3 | CASTOR with regime labels| 135      | 24.1    | 135     | 29.6   | x       |\n| standardized | 3 | CASTOR with regime labels| 37       | 0.0     | 48      | 0.0    | x       |\n| Raw          | 3 | FANTOM | 3        | 96.1    | 10      | 91.5   | 92.2    |\n| standardized | 3 | FANTOM | 2        | 97.3    | 9       | 92.3   | 92.2    |\n\n>In Eq. (9) (Appendix C) should small $k$ be $\\tau$?\n\nWe thank the reviewer for pointing out this typo. In Eq. (9), τ should replace k. We will correct this in the revised manuscript.\n\n>What about L>1? Line 323 references Appendix F.1 for such cases, but I could not find results in that section.\n\nWe want to thank the reviewer for this remark, we acknowledge that due to the extensive experiments we conducted: *ablation studies of the hyperparameters, 6 tables for non-Gaussian and 6 tables for heteroscedastic settings.* We missed adding the results for L>1. We apologize for this oversight and will add them in the revised version. \n\nWe also want to clarify that in real world scenarios we used FANTOM with L=8 for causal chambers and Epilepsy experiments. Here are the results on synthetic data with heteroscedastic noise and L=2:\n\n**Table 2: Performance of FANTOM on $𝐿=2$ time series under varying node counts $d$ 10,20,40 and regime settings $K$ 2 and 3.**\n| K | d  | Inst SHD | Inst F1 | Lag SHD | Lag F1 | Reg Acc |\n|---|----|-----------|----------|---------|--------|----------|\n| 2 | 10 | 1         | 97.7     | 0       | 100.0  | 94.3     |\n| 2 | 20 | 14        | 87.0     | 4       | 98.7   | 99.0     |\n| 2 | 40 | 15        | 94.5     | 3       | 99.5   | 98.9     |\n| 3 | 10 | 1         | 98.5     | 1       | 99.5   | 91.5     |\n| 3 | 20 | 5         | 97.5     | 11      | 97.4   | 97.9     |\n| 3 | 40 | 43        | 88.6     | 3       | 99.0   | 99.0     |\n\nFANTOM maintains similar performance with a larger lag ($L=2$). It achieves an F1 score above 87% on instantaneous links for 10, 20, and 40 nodes in MTS with 2 or 3 regimes. For time‑lagged links, we evaluate using the summary graph of lagged relations; the table shows that FANTOM detects them effectively."}, {"review_id": "sKruuHZ2hL", "reviewer": "Reviewer_FQMP", "summary": "This paper introduces FANTOM, a Bayesian framework for causal discovery in multivariate time series that exhibit nonstationarity, regime shifts, and complex noise structures such as heteroscedasticity. Unlike existing methods that assume stationarity and homoscedastic Gaussian noise, FANTOM jointly learns the number of latent regimes, their temporal boundaries, and a separate causal graph (DAG) for each regime. It does so using a Bayesian EM algorithm that integrates conditional normalizing flows to model input-dependent noise. The authors establish theoretical identifiability guarantees under heteroscedastic (Gaussian and non Gaussian) noise settings and for mixtures of identifiable causal models, ensuring the recoverability of the true graph structure. Experimental results on synthetic and real-world data show that FANTOM consistently outperforms prior approaches, especially in challenging nonstationary or heteroscedastic conditions.", "questions": "1.\tInitialization sensitivity and “regime-collapse” heuristics. The E-step begins with an Nw >K equal-window initialization and then prunes or merges regimes via heuristic rules (e.g., discarding windows with fewer than ζ points and re-assigning their samples). While the qualitative intuition is clear, the paper does not quantify how sensitive FANTOM’s performance is to the initial choice of Nw, window length, or the pruning threshold. Canyou provide insights on this? \n\n2.\tClarification of ELBO derivation from Equation (2) to Equation (6) The derivation from the general ELBO expression in Equation (2) to the explicit M-step optimization objective in Equation (6) is not fully detailed in the main text. Several substitutions appear to be made (e.g., replacing the latent regime distribution p(zt∣xt,x<t) with soft assignments βtr, separating the graph prior and entropy terms, and incorporating the time-varying regime prior πt(ωr)but these steps are not explicitly shown. Can you provide details of derivation? \n\n3.\tClarification on neural network design for fi,r and the functional constraint. In Section 3, the authors state that the structural function fi,r is designed to satisfy the condition ∂ fi,r / ∂xt−τ j= 0 and they propose a flexible neural architecture in Equation (4) to enforce this. However, the design and enforcement of this partial derivative constraint within the neural network is not clearly explained. Could the authors clarify how the architecture in Equation (4) ensures that with respect to a non-parent input is exactly zero?", "limitations": "No, the authors have not adequately addressed the limitations of their work nor any potential societal impact. They did mention Theoretical guarantees (e.g., identifiability) are currently restricted to stationary settings, while empirical results are shown under nonstationary conditions. While the method is general-purpose, it would be valuable to include a short discussion on potential societal implications—especially when applying causal discovery to sensitive domains (e.g., public policy, healthcare, or financial modeling)—and how to responsibly interpret or validate discovered causal structures.", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\n\n•\tQuality: The paper demonstrates a high level of technical rigor, offering a unified probabilistic model that addresses both structure learning and temporal regime segmentation. The use of conditional normalizing flows to capture heteroscedasticity is well-justified, and the variational EM formulation is grounded in a clearly articulated ELBO. The proofs of identifiability under heteroscedastic noise and their extension to mixtures of regimes further strengthen the theoretical foundations.\n\n•\tOriginality: The work is novel in its comprehensive treatment of causal discovery for multivariate time series under both nonstationary and heteroscedastic conditions. The idea of integrating regime inference, graph learning, and flow-based noise modeling in a single Bayesian framework is a substantial innovation compared to prior works that treat these steps separately or rely on simplistic assumptions.\n\n•\tSignificance: The proposed method addresses a major limitation in the causal discovery literature—namely, its sensitivity to regime shifts and complex noise. This makes the contribution relevant to real-world domains such as finance, climate, and neuroscience, where such assumptions are frequently violated. The empirical results demonstrate clear gains over strong baselines.\n\n•\tClarity: Overall, the paper is decently written and structured, with clear motivation and appropriate mathematical detail. \n \nWeaknesses:\nClarity: the notations are not well explained. For example some aspects of the ELBO derivation from eqn (2) to (6) are not explained,  and the main Figure 2 can be more clearly explained- the terms in eqn 4. \n\nTechnical: 1) The initialization strategy and merging procedure in the E-step rely on heuristic pruning and regime collapse rules, which, while practical, introduce some ad hoc elements. The sensitivity of the results to this initialization (e.g., number of initial regimes $N_w$) is not deeply explored.\n\n2)\tWhile the identifiability results are a strong theoretical asset, they are derived under stationary settings. The implications for nonstationary time series—where regime transitions are latent and evolving—are not theoretically addressed and are instead treated empirically.", "was_revised": true, "final_justification": "This paper presents FANTOM, a novel Bayesian framework for causal discovery in multivariate time series with nonstationarity and heteroscedastic noise. The work is technically rigorous, integrates regime segmentation and graph learning in a unified probabilistic model, and is supported by identifiability guarantees and strong empirical results. While some aspects (e.g., initialization heuristics, clarity of derivations) could be improved, the rebuttal addressed my main concerns satisfactorily. Overall, I find this a solid and impactful contribution and recommend acceptance.", "rebuttal": "We thank the reviewer for their valuable feedback and suggestions. We are pleased they find that our paper is novel, demonstrates a high level of rigor, and our choices are well-justified. The reviewer finds that our method addresses a major limitation in causal discovery and is supported by theoretical results for the single and multi regime settings. Below, we address their specific comments and suggestions. We hope that our responses are convincing.\n\n>Clarification of ELBO derivation from Equation (2) to Equation (6) \n\nWe thank the reviewer for pointing out this question, we will add a section in our appendices in which we will clarify the mathematical derivation Eq.(2) to Eq.(6) and the intuition behind it.\n\nWe perform the maximization of ELBO presented in proposition 3.1, using a BEM procedure, where we alternate between E-step (updating posterior probabilities while fixing all the parameters $\\Theta = \\\\{\\theta^r,\\phi^r, \\omega^r\\\\}\\_{r=1}\\^{N_w}$) and M step (Updating the parameters while using the posteriors learned in the E-step), we can summarize the process as follows:\n\n* In the Estep: we learn $\\beta\\_t^r = p(z_t|x_t,x_{<t}, \\Theta\\_{old})$, where $\\Theta\\_{old}$ is $\\Theta$ of the previous iteration.\n* In the M-step: we fix the learned posterior probabilities to the values $\\beta\\_t^r $ and we update $\\Theta$.\n* By fixing the posterior probabilities in the M-step, the entropy of these probabilities $H\\left(p(z_t \\mid x_t, x_{<t},\\Theta_{\\mathrm{old}})\\right)$ is a constant, which allow us to discard it.\n\nThe detail derivation from Eq.(2) to Eq.(6) is the following:\n\n$\\\\mathrm{ELBO}(\\Theta) = \\\\sum_{t=1}^{|T|}\\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})} \\Big[\\\\mathbb{E}\\_{p(z_t \\mid x_t, x_{<t})}\\big(\\log p_{\\theta_{z_t}}(x_t \\mid x_{<t}, \\mathcal{G}^{z_t})+ \\log p(z_t) \\big)+ H\\big(p(z_t \\mid x_t, x_{<t})\\big)\\Big]+ \\sum\\_{r=1}^{N_w}\n\\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q_{\\phi_r}(\\mathcal{G}^r)\\right)$\n \n$= \\sum_{t=1}^{|T|} \\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})}\\Big[\\sum_{z_t} p(z_t \\mid x_t, x_{<t}, \\Theta_{\\mathrm{old}})\\big(\\log p_{\\theta_{z_t}}(x_t \\mid x_{<t}, \\mathcal{G}^{z_t}) + \\log p(z_t)\\big) \\Big] + \\\\underbrace{H\\left(p(z_t \\mid x_t, x_{<t},\\Theta_{\\mathrm{old}})\\right)}_{\\text{Cte}} + \\\\sum\\_{r=1}^{N_w} \\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q\\_{\\phi_r}(\\mathcal{G}^r)\\right)$\n\nWe replace $p(z_t \\mid x_t, x_{<t}, \\Theta_{\\mathrm{old}})$ by $\\beta_t^{r}$, $p(z_t)$ by our prior choice $\\pi_t(\\omega^r)$, and we discard the entropy of the posterior because it is a constant in these steps. Hence, we got the Eq (6):\n\n$=\\\\mathbb{E}\\_{q\\_{\\phi}(\\mathcal{G})}\\Big[\\\\sum\\_{t=1}^{|T|}\\sum_{r=1}^{N_w}\\beta_t^r \\big(\\log p_{\\theta_r}(x_t \\mid x_{<t},\\mathcal{G}^r) + \\log \\pi_t(\\omega^r)\\big)\\Big]  + \\sum_{r=1}^{N_w} \\\\mathbb{E}\\_{q\\_{\\phi_r}(\\mathcal{G}^r)}\\left[\\log p(\\mathcal{G}^r)\\right]+ H\\left(q_{\\phi_r}(\\mathcal{G}^r)\\right).$ Eq(6)\n\nAfter the maximization of Eq.(6), FANTOM update the posteriors $p(z_t \\mid x_t, x_{<t}, \\theta_{\\mathrm{old}})$ following Eq.(7).\n\n>The initialization strategy and merging procedure in the E-step rely on heuristic pruning and regime collapse rules, which, while practical, introduce some ad hoc elements. The sensitivity of the results to this initialization (e.g., number of initial regimes Nw) is not deeply explored.\nWhile the qualitative intuition is clear, the paper does not quantify how sensitive FANTOM’s performance is to the initial choice of Nw, window length, or the pruning threshold. Can you provide insights on this?\n\nWe thank the reviewer for this suggestion. FANTOM uses two hyperparameters: (i) the window initialization that gives rise to initial regimes, and (ii) ζ which cancels short regimes during the BEM training loop. The choice of window length determines the number of initial regimes, starting with small windows will give rise to a bigger number of regimes. For example, in the presence of two unknown groundtruth regimes [3000,2500] a window length = 1000 will give rise to **five initial regimes, hence N_w=5** [1000,1000,1000,1000,1500], then FANTOM will smoothly converges to the two ground-truth regimes.\n\nTo quantify the sensitivity of FANTOM to the choice of window length and ζ, we conducted new experiments that shows the robustness of FANTOM towards the choices of these hyperparamerters. To show efficiency under long-horizon conditions, we used regimes with lengths [2000, 3000, 2000, 2500], and we also show that FANTOM can scale to more than 3 regimes with different initializations and for different number of nodes, the results are presented in these tables:\n\n**Table 1: Performance of FANTOM with varying hyperparameters (window size and ζ) on 10‑node graphs with 2, 3, and 4 regimes.**\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Numb | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   10  | 4    | 17'8s   |           5           | 91.2 |  13 | 85.7 |  96.5  |\n|    (1500,1200)    | 2 |   10  | 4    | 10'12s  |           1           | 97.9 |  9  | 89.8 |  95.1  |\n|     (1000,900)    | 3 |   10  | 4    | 30'33s  |           5           | 94.2 |  8  | 94.4 |  96.4  |\n|    (1500,1200)    | 3 |   10  | 4    | 16'2s   |           3           | 96.7 |  10 | 91.9 |  91.8  |\n|     (1000,900)    | 4 |   10  | 4    | 40'2s   |           5           | 95.8 |  12 | 92.7 |  92.9  |\n|    (1500,1200)    | 4 |   10  | 4    | 25'38s  |           3           | 97.5 |  16 | 90.5 |  91.5  |\n\n**Table 2: Performance of FANTOM with varying hyperparameters (window size and ζ) on 40‑node graphs with 2, 3, and 4 regimes.**\n|                   |   |       |      |         | Heteroscedastic noise |      |     |      |        |\n|:-----------------:|---|:-----:|------|---------|:---------------------:|:----:|:---:|:----:|:------:|\n| (window, $\\zeta$) | K | Nodes | Numb | Running |         Inst.         |      | Lag |      | Regime |\n|                   |   |       | iter | Time    |          SHD          |  F1  | SHD |  F1  |  Acc.  |\n|     (1000,900)    | 2 |   40  | 4    | 22'15s  |           33          | 81.9 |  27 | 91.9 |   100  |\n|    (1500,1200)    | 2 |   40  | 4    | 10'20s  |           35          | 83.8 |  25 | 91.7 |   100  |\n|     (1000,900)    | 3 |   40  | 4    | 34'25s  |           49          | 85.6 |  58 | 88.6 |  99.8  |\n|    (1500,1200)    | 3 |   40  | 4    | 15'52s  |           43          | 87.3 |  72 | 87.3 |  99.9  |\n|     (1000,900)    | 4 |   40  | 5    | 66'30s  |           34          | 86.5 |  53 | 91.9 |  99.7  |\n|    (1500,1200)    | 4 |   40  | 4    | 21'38s  |           32          | 87.1 |  54 | 91.8 |  99.6  |\n\nFrom the table, FANTOM scales to 40 nodes and converges to four regimes under different initializations. **With initialization = 1500 and ζ =1200**, it starts with **six initial regimes** [1500,1500,1500,1500,1500,2000] and converges smoothly to four in 21 min 38 s. With **initialization = 1000 and ζ =900**, it starts with **nine initial regimes** [1000,1000,1000,1000,1000,1000,1000,1000,1500] and converges to four with 99.7% accuracy in 66 min 30 s. The choices of hyperparameters does not affect FANTOM performances, it only affects the runtime of the method. Also, in our real world experiment for **Epilepsy detection**, FANTOM is initialized with **eight initial regimes** and reliably converges to the **two ground-truth states [non seizure, seizure]**.\n\n>Clarification on neural network design for fi,r and the functional constraint.\n\nWe thank the reviewer for this question. We propose flexible functional designs for $f^{i, r}$, which must respect the relations encapsulated in $\\mathcal{G}^r$. Namely, if $x\\_{t-\\tau}^j \\notin \\\\mathbf{P a}\\_{\\\\mathcal{G}^r}^i(<t) \\cup \\\\mathbf{Pa}\\_{\\mathcal{G}^r}^i(t)$, then $\\partial f^{i, r} / \\partial x_{t-\\tau}^j=0$. \n\nWe instantiate $f\\_{\\theta^r}^{i, r}$ as:\n\n$f\\_{\\theta^r}^{i, r}\\left(\\mathbf{P a}\\_{\\mathcal{G}^r}^i(<t), \\mathbf{P a}\\_{\\mathcal{G}^r}^i(t)\\right)=\\psi^r\\left(\\sum\\_{\\tau=0}^L \\sum\\_{j=1}^d G\\_{\\tau, j i}^r \\cdot \\vartheta^r\\left(x\\_{t-\\tau}^j, e\\_{\\tau, j}^r\\right), e\\_{0, i}^r\\right)$,\n\nwhere $\\theta^r$ includes the learnable parameters of the neural network $\\vartheta^r$ and $\\psi^r$ and the learnable embedding $e\\_{\\tau, j}^r$ and $e\\_{0, j}^r$. \n\n$G\\_{\\tau, j i}^r$ encodes whether the edge $j \\rightarrow i$ is present at lag $\\tau$ in regime $r$. When $x\\_{t-\\tau}^j \\notin \\\\mathbf{P a}\\_{\\\\mathcal{G}^r}^i(<t) \\cup \\\\mathbf{Pa}\\_{\\mathcal{G}^r}^i(t)$, we have $G\\_{\\tau, j i}^r =0$; consequently, variations in $x\\_{t-\\tau}^j$ do not impact the value of $f\\_{\\theta^r}^{i, r}$, hence $\\partial f^{i, r} / \\partial x_{t-\\tau}^j=0$.\n\n>No, the authors have not adequately addressed the limitations of their work nor any potential societal impact.\n\nWe thank the reviewer for bringing to our attention that we did not mention that a limitation section is introduced in our appendices. In our revised version, we will explicitly reference the Limitations section in the main text, and we will add a paragraph on societal impacts."}]

Author rebuttal:

Meta-review: The authors propose a causal discovery framework in the realistic setting of non-stationary, multivariate time series data. Overall, the contribution is appreciated by the reviewers and they all recommend acceptance of the paper.  The novelty, experiments, and the presentation were all found to be adequate.
Decision: Accept (poster)