GDBN: a Graph Neural Network Approach to Dynamic Bayesian Network
Identifying causal relations among multi-variate time series is one of the
most important elements towards understanding the complex mechanisms underlying
the dynamic system. It provides critical tools for forecasting, simulations and
interventions in science and business analytics. In this paper, we proposed a
graph neural network approach with score-based method aiming at learning a
sparse DAG that captures the causal dependencies in a discretized time temporal
graph. We demonstrate methods with graph neural network significantly
outperformed other state-of-the-art methods with dynamic bayesian networking
inference. In addition, from the experiments, the structural causal model can
be more accurate than a linear SCM discovered by the methods such as Notears.
Yang Sun, Yifan Xie
# GDBN: a Graph Neural Network Approach to Dynamic Bayesian Network

###### Abstract

Identifying causal relations among multi-variate time series is one of the most important elements towards understanding the complex mechanisms underlying the dynamic system. It provides critical tools for forecasting, simulations and interventions in science and business analytics. In this paper, we proposed a graph neural network approach with score-based method aiming at learning a sparse DAG that captures the causal dependencies in a discretized time temporal graph. We demonstrate methods with graph neural network significantly outperformed other state-of-the-art methods with dynamic bayesian networking inference. In addition, from the experiments, the structural causal model can be more accurate than a linear SCM discovered by the methods such as Notears.

## 1 Introduction

Multi-variate time series data has been carefully studied for many years. Sensors have been recording large amount of collected time series data including traffics, electricity consumption and so on. Interestingly, these time series variables are subtly internally interlinked. These links with causal dependencies carry critical information for decision makers.

Causality has been formalized as a theory for causal and counterfactual inference [19] from a structural model. Once knowledge about the causality among time series variables is not a priori, we have to discover these relational structure. Moreover, the dynamics described from a structural equation model can articulate the causal dependencies via a directed acyclic graph (DAG). Causal models with explicitly solved DAG enables computation of posterior distribution with treatment[26, 25] and provides us with actionable insights. A careful treatment of causal relation will also avoid misleading judgement from misidentifying spurious correlations as illustrated from Figure 1.

Recently, research done on time series data mainly focuses on the summary graph [15][21] for multi-variate time series. We are more interested with latent structure that governs the multi-variate time series data. The other closedly related topic is graph neural network. Deep neural network has been known as universal function approximator for data distribution and a powerful tool for representation learning. Graph neural network is known as the specific class of deep neural network that is suitable for graph-structured data. With a clear definition of interactions (with directions) among time series data with lags, we recast our problem of Bayesian Dynamics Networks as a link prediction problem via Graph Neural Network (GNN) [30]. The adjacency matrix of the DAG is provided as an input to the graph neural network which will be the objective of our inference problem. One of the most closely related work is [18] based on the linear version of Notears. However the generation process of real world time series data is highly complex and it is difficult to describe them as the usual linear relations.

The main contributions of ours can be summarized in the following.

* By understanding the intrinsic property of the time series data, we define the causal temporal graph and its associated adjacency matrix. We introduce the fundamental elements used in our causal discovery framework.
* We synthesize a wide range of linear dataset and nonlinear dataset from a generation process of VAR(p) with stationary property for causal discovery.
* We further designed a graph-based neural network with a neural message passing scheme that caters for the intrinsic property of multi-variate time series. The graph structure is optimized based on a variational autoencoder structure which generalizes better on understanding the distribution of time series data. We formulate
the optimization objective with variational lower bound and sparsity constraint to solve for causal temporal graph and the SEM.
* We carried out experiments over our synthesized data of different settings and datasets published on [14]. Experiments show that our approach significantly outperformed methods from linear method as Notears [32] and conditional independence framework [24].

## 2 Related Work

