Rhino: Deep Causal Temporal Relationship Learning with History-dependent Noise
Authors: Wenbo Gong
Joel Jennings
Cheng Zhang
Nick Pawlowski
Venue: ICLR 2023 notable top 25% (2023)
Keywords: Structure learning
Causal discovery
Time series
Structure equation model
deep generative model

Abstract: Discovering causal relationships between different variables from time series data has been a long-standing challenge for many domains. For example, in stock markets, the announcement of acquisitions from leading companies may have immediate effects on stock prices and increase the uncertainty of the future market due to this past action. To discover causal relations in such case, the model needs to consider non-linear relations between variables, instantaneous effect and the change of noise distribution due to past actions. We name the latter as history-dependent noise. However, previous works do not offer a solution addressing all these problems together. In this paper, we propose a structural equation model, called Rhino, which combines vector auto-regression, deep learning and variational inference to model non-linear relationships with instantaneous effects while allowing the noise distribution to be modulated by history observations. Theoretically, we prove the structural identifiability of Rhino. Our empirical results from extensive synthetic experiments and two real-world benchmarks demonstrate better discovery performance compared to relevant baselines, with ablation studies revealing its robustness under model misspecification.

RHINO: DEEP CAUSAL TEMPORAL RELATIONSHIP LEARNING WITH HISTORY-DEPENDENT NOISE
Wenbo Gong, Joel Jennings, Cheng Zhang & Nick Pawlowski
Microsoft Research Cambridge, UK
(wenbogong, jennings, .hang, nick. Hick.pawlowske @mi emicroseft.cm com
ABSTRACT
Discovering causal relationships between different variables from time series data has been long-standing challenge for many domains such climate science, finance, and healthcare. Given the complexity of real-world relationships and the nature of observations in discrete time, causal discovery methods need con- sider non-linear relations between variables, instantaneous effects and history- dependent noise (the change of noise distribution due to past actions). However, previous works do not offer a solution addressing all these problems together. In this paper, propose novel causal relationship learning framework for time- series data, called Rhino, which combines vector auto-regression. deep learning and variational inference to model non-linear relationships with instantancous ef- fects while allowing the noise distribution to be modulated by historical observa- tions. Theoretically, we prove the structural identifiability of Rhino. Our empir- ical results from extensive synthetic experiments and two real-world benchmarks demonstrate better discovery performance compared to relevant baselines, with ablation studies revealing its robustness under model misspecification.
INTRODUCTION
Time series data is collection of data points recorded at different timestamps describing a pattern of chronological change. Identifying the causal relations hetween different variables and their
teractions through time (Spirtes et al.] 2000; Berzuini et al.] 2012 Guo et al. 2020: Peters et al.) 2017) is essential for many applications e.g. climate science, health care, etc. Randomized control trials are the gold standard for discovering such relationships, but may be unavailable due to cost and ethical constraints. Therefore, causal discovery with just observational data is important and fundamental to many real-world applications (Löwe al. 2022 Bussmann et al., [021]] Moraffah
et al. 2021 2020
Wu et al. 2020 Runge 2018: Tank et al. 2018 Hyvärinen et al.] 2010 Pamfil et al.
The task of temporal causal discovery can be challenging for several reasons: (1) relations between variables can be non-linear in the real world; (2) with a slow sampling interval, everything happens in between will be aggregated into the same timestamp, i.e. instantancous effect; (3) the noise may be non-stationary (its distribution depends on the past observations), i.e. history-dependent noise. For example, in stock markets, the announcements of some decisions from a leading company after the market closes may have complex effects (i.e. non-linearity) on its stock price immediately after the market opening (i.e. slow sampling interval and instantaneous effect) and its price volatility may also be changed (i.e. history-dependent noise). Similarly, education, students that recently earned good marks algebra tests should also score well on upcoming algebra exam with little variation (i.e. history-dependent noise).
To the best of our knowledge, existing frameworks' performances suffer in many real-world sce- narios as they cannot address these aspects in a satisfactory way. Especially, history-dependenn noise has been rarely considered in past. A large category of the preceding works, called Granger causality (Granger, 1969). is based on the fact that cause-effect relationships can never go against
time. Despite many recent advances (Wu et al. 2020 Shojaie & Michailidis 2010 Siggiridou &

---

Kugiumtzis 2015; Amornbunchornvey al. 2019 Lowe al. 2022 Tank et al. 2018 Bussmann et al. 2021 Dang et al. 2018: Xu et al. 2019). they all rely on the absence of instantaneous effects with fixed noise distribution. Constraint-bas based methods have also been extended for time series causal discovery (Runge, 2018: 2020), which is commonly applied by folding the time-series. This introduced new assumptions and translated the aforementioned requirements challenges in condi- tional independence testing (Shah & Peters 2020).Addiite Additiinnlly, they require stronger faithfulness assumption and can only identify the causal graph up to Markov equivalence class without detailed functional relationships.
An alternative line of research leverages the development of causal discovery with functional causal models (Hyvärinen et al.] 2010 Pamfil et al., 2020) Peters et al. (013). They can model both instan- taneous and lagged effects as long as they have theoretically guaranteed structural identifaaiilit. Unfortunately, they do not consider history-depende
noise. One central challenge of modelling this dependency is that noise depending on the lagged parents may break the model structural iden- tifiability. For static data, Khemakhem et al. (2021) proves the structural identifiability only when this dependency is restricted to simple functional form. Thus, the key research question is whether the identifiability be preserved with complex historical dependencies in the temporal setting.
Motivated by these requirements, we propose a novel temporal discovery framework called Rhino (deep causal temporal relationship learning with history dependent noise), which can model non- linear lagged and instantaneous effects with flexible history-dependen noise. Our contributions are:
A novel causal discovery framework called Rhino, Revision((2)- hiTaccnsiitt- ing of novel functional form of its SEMs and variational training framework, where the proposed form of its SEM combines vector auto-regression and deep learning to model non-linear lagged and instantaneous effects with history
noise.
We prove that Rhino SEMs with the proposed form are structurally identifiable. To achieve this, we provide general conditions for structural identifiability with history-dependent noise, of which the form Rhino SEMs is a special case. Furthermore, we clarify re- lations to several previous works.
We conduct extensive synthetic experiments with ablation studies to demonstrate the ad- vantages of Rhino and its robustness across different settings. Additionally, we compare its performance to a wide range of baselines in two real-world discovery benchmarks.
2 BACKGROUND
In this section, we briefly introduce necessary prerequisite knowledge. In particular, we focus on structural equation models, Granger causality (Granger 1969) and vector auto-regression. For view of more recent related work, please refer to Section 5
Structural Equation Models (SEMs) Consider X IRD with D) variables, SEM describes the causal relationships between them given a causal graph G:
X f.
where Pat: are the parents of node i and e' are mutually independent noise variables. Under the context multivariate time series, X, = (xx) kev where V is a set of nodes with size D, the corresponding SEM given a temporal causal graph G is
X ft, (Pab(< t), (t),
(2)
where Pa (< t) contains the parent values specified by G in previous time (lagged parents); (ta are the parents at the current time (instantuneous parents). The above SEM induces joint distri- bution over the stationary time series ) (see Assumppttiiio 1 in Appendix B for the definition). However, functional causal models with the above general form cannot be directly used for causal discovery due to the structural unidentifiability (Lemma 1, ZZZang et al. (2015) One way to solve this sacrificing the ficxibility by restricting the functional class. For example, additive noise models (ANM), (Hoyer et al. 2008)
X f. f.Paa((X)) + +
which have recently been used for causal reasoning with non-temporal data (Geffner et al. 2022
2

---

Granger Causality Granger causality (Granger 1969) has been extensively used for temporal causal discovery. is based on the idea that the series X does not Granger cause X' the history, X 20 does not help the prediction of X; for some given the past of all other time series X for f
Definition 2.1 (Granger Causality (Tank et al. 2018 Löwe et al. 2022 Given a multivariate stationary time series k and SEM Ji, defined as
x = ft.. (P
X Granger causes X if E [1 such that x E ( t) and fi, depends on X
Granger causality is equivalent to causal relations for directed acyclic graph (DAG) if there are no latent confounders and instantancous effects (Peters et al. 2013 2017). Apart from the lack of instantaneous effects, it also ignore the history-dependent noise with independent eç.
Vector Auto-regressive Model Another line research focuses on directly fitting the identifiable SEM to the observational data with instantancous effects. One commonly-usee approach is called vector regressgrnsi (Hyvärinen et al. 2010
Pamfil et al. 2020):
K
D
x = 8 + +E>
B
BrsoX/ +
(5)
T=0 T=Uj=1 j=1
where Be is the offset, is the model lag, B, DxD is the weighted adjacency matrix specifying the connections at time T (i.e. if B+j means no connection from X to X;) and ej is the independent noise. Under these assumptions, the above linear SEM is structurally identifiable, which is a necessary condition for recovering the ground truth graph (Hyvärinen et al.) 20100) Peters er al. 2013] Pamfil et al., 2020). However, the above linear SEM with independent noise variables too restrictive fulfil the requirements described in [Section T. Therefore, the research question how to design structurally identifiable non-linear SEM with fiexible history-dependent noise.
3 RHINO: RELATIONSHIP LEARNING WITH HISTORY DEPENDENT NOISE
This section introduces Rhino: Section 3.1 describes the novel functional form for Rhino SEMs, allowing for history-dependent noise. Section details the associated variational inference frame- work for causal discovery.
3.1 MODEL FORMULATION
For multivariate stationary time series (X,) we assume that their causal relations follow temporal adjacency matrix Go:K with maximum lag K where Grę(1.8] specifies the lagged effects
;. Go specifies the instantaneous parents. We define Gr,ij if X X and otherwise. | We propose a novel functional form for Rhino's SEM that incorporates non-linear relations, instantancous effects, and flexible history-depente
between X, and X,,
noise:
X fi (Pa' t), (t) + gi (Pac t),
(6)
where fi general differentiable non-linear function, and 9; is differentiable transform s.t. the transformed noise has proper density. Despite of an additive structure. formulation of- fers much more flexibility in both functional relations and noise distributions compared to previous works (Pamfil et al. 2020; Peters et al. 2013). By placing few restrictions on fo Is. it can capture functional non- lineaiiyyy through fi and transform ej through a flexible function gi. depending on Pa t), to capture the history dependency of the additive noise.
Next, we propose flexible functional designs for fa, 9t. which must respect the relations encapsulated in G. Namely, if X & Pa t) ()) then ff//ax/ 0 and similarly for g:. We design
A(Pa((< Pala =
ΣΣ 7=0
Griiees (xL_.)
In the following, interchange the usage of the notation and Go:K for brevity.
3

---

where {a and e. (i D] and [0, K]) are neural networks. For efficient computation, weight sharing across nodes and lags: (x(-) = ((,, uo and (T(() = E(', Ur,j), where Ur is the trainable embedding for node i at time
The design of 9j needs to properly balance the flexibility and tractability of the transformed noise density. We choose conditional normalizing flow, called conditional spline flow (Trippe & Turner 2018 Durkan et al. 2019 Pawlowski et al. 2020), with a fixed Gaussian noise cj for all t, i. The spline parameters are predicted using a hyper-network with similar form [Eq. (7) incorporate history dependency. The only difference is now r is summed over [1, K] to remove the instantancous parents, Due to the invertibility of gs. the noise likelihood conditioned on lagged parents is
0g, Oej
Ps. (gi(E)) Pac(< ()) (e(
(8)
3.2 VARIATIONAL INFERENCE CAUSAL DISCOVERY
Rhino adopts a Bayesian view of causal discovery (Heckerman et al., 2006). which aims to learn graph posterior distribution instead of inferring a single graph. For N observed multivariate time series 0.7
(N) the joint likelihood with model parameter 0 is 0.T X
P X O.
G) P(G) O.T
II 1==1
(m) Pa 0:7
G).
Graph Prior When designing the graph prior, we combine three components: (1) DAG
straint; (2) graph sparseness prior; (3) Nomain-specifit prior knowledge (optional). Inspired by the
NOTEARS (Zheng al. 2018) Geffner et al., 2002 Morales-Alvarez et al.) (02T), we propose the
following unnormalised prior
PGG) exp (->x
ph? (Go) ah((G)) >p Go
xxx11)
(10)
where h h(G)- GoC)- D is the DAG penalty proposed in (Zheng et al. 2018 and if and only if G is a DAG; ④ is the Hadamard product; GO is an optional domain-specific prior graph, which can be used when partial domain knowledge available; A, Ap specify the strength of the graph sparseness and domain-specific prior terms respectively; a, P characterize the strength of the DAG penalty, Since the lagged connections specified G 1: can only follow the direction of time, only the instantancous part, Go, can contain cycles. Thus, the DAG penalty is only applied to Go-
Variational Objective Unfortunately, the exact graph posterior p O........... Xo-T) is tractable due to the large combinatorial space of DAGs. To overcome this challenge. we adopt variational inference (Blei et al., 2017]] Zhang et al. 2018), which uses a variational distribution (G) to upproximate the true posterior. We choose Io distributions (refer to (Appendix El for details). The corresponding evidence lower bound (ELBO) is
( (G) to be product of independent Bernoulli
log po (
(1
X 0:T
2 (G)
Σ log Po
G + logp(G) + H(gs (G))
(11)
0.1
ELBO(8.6)
where Hgg (G)) is the entropy (G). From the causal Markov assumption and auto-regressive nature, we can further simplify
T D
log Po x
G)
log Po (X
(n) Pa t), Pa)(())
(12)
1=0
1=1
and from Rhino's functional form
Eq. (6)
proposed
Section 3.1
log leg (X (n)
Pa (<<)) t), (t)) log Ps.
(n))) Pa (< t)
(13)
where 2m((m) = (())) fo (Pa (< Pa and Py. is defined in Eq. (8) Appendix A for details).
The parameters 0, learned by maximizing the ELBO, where the Gumbel-soltmax gradient
4

---

estimator is used (Jang al. 2016 Maddison et al. 2016). We also leverage augmented Lagrangian training (Hestenes) 1969; Andreani et al. 2008), similar as Geffner et al. (2022), to anneal a, p to make sure Rhino only produces DAGs (refer to Appendix B.1 in Geffner et al. (2022).. Once Rhino
is trained, the temporal causal graph can be inferred by G ~ (G(G).
Treatment effect estimation As Rhino learns the causal graph and the functional relationship simultaneously, it can be extended for causal inference tasks such as treatment effect estimation (Geffner al. 2022). See Appendix D for details.
4 THEORETICAL CONSIDERATIONS
Here, we show the theoretical guarantees of Rhino including (1) the structural identifiability of Rhino SEMs and (2) soundness of the proposed variational inference framework. Together, they guarantee the validity of Rhino as causal discovery method. In the end, we clarify relations to existing works.
4.1 STRUCTURAL IDENTIFIABILITY
One of the key challenges for causal discovery with a flexible functional relationship to show structural identifiability. Namely, cannot find two different graphs that induce the same joint likelihood from the proposed functional causal model. In the following, we present a theorem for Rhino SEMs that summarizes main theoretical contribution.
Theorem 1 (Identifiability of Rhino SEMs). Assuming Rhino SEMs satisfy the causal Markov property, minimality, sufficiency, DAGness and the induced joint likelihood has proper density (see Appendix for details), and (1) all functions and induced distributions are third-order differentiable; (2) function ft is non-linear and not invertible w.r.t. any nodes Pa (t); (3) the double derivative (logPa Po. (g (<i) Palg << t)))" w.r.t Ef is zero at most at some discrete points, then the SEM with the form defined in Eq. (6) is structural identifiable for both bivariate and multivariate case.
Sketch of proof. This theorem is a summary of a collection of theorems proved in Appendix The strategy is instead of directly proving the identifiability of the model, we provide identifiability conditions for a general temporal SEMs, followed by showing a generalization of Rhino SEMs satisfies these conditions. The identifiability of Rhino SEMs directly follows from it.
The first step is to prove the bivariate identifiability conditions that general temporal SEMs (Eq. (2) should satisfy (refer to Theorem in Appendix B.I). Inspired by the techniques from Peters et al. (2013). we proved that temporal SEMs are bivariate identifiable if (1) the model for iiiiill conditions is identifiable; (2) the model is identifiable w.r.t. instantaneous parents. Compared to Petees et al. (2013), we relaxed the identifiable model class condition that it can applied to history-dependent noise. In particular, (2) implies we only need to pay attention to instantaneous parents, rather than the entire parents (Peters al. 2013), and opens the door for flexible lagged dependency.
Prove bivariate identifiability conditions for general temporal SEMs
Identifiability of history-dependen post non-linear model Next, we propose novel gener- alization of Rhino SEMs, called history-dependent PNL. [Teorem 4 and Corollary 4.1 in Ap- pendix
B.2 prove it is bivariate identifiable w.r.t. instantaneous parents (i.e. satisty the conditions
in Theorem 3) with additional assumptions (1), (2) and (3) in Theorem The functional form of destdryd PNL is defined as
X; Vie (fi (Pab(<< (< t), (t)( + gst (Pai eu) Pa t))
where is invertible w.n.t. the first argument. Due to its similarity to PNL (Zhang & Hyvarinen, 2012) and ANM (Hoyer et al. 2008), we combine their techniques to prove our results. The bivariate identifiability of Rhino SEMs directly follows from this with being the identity mapping.
Generalization to multivariate case In the end, we prove the above bivariate identifiability can be generalized to the multivariate case by adapting the techniques from Peters al.) (2012) and combining it with the proof strategy from step 1. Refer to Theorem Ain/Appediix B.3 for details.
5

---

4.2 VALIDITY OF VARIATIONAL OBJECTIVE AND RELATIONS OTHER METHODS
Next, we show the validity of the variational objective (Eg. (1I) in the sense that optimizing it can lead the ground truth graph. Theorem I in Geffner al. (2022) justifies the validity of the variational objective under the same set of assumptions as Rhino.
Theorem 2 (Validity of variational objective (Geffner et al., 2022 Assuming the conditions in Theorem I are satisfied, and further assume no model misspecification (see DDefinition B.I). then the solution (0', L/s (G)) from optimizing Eq. (11) with infinite data satisfies 4o (G) 0 G.), where G' is unique graph. In particular, G' = G' and per (Xo:T; G¹) p o(o:;; ; G*), where G is the ground truth graph and p( (oo;:; G*) is the true data generating distribution.
G
We emphasize that Theorem 2 guarantees the correctness with the global optimum, rather than characterizing the convergence during optimizing Eq. (IT) We assume the Rhino only searches in the DAG space. In practice, with augmented Lagrangian (Wei al., 2020 Ng al. 2022). we can only obtain an approximate posterior over DAGS, converging to a local optimum.
Relation to other methods Many previous works of using functional causal model for causal time series discovery (Hyvärinen et al. 2010 Pamfil et al.] 2020 Tank et al., 2018 Peters et al., (013)) closely related to Rhino. Since Rhino SEMs incorporate history-dependent noise with flexible function non-linearity, it is the most flexible member of this family. Refer to AAppendix C Cfor details.
RELATED WORK
Assaad al. 2022) provides comprehensive overview of causal discovery method for time series.
The first is Granger causality, which can be further split into (1) vector auto-regression (Wu et al.
2020 Shojaic & Michailidis 2010 Siggiridou Kugiumtzis 2013 2015| Amonbunchornvv al. 2019) and (2) deep learning (Löwe et al. [0222] Tank et al.) 2018 Bussmann et al. [021] Dang et al,
2018 Xu et al., 2019). Granger causality methods cannot handle instantaneous effects, which can
be observed in slow-sampling system. Additionally, they also assume fixed noise distribution.
Using funcitonal causal models can mitigate the aforementioned two problems. VARLiNGaM (Hyvärinen et al. 2010) extends the identifiability theory of linear non-Gaussian ANM (Shimizu ct al. 2006) to vector auto-regression for modelling time series data. DYNOTEARS (Pamfil et al., 2020 leverages the recently proposed NOTEARS framework (Zheng et al. 2018) to continuously relax the DAG constraints for fully differentiable DAG structure learning. However, the above ap- proach is still limited to linear functional forms. TiMINo (Peters al. 2013) provides general theoretical framework for temporal causal discovery with functional causal models. Our theory leverages some of their proof techniques. Unfortunately, all the aforementioned methods assume no history dependency for the noise. On the other hand, Rhino can model (1) non-linear function relations; (2) instantancous effect; (3) and distory-dependee
noise at the same time.
The third category is constraint history-dependent noise. PCMCI (Runge et al., 2019) combines PC (Spirtes et al. (2000 and con- ditional independence test for discovery from time series. PCMCI (Runge) 2018: 2020) further extends PCMCI infer both lagged and instantaneous effects. CD-NOD (Huang al. 2020) has recently been proposed to handle non-stationary heterogeneous data, where the data distribution can shift across time. Despite their generality, they can only infer MECs; cannot learn the explicit func- tional forms between variables; and require a stronger assumption than minimality (i.e. faithfulness).
approaches. Due to its non-parametric nature, it can handle
6 EXPERIMENTS
We release the code of Rhino for reproducing the following experiments
6.1 SYNTHETIC DATA
We evaluate our method on a large set of synthetically generated datasets with known causal graphs. We use the main body of this paper to present the overall performance our method compared to rel-

---

evant baselines and one ablation study on the robustness lag mismatch. In Appenddi F.3 duct extensive analysis, including (1) on different graph type; (2) ablation on history-dependenny; (3) ablation study instantaneous effect. This set of datasets are generated by various settings (e.g. type of graphs, instantaneous/no instantaneous effect, etc.). datasets are generated for each combination of settings with different seeds, yielding 160 datasets in total. In order to compre- hensively test Rhino's ' robustness, we deliberately generated 75% of the datasets that mismatch the Rhino configurations. Details of the data generation can be found in Appendix F.I
We compare Rhino to wide range of baselines, including VARLiNGaM (Hyvärinen et al. 2010), PCMCI* (Runge 2020) and DYNOTEARS (Pamfil et al., 2020). PCMCI only outputs Markov equivalence classes (MECs). We resolve this by enumerating all DAGs in the MEC. For details on the methods, see Appendix F.2 Additionally, we include two variants of Rhino: (1) Rhino+g. where an independent Gaussian noise is used; (2) Rhino+s, where Gaussian 6j is transformed by an independent spline.
Figure 1 presents the F score of the lagged, instantaneous and temporal adjacency matrix of all methods aggregated over all dataset denoted as 'Lag', 'Inst.' and 'Temporal", respectively. Rhino achieves overall competitive or the best performance in terms of the full temporal adjacency matrix across all possible datasets, especially for lower dimensions. Comparing Rhino's 's lagged discovery to its two variants, the better score indicates the -eseede-dependent noise is useful to the lagged graph discovery, contributing to the better overall F1 performance (Appendix F.3 for ablation).
Despite of the strong performance from PCMCI+ it can only identify the graph up MECs without explicit functional relations. Computationally, PCMCI* exceeds the maximum training time week on 40 nodes (see Appendix F.3). suggesting its computation bottleneck high dimensions.
DYNOTEARS achieves inferior results in general due to limited modelling power from the linear nature. This is much clearer in high dimensions due the increasing difficulty of the problem.
Dim 5
Dim 10
1.0
d 0.5
Model Name
Rhino (L=2)
0.0
Rhino+ (L=2)
Dim 20
Dim = 40
Rhino+s (L=2)
1.0
DYNOTEARS
VAR-LiNGaM
0.5
PCMCI*
0.0
Inst.
Lag
Temporal
Inst.
Lag
Temporal
Figure I: Fi-scores of Rhino (light yellow) compared to all baseline methods. The different subplots show the performance for dataset with different number of nodes. 'L=2' refers to models with lag 2.
We explore the behaviour of Rhino with different lag parameters other than the ground truth lag 2. From Table worse training ligeliielssoods suggest that Rhino with insufficient history (lag 1) unable to correctly model the data and this leads to decrease in Fi scores. Interestingly, Rhino also robust with longer lags. Despite of the slightly better likelihood (lag = 3), achieves comparable performance the model with the correct lag. Also, from their similar Fj Lag score, suggests the extra adjacency matrix is nearly empty.
6.2 DREAM3 GENE NETWORK
In this section. we evaluate Rhino performance with a real-world biology benchmark called DREAM3 (Prill et al., 2010 Marbach al. al. 2009). These datasets often used to evaluate Granger
causality (Khanna & Tan 2019 Tank et al. 2018 Nauta et al.) 2019 Bussmann et al.] 202T but
We note that each method 40 different dataset settings for all possible numbers of nodes.

---

Dim
Rhino (I==1)
Rhino (I.=2)
Rhino (1.=3)
5
F1 Inst.
0.11 = 0.17
0.20 0.22
0.21 0.23
Fj Lag
0.28 + 0.13
0.59 + 0.22
0.57 + 0.24
F1 Temporal
0.34 + 0.12
0.59 + 0.20
0.56 0.22
LL
4.14 1.63
3.83 + 1.62
-3.75 1.64
10
Inst.
0.13 0.17
0.19 + 0.23
0.19
Fí Lag
0.26 0.08
0.51 + 0.17
Fi Temporal
0.28 0.11
0.49 0.18
0.45 0.20
LL
77.77 + 2.09
-7.21 + 2.22
-7.01 1.91
20
F1 Inst.
0.15 + 0.16
0.16 0.17
0.18 0.19
F1 Lag
0.24 + 0.12
0.42 + 0.22
0.40+
F1 Temporal
0.25 0.13
0.39 0.22
LL
15.62 3.16
14.70 + 2.87
14.72 4
40
F1 Inst.
0.13 + 0.16
0.22 + 0.23
0.18 + 0.21
F1 Lag
0.20 = 0.18
0.40 0.31
0.34 + 0.30
F1 Temporal
0.20 0.18
0.37 + 0.30
0.32 + 0.29
LL
31.44 5.16
-30.10 4.71 300.20 = 4.74
Table 1: Comparison of the causal discovery performance of Rhino with different hag-parameters (L E (]. Apart from the 3 F1 scores, LL shows the log-likelihood of the training data.
recently adopted for SEM-based method (Pamfil et al. 2020). The dataset consists in silico mea- surements of gene expression levels for 5 different networks. Each network contains d 100 genes. Each time series represents a perturbation trajectory with time length T = 21. For each network, recorded. The goal is to infer the causal structure of each network. We use the area under the ROC curve (AUROC) as the performance metric. We consider the same baselines in the synthetic experiments (i.e. DYNOTEARS and PCMCI+) without VARLiNGaM since its default implementation fails when the number of variables (d = 100) is greater than the series length (I' 21). Additionally, we also consider relevant Granger causality methods, includ- ing cMLP, cLSTM TTank et al., 2018); TCDF
46 perturbation trajectories
et al.) 2019): SRU and eSRU (Khanna & Tan 2019). Their corresponding results are directly cited from Khanna & Tan (2019). Apppendix G. specifies Rhino hyperparameters. Since the ground truth graph is summary graph (see Defini- tion G.T in Appendix G.2). [Appendix G.2 details about the post-processing step on aggregating
Temporal graph to summary graph for Rhino, DYNOTEARS and PCMCI+.
Method
E.Coli
E.Coli
Yeast
Yeast 2
Yeast
cMLP
0,644
0.568
0.585
0.506
0.528
cLSTM
0.629
0.609
0.579
0.519
0.555
TCDF
0,614
0.647
0.581
0.556
0.557
SRU
0.657
0.666
0.617
0.575
0.55
eSRU
0.66
0.629
0.627
0.557
0.55
DYNO.
0.590
0.547
0.527
0.526
0.510
PCMCI+
0.530 + 0.002
0.519 + 0.002
0.530 + 0.003
0.510 + 0.001
0.512 +:
Rhinotg
0.6731:0.013
0.665+0.009
0.659-10.005
0.59810.004
0.588+10.005
Rhino
0.685+000033
0.680=0.007
0.66410.006
0.585+0.004
0.567+0.003
Table 2: The AUROC of the aggregated adjacency matrix for 5 DREAM3 datasets without self- connections. DYNO. means DYNOTEARS, For Rhino and PCMCI the results are reported by averaging over 5 runs. Khanna & Tan 2019) only reported the single-run results for baselines.
Table 2 demonstrates the AUROC of the summary graph inferred after training. is clear that Rhino and its variant outperform all other methods. Although Rhino is not formulated solve the summary graph discovery, it shows a clear advantage compared to the state-of-the ofth-of-the-ar Granger causality. Thus, Rhino be used to infer either temporal summary graph depending users' needs.
8

---

By inspecting the hyperparameters of Rhino in Appendix G.1 no obvious help for in these datasets. It suggests the recording intervals are fast enough to avoid any aggregation effect. This explains why the Granger causality can also perform reasonably well.
instantaneous effects to provide
Unlike the strong performances of DYNOTEARS and PCMCI+ in synthetic experiments, they per- form poorly in DREAM3. The linear nature of DYNOTEARS seems harm its performance dras- tically. PCMCI+ suffers from the low independence test power under small training data.
Another interesting ablation is to compare with Rhinotg, which performs par with Rhino and achieves better scores on 2 out of 5 datasets. Although we have no access to the true noise mech- anism, we suspect that the added noise is not history-dependent and highly likely be Gaussian. Despite the model mismatch, Rhino is still one of the best methods for this problem. This further strengthens our belief in the robustness of our model under different setups.
6.3 NETSIM BRAIN CONNECTIVITY
In this section, we evaluate Rhino using fMRI imaging data, which has also been used as a bench- mark for temporal causal discovery (Lowe al. 2022 Khanna & Tan 2019 Assaad et al. 2022)
Each time series represents the signal simulated for a human subject, which describes d = 15 dif- ferent regions in the brain with T 200 timestamps. The goal is to infer the connectivity between different brain regions. We assume that different human subjects share the same connectivity. We only use the data from human subject 6 in Sim-3.mar from https: / /www. fmr ib. ox aa.
hik/dataeets////////// /net im/ index. html with self-connections. We use the same set of baselines as DREAM3 (Section 6.2) plus VARLiNGaM, Appendix G.4 describes hyperparameter settings.
Table 3 shows the AUROCs for different meth- ods. Remarkably, the proposed Rhino and its variants achieve significantly better AUROC compared to the baselines. Especially, Rhino obtains nearly optimal AUROC, demonstrating its robustness to the small dataset and good balances between true and false positive rates (see [Appendix H). By comparing Rhino and Rhino+Nolnst,, we conclude that modelling in- stantaneous effects is important in real applica- tion, indicating the sampling interval is not fre- quent enough to explain everything as lagged effects. This can double confirmed by com- paring Rhino+NoInst with Granger causality, where it performs on par with the state-of- the-art baseline when disabling the instanta- neous effect. Last but not least, by compar- ing Rhino+g with Rhino, find that history- dependent noise is also helpful in this dataset.
Method
AUROC
cMLP
0.93
cLSTM
0.83
TCDF
0,91
SRU
0,80
eSRU
0.88
DYNO.
0.90
PCMCI*
0.83 = 0
VARLiNGaM
0.84 + 0
Rhinotg
0.002
Rhino+Nolnst..
0.93 + 0.006
Rhino
0.99 + 0.001
Table 3: AUROCs of the summary graph. Rhino+NoInst is Rhino without instantaneous ef- fects. For Rhino, VARLiNGaM, PCMCI* results are obtained by averaging over different runs.
7 CONCLUSION
Inferring temporal causal graphs from observational time series is an important task in many sci- entific fields. Especially, some applications (e.g. education, climate science, etc.) require the mod- elling of non-linear relationships: instantancous effects and history-dependent noise distributions at the same time. Previous works fail to offer an appropriate solution for all three requirements. Mo- tivated by this, we propose Rhino, which combines vector auto-regression with deep learning and variational inference to perform causal temporal relationship learning with all three requirements. Theoretically, we prove the structural identifiability of Rhino with flexible history-dependent noise, and clarify its relations to existing works. Empirical evaluations demonstrate its superior perfor- mance and robustness when Rhino is misspecified, and the advantages of -deeedeet noise mechanisms. This opens an exciting route of extending Rhino to handle non-so non-stationary time-series and unobserved confounders in future work.
9

---

8 REPRODUCIBILITY STATEMENT
Theoretical Contributions The main theoretical contribution summarized in Theorem 1 This theorem is the result from a collection of theorems proved in Appendix In Appendix B, de- tailed the fundamental assumptions (Assumption 1 Assumption 5) required for the all theorems. The -heorem-specific assumptions are mentioned in the statement of the theorem. To ease the under- standing of the proof, we also provide the skecth of proof in Theorem I] Since JinceTheorem 2 directly cited from (Geffner al. 2022) without major modification, the proof can be found Appendix A
in Geffner et al. 2022).
Empirical Evaluations For synthetic, DREAM3 and Netsim experiments, we listed the hyperpa- rameters in Appendix F.2 Appendix G.1 and Appendix G.4 respectively, [Appendix F. explaains the synthetic data generation. For DREAM and Netsim, the dataset can be found in the public github repo / iigggiithuhu thub. com/ s sakhanna/ SRU_ for_roroooooo for__oor////ofaarrrrrrrrrrrrfrrr /CT/GCT//rr ree mastterdaattt The post processing steps for DREAM3 and Netsim evaluations are described in Appendix G.2
REFERENCES
Chainarong Amornbunchornvej Elena Zheleva, and Tanya Y Berger- Wolf. Variable-lag granger causality for time series analysis. In 2019 IEEE International Conference on Data Science and Advanced Analytics (DSAA), 21-30. IEEE, 2019.
Roberto Andreani, Ernesto G Birgin, José Mario Martínez, and María Laura Schuverdt. On aug- mented lagrangian methods with general lower-level constraints. SIAM Journal on Optimization, 18(4):1286-1300) 2008.
Charles K Assaad, Emilie Devijver, and Eric Gaussier. Survey and evaluation of causal discovery methods for time series. Journal of Artificial Intelligence Research, 73:767-819, 2022.
Carlo Berzuini, Philip Dawid, and Luisa Bernardinell. Causality: Statistical perspectives and ap- plications. John Wiley & Sons, 2012.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe. Variational inference: review for statisti- cians. Journal of the American statistical Association,
877, 2017.
Bart Bussmann, Jannes Nys, and Steven Latré. Neural additive vector autoregression models for causal discovery in time series. International Conference on Discovery Science, pp. 446-460. Springer, 2021.
Xuan-Hong Dang, Syed Yousaf Shah, and Petros Zerfos. seq2graph: Discovering dynamic depen- dencies from multivariate time series with multi-level attention. arXiv preprint 2018.
Hadi Mohaghegh Dolatabadi, Sarah Erfani, and Christopher Leckie. Invertible generative modeling using linear rational splines. In International Conference Artificial Intelligence and Statistics, pp. 4236 4246. PMLR, 2020.
Conor Durkan, Artur Bekasov, lain Murray, and George Papamakarios. Neural spline flows. Ad- vances in neural information processing systems, 32, 2019.
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.02111 02515, 2022.
Clive WJ Granger. Investigating causal relations by econometric models and cross-spectral methods. Econometrica: journal of the Econometric Society, pp. 424 438, 1969.
Ruocheng Guo, Lu Cheng. Jundong Li, P Richard Hahn, and Huan Liu. A survey of learning causality with data: Problems and methods. ACM Computing Surveys (CSUR), 53(4):1-37, 2020.
David Heckerman, Christopher Meek, and Gregory Cooper. A bayesian approach to causal discov- ery. In Innovations in Machine Learning, pp. 1-28. Springer, 2006.
10

---

Magnus R Hestenes. Multiplier and gradient methods. Journal of optimization theory and applica- tions, 4(5):303-320 1969.
Patrik Hoyer, Dominik Janzing, Joris M Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. Advances in neural information processing systems, 21, 2008.
Biwei Huang, Kun Zhang, Jiji Zhang, Joseph D Ramsey, Ruben Sanchez-Romero, Clark Glymour, and Bernhard Schölkopf. Causal discovery from heterrgeneouss/nonttttionnnn data. J. Mach. Learn. Res., 21(89):1-53, 2020.
Aapo Hyvärinen, Kun Zhang. Shohei Shimizu, and Patrik O Hoyer. Estimation of a structural vector autoregression model using non-gaussianity. Journal of Machine Learning Research, 11(5), 2010.
Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparamcterization with gumbel-softmax, arXiv preprint arXiv: /6/ 1.01 144, 2016.
Saurabh Khanna and Vincent YF Tan. Economy statistical recurrent units for inferring nonlinear granger causality. arXiv preprint arXiv:1911.09879, 2019.
Hlyes Khemakhem, Ricardo Monti, Robert Leech, and Aapo Hyvarinen. Causal autoregressive flows. In International conference on artificial intelligence and statistics, pp. 3520-3528.. PMLR, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6988, 2014.
Steffen Lauritzen. Graphical models, volume 17. Clarendon Press, 1996.
Sindy Löwe, David Madras, Richard Zemel, and Max Welling. Amortized causal discovery: Learn- ing to infer causal graphs from time-series data. In Conference on Causal Learning and Reason- ing. pp. 509-525. PMLR, 2022.
Chao Ma and Cheng Zhang. Identifiable generative models for missing not at random data imputa- tion. Advances in Neural Information Processing Systems, 34:27645-27658, 2021.
Chris J Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous
relaxation of discrete random variables. arXiv preprint arXiv:1611.00 1.00712, 2016.
Daniel Marbach, Thomas Schaffter, Claudio Mattiussi, and Dario Floreano. Generating realistic in silico gene networks for performance assessment of reverse engineering methods. Journal of computational biology, 16(2):229-22) 239, 2009.
Raha Moraffah, Paras Sheth, Mansooreh Karami, Anchit Bhattacharya, Qianru Wang, Anique Tahir, Adrienne Raglin, and Huan Liu. Causal inference for time series analysis: Problems, methods and evaluation. Knowledge and Information Systems, pp. 1-45, 2021.
Pablo Aolales-Alvarez, Angus Lamb, Simon Woodhead, Simon Peyton Jones, Miltiadis Allamanis, and Cheng Zhang. Vicause: Simultaneous missing value imputation and causal discovery with groups. arXiv preprint 20X.0:21 2021.
Meike Nauta, Doina Bucur, and Christin Seifert. Causal discovery with attention-based convolu- tional neural networks. Machine Learning and Knowledge Extraction, 1(1):19, 2019.
Ignavier Ng, Sébastien Lachapelle, Nan Rosemary Ke, Simon Lacoste-Julien, and Kun Zhang. On the convergence of continuous constrained optimization for structure learning. In International Conference on Artificial Intelligence and Statistics, pp. 8176-8198 8198. PMLR, 2022.
Roxana Pamfil, Nisara Stiwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Geor- gatzis, Paul Beaumont, and Bryon Aragam. Dynotears: Structure learning from time-series, data. In International Conference on Artificial Intelligence and Statistics, pp. 1595-1605. PMLR, 2020.
Nick Pawlowski, Daniel Coelho de Castro, and Ben Glocker. Deep structural causal models for tractable counterfactual inference. Advances in Neural Information Processing Systems, 33:857- 869, 2020.
11

---

Jonas Peters, Joris Mooij, Dominik Janzing, and Bernhard Schölkopf. Identifiability of causal graphs using functional models. arXiv preprint arXiv: 2012.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Causal inference on time series using restricted structural equation models. Advances Neural Information Processing Systems, 26, 2013.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of causal inference: foundations and learning algorithms. The MIT Press, 2017.
Robert Prill, Daniel Marbach, Julio BadiiRzd,,, Peter K Sorger, Leonidas G Alexopoulos, Xi- aowei Xue, Neil D Clarke, Gregoire Altan-Bonnet, and Gustavo Stolovitzky. Towards a rigorous assessment of systems biology models: the dream3 challenges. PloS one, 5(2):e9202, 2010.
Jakob Runge. Causal network reconstruction from time series: From theoretical assumptions to practical estimation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(7):053100, 2018.
Jakob Runge. Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In Conference on Uncertainty in Artificial Intelligence, pp. 1388-1397, PMLR, 2020.
Jakob Runge, Peer Nowack, Marlene Kretschmer, Seth Flaxman, and Dino Sejdinovic. Detecting and quantifying causal associations in large nonlinear time series datasets. Science advances, r11)aaau4996, 2019.
Rajen D Shah and Jonas Peters. The hardness of conditional independence testing and the gener- alised covariance measure. The Annals of Statistics, 48(3):1514-1553 2020.
Shohei Shimizu, Patrik Hoyer, Aapo Hyvärinen, Antti Kerminen, and Michael Jordan. A linear non-gaussian acyclic model for causal discovery. Journal of Machine Learning Research, 7(10), 2006.
Ali Shojaie and George Michailidis. Discovering graphical granger causality using the truncating lasso penalty. Bioinformatics, 26(18):1517-i22 2010.
Elsa Siggiridou and Dimitris Kugiumtzis. Granger causality in multivariate time series using a time- ordered restricted vector autoregressive model. IEEE Transactions Signal Processing, 64(7): 1759-1773, 2015.
Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman. Causation, prediction, and search. MIT press, 2000.
Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaic, and Emily Fox. Neural granger causality for nonlinear time series. stat, 1050:16, 2018.
Brian I. Trippe and Richard Turner. Conditional density estimation with bayesian normalising flows. arXiv preprint arXiv: 180::8002000 04908, 2018.
Dennis Wei, Tian Gao, and Yue Yu. Dags with no fears: closer look at continuous optimization for learning bayesian networks. Advances in Neural Information Processing Systems, 33:3895-3906, 2020.
Tailin Wu, Thomas Breuel, Michael Skuhersky, and Jan Kautz. Discovering nonlinear relations with minimum predictive information regularization. arXiv preprint arXiv:2001.0 2020.
Chenxiao Xu, Hao Huang, and Shinjae Yoo. Scalable causal graph learning through a deep neural network. In Proceedingy of the 28th ACM international conference on information and knowledge management, pp. 1853-1862, 2019.
Cheng Zhang, Judith Bütepage, Hedvig Kjellström, and Stephan Mandt. Advances in variational inference. IEEE transactions on pattern analysis and machine intelligence, (0)):2008-0226 2018.
12

---

Kun Zhang and Aapo Hyvarinen. On the identifiability of the post-nonlinear causal model. arXiv preprint arXiv:1205.. 2599, 2012.
Kun Zhang, Zhikun Wang, Jiji Zhang, and Bernhard Schölkopf. On estimation of functional causal models: general results and application to the post-nonlinear causal model. ACM Transactions on Intelligent Systems and Technology (TIST), 7(2):1-22, 2015.
Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing. Dags with no tears: Continuous optimization for structure learning. Advances Neural Information Processing Systems, 31, 2018.
13

---

A ELBO AND LIKELIHOOD DERIVATION
The goal is to derive a lower bound for the joint likelihood po (x 0:7
x
0:2
(1)
X
(NN
Po
0:7
0:7
= log
po
(1)
0:7
(N) 0.7
G
p(G)dG
log
4s (G) P6
2o GG)
(X
IG) P(G)&G
2 qo (G) log
(X x (o)T)))
p(G)dG+ ((G)))G+((((GG)) H())(G))
(14)
(n) (G) log po X
G) log
H(qg(G))
where Eq. (14) is obtained by using Jensen's 's inequality.
We can further simplify the likelihood pol X
(n)
G¹):
T
Pog X (n) log
0:7
G)
pa
X
(z)
X (n)
G)
X ( G)
(n)
logpo Pa
T
log (
Pa t), (t)
(15)
t=
where Eq. (15) obtained through Markov factorization (Lauritzen 1996).
B STRUCTURAL IDENTIFIABILITY
In this section, we will focus proving the structural identifiability of Rhino SEMs. Before diving into the details, let us clarify the required assumptions.
Definition B.1 (Correctly specified model). For true data generating mechanism X = fī (Pa', e*) for i = 1, .., D, we say a model with functional space J is correctly specified if there exists function f: E F s.t. ft = f for all i == =.....
Here, we emphasize that the above definition of model specification does not require identifiability in parameters space, but in function space instead. Namely, it allows multiple sets of parameters that correspond to the same function. Our definition more general than some of the previous work, which enforces parameter identifiability (e.g. Def 2.1 in Ma & Zhang (202I)
Assumption (Causal Stationarity (Runge 2018).. The time series process X; with graph G is called causally stationary over time index set T if and only if for all links X + X in the graph
XL_ A X/ holds for E ET
This characterizes the nature of the time-series data generating mechanism, which validates the choice of the auto-regressive model.
Assumption (Causal Markov Property (Peters et al. 2017).. Given DAG G and joint distribu- tion p. this distribution is said to satisfy causal Markov property w.r.t. the DAG G each variable is independent of its non-descendants given its parents.
This is a common assumptions for the distribution induced by an SEM. With this assumption, one deduce conditional independence between variables from the graph.
14

---

Assumption 3 (Causal Minimality). Consider distribution and a DAG G, say this distribution satisfies causal minimality w.r.t. G it is Markovian w.r.t. G' but not to any proper subgraph of G.
Minimality is also common assumption for SEMs (Hoyer et al., 2008 Zhang & Hyvarinen 2012 Peters et al. 2012 which can be regarded as a weaker version of faithfulness (Peters et al. 2017).
Assumption (Causal Sufficiency). A set of observed variables V is causally sufficient for pro- cess e and only if in the process every common cause of any two or more variables V is in V or has the same value for all units in the population.
This assumption implies there are no latent confounders present in the time-series data. Assumption (Well-defined Density). We assume the joint likelihood induced by the Rhino SEM (Eq. (6) is absolutely continuous w.r.. Lebesgue counting measure and log p ((.oT;; G)| < co for all possible G.
This assumption to make sure the induced distribution has well-defined probability density function. It is also required for the equivalence of the global, local Markov property and Markov factorization property (Theorem 6.22 from Peters et al. (2017))
Assumption 6 (Rhino in DAG space). We assume the Rhino framework can only return the so- lutions from DAG space. Namely, the posterior distribution from Rhino can only put weights on DAGs.
This assumption regularizes the search space of Rhino to be DAG space, which aligns with assump- tion 2 and 3 on causal graphs.
In the following, will structure the entire proof into three steps:
1. Prove general conditions that the bivariate time series model needs to satisfy for structural identifiability. This adapts from the theorem in Peters et al. (2013).
2. Prove that a generalized form of SEM, modified from the post non-linear (PNL) model (Zhang & Hyvarinen| 2012), satisfies the conditions mentioned in step 1. The proposed Rhino Eq. (6) special case of the above SEM.
3. In the end, generalize the above indentifiability the multivariate case.
B.1 GENERAL IDENTIFIABILITY CONDITIONS
First, we derive the conditions required for identifiability for general bivariate time series SEM, defined as
X = fú, (Pa 1), (t))
(16)
We call the above SEM transition model, since it only defines the transition behavior rather than the initial conditions. We also need to incorporate a source model, which characterizes the initial conditions:
X ft. (Pa
(17)
for [0, S). where S is the length for the initial conditions and Pat contains the parents for node i. We define p. (Xo.s) ) as the induced joint distribution for the initial conditions.
Now, we prove the following theorem.
Theorem (Identifiability conditions for bivariate time series). Assuming Assumption 1 are sat- isfied, given a bivariate temporal process Xo.T and Yo:T that are governed by the above SEM Bqq (16) with source model (Eq. (17)), then the above SEM for the bivariate temporal process is structural identifiable if the following conditions are true:
1. Source model fr.e is structural identifiable for all = and s [0.S].
2. The transition model (Eq. (16) is bivariate identifiable w.r. the instantaneous parents. Namely. if graph G induced conditional distributions p(X,, Y, Pa
X.Y (< t)). then IG'
X.Y such that * G' and the induced conditional p(X,, Y, Pag (< t)) = p for all t e [S 11,,, 1, 1.
15

---

X.Y where Pa X.Y (< is the union of the lagged parents of X, and Y under G, and Pagi (< t) the union of parents under G'.
Proof. We prove this by contradiction. Assume we have an induced joint distribution p( (ooT, Yo:) under G. and corresponding p under G". We further assume the above two conditions in the theorem are met and p = p but G + G'.
Thus, have DRL.pplp = O. Due to the temporal nature of the model, we can further decompose as the following:
DKL [plp]]
- -/ pl p(Xo:T, Yo.T) log
(oTT, Yo.T),
d oo.Tdoo
p (XTT, Yo.T)
-KL (Xuus, Yu.s) p (Xuss Yo:s)] + p(Xu.s, Yo.s Yuss))k Ys41. / Yu:s)
Do
((X++IT YSHIT Xo.s, Yo.s) Yo.s)]dXo:sdYGo.
-DKL.PP.lp.]
Σ S=1
Ep(Xor- Yo.- (KL (p(X,, Xo Yo.- Yi Yo--))
=0.
This means have DKI.(pslie [lp]]
and DKL. [p(X,, Y, Xox You >lp (x, X Yo--1))
0 almost everywhere. Inspired by the strategy used in (Peters et al. 2013). We consider the following three scenarios:
Disagree initial conditions We assume G and G' disagree the initial conditions. From the condition 1, we know the source model fi,s identifiable. Namely, we cannot find G / G' with disagreement on initial conditions such that DKL (Ps [p.]] 0. This is contradiction, meaning that G' and G' must agree on the connections between initial set of nodes.
Disagree on lagged parents only This means for all E [S 1, 7]. the instantaneous connections for G and G' are the same, and št 1,7) such that Pa Pax.Y (< t) + Pa can use a similar argument as the theorem in Peters et al. (203)). W.l.o.g, we assume under G, we have X,-- - Y, and there no connections between them under G'. Thus, from Markov conditions, we have
XXY (< t). We
Y, I. X Xox-11 U Yo..
! (Y,,
under G', where NDİ are the non-descendants of node Y, at some time t. However, from the causal minimality and proposition 6.16 in Peters et al. (017), we have
Y, L X,_, Xo.-1 YY.-IUPP UND \Y,,
under G. This means under this case, DKL. [p(X, (X,, V, Xox-1, (ex-1)ll(( Y1|X Yo:- # 0, which is a contradiction.
Disagree also on instantaneous parents This scenarior means 31 E [S + 1. such that they disagree on instantaneous parents. W.l.o.g. we assume X, Y, under and Y, X, under G'.
Let's define X Yu:-1 h. h/ C h contains the values of Pay(<< t) under G. h C h contains the parent values under G', and h.y. hxx accordingly. Thus, the induced conditional distributions from SEM Eg. (16) with G, G' are
p(X,, Yi U kt) and K(X,YYlh u
From the Markov conditions, have
p( (x,, Y1 Xo Yo. p(X. Y, Pa X.Y (<1))
16

---

Therefore, have
DKL [p(,,, Y:lh)) Y:(k)]
=0
=DKLlp(Xe Yil U ( u
for arbitrary h, which contradicts the strucutral identifiability w.r.t. the instantaneous parents.
In summary, with the two conditions, we cannot find G * G' such that the induced joint p( ((.u. Yo.T) p( Yor)). meaning that the SEMs defined as Eq. (16) and Eq. (17)
identifiable w.r.t bivariate time series.
Since can any identifiable static models characterize the initial behavior of the time series, we will focus on condition 2 for the transition model. In the following, we will show that generalization of PNL, called distory- PNL, satisfies condition under assumptions.
B.2 IDENTIFIABILITY OF HISTORY-DEPENDENT PNL.
First, we propose generalization of PNI. (Zhang & Hyvarinen, 2012 so that can be history- dependent. For multivariate temporal process Xo.r, we propose sestorydeepndder PNL as
x = Vic (fa (Pac(< t), (t) gar (Pab(<< t), cxe). Pa'((< t))
(18)
where Vu an invertible transformation w.r.t. the first argument. The main differences of the above SEM compared to typical PNL are (1) the invertible transformation Vie can be history dependent; (2) the inner noise distribution can also be history-dependent.
Next, we show the main theorem about its bivariate identifiability w.r.t. its instantaneous parents.
Theorem 4 (History-dependee) PNL Bivariate Identifiability). Assume Assumption 15 are satis- fied, all transformations in Eq. (18) and corresponding induced distributions -rddor differen- tiable. Given a bivariate temporal process Xo.r, Yo.r, then the history- PNL defined as Eq. (18) bivariate identifiable w.r.t its instantaneous parents (i.e. satisfy condition 2 Theorem 3), except for some special cases.
Proof. W.l.o.g. time E 1, we assume X; + Y, for instantaneous connection under G and Y + X2 under G'. We fix value h for their entire history Xo.-1 Yo:-1 h. With h, we further define their lagged parents as Pa (< t) = hg C h, Pa t) = h C h under G and Pa (< = K C h. Pa (< t) h under G'. Therefore, the SEM at time t can be written as
Y, (f Xi) g(hg ey), .1/)
(19)
and
X, = (Jhhz Y)) + gh ex) hax)
(20)
under G and G', respectively. Let's assume that their induced conditional distributions at time are equal (i.e. violating the identifiable condition (2) in Theorem 3):
log pXX,, Yo kb) Jogp(XX, Yh U
under G
under G'
From the Markov properties, the above equation is equivalent to logp(X,, Ya log p(X,, Y/
with fixed value h of the entire history.
Now, let's define
0x = and 3, =-=--(())
17

---

where we omits the dependence of 5- to k and vー to hyg is easy observe that we have an invertible mapping between (X,, Y2) and (t), Thus, from the change of variable formula, we have
log p( Xt, Yi hh) = Po (cr, B1(h) + log |J]
and
log logp(X, (X,, Y, hh) = log Po Br (h) + log [J]
where J is the Jacobian matrix of the transformation. Thus, the equivalence of log p and log p in the (X,, Yı) space can be translated (ar, space.
Thus, from Eq. (19) we have
= ((a)) + cyh
(21)
under G. And from Eq. (20)
we have
ar = (3:) + g( ex)
(22)
under G'. This forms an additive noise model between ar, B, with history-dependent noise. Next, we can use similar proof techniques in Hoyer al. (2008). Here, .b-) = f(h 1 gr,
and (-) = f vh). We further define
7/1 (a)) = log p(au(h)
m (3,) log h))
(g Ey)) log Py (g (g kh)
(2 ex)) log Py (g cx) h)
Thus, under G (i.e. Eq. (21)). have
log Be(h) = logg(BB (Br h) + )) =12 2(a:)) + (at)
(23)
Similarly, under G' (i.e. Eq. (22). have
log pl Be) 72(00 (T())) (3.)
(24)
Based on Eq. (24)
we have
8 log 03,03
0 log p da,
V and
7/2
Thus, we have
0
8 log logg//a,,03,
Oot
82 log togp/o Do,2
Due to the equivalence of log and log we apply the above operations to Eq. (23) algebraic manipulation, obtained the following differential equations for 1/2 40 0:
After some
171
7/1 "b"
1/2222 21722
4"
1/2 4 4'Th' T
7/2
nh
($")2
(25)
T/
Interestingly, this is exactly equivalent (q. in Zhang & Hyvarinen (012))
The main difference
the definition of variables and transformations in here are all destodedd
Further, also observe that
Br 1 oh ey ) Xor- U You-1 h.
Since 3 +(a)) + ghi,, ey) and g( ex) a₂ (3t,, it is trivial to show the determinant of the Jacobian of the transformation (au, g) (Be, g) is Thus, by similar argument in theorem 1 from Zhang & Hyvarinen (2012), we can derive
" 1/2 (b))² 1/2 1"
1/2 T
for 4' * 0.
Thus, the above two differential equations has the same form as theorem 1 in Zhang Hyvarinen 2012) where the main difference is that all distributions and transformations involved in depend on history h.
Therefore,
directly cite the theorem 8 from Zhang & Hyvarinen (012). which proves that the above differential equations hold true only for 5 types of special cases. One can refer to Table in Zhang & Hyvarinen (2012) for details.
n
18

---

Corollary 10 from Zhang & Hyvarinen (2012) validates the choice of using nueral network for the transformation f. For completeness, we include it here with slight modification:
Corollary 4.1 (Identifiability with neural netowrk f). Assuming the assumptions Theorem 4 are true. and the double derivative (log Py (g (Pa) (<t), ey
w.r. ey is zero at most at some discrete points. If function is not invertible w.kf. the instantaneous parents, then, the history-dependent PNL defined as EEq. (18) is bivariate identifiable w.r.. the instantaneous parents (i.e. satisfy condition 2 in Theorem 3).
clear to see that the form of Rhino SEMs Eg. (6) is a special case of the history-dependent PNL (Eq. (18)), where the outer history-dependent invertible transformation is the identity map- ping. Thus, we can directly leverage [Theorem ttogether with [Theormm 4jt show Rhino SEMs are identifiable w.r.t bivariate time series, and Corollary 4.1 to validate our design choice (Eq. (7)).
B.3 GENERALIZING TO MULTIVARIATE TIME SERIES
Previously, we prove the identifiability conditions for bivariate time series. In this section, we will generalize the multivariate case.
Theorem 5 (Generalization to multivariate time series). Assuming the assumptions in Theorem 4 are satisfied, we further assume that the multivariate SEM defined in Eq. (18) for each pair of node i E V. the SEM
X3 = va fa
Pa)(tt))))
+ (Pab(< (), t), Pab(<
bivariate identifiable w.r.t. the input, and an identifiable source model is adopted. Then, the history-dependent PNL is identifiable except for some special cases.
Proof. For this proof, we can follow the strategy used in Theorem and Peters et al. (2013). We categorize the difference of the graph G and G' into three types. Following the same analysis of the KL divergence of the two induced joint distributions, we can see that (1) DKL pps pp. = 0 and DKL (p( X, X (ilett1))p((,,,,, xxx) O.
Disagree on initial conditions diets [KL [p..)
Since we assume that the source model is identifiable, this contra-
Disagree lagged parents only We notice that the analysis used in Theorem for this disagree- ment can be directly translated to multivariate case. The only difference is that the notation Y,, X₁ changed accordingly.
Disagree also on instantaneous parents For this case, with a fixed history value f Xoa-11 the aim is compare the conditionals DkL[p( X₂ X = h)) p))p XX:t-1 h)]. Thus, the problem becomes to how to generalize the bivariate identifiability for instantaneous parents to the multivariate case. We leverage the theorem 2 from Peters et al. (2022), which proves the multivariate identifiability for any models that belongs to IFMOC. It is easy to see that if the assumptions in Theorem 3) are met, the history-dependent PNL. belongs to IFMOC w.r.t. the instantaneous parents. It should be noted that the entire history deeeede-de PNL DOES NOT belong to IFMOC, but this does not affect our results since we only care about the instantancous parents under this case.
C RELATION TO OTHER METHODS
VARLINGaM (Hyvärinen et al. 2010) VARLiNGaM (Hyvärinen et al., 2010) is a causal dis- covery method for time series data based on the linear vector auto-regression, which can model both lagged and instantaneous effects. Its SEM is defined as Eq. (5) where the noise ei is an independent non-Gaussian noise. It is easy to observe that this is special fi as the matrix multiplication of the weighted adjacency Go.k with the nodes, and 9i as the iden- tity mapping. For the training objective, VARLiNGaM adopted two stage training to sidestep the
of Rhino (Eq. (6) by setting
19

---

difficulty of directly optimizing the log likelihood. From the Theorem for Rhino, we note that the solution from optimizing the variational objective equivalent to maximizing the log likelihood under infinite data limit. Therefore, by setting large enough DAGness penalty coefficient a, p. the inferred graph from both methods should be equivalent.
DYNOTEARS (Pamfil al., 2020) The formulation of DYNOTEARS is the same VAR- LiNGaM, which is based on linear vector auto-regression. The main novelty is the usage of the DAGness penalty h (G),, which continuously relaxes the DAG constraint. The training objective is the mean square eтror with augmented Lagrange scheme for DAGness penalty, Thus. it is obvious that DYNOTEARS is a special case of Rhino with linear transformations and identity g.. Similarly, Theorem 2 shows the connections between the variational objective and maximum likelihood, which is equivalent to mean square error if the noise distribution Gaussian with equal variances.
cMLP cMLP Tank al. 2018 combines Granger causality with deep neural networks. The model formulation is
X fi(X.:- x +
where Ji is a function based on MLP. Although the input is the entire history. the one that matters the node that has the connection X (i.e. lagged parents). Therefore, it is easy see they closely related Rhino without instantaneous parents (ta and history-dependent noise. Since the training objective of cMLP is based on the mean square error with sparseness constraint, by the same argument as before, the variational objective is equivalent to mean square error with equal variance Gaussian noise and large training data.
TiMINo (Peters et al. 2013) TiMINO is most similar to our work among all the aforementioned methods in terms of model formulation. TiMINo proposed a very general formulation based on IF- MOC (Peters al. 2012) and showed the conditions for structural identifiability. Rhino generalizes the TiMINO in a way such that noise history dependency can be incorporated. Thus, Rhino only belongs to IFMOC w.r.t. the instantaneous parents. Therefore, Rhino without the history-dependen noise is TiMINo model. The training objective of TiMINo is based on the dependence minimiza- tion between the noise residuals and causes, and can only infer summary graph instead of temporal causal graph. /Zhan al. (2015) proved the equivalence of the mutual information minimization maximum likelihood, which is equivalent to our variational objective under infinite data.
D TREATMENT EFFECT ESTIMATION
We now show how to leverage the fitted Rhino for estimating the conditional average treatment effect (CATE). For simplicity, only consider a special case of CATE defined
CATE(a, b) = (g)(G) ( - (Xx,x
(X)
(26)
We assume the conditioning variable can only be Xee (i.e. the entire history before (), and the intervention and target variable can only be either at current time or sometime in the future f+r. We emphasize that this formulation is for simplicity, and Rhino can be casily generalized to more cases as Geffner et al. 2022). Once fitted, the idea distribution p( X X X Kee do)))) G) for each graph sample G (G) Then, unbiased Monte Carlo estimation can be used to compute CATE. For sampling from the interventional distribution, we can use the "multilated" graph G to replace G, where all incoming edges to X are removed. The intervention samples can be obtained by simulating the Rhino with history X X = a or b and Gdo(X
draw target samples X from the interventional
D.1 CAUSAL INFERENCE RESULTS
Here, provide the preliminary results for CATE performance of Rhino by calculating the RM- SEs of the estimated CATES comparing to the true CATE from the interventional samples (lower is better). We present boxplots of the performance in Fig. 2 All Rhino-based method perform sim- ilarly. Surprisingly, the CATE performance seems to have little correlation to the causal discovery performance and warrants further study in the future.
20

---

Dim 10
RMSE
Model Name
Rhino (L=2)
Dim 20
Dim 40
Rhino+ (L=2)
Rhino+s (L=2)
RMSE
E[ATE]
ML ATE
E[ATE]
ML ATE
Figure 2: Comparison of the RMSE of the average treatment effects (CATEs) the different in- stantiations of Rhino depending the dimensionality. E[CATE] refers RMSE of the expected CATE over the posterior graph distribution (i.e. G ~ (G) ML ATE uses the most likely graph to calculate the ATE. These results are obtained averaging 160 datasets, similar to the discovery setup.
E VARIATIONAL DISTRIBUTION FORMULATION
Here provide the detailed formulation of the independent Bernoulli distribution 4o (G). Since this distribution is responsible for modelling the temporal adjacency matrix Ga.K. we use Σg to represents the edge probability in GR.- We further split the edge probability matrices into the instan- taneous part Σg and lagged parts E1:K-
To avoid the constrained optimization of E (i.e. the value needs to be within [0, 1)). we adopt the following formulation:
exp(u).is) Th,ij (xp(u + exp("k.3))
(27)
where Uk.6j UR Uk.y) V4 and Uk. Vx E DxD for all I,..., K. Since we do not require lagged adjacency matrix to be DAG, Us. V₁ has constraints during optimization.
On the other hand, Go needs to DAG for instantaneous effect. By smart formulation, we can get rid of the length-1 cycles. The intuition is that for pair of node i, only three mutually exclusive possibilities can exist: (1) i + j; (2) j + i; (3) no edge between them. Thus, instead of using a full probability matrix Eg. we three lower triangular matrices Un, Vo and E, to characterise the above three scenarios. For node > j.
exp(ur)
p(i + +1) exp(u;)) + exp((is) exp(eii)
exp(Uig)
p(j + i) = exp(u)) + exp(U13) (1ij) <xx((ci)
exp(e:))
p(no edge) exp(u.)) exp((;;) exp(ei)
Thus, by this formulation, the corresponding instantaneous adjacency matrix will not contain length- 1 cycles.
21

---

F SYNTHETIC EXPERIMENTS
F.1 DATA GENERATION
We create the synthetic datasets in a four step process: 1) generate random Brdös-Réénl) (ER) or scale-free (SF) graphs that specify the lagged and instantaneous causal relationships; 2) drawing random MLPs for the functional relationships as well as random conditional spline transformation modulate the scale of the Gaussian noise variables e; 3) sample initial starting conditions and follow Eq. (2) with the additive noise to simulate the temporal progression; 4) removing the burn-in period and return stable timeseries. We consider four different axes of variation for the data genera- tion: number of nodes Nnodes E [5, 10, 20, 40): ER or SF graphs: instantaneous or no instantaneous effects; and history-dependent or history-independent noise (i.e. Gaussian noise). All combinations generated with 5 different seeds, yielding 160 different datasets. Datasets with instantancous effects have X Nodes edges in the instantaneous adjacency matrix. All datasets have 2 x Nnodes connections in the lagged adjacency matrices. The MLPs for the functional relationships are fully- connected with two hidden layers,64 units and ReLU activation. In case of history- -istery-ideeeedeent noise, we are using Gaussian as the base distribution. The history dependency is modelled product of a scale variable obtained by the transformation of the averaged lagged parental values through a random-samples quadratic spline, and Gaussian noise variable.
The datasets with 40 nodes are generated with a series length of 400 steps, a burn-in period of 100 steps, and 100 training series. All other datasets are generated with a time-series length of 200, burn- in period of 50 steps and 50 training series. We generate random interventions for all the datasets by setting the treatment variable 10 for intervention and -10 for reference. 5000 ground-truth intervention samples are used to estimate the true treatment effect.
F.2 METHODS
All benchmarks for the synthetic experiments are run by using publicly available libraries: VARLiNGaM Hyvärinen et al. (2010) is implemenented in the lingan python package. 5 We use the implementation in et al. 2020 We use the default parameters for all
PCMCI* (elcrt(Rungee 2020) is implemented in Ti Tigrami
causalnex to run DYNOTEARS these baselines. For PCMCI+ enumerate all graphs in the Markov equivalence class to evaluate the causal discovery performance (see Appendix G.2 for details).
For Rhino and its variants, we use the same set of hyper-parameters for all 160 datasets demon- strates our robustness. By default, we allow Rhino and its variants to model instantaneous effect; set the model lag be the ground truth except for ablation study; the ga (G) is initialized to favour sparse graphs (edge probability< < 0.5); quadratic spline flow is used to for history-dependent noise. For the model formulation, we use 2 layer fully connected MLPs with 64 (5 and 10 nodes), 80 (10 nodes) and 160 (40 nodes) for all neural networks in Rhino-based methods. We also apply layer normalization and residual connections to each layer of the MLPs. For the gradient estimator, we use the Gumbel softmax method with a hard forward pass and a soft backward pass with tempera- ture of 0.25. All spline flows uses 8 bins. The embedding sizes for transformation (i.e. Eq. and conditional spline flow) is equal to the node number.
For the sparseness penalty A in Eq. (10) we use for Rhino and Rhino+s, and 5 for Rhino+g. We set p 1 and a 0 for all Rhino-based methods. For optimization, we use Adam (Kingma & Ba, 2014) with learning rate 0.01. The training procedure follows from Appendix B.1 in Geffner et al. 2022
F.3 ADDITIONAL CAUSAL DISCOVERY RESULTS
different type of graphs The first study is test our model robustness to different types of graphs. Fig. 3/shows the discovery performance over ER or SF graph averaged over all other possible data setting combinations. Most methods perform better on ER graphs than on SF graphs,
Ablation:
ssee ht
see
ht htpps///kkoobuuuggg.........bbbb ttt.tiiiiiiiiiiggiii///e/
'see
htpa://causalrex.eeaaa www.aaeeeeeeeeeeeeeeeeeeeeee httptthheeoo..ooo///aaeeeeee
22

---

Type
5| Type
1.0
0.5
0.0
Dim 10 Type ER
Dim = 10 Type = SF
1.0
0.5
Model Name
Rhino (L=2)
Rhino+g (L=2)
Dim 20 Type ER
Dim 20 Type
SF
Rhino+ss (L=2)
1.0
DYNOTEARS
VAR-LiNGaM
[ 0.5
PCMCI
0.0
Dim = 40 Type = ER
Dim = 40 Type = SF
1.0
0.5
0.0
Inst.
Lag
Temporal
Inst.
Lag
Temporal
Figure Comparison of the F, score of the different baseline methods as well as Rhino (light yellow) depending on the dimensionality and the graph type. Inst. refers to the performance on the instantancous adjacency matrix, Lag refers to the lagged adjancency matrices and temporal considers the full temporal matrix.
with only DYNOTEARS (Pamfil et al. 2020) as an exception. We note that the PCMCI* runs on SF graphs with 40 nodes exceed our maximum run time of week, showing its computational limitation in high dimensions. Nevertheless, Rhino achieves consistent performance throughout all graph settings.
Ablation: history dependency Figure 4 explores the performance difference of all methods on data generated with/without history-dependent noise. Interestingly, most methods perform better on the history-dependent datasets than the history-independent ones. The possible reasons are (1) the difficulty of the discovery also depends on the randomly sampled functions; (2) the default hyperparameters of all methods are initially chosen to favor the datasets with history-dependent noise and instantaneous effects. We find that PCMCI+ is the most robust across both settings, followed by Rhino and DYNOTEARS. On the other hand, the two variants of Rhino seems to be less robust. When the Rhino is correctly specified, it achieves the best performance. summary, Rhino demonstrates reasonable robustness to history- deporyyde mismatch and achieves the best when correctly specified.
Ablation: instantancous effect We investigate the impact of instantaneous effects in the data. Figure 5]shows the F2 score averaged over all possible setting combinations other than instantaneous effect. All methods seem to be robust across both settings with PCMCI+ and Rhino performing the best. The score of the instantancous adjacency matrix when instantaneous effects are disabled not defined and therefore not plotted.
23

---

Dim Hist.-Dep.
Dim 5 Hist.-Dep.
1.0
0.5
0.0
Dim 10 Hist.-Dep. False
Dim 10 Hist.-Dep. =
True
1.0
0.5
Model Name
Rhino (L=2)
Rhino+g (L=2)
Dim Hist.-Dep. False
Dim 20 Hist.-Dep. True
Rhino+ss (L=2)
1.0
DYNOTEARS
VAR-LiNGaM
[ 0.5
PCMCI
0.0
Dim 40 Hist.-Dep. False
Dim 40 Hist.-Dep. = True
1.0
0.5
0.0
Inst.
Lag
Temporal
Inst.
Lag
Temporal
Figure 4: Comparison of the F, score of the different baseline methods as well as Rhino (light yellow) depending on the dimensionality and whether the data is generated with history-depence or not. Inst. refers to the performance on the instantaneous adjacency matrix, Lag refers to the lagged adjancency matrices and temporal considers the full temporal matrix.
G REAL-WORLD EXPERIMENT DETAILS
G.1 DREAM3 HYPERPARAMETER SETTING
For tuning the hyper-parameter parameters of Rhino, its variants and DYNOTEARS, split each of the datasets into 80% 20% training/validation. We tune Rhino and its variants based on the validation likelihoods, and DYNOTEARS based the validation RMSE emror. For PCMCI+ fault settings recommended in the Tigramite package (https: : /github...../ com/ jatp://jkobunggee
the de-
tigramite) For other Granger causality baselines, refer Table 7-11 in Khanna Tan (011)).
Other than the hyper-parameters reported in Table 4 we use 1-layer MLPs with 10 hidden units for both ( <; Eq. (7) and the hyper-network for conditional spline flow bins). All the MLPs use residual connections and layer-norm at every hidden layer. We use linear conditional spline flow (Dolatabadi et al. 2020) instead of the original quadratic version (Durkan al.) (019) for better training stability. We also initialise the Bernoulli probability (G) to favour dense graphs (i.e. and a For the gradient estimator, we use the Gumbel softmax method with hard forward pass and a soft backward pass with temperature 0.25, We use batch size 64, learning rate 0,001 with Adam optimizer (Kingmal & Ba 2014). The training procedure follows from Appendix B.1 in Geffner et al. 2022
edge probability 0.5). For prior p( p(G), we set the initial value
Table contains the parammersssserss setup for DYNOTEARS. We set the maximum training iter- ations be 1000 with DAGness tolerance 10 - The threshold value for the weighted adjacency matrix is 0.05. For PCMCI* the maximum lag is set to The conditional independence test is set parcorr, which based linear ordinary least square (OLS). A powerful choice can be
24

---

5 Inst. Effect False
Dim Inst. Effect True
1.0
0.5
0.0
Dim 10 Inst. Effect = Falss 1.0
Dim 10 Inst. Effect True
0.5
Model Name
Rhino (L=2)
Rhino+g (L=2)
Dim 20 Inst. Effect
Fallse
Dim Inst. Effect True
Rhino+ss (L=2)
1.0
DYNOTEARS
VAR-LiNGaM
[ 0.5
PCMCI
0.0
Dim 40 Inst. Effect = False 1.0
Dim 40 Inst. Effect True
0.5
0.0
Inst.
Lag
Temporal
Inst.
Lag
Temporal
Figure 5: Comparison of F₁ score different baseline methods as well as Rhino (light yel- low) depending on the dimensionality and whether the data is generated with instantaneous effects or not. Inst. refers to the performance on the instantancous adjacency matrix, Lag refers to the lagged adjancency matrices and temporal considers the full temporal matrix.
Hyperparams
Node Embedding
Instantaneous eff.
Node Embed. (flow)
lag
Auglag
Rhino (Ecolil)
16
False
16
2
19
30
Rhino (Ecoli2)
16
False
100
25
80
Rhino Yeast
32
False
100
2
25
10
Rhino (Yaast2)
False
100
25
80
Rhino (Yeast3)
32
False
16
25
5
Rhino+g (Ecolil)
100
False
N/A
15
60
Rhinotg (Ecoli2)
100
False
N/A
2
25
25
Rhinotg (Yeastl)
100
False
N/A
15
5
Rhino+g (Yaast2)
100
False
N/A
19
125
Rhino+g (Yeast3)
100
False
N/A
9 10
Table The hyperparameter setup for Rhino. Node embedding the dimensionality of u, be- low Eq. (7), Instantaneor eff. specifies whether it models the instantaneous effect or not; Node Embed. (flow) represents the dimensionality of the node embedding for the hyper- network used for conditional spline flow yi since the hyper-network shares the similar structure as Eq. (7) lag defines the model lag order; and X, is the sparseness penalty in the prior (Fq. (10); Auglag is the number of augmented Lagrangian steps, each step consists of 2000 training itera- tions.
a nonlinear independence test based on GP, called GPDC, However, PCMCI with GPDC is too slow finish the training.
25

---

Hyperparams
lag
Xa
Ecolil
0.01
0.5
Ecoli2
2
0.1
0.01
Yeastl
2
0.005
0.1
Yeast2
0.01
0.01
Yeast3
2
0.01
0.005
Table The hyperparameter setup for DYNOTEARS.
Metries
Ecolil
Ecoli2
Yeast
Yeast2
Yeast3
Rhino
0.183= 0.012
0.2141 0.009
0.261 0.006
0.1361: 0.005
0.1133
0.005
Rhino+g
0.169- 0.010
0.211+0.0010
0.254 0.006
0.148 0.007
0.126 0.003
Dynotears
0.120
0.066
0.059
0.092
0.045
PCMCI
0.051=0.006
0.049 0.003
0.068.10.005
0.046 0.002
0,060
Table 6: Orientation F1 score of DREAM3 datasets.
G.2 POST-PROCESSING TEMPORAL ADIACENCY MATRIX
The ground truth graphs for DREAM3 and Netsim datasets are summary graph, which is essentially the temporal graph aggregated over time. We provide a formal definition of summary graph:
Let X; be a multivariate temporal process, and G = (V, EE) be summary graph. The edge + exists if and only if there exists some time t and some lag r such that X causes X at time t with lag 0 i for p + q and with a time lag of 0 < i for q.
Definition G.1 (Causal summary graph (Assaad et al. 2022
Unlike the some of the Granger causality baselines, Rhino (and its variants), DYNOTEARS, VAR- LiNGaM produces the temporal adjacency matrix after training. For DREAM3 and Netsim datasets, this creates the incompatbiliity during evaluation. Thus, we need to aggregate the temporal graph into a summary graph before comparing to the ground truth. For binary adjacency matrix, we sum over the time steps followed by step function, i.e. step Σs G3). Thus, there will be edge in summary graph as long as there is a connection from to j at any timestamp. For the Bernoulli probability matrix from Rhino and its variants, we take max(-) over the timestamp to generate the probability matrix for the summary graph.
An exception is PCMCI+ which only produce MECs for the instantaneous adjacency matrix. In such case, we will enumerate up to 10000 possible instantaneous DAGs from the MECs. Together with the lagged adjacency matrix, we will perform the above post-processing step to generate the corresponding aggregated adjacency matrix. We also estimate the corresponding edge probabilities by taking the average over all possible DAGs.
For DREAM3 experiments, we ignore the self-connections by setting the diagonal of the aggregated adjacency matrix be
For Netsim, self-comnections are not ignored, following the same settings as Khanna & Tan 2019
G.3 ADDITIONAL DREAM3 RESULTS
shows the additional ROC curve plots for all 5 datasets in DREAM3. For the visual- ization purpose, we only select a single run for Rhino and this will not affect the curve much due to small standard error in Table 2
Here, Fig.
In addition, we provides the additional metrics (Orientation Fl and SHD) for DREAM3 datasets in Table 6 and Table 7 These results are obtained by using the same hyperparameters mentioned in Appendix G.I In particular, we tune the threshold for rounding the continuous- valued adjacency matrix to binary adjacency matrix for both Rhino and baselines. It be observed that the Fl and SHD agree with the trend of AUROC, where Rhino and Rhinotg achieve the best performance compared to baselines. This further supports the advantages of our proposed methods.
26

---

eSRU
TCDF
cMLP
DYNOTEARS
Rhino +
SRU
CLSTM
Rhino
pcmci+
Ecolil 100
Ecol12 100
1.0
1.0
0.8
0.8
0.6
0.6
TPR
0.4
0.4
0.2
0.2
0.0
0.0
0.5
FPR
FPR
Yeast1_100
Yeast2 100
1.0
0.8
0.2
0.2
0.0
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
FPR
FPR
1.0
0.8
0.6
0.4
0.2
0.0
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
FPR
Figure 6: The ROC curve plots Rhino and other basclines for DREAM3 datasets. For illustration purpose, we only select single run of Rhino, Rhino+g, DYNOTEARS and PCMCI to plot ROC curve. Since the standard error reported in Ta/Table i relatively small, the plot should not vary much for other runs. The ROC curve of other baselines are directly taken from figure in Khanna & Tan (201
27

---

Metrics
Ecolil
Hcoli2
Yeast
Yeast2
Yeast3
Rhino
157.8=13.23
120.4 1.19
168
6.7.8
627.6+15.47
Rhino+g
162.8=33113
122
182
401.4+3.85
760.6+66.1 18
Dynotears
372
422
397 288.6.+-2.24
630
838
PCMCI
266.6-12.01
273.2++1.48
50..4+00.2
704-100
Table 7: SHD of DREAM3 datasets.
G.4 NETSIM HYPERPARAMETER SETTING
For the Netsim experiment, we extract subject 2-6 in Sim-3.mat to form the training data and use subject 7-8 validation dataset. Following the same settings DREAM3 (Appendix G.T., we tune the hyperparameters of Rhino and its variants based on the validation log likelihood; DYNOTEARS with MSE on validation dataset; and use default settings of PCMCI* from Tigramite package.
It is worth noting that unlike DREAM3 experiment, where the results and hyperparameters of Granger causality baselines can be directly taken from Khanna Tan (2019). Their setup of Netsim experiment is different from ours, where they train the baselines using single subject and compute the corresponding AUROC, followed by averaging over subjects 2-6. Our setup is to train all meth- ods using the entire data from subject 2-6 before computing AUROC. Thus, the hyperparameters for Granger causality are slightly different, and the AUROC increases for the baselines compared to those reported in Khanna & Tan (2019).
Rhino The hyperparameters are the same as DREAM3, except for the following: we initialise the Bernoulli probability of qa (G) have preference (i.e. edge probability= 0.5); the A 25; we use 2 layer MLPs with 64 hidden units for both functional model (Eq. (7)) and hyper-network with embedding size 15; the augmented Lagrangian step is 5, For Rhino variants, we use the above settings as well.
DYNOTEARS, PCMCI* and VARLINGaM For DYNOTEARS, we set lag to be 2, Aa = 0.5 and An 0.5. For PCMCI* we use parcorr independence test with lag 3. For VARLiNGaM, we use lag with default settings as
readthedocs io/en/////iiie //ttess/
Granger Causality For computing AUROC, we follow the same method as Khanna & Tan (019); Tank et al. (018) by sweeping through range of hyperparameters. Specifically, we use the same hyperparameters for SRU and eSRU as (Khanna & Tan 2019). For cMLP, choose the ridge penalty as 0.43 and sweep through the group sparse penalty in range [0.1, 1]. For cLSTM, we set the ridge penalty be 0.045, and sweep the group sparse penalty in range [0.1, Fo. TCDF, we sweep through the threshold in range [-1,2] 2] for the attention scores. Other than the above hyperparameters, everything else follows the setup as in Khanna & Tan (2019)
G.5 ADDITIONAL NETSIM RESULTS
Figure 7 shows the ROC curve plot for Rhino and other baselines. It is clear that Rhino achieves significantly better TPR-FPR trade-offs compared to others. [Table 8 shows additional discovery metrics of Rhino and baselines for Netsim dataset. We can observe that F1 score and SHD in general agree with AUROC reported in Table where Rhino-based methods outperform the baselines, apart from Rhino+Nolnst. This again confirms the necessity of modelling instantaneous effect for real- world challenges. Rhino outperforms Rhinotg on two out of three metrics (including AUROC), which shows the advantage of history-dependent noise.
H AUROC METRIC
AUROC metric is one of the standard metrics for evaluating the causal discovery, which measures the trade-off between the true positive rate (TPR) and false positive rate (FPR). However, during the experiments, we found out that AUROC does not necessarily correlate well with other discovery metrics. From Fig. 8, is clear that the F, score continues to increase whereas AUROC and vali- dation likelihood starts decrease after few steps. Since the dataset of Netsim is relatively small,
28

---

Method
Ori. F1
SHD
DYNO.
0.341
PCMCI
0.41
18
VarLiNGaM
0.44
18
Rhino+g
0,539 + 0.036
10.4 + 1.08
Rhino+NoInst
0.212 +: 0.014
29 + 1.5
Rhino
0.551 + 0.048
13.8 + 1.5
Table 8: Orientation F1 and SHD of Rhino and baselines for Netsim dataset.
Netsim Roc Curve
0.8
Rhino
Rhino+g
DYNOTEAR
TPR
VARL NGaM
cMLP
0.4
cLS TM
TCDF
0.2-
0.0
0.0
0.2
0.4
0.6
0.8
1.0
FPR
Figure 7: The ROC curve plots of Rhino and other baselines for Netsim dataset. Similar to Fig. 6 only select run out of 5 for Rhino, Rhino+g, DYNOTEARS, PCMCI- for illustration purpose.
29

---

AUROC/val likellhood v.5. Score
-30
40
0.8
AUROC
Orientation F1 Val. likelihood
-50
va likellhood
FL/AUROC
60
0.4
-70
0.2
-80
10
30
40
Steps
Figure 8: The curves of orientation Fi. AUROC and validation likelihood during training. Each curve is obtained by averaging over 5 random seeds. The validation curve agrees well with the AUROC curve, but shows opposite trends as F1 curve. This potentially indicates model overfitting in the later stage of training.
this indicates the possible overfitting. This disagreement originates from the different aspects these metries care about. For AUROC, cares about the trade-off between TPR and FPR with various decision thresholds, and it penalizes the wrong decisions with certainty harshly. On the other hand, Fi score cares about the final inferred binary adjacency matrix with fixed decision threshold. For example, if multiply the Bernoulli probability matrix by small factor (e.g. 10-5), the AUROC score will remain the same but the F1 score will tends to with the default decision threshold 0.5.
Thus, model overfitting tends drive the edge probabilities towards score but these extreme decisions can result in a large decrease in the AUROC score. Thus, for small dataset, we believe AUROC is better metric than F,, which also agrees with validation likelihood.
O, which may help the Fj
In addition, the Bayesian setup of Rhino may also help with better AUROC for small dataset. From the same figure, even the large decrease of validation likelihood suggests potential model overfitting, the AUROC still maintains a reasonable value. This may be due to the Bayesian view of the causal graph, where the posterior edge probability does not converge to extreme values.
30


Reviews:
[{"review_id": "ouji7c-Yn1", "reviewer": "Reviewer_4WVg", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "In this paper, the authors investigate the problem of temporal causal discovery. Existing work relies on the absence of instantaneous effects with fixed noise distributions, whereas constraint-based methods require stronger faithfulness assumptions. To alleviate this problem, the authors propose Rhino, which combines vector auto-regression and deep learning, and is trained by variational inference. They further demonstrate that Rhino is structurally identifiable under similar assumptions to ANMs. Experiments are conducted on two discovery benchmarks.\n", "summary_of_the_review": "The paper is well written. The proposed method is novel. ", "strength_and_weaknesses": "Strength: \n1. The authors propose to combine autoregression and deep learning, and train the network using variational inference.\n2. The authors demonstrate that the proposed Rhino is structurally identifiable under similar assumptions to ANMs.\n3. The paper is well written. \n \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "m1fdVIsA94n", "reviewer": "Reviewer_wW3u", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The author proposed a new structural equation model, named Rhino, that can be used to model causal relations between variables in temporal data. Compared to ANM, the proposed model is more flexible, but there is no sacrificing in identifiable properties. The proposed Rhino model is actually identifiable under some common assumptions. The experiments show that the proposed approach can be applied in various different applications to produce satisfactory results. ", "summary_of_the_review": "The authors proposed a new SEM, named Rhino, for temporal data, and proved the identifiable of the proposed model. Experiments on data from different domains show the effectiveness of the proposed model.", "strength_and_weaknesses": "Strength:\n\n1. The motivation of the paper is clearly stated, and the writing of the paper is good. It is easy to follow. Even for the hardest proof part, the sketch of the proof is easy to follow.\n\n2. The experiment part is solid. The experiments are based on datasets from different domains, and all of them show the effectiveness of the proposed model.\n\n3. The mathematic derivation of the paper is solid.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "OpAQLoHw2j", "reviewer": "Reviewer_w2Xd", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper is trying to deal with non-linear relationships in time series data and proposes a structural equation model called Rhino, which combines vector auto-regression, deep learning as well as variational inference techniques.", "summary_of_the_review": "Overall, I think the paper is of good quality, with some improvements required on the experiments and writing. But the overall the methodology and theory part is good.", "strength_and_weaknesses": "Strength\n1. The problem causal discovery the paper is trying to tackle is quite important and interesting.\n2. The method is supported by the structural identifiability theory. \n\nWeakness\n1. I do not understand why related works is not Section 2, why Section 5. This is quite confusing.\n2. In the experiments section, I was wondering if we should only use AUC as the evaluation metrics since one major goal seems to be discover the causal relationships.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "DuZUQX5M7j", "reviewer": "Reviewer_hiTa", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper provides a new structure identifiability result for temporal causal discovery, specifically time-series structural equation model with nonlinear lagged and instantaneous effects. The key advancement is that the identifiability result allows for history-dependent noise, while existing works typically assume independent noise. The paper further proposes a structure identification method based on continuous optimization, variational inference, and deep learning. Empirical results for synthetic and real world data are further reported.", "summary_of_the_review": "The setting considered with history-dependent noise is well-motivated. Theoretical results are provided to justify the proposed method. However, there are some concerns with Theorem 2.", "strength_and_weaknesses": "Strengths\n- The identifiability result and identification method both seem to be sound.\n- Theoretical and empirical results are provided.\n- The motivation of using history-dependent noise is convincing as far as I know.\n\nWeaknesses/Main Points\n- Several proof techniques are adapted from prior works. It would be clearer to explicitly acknowledge/summarize in the main paper (Section 4.1) if any of the important proof technique are adapted from prior works, instead of only mentioning in the Appendix.\n- Some of the writing may be confusing. Is Rhino a type of structural equation model (SEM), or a structural identification method? In Section 3.1 and Theorem 1, it seems that Rhino is defined as a certain type of structural equation model Eq. (5). However, in Section 3.2 and many other parts, it seems that Rhino is a structure identification method that refers to Eq. (10). I think it would be clearer to make a distinction between SEM/identifiability and identification. E.g., we often use LINGAM to refer to the specific structural equation model, and ICA-LINGAM/Direct-LINGAM to refer to the specific method.\n- Section 3.2 mentions that augmented Lagrangian training is used to make sure Rhino only produces DAGs. Does Theorem 2 require the augmented Lagrangian training or the DAG assumption of Rhino's solution?\n    - If yes:\n        - The authors should make this clear in the theorem statement. The current statement seems like any value of $\\rho$ and $\\alpha$ suffices, and the solution (or search space) need not be DAG.\n        - Is there any theoretical guarantee that Rhino only produces DAGs with such training? What assumptions are needed? Such assumptions should also be clearly stated in the theorem statement. \n    - If no:\n        - Why do we still need such DAG penalty and augmented Lagrangian training?\n- What values of $\\lambda$ in Eq. (10) is needed for Theorem 2 to hold? For infinite data, do we still need sparseness prior? (In many statistics/ML settings, sparseness is often used for finite samples.)\n- Theorem 2: What does it specifically mean to be no model misspeficiation? Could the authors provide a precise technical definition?\n- The paper should make it very clear that Theorem 2 is based on global solution of Eq. (10). In practice, it is impossible and we often end up with a local solution.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: notable-top-25%