The causal inference problem and causal feature selection for time series data has been an issue under different contexts, including econometric, biology, physics and etc. Granger Causality [8] based on prediction method has been proposed as a corner stone methodology for time series causal analysis. It used a variation of the VAR(p) model and the lagging causal effect is modeled with different matrices with index p. Despite its advantages, it suffers from the often lack of the causal sufficiency condition in real scenario. Lower sampling rate on the time dimension could also lead to unidentifiability of the Granger Causality, although some exceptional cases are studied [7]. An alternative is the PC algorithm that tests the conditional independence along with the autocorrelation [22][24]. Despite its capacity of accurate causal discovery, this type of methods usually scale with number of variables and \(\tau_{max}\) in polynomial time. For example [24] shows the time complexity is \(N^{3}\tau_{max}^{2}\) where \(N\) is the number of variables while the complexity of our inference method is of \(N^{2}\tau_{max}\). [17] gives the necessary and sufficient condition to select time series random variable with conditional dependencies based on certain assumptions for the time series data. [2] introduces spectral method for identifying causal effect with robustness to down-sampling.

Score-based methods formulate the causal graph inference problem with an objective function which regularizes the adjacency matrix to satisfy the constraint with respect to being a DAG [11]. NO TEARS[32] recasts the combinatoric graph search problem as a continuous optimization problem, while DYNATEARS[18] provides a vector-autoregressive extension to multivariate time series data. More recently, [31] derived a first principle relation between Graph Neural Network and SCM. Our approach is based on [30] where the authors treat the problem of DAG inference as a link prediction problem.

Neural ODE [4] has started a new chapter for understanding the dynamics of multi-variate time series data in continuous time domain [3] had revealed a fundamental situation where ordinary differential equation can be discovered and [1] approximates a vector fields parametrized by a neural network for modelling the continuous dynamical system. [5] further develops the techniques for conditional treatment effect estimation.

## 3 Problem Formulation

### Assumptions

We assume Markov property and faithfulness holds, and that there are no instantaneous effects or hidden common causes. Then identifiabilty is guaranteed [16, 20]. Our causal discovery work is based on the assumption used in [17]. Besides, the causal relations are invariant under a joint time shift of all variables.

We represent the observations of the multi-variate time series \(X\) as \(\left(X_{1}^{t},...,X_{d}^{t}\right)\) with discrete timestamps.

Following Mastakouri et al.[17], we further assume, in the process of data generation,

**Hypothesis 1** (**H1**): _Lag is one for univariate time series: \(x_{i}^{t-\tau}\) influence \(x_{i}^{t}\) if and only if \(\tau=1,\forall i\in\mathbb{N}^{*}\)._

**Hypothesis 2** (**H2**): _Single-lag:There is only one lag for each pair of time series \((x_{i},x_{j}),\forall i,j\in\mathbb{N}^{*}\)._

where \(\tau\) denotes the time lag between two temporal nodes.

### Definitions

**Definition 3.1**: _causal temporal graph: Causal temporal Graph \(G=\{V,E\}\) is an equivalence class up to time translational invariance. It is composed of vertices \(V=\{X^{t},X^{t-1},\cdots,X^{t-\tau_{max}}\}\) and edges \(E=\{\text{Pa}_{X^{t}}\to X^{t}\}\), \(\text{Pa}_{X^{t}}\) denotes the parent nodes of \(X^{t}\) and \(\tau_{max}\) denotes the size of the observation window size. \(\tau_{max}\) is a hyperparameter in the model training [18]._

**Definition 3.2**: _causal temporal adjacency matrix: Causal

Figure 1: An example of time series (a) with the corresponding full time graph (b) and summary graph (c) under our assumptions. The bold nodes and arrays in (b) constitute a causal temporal graph (Definition 3.1), which highlights the connections between variables at a time slice \(t\) with their past. The maximum time lag \(\tau_{\text{max}}\) is 3. \(x_{2}\) influence \(x_{3}\) with lag=1, which is denoted as \(a_{32}^{1}\). There exist spurious associations represented by gray dashed arrows in (c), which could lead to errors of misidentification of causal relations.
Temporal Adjacency Matrix (TAM) \(\mathcal{A}\) is the adjacency matrix defined over the causal temporal graph G.

The ideal situation would be \(\tau_{max}\) exceeds the maximum lag between time series variables. These definitions enables our proposed graph neural network. Note that the temporal causal graph is a sub-component of the full time graph. Since the full time is translational invariant, it can be generated from the causal temporal graph.

### Causal Discovery for Multi-variate Time Series Data

Consider a \(d\)-variate time series \((X^{t})_{t\in\mathbb{Z}}=(x^{t}_{i},\cdots,x^{t}_{d})\), where the variables influence each other in a time-lagged manner. A general structural causal model (SCM) for the stochastic process can be described as

\[x^{t}_{i}=f_{i}(\text{Pa}^{\tau_{max}}_{i},\cdots,\text{Pa}^{1}_{i},z^{t}_{i}) \tag{1}\]

where \(\text{Pa}^{\tau}_{i}\) denotes a set of variables that influence \(x^{t}_{i}\) at \(t-\tau\). \(z^{t}_{i}\) are jointly independent noise terms, \(i=1,\cdots,m,t\in\mathbb{Z},\tau=1,\cdots,\tau_{max}\). \(\tau_{max}\) is the maximal lag of the causal temporal graph.

#### 3.3.1 General SCM for Multi-variate Time Series Data

A special case is that the model is a class of VAR model with additive noise. Given \(m\) variables \(\mathbf{x}\in\mathbb{R}^{d}\). For a time series of length \(T\), we use \(a^{\tau}_{ij}\) to denote the impact of \(\mathbf{x}_{j}\) on \(\mathbf{x}_{i}\) with a time lag of \(\tau\).

\[\mathbf{x}^{t}_{i}=\sum_{j}\sum_{\tau}a^{\tau}_{ij}\mathbf{x}^{t-\tau}_{j}+\mathbf{z}^{t}_ {i} \tag{2}\]

The variable \(\mathbf{x}\) can be a \(d\)-dimensional vector, or reduced to a scalar when \(d=1\). \(\mathbf{z}\) is the random noise of the same shape with \(\mathbf{x}\), typically independent Gaussian exogenous factors. By taking the matrix form,

\[X^{t}=\sum_{\tau}A^{\tau}X^{t-\tau}+Z^{t} \tag{3}\]

where \(X^{t}=[\mathbf{x}^{t}_{1}\mid\cdots\mid\mathbf{x}^{t}_{m}]\in\mathbb{R}^{m\times d}\), \(Z^{t}\in\mathbb{R}^{m\times d}\) is the noise matrix, \(A^{\tau}=(a^{\tau}_{ij})\in\mathbb{R}^{m\times m}\) is introduced to describe the causal relation between variables, and the compact form can be derived as

\[X^{t}=\mathcal{A}X+Z^{t} \tag{4}\]

where \(\mathcal{X}=[X^{t-1\tau}\mid\cdots\mid X^{t-p\tau}]^{\top}\in\mathbb{R}^{pm \times d},\mathcal{A}=[A^{1}\mid\cdots\mid A^{p}]\in\mathbb{R}^{m\times pm}\) corresponds to the TAM.

For a generalized case we have,

\[X^{t}_{j}=f(\text{Pa}(X^{t}_{j}),Z^{t}_{j}) \tag{5}\]

\(f\) is an arbitrary function on parents node \(X^{t}_{j}\) and independent noise. We only formulate the evolution of the time series with a specific time slice from its past observation window. The reason for this simplification comes from the fact that normally the time series is lengthy and sliding window is a reasonable approach for this type of problem.

#### 3.3.2 GNN-based SCM

Graph neural network is known for its phenomenal power in representation learning for graph. The nodes in the time series contains abundant information from the past, and hence we might utilize the message passing scheme of GNN and approximate complex functions for the dynamics. We parametrize the neural network for SCM along with a graph neural networks,

\[X^{t}=f_{\mathcal{A}}(\mathcal{X},\,Z^{t}) \tag{6}\]

where \(f_{\mathcal{A}}\) is a graph neural network and \(\mathcal{A}\) parametrizes the temporal adjacency matrix which encodes the causal graph, where the non-zero elements suggest the existence of causal relation between the corresponding variables.

## 4 Framework of Gdbn

We propose a bayesian deep generative model with a graph neural network to capture the inter connections among time series nodes. The dependency relations can be described with message passing mechanism in a graph neural network [12][15]. Instead of treating the graph edges as hidden variables, we allow straightforward gradient descent through elements of the temporal adjacency matrix.

### Time Series Window Construction

Now we build the sliding window in the time series, which consists of two parts, one for observation followed by the other one for prediction. By setting the observation window size \(s_{o}\) and the prediction window size \(s_{p}\), we assume that the current effect to the future in time series \(\mathbf{x}\) lasts at most \(s_{o}\) steps, and rebuild the following time series of \(s_{p}\) steps is reasonable. If \(s_{o}>p\), we expect all learned \(A^{p+1},\cdots,A^{s_{o}}\) to be zero matrices. Please see Figure 2 for detail.

### Model Architecture

Here we propose a neural network for inference of the causal relation among multi-variate time series data. It has two main components. The first one is a proposed graph neural network with message passing scheme following a strict causal relational path. Nodes and edges form the representation of values of the time series and their time difference, e.g. lags. The second component is variational inference network
that produce a variational objective for the optimization of the temporal causal graph and structural equation end to end, namely, our framework, is not only capable of achieving causal discovery in a computationally friendly manner, but also provision a tool for prediction purpose of the time series data with strong explainability.

#### 4.2.1 Proposed Graph Neural Network

Graph Neural Network (GNN) is a powerful machinery to introduce relational inductive bias. The neural networks learns the node representation by local message passing mechanism via the graph structure. All variants of GNN [13][9][27][29] have been shown effective on graph structure data related tasks such as graph classification, node classification, link prediction and etc.. Following [6], the local message passing mechanism between vertices (v) and edges (e) can be summarized in the following equations:

\[v\to e: h^{l}_{(i,j)}=f^{l}_{e}(\{h^{l}_{i},h^{l}_{j},x_{(i,j)}\}) \tag{7}\] \[e\to v: h^{l+1}_{j}=f^{l}_{v}(\{\sum_{i\in N_{j}}h^{l}_{(i,j)},x_{j}\}) \tag{8}\]

\(i\) and \(j\) denote the node in the graph. \(x\) represents the raw feature, while \(h^{l}\) represents the embedding in layer \(l\). For example, \(x_{(i,j)}\) is the raw feature of the edge \(i\to j\) and \(h^{l}_{j}\) is the \(l\)th embedding layer of the node \(j\). \(\mathcal{N}_{j}\) represents the set of neighboring nodes for node \(j\). The square bracket denotes the concatenation of the features. The skip connection[10] with raw feature of \(x_{(i,j)}\) in (7) or \(x_{j}\) in (8) is optional. \(f_{e}\) and \(f_{v}\) can be modeled with neural networks. [13] and [28], including the GNN we had proposed can be treated as special case of (8).

_Proposed Graph Neural Network_ (\(\mathcal{G}\))

* Node Embedding (NE): NE maps the node feature to a new embedding space.
* Edge Neural Network (ENN): edge feature can be attained with a edge neural network of a non-linearized transformation. For example, we can choose to encode the lag \(\tau\) as edge feature. We aggregate the node embedding and edge embedding by concatenation to update the embedding of edge, followed by an element-wise multiplication operation with the temporal adjacency matrix \(\mathcal{A}\).
* Node Neural Network (NNN): NNN maps the aggregated edge features back to the node again.

Notice that the propagation of the messages follows a strict topological order which is faithful to the temporal order. Considering the causal effect from the node \(j\) at \(t^{\prime}\) to node \(i\) at \(t\), \(\tau=t-t^{\prime}=1,\cdots,p\), \(\mathcal{G}\) can be formulated as

\[\text{NE}: h_{it}=f_{emb}(x_{it})\] \[\text{ENN}: h^{l}_{(i,j\prime t^{\prime})}=A^{T}_{ij}f^{l}_{e}(\{h^{l}_{it},f_{\tau}(\tau)\})\] \[\text{NNN}: h^{l+1}_{it}=f^{l}_{v}(\sum_{j,t^{\prime}}h^{l}_{(it,j\prime t^{ \prime})})\]

\(f_{emb},f_{\tau},f_{e},f_{v}\) are parameterized by the neural networks, which take the form of Multi-Layer Perceptron (MLP) in this work. See Algorithm 1 for the module of computing GDBN.

#### 4.2.2 Variational Autoencoder

For each step of reconstruction of the time window as in Figure 2, the noise can be generalized from linear Granger model in (6) with our proposed GNN \(\mathcal{Y}\). We assume they yield independent multi-variate Gaussian distribution. At the same time, prediction window is computed with a decoding

Figure 2: Left: Construction of the Causal Temporal Graph with maximal time length of \(S_{o}\). Each slice of the tensor \(A\) represents the connection of \(i\)th node and \(j\)th node with a certain lag. Right: Reconstruction of the time window. Given the observation window of length \(s_{o}\), at every step we predict next value of the time series based on recent past \(s_{o}\) values. After \(s_{p}\) steps, the reconstruction is completed for the current time window. The dotted squares are rolling observation windows from left to right which denotes the direction for recurrent prediction.
procedure faithful to the causal graph. The encoder and decoder are

\[[\mu_{Z^{t}}|\log\sigma_{Z^{t}}] =F_{1}\big{(}F_{2}(X^{t})-\mathscr{C}(X)\big{)} \tag{9}\] \[[\mu_{X^{t}}|\log\sigma_{X^{t}}] =F_{3}(\mathscr{C}(X)+F_{4}(Z^{t})) \tag{10}\]

respectively, where \(Z^{t}\) is the hidden variable in the VAE model, \(F_{i},i=1,\cdots,4\) can be either an identity map or an MLP layer that acts upon the feature dimension.

We used a recurrent decoder that can forecast for multi steps based on one time decoder described above while the reconstruction loss and the KL-divergence for multi-steps are jointly trained. During the training, we model \(P_{\mathscr{C}}(X_{t+1}|X^{t}_{t},X_{t-1}..X)\), where \(X^{t}_{t}\) as predicted from (10). The observation and prediction window described in Figure 2 slides forward along with the causal path. We found in the training process, it is a good choice to use only the expectation value \(\mu_{X^{t}_{t}}\) instead of the sampled value for successive recurrent prediction. This has been proved to be effective for the causal learning as shown in experiments section. For the detail of GDBN network, see a summarized Algorithm 2

```
0: Number of observed variables: \(m\)  Size of the observation window: \(s_{o}\)  Size of the prediction window: \(s_{p}\)  Time series of length \(s_{o}+s_{p}\): \(X^{1:s_{o}+s_{p}}_{1:m}\)  EMB: MLP on feature dimension Output:  TAM \(\mathscr{A}\), Reconstructed Series \(X^{s_{o}+1:s_{o}+s_{p}}_{1:m}\) \(Z^{s_{o}+1:s_{o}+s_{p}}_{1:m}\) for\(i\gets 1\) to \(s_{p}\)do  ENCODE:  logits \(\leftarrow\) EMB(\(\text{EMB}(X^{i+s_{o}}_{1:m})\)-\(\text{GENC}(X^{i+s_{o}-1}_{1:m})\)) \(M\), \(\log\Sigma\leftarrow\text{EMB}(logits)\) \(Z^{s_{o}+s_{p}}_{1:m}\leftarrow\text{SAMPLE}(M,\log\Sigma)\)  DECODE: \(X^{i+s_{o}}_{1:m}\leftarrow\text{EMB}(\text{GENC}(X^{i+s_{o}-1}_{1:m})+\text{ EMB}\ (Z^{i+s_{o}}_{1:m}))\)  REPLACE: \(X^{t+s_{o}}_{1:m}\gets X^{i+s_{o}}_{1:m}\) endfor return\(\mathscr{A}\), \(X^{{}^{\prime}_{o}+1:s_{o}+s_{p}}_{1:m}\), \(M\), \(\log\Sigma\)
```

**Algorithm 2** Forward GDBN Network Module (**GDBN**)

#### 4.2.3 Optimization Objective

**Evidence Lower Bound**

Since the exact inference is intractable, the Evidence Lower

Figure 3: Model architecture of GDBN. Time series of length \(s_{o}+1\) is fed into the model and the output is the prediction \(X^{\prime}_{p}\). At time \(t\), \(X_{p},X_{o}\) are \(X^{t},\mathcal{X}\) respectively. The rounded rectangle is a repeated unit for \(s_{p}\) times. \(h_{\tau}\) is the embedding of raw edge feature \(\tau\).
Bound (ELBO) is optimized in order to approximate the posterior distribution,

\[\text{ELBO}(\theta,\phi;X)= \mathbb{E}_{q_{\phi}(Z|X)}\left[\log p_{\theta}(X|Z)\right]\] \[-KL(q_{\phi}(Z|X)||p_{\theta}(Z)) \tag{11}\]

where \(p_{\theta}(Z)\) is the prior distribution of the latent variables, \(p_{\theta}(X|Z)\) is the conditional distribution of observed \(X\) given \(Z\) and \(\theta\) are the true parameters. \(q_{\phi}\) is the variational posterior parameterized by \(\phi\).

For the inference network, the variational posterior is a factored Gaussian with mean \(M_{Z}\in\mathbb{R}^{m\times d}\) and variance \(\Sigma_{Z}\in\mathbb{R}^{m\times d}\). For the generative network, the likelihood \(p(X|Z)\) is a factored Gaussian with mean \(M_{X}\in\mathbb{R}^{m\times d}\) and variance \(\Sigma_{X}\in\mathbb{R}^{m\times d}\). Then the complexity loss of the latent representation is derived as

\[-KL(q_{\phi}(Z|X)||p_{\theta}(Z))\] \[= \frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{d}\left(1+\log(\Sigma_{Z})_{ ij}-(M_{Z})_{ij}^{2}-(\Sigma_{Z})_{ij}\right) \tag{12}\]

and the reconstruction loss can be computed with Monte Carlo approximation:

\[\mathbb{E}_{q_{\phi}(Z)}\left[\log p_{\theta}(X|Z)\right]\approx\] \[-\frac{1}{L}\sum_{l=1}^{L}\sum_{i=1}^{m}\sum_{j=1}^{d}\left( \frac{(X_{ij}-(M_{x}^{l})_{ij})^{2}}{2(\Sigma_{X}^{l})_{ij}}+\log(\Sigma_{X}^{ l})_{ij}+c\right) \tag{13}\]

where \(L\) denotes the number of Monte Carlo samples, \(c\) is a constant.

### Sparsity and Acyclicity Constraint

We add \(L1\)-norm penalty for the temporal adjacency matrix as a standard approach to introduce sparsity in the optimization objective. Since, we have ignored inter-slice dependencies on the temporal graph, causal relations suffices to have strict topological order and hence no acyclicity constraint needs to be enforced. However, our method can be naturally adapted to inter-slice relations by adding acyclicity constraint in the inter-slice temporal adjacency graph.

### Training

Now that we have all the components for optimization, the training process is composed of the following,

* Given training example \(X=X^{1:s_{0}+s_{P}}\), we first run the encoder and compute \(q(Z|X)\), then we sample \(Z\) and utilize the decoder to compute the prediction values for \(X^{s_{0}+1},\cdots,X^{s_{0}+s_{P}}\).
* We compute the ELBO objective with estimated KL loss and reconstruction error in (12) and (13). We add the regularized term for sparsity constraint and the optimization objective is estimated, \[\min_{\Phi}\min_{\mathcal{A}\in\mathbb{R}^{m\times p}}\frac{1}{n}\sum_{n}I_{ \mathcal{A},\Phi}(X)+\lambda||\mathcal{A}||_{1}\] (14) where \(n\) is the number of samples, \(\lambda\) is the coefficient of sparsity and \(l\) is the opposite of ELBO defined in (11). \(\Phi\) denotes the parameters of the neural network.
* The code is implemented in Pytorch. We perform gradient descents through the neural network parameters and \(\mathcal{A}\). When the optimization process finishes, a hard threshold \(\omega\) on \(\mathcal{A}\) is set. Any weight smaller than \(\omega\) in absolute value is treated a vanishing matrix element, which provably reduces the number of false discoveries[34].

Training process of the network is given in 3,

```
0: Number of observed variables: \(m\)  Size of the observation window: \(s_{o}\)  Size of the prediction window: \(s_{p}\)  Addency Matrix \(\mathcal{A}\) repeat  Sample time series of length \(s_{o}+s_{p}\): \(X^{1:s_{0}+s_{p}}_{1:m}\) from the whole time series \(X\)  Compute \(\mathcal{A}\), \(X^{s_{o}+1:s_{o}+s_{p}}_{1:m}\), \(M\), \(\log\Sigma=\text{GDBN}(X^{1:s_{0}+s_{p}}_{1:m}\), \(\mathcal{A})\)  Compute ELBO(\(\theta,\phi\).\(\mathcal{A}\)) with \(M\), \(\log\Sigma\), \(X^{1:s_{0}+s_{p}}_{1:m}\) from (12) and (13).  Compute the sparsity term \(||\mathcal{A}||_{1}\)  Compute the gradients \(g\leftarrow\nabla ELBO+\lambda||\mathcal{A}||_{1}\) \(\theta\), \(\phi\), \(\mathcal{A}\)\(\leftarrow\) Update parameters using gradients until Convergence of \(\theta,\phi,\mathcal{A}\) return\(\mathcal{A}\)
```

**Algorithm 3** Causal Inference for causal temporal graph

## 5 Experiments

### Datasets, Baselines and Performance Metrics

Since we are evaluating the performance of our inference network for causal discovery. We focus on two datasets that has validation ground truth of causal links.

#### 5.1.1 Dataset

Our own synthetic DatasetWe build our own synthetic dataset following the two steps:

_Temporal Adjacency Matrix Generation_
The temporal adjacency matrix \(\mathcal{A}=[A^{1}\mid\cdots\mid A^{p}]\in\mathbb{R}^{m\times p\mid}\) suffices to generate the full time series data given initial condition. We first assign the positions of nonzero elements. In order for the causal temporal graph to meet the assumptions mentioned in section 3.1, \(\mathcal{A}\) has to satisfy Hypothesis 1 and Hypothesis 2. Hence among all possible positions for nonzero elements, several are assigned nonzero values in a random way. We set the ratio of zero elements in the temporal adjacency matrix \(r\in[0,1]\), to control the sparsity of the causal graph. Besides, about half of the final nonzero positions are assigned with negative values. Then we sample the weights of \(\mathcal{A}\) uniformly with absolute values in the range \([0.7,0.95]\), as suggested by Mastakouri et al.[17].

#### Linear Dataset and Nonlinear Dataset Generation

Now we simulate the time series with TAM \(\mathcal{A}\) based on (4). For the linear case, we set a time series of the \(m\)-variable vector \(\mathbf{x}\), which corresponds to a \(p\)th order Vector Autoregressive (VAR) model:

\[\mathbf{x}^{t}=A_{1}\mathbf{x}^{t-1}+\cdots+A_{p}\mathbf{x}^{t-p}+\mathbf{\epsilon}^{t} \tag{15}\]

where \(A\) is an \(m\)th-order parameters matrix, \(\mathbf{\epsilon}^{t}\) is an independent noise, \(\mathbf{\epsilon}\sim\mathcal{N}(\mathbf{0},\sigma I)\). To guarantee the stability, all roots of \(\det(I-A_{1}z-\cdots-A_{p}z^{p})\) are limited to lying outside the unit circle. See (1) for details. Notice that the variables are scalar-valued for simplicity but can be easily generalized to the vector-valued ones.

For the nonlinear case, we consider two ways:

\[X^{t}=f(\mathcal{A}\mathcal{X})+Z \tag{16}\] \[X^{t}=\mathcal{A}g(\mathcal{X})+Z \tag{17}\]

where the nonlinear function \(f\), \(g\) perform elementwise. (16) is a widely adopted form considering nonlinearity. We also notice that (17) will not change the causal graph structurally, as Yu et al.[30] reason its property with Taylor expansion. Here we consider \(f(\cdot)=g(\cdot)=\sin(\cdot)\) as the nonlinear form for convenience.

With \(X^{1},\cdots,X^{p}\) initialized independently in a random way, the time series of length \(T\) can be computed by generating \(X^{t}\) at each time step \(t\) repeatedly for \(t=p,\cdots,T\). We construct time series samples of length \(s\) with sliding windows.

Other Benchmark DatasetTo further test our proposed model, we also adopt the dataset generated by [14], which provides a system to generate time series data with flexible configurations for causal discovery. Notice the data is different from ours in that it allows for a varying likelihood of nonlinear functions given a TAM. We set the maximum lag \(p=2\), the number of variables \(d=5\) and generate TAM under out assumptions, then simulate time series with gaussian noise. Here we adopt the concept of complexity, which is set with different default levels in [14]. For level 0, the time series is generated only in a linear way. Increasing level stands for increasing percentage and variety of non-linear dependencies, including piecewise linear, monotonic and periodic functions.

#### 5.1.2 Baselines

To validate the proposed model GDBN, we compare its performance with other popular approaches for bayesian networks inference. Our baseline is PCMCI [24], PCMCI+[23] and DYNAPARS [18] without intra-slice acyclicity constraint according to our assumptions about the time series data. PCMCI is strong baseline based on conditional independence test. It is compatible with our assumption of Causal stationarity, no contemporaneous causal links and no hidden variables. DYNAPARS is similar to ours that the optimization objective is the same, i.e., the adjacency matrix of DAG.

#### 5.1.3 Performance Metrics

The inference results are evaluated on some common graph metrics: 1) False discovery rate (FDR), 2) True positive rate (TPR). 3) F1 computed as the harmonic mean of (1-FDR) and TPR, 4) Structural Hamming Distance (SHD), which is the number of required changes to the graph to match the ground truth and computed as the sum of missing edges and extra edges. In view of the time order, the direction of all edges are already determined according to our assumptions. True (T) is the set of true edges, False (F) is the set of non-edges in the ground truth graph. Positive (P) is the set of estimated edges, where True Positive (TP) is an edge that is in the true graph, while False Positive (FP) is not.

* FDR = # FP / # P
* TPR = # TP / # T
* FDR) \(\times\) TPR / (1
- FDR + TPR)
* 2 \(\times\) # TP

Considering these metrics relies on the threshold set manually(4.3), we also introduce AUROC.

### Results and Hyperparameters

We set the maximum lag \(p=5\), the number of variables \(d=10\) to simulate the nonlinear time series of length \(T=600\). We set the coefficient for L-1 norm \(\lambda=0.01\) for GDBN and DYNAPARS, regularization parameter \(\alpha_{\text{PC}}=0.1\) for PCMCI, and threshold are tuned for each configuration to achieve a desirable F1 score. For GDBN, there are more hyperparameters. we set the observation window size \(s_{o}=10\), the prediction window size \(s_{p}=3\), the feature dimension
of the latent variables \(d_{Z}=8\), and the number of hidden units for MLP \(h=32\). The performances are shown in Table 1 for the nonlinear data (16), and Table 2 for the nonlinear data (17). Note that we run 5 repeated random experiments and take an average to report the performances of GDBN and DYNAGEARS, while for PCMCI, which is significantly time-consuming, we only run once.

GDBN outperforms other methods, especially in the detection of complicated nonlinear causal relations. And GDBN has a particular advantage when applied to relatively large datasets with high accuracy and efficiency.

For nonlinear data, GDBN gives an outstanding performance, especially decreasing the number of false discoveries within range of the maximum lag \(p\) when compared with DYNAGEARS, shown in Figure 4 and Figure 5. This shows the effectiveness of GDBN on nonlinear data and DYNAGEARS is generally not sufficient for nonlinear scenario.

We show the GDBN performance on dataset in [14] for different level of complexity in Table 3. In case of large number of edges, GDBN reaches comparable performance with saving a signifcant ammount of time shown in Table 4. The experiment is carried on Intel(R) Core(TM) i7-8557U CPU @ 1.70GHz. The time recorded for GDBN is for 100 epochs so that GDBN reaches stable AUC performance.

## 6 Conclusion and Discussion

In this work, we introduce GDBN, a method that can infer causal dependencies for a full time graph in a discrete multi-variate temporal graph. As a preliminary, we define a fundamental structure of temporal graph and the relevant adjacency graph matrix. We formalize the inference network to causal discovery and structural causal model, and designed the way to train the network jointly from end to end. We conducted experiments and benchmark our methods on two main different types of time series data with the known ground truth. Our approach reaches significantly higher precision than popular benchmarks such as Dynotears[18] and PCMCI [24]. Our method can be utilized widely in the real world problems, as understanding causal relational knowledge is becoming a pressing need for decision makers.

In despite of the success of our method, the graph neural networks of higher expressive power is expected. The question remains regarding what type of graph neural network suites the causal discovery is suitable for general SCM defined in [33]. The other important problem is what if our time series data is downsampled or irregularly sampled from the sufficient time series data and how our framework can be adjusted to adapt to such cases.
