11institutetext: Department of Information Systems, University of Maryland,
Baltimore County, Baltimore, MD, United States
11email: {omarfaruque, sali9<EMAIL_ADDRESS>
https://bdal.umbc.edu 22institutetext: Climate Science Section, Lawrence
Livermore National Laboratory, Livermore, CA, United States
22email<EMAIL_ADDRESS>
# TS-CausalNN: Learning Temporal Causal Relations from Non-linear Non-
stationary Time Series Data

Omar Faruque 11 0009-0006-8650-4366 Sahara Ali 11 Xue Zheng 22 Jianwu Wang 11
0000-0002-9933-1170

###### Abstract

The growing availability and importance of time series data across various
domains, including environmental science, epidemiology, and economics, has led
to an increasing need for time-series causal discovery methods that can
identify the intricate relationships in the non-stationary, non-linear, and
often noisy real world data. However, the majority of current time series
causal discovery methods assume stationarity and linear relations in data,
making them infeasible for the task. Further, the recent deep learning-based
methods rely on the traditional causal structure learning approaches making
them computationally expensive. In this paper, we propose a Time-Series Causal
Neural Network (TS-CausalNN) - a deep learning technique to discover
contemporaneous and lagged causal relations simultaneously. Our proposed
architecture comprises (i) convolutional blocks comprising parallel custom
causal layers, (ii) acyclicity constraint, and (iii) optimization techniques
using the augmented Lagrangian approach. In addition to the simple parallel
design, an advantage of the proposed model is that it naturally handles the
non-stationarity and non-linearity of the data. Through experiments on
multiple synthetic and real world datasets, we demonstrate the empirical
proficiency of our proposed approach as compared to several state-of-the-art
methods. The inferred graphs for the real world dataset are in good agreement
with the domain understanding.

###### Keywords:

Causal Discovery Time-Series Data Non-Stationarity Non-Linearity, Neural
Network

##  1 Introduction

Time series data generated by different natural systems, such as climate and
environment, can possess different features like non-linearity, non-
stationarity, presence of different noise categories, and autocorrelation,
which are mostly multivariate [34]. These intricate characteristics of time
series data ingraft complex challenges in understanding the dependencies of
different components of these natural systems. One popular approach to
simplify the understanding of large multivariate time series datasets is to
graphically represent the data generation model using directed acyclic graphs
(DAGs), which is a very convenient way to express complex systems in a highly
interpretable manner and also provide causal insights into the underlying
processes of the system [29]. DAG representation of a system plays a vital
role in decision-making and future condition prediction in different
applications like causal inference [28, 43], neuroscience [31], medicine [16],
economics [2], finance [36], and machine learning [22]. Learning DAG through
causal discovery from observational time series data is very challenging when
controlled experiments with different population sub-groups are not possible
or unethical [43, 30].

Several state-of-the-art methods have been developed for causal discovery from
temporal data based on constraint-based and score-based methodologies.
Constraint-based causal discovery methods [35, 32, 11, 8, 18] learn the
conditional independencies from the dataset using different tests and
construct the DAG to reflect this independence knowledge. One major limitation
of the conditional independence test is that it requires a large number of
observation samples to generate reliable test scores [38]. Score-based causal
discovery methods use a score function to quantify the predicted causal graph
based on the adjacency matrix and try to optimize the score function gradually
by enforcing the acyclicity constraint of the graph. The large search domain
of the score function based on the adjacency matrix combination makes the
optimization process very challenging and sometimes requires additional
knowledge of the DAG. The combinatorial characteristic of the score
optimization transformed into a continuous optimization problem by Zheng et
al. [49] formulating an equivalent acyclicity constraint using the trace
exponential of the predicted adjacency matrix. Through this continuous
formulation of the acyclicity constraint now causal graphs can be optimized
using several widely used gradient-based optimizers. This drastically changes
the domain of multivariate time series causal discovery using neural network-
based approaches and several methods have been proposed in recent years.

Time series causal discovery leveraging the power of neural networks has
gained much attention as an active area of research in several fields.
However, the majority of existing methodologies are designed under the
assumption of stationarity, a condition often violated in real world scenarios
where dynamic systems evolve periodically [15, 33]. Some existing approaches
also require prior knowledge of the domain model like linearity, noise
distribution, and parametric information. Motivated by the tremendous success
of neural networks [13], in this paper, we propose Time-Series Causal Neural
Network (TS-CausalNN) - a novel causal discovery method from temporal data
using a convolution neural network. Where the causal relationship between each
child to its temporal parents is learned by using a custom 2D convolution
layer. Our proposed model is capable of capturing the causal structure from
multivariate temporal data without any noise and data distribution
assumptions. Our method is designed to be applicable across diverse domains
that scale gracefully for multivariate time series data. The contributions of
this paper are three-fold. First, we propose a 2D convolutional neural network
layer to learn the causal relationships from the multivariate time series
dataset. By utilizing the power of the CNN, our proposed model can handle both
the stationary and non-stationary data for linear and non-linear structural
models with the presence of different noise distributions. The proposed
parallel network architecture has not been used earlier. Second, the
simplified optimization routine of the proposed model can identify lagged and
contemporaneous causal links simultaneously. The integration of the acyclicity
constraint and sparsity penalty into the optimization process helps to learn
better causal graphs. Finally, we conduct extensive evaluations of the
proposed model with state-of-the-art methods using synthetic and real world
datasets. The proposed model achieves better evaluation scores for generated
causal discovery graphs for most of the cases, making TS-CausalNN a strong
contender for time-series causal discovery.

The remainder of this paper is organized as follows. Section 2 describes the
preliminaries of temporal casual discovery and provides an overview of
existing literature highlighting key methods and their limitations. In section
3 we describe the proposed method in detail, emphasizing its key components
and innovations. Section 4 presents the datasets, preprocessing steps, and
experimental configurations used to evaluate the proposed method. The
quantitative and qualitative comparison of the performance of our method with
baseline approaches and real world applications are discussed in section 5.
Section 6 summarizes the main contributions of the paper, acknowledges
limitations, and suggests directions for future work.

##  2 Preliminaries

Temporal causal discovery learns directed acyclic graphs (DAG) from time
series data. Each directed edge of the DAG represents the influence of the
cause variable on the target variable. As the order of data is strictly
maintained in the time series case, influences from cause variables to the
target variable can only come from the same and previous timesteps, also
called temporal precedence (the right side of Figure 1). This feature makes
causal graph learning from time series data more challenging compared to
independent and identically distributed (IID) data.

Let’s consider a multivariate time series dataset
$X=\\{x^{1},x^{2},x^{3},…,x^{n}\\}$ consisting of $n$ variables, and each
variable is measured for $T$ timesteps. Each variable $x^{i}$ at a specific
time point $t$ can be caused by other variables in the same time point ($t$)
and all variables from previous timesteps ($0$ to $t-1$). The effect from
previous timesteps, also called lagged effects, can propagate from infinite
earlier time points, but for DAG learning purposes we will consider a maximum
time lag, $l_{max}$. So the set of possible cause variables of each time
series $x^{i}$ at time $t$ is
$PA_{x^{i}}\in[\\{X_{(t-l_{max})},X_{(t-l_{max}-1)},...,X_{(t-1)},X_{t}\\}-x^{i}]$.
The goal is to learn a causal graph $G(V,E)$ from the provided dataset such
that its vertices resemble time-lagged and current time variables in the time
series and its directed edges express causal links from parents to target
variables. So the vertices and edges are denoted as
$V=\\{X_{(t-l_{max})},X_{(t-l_{max}-1)},...,X_{(t-1)},X_{t}\\}$,
$E=\\{(V_{i},V_{j}):V_{i},V_{j}\in\\{X_{(t-l_{max})},X_{(t-l_{max}-1)},...,X_{(t-1)},X_{t}\\}\\}$,
respectively. Let the weighted adjacency matrix of causal graph G be denoted
by $W\in R^{(n\times(l_{max}+1))\times n}$, which contains both the time-
lagged and instantaneous part of the causal links. The structural equation
model of the time series can be defined as:

$\begin{split}X_{t}=f_{W}(X_{(t-l_{max})},X_{(t-l_{max}-1)},...,X_{(t-1)},X_{t})+e_{t}\end{split}$
(1)

, where the noise term $e_{t}$ can be of any type, independent of cause and
effect, and the structural function $f_{W}()$ can be any linear or non-linear
data generation process.

Figure 1: Temporal causal graph learned (right) from multivariate time series
data (left) with causal links from the same and previous timesteps. Each node
in the graph represents one variable at a specific timestep ( $t$ is present
time and $t-i$ is a previous timestep). A directed edge denotes a causal
relationship between cause and effect.

In the time-lagged part of the adjacency matrix $W$, the nodes from previous
$l_{max}$ timesteps of $t$ will always be the source/cause node of the lagged
causal links. So, we do not need to explicitly focus on the acyclicity of the
time-lagged part of $W$. However, for the contemporary part of the $W$ at $t$,
each node can serve as both the source and the target of the causal links
where we have to maintain the acyclicity of the DAG.

Several methods can be employed to recover the adjacency matrix $W$ of the
causal graph from time series data, like search-based, constraint-based,
score-based, graph-based methods, etc. To formulate this combinatorial search
problem into a continuous optimization nature Zheng et al. [49] proposed an
algebraic representation of the acyclicity constraint of DAG with the
adjacency matrix exponential. This new representation of acyclicity opens a
window to learning the DAG using generic continuous optimization routines like
neural networks by minimizing the score. However, simultaneously learning the
lagged and contemporaneous part of the adjacency matrix is very challenging
for complex datasets. Since any variable might be the cause of another effect
variable, cycles can occur in the contemporaneous part of the adjacency
matrix. Consequently, these two parts of the adjacency matrix require a
different set of optimization criteria.

##  3 Proposed Methodology

We can consider the causal graph generating task as an unsupervised learning
process of the adjacency matrix $W$ given the multivariate time series data
$X=\\{x^{1},x^{2},x^{3},…,x^{n}\\}$ of $T$ observations. To learn the directed
adjacency matrix $W$ of the temporal causal graph $G$, we propose a neural
network-based unsupervised model. The proposed time series causal neural
network model will learn the instantaneous $(X_{t}\to X_{t})$ and time-lagged
$(\\{X_{(t-l_{max})},X_{(t-l_{max}-1)},...,X_{(t-1)}\\}\to X_{t})$ causal
links of $W$ for maximum time lag ($l_{max}>0$) in the same gradient
propagation using a single neural network. To perform the temporal causal
learning, we propose a causal convolution 2D layer.

###  3.1 Time Series Causal Neural Network

As illustrated in Figure 2, our proposed TS-CausalNN model consists of two
separate blocks of 2-dimensional convolution layers, which helps its
simplicity and efficiency. The first 2D convolution layer of the model
transforms input data into a latent representation, which helps the model to
learn any non-linear relationship present in the data generation process. The
Causal Conv2D block contains $n$ number of parallel custom causal layers to
learn the time-lagged and instantaneous causal relationships of each input
variable to its parent variables. Following a similar analogy used by Zheng et
al. [50], DAG-GNN [48] and NOTEARS-MLP [49], we will learn the causal links
from the parameters of the Conv2D block. Here the links learned by this causal
block are always unidirectional and the parallel blocks help to learn the
causal links for each input variable independently. Finally, the results of
each parallel causal layer are aggregated to generate the model output, which
will be used in the optimization process of the learned causal graph.

Figure 2: Proposed TS-CausalNN model architecture to learn full temporal
causal graph. The graph can be learned from the parallel Causal Conv2d layers.

Our proposed Causal Convolution 2D (Causal Conv2D) layer takes input in a
similar structure as shown in the full causal graph on the right side of
Figure 1, with lagged data followed by the current time point’s data. Each
Causal Conv2D layer is designed to learn the causal links of an input variable
for example $x^{1}$ from all possible parents
$PA_{x^{1}}\in\\{x^{1}_{(t-l_{max})},x^{1}_{(t-l_{max}-1)},...,x^{1}_{(t-1)},x^{2}_{(t-l_{max})},x^{2}_{(t-l_{max}-1)},...,x^{2}_{(t-1)},x^{2}_{t},$
$...,x^{n}_{(t-l_{max})},x^{n}_{(t-l_{max}-1)},...,x^{n}_{(t-1)},x^{n}_{t}\\}$
of that variable. The variable itself cannot be included in the set of its
parent variables. Let’s assume we have a time series dataset with 4 variables
$X=\\{x^{1},x^{2},x^{3},x^{4}\\}$, and for lagged effects consider the maximum
time lag $l_{max}=4$. So the input data will be a matrix of size $(4\times 5)$
one row for each variable and $l_{max}+1=5$ column for lagged and
contemporaneous variables. To learn the temporal causal graph for these 4
variables, as shown in Figure 3, we have to employ 4 parallel Causal Conv2D
layers one layer for each variable. Each of these layers predicts the
expectation of the target variable at the current timestep $t$ given all
lagged and instantaneous parents (Equation 2). To exclude the target variable
from the parents list, its weight is set to zero.

Figure 3: Proposed custom Causal Conv2D Layers.

$E[x^{i}|PA_{x^{i}}]=f_{W^{x^{i}}}(PA_{x^{i}})$ (2)

Here $f_{W^{x^{i}}}()$ is the function learned at the Causal Conv2D layer
considering $x^{i}$ as the target variable and $W^{x^{i}}$ is the set of
weight parameters of that layer. Motivated by NOTEARS-MLP [49] and DAG-GNN
[48], we derive the adjacency matrix of causal DAG from the weight parameters
of the Causal Conv2d layers. The value of the weight parameter of the layer
for a target variable represents the strength of the causal link from its
parent. The weight parameter $W^{x^{k}}_{ij}=0$ means the target variable
$x^{k}$ is independent of the cause variable $x^{i}$ at timestep $j$ in the
$PA_{x^{1}}$ list. If the $W^{x^{k}}_{ij}>0$ then the variable $x^{k}$ has a
causal edge from the parent variable $x^{i}$ at a specific time lag $j$. After
learning the final weight parameters of the target variable, we apply a
thresholding operation to prune causal links with weak dependency strength,
$W^{x^{k}}_{\omega}=(W^{x^{k}}>\omega)$, where $\omega$ is the threshold
value. Finally, the weight parameters of all variables after thresholding are
concatenated to generate the combined adjacency matrix of the final causal
graph.

###  3.2 Acyclicity Constraint

Maintaining the acyclic property of the adjacency matrix makes it very
difficult to learn causal graphs from time series data. Using a continuous
optimization process like a neural network to learn the causal graph does not
guarantee the acyclicity of the learned adjacency matrix of the DAG. To impose
the acyclicity restriction in the adjacency matrix of the learned causal DAG,
we will use a similar equality constraint $h(W)=0$ as Zheng et al. [49]. The
function $h(W)$ is defined using the trace exponential of the elementwise
product of the adjacency matrix with itself.

$h(W)=tr\>e^{WoW}-n$ (3)

Here $n$ is the number of variables in the dataset. We cannot use the learned
adjacency matrix $W$ of the proposed TS-CausalNN method directly in this
equality function because $W$ contains both the time-lagged
$(t-l_{max},t-l_{max}-1,…,t-1)$ and contemporaneous $(t)$ edges of the causal
graph. The time-lagged edges of the causal graph will always go forward in
time from the previous timesteps to the current timestep $t$, therefore we do
not need to compute acyclicity for this part of the graph. So we have to apply
the acyclicity constraint on the contemporaneous part of the adjacency matrix,
$W^{t}$. Hence the function can be stated below and the equality will be
satisfied if and only if the contemporaneous part of the adjacency matrix
($W^{t}$) is acyclic.

$h(W^{t})=tr\>e^{W^{t}oW^{t}}-n=0$ (4)

###  3.3 Optimization

The goal of the proposed causal discovery method is to learn the adjacency
matrix $W$ of the DAG from the given time-series dataset $X$. The adjacency
matrix $W$ contains the edges from the time-lagged variables together with the
contemporaneous variables to the target variables. So the proposed TS-CausalNN
model has to estimate the target variables using all possible parents and find
the causal influence ($W$) of each parent on the target child, where the
contemporaneous part of $W^{t}$ must satisfy the acyclicity requirement. So
the optimization objective of the model is to minimize the least-square loss
($L$) of the target variables with the acyclicity constraint on $W^{t}$.

$\min_{W}L(W)\;\;\mathrm{subject\;to}\;\;h(W^{t})=0\;\;\mathrm{for\;acyclicity}$
(5) $\mathrm{where}\;\;L(W)=\frac{1}{T}||X-WX||^{2}_{F}$ (6)
$h(W^{t})=tr(e^{W^{t}oW^{t}})-n=0$ (7)

The function $h(W^{t})$ will be equal to $0$ if and only if the corresponding
sub-graph graph of matrix $W^{t}$ does not have any cycle. We cannot directly
integrate the equality-constraint problem of the acyclicity into the
continuous optimization problem. But the equality constraint $h(W^{t})=0$ can
be solved using continuous optimization after converting this into an
unconstraint subproblem [49]. Therefore we will use the augmented Lagrangian
method to solve the equality constraint problem. An additional penalty is
incorporated with the objective function to enforce the sparsity of the
adjacency matrix using the $L1$ norm of $W$.

The least-square loss and the acyclicity constraint (Equations 6 and 7) will
try to increase the weight parameter to minimize the loss. On the contrary,
the L1 norm of W will try to reduce weight values to zero to keep a minimum
number of non-zero entries. By working contrariwise to each other the loss
function will be optimized in an equilibrium way. Hence the final
unconstrained objective function is:

$\min_{W}\left[L(W)+\frac{\rho}{2}|h(W^{t})|^{2}+\alpha
h(W^{t})+\lambda||W||_{1}\right]$ (8)

, where 2nd and 3rd terms are the augmented Lagrangian for the acyclicity
constraint, $\alpha$ is the Lagrange multiplier, $\rho>0$ is the penalty
parameter of the augmented Lagrangian, and $\lambda$ is the sparsity penalty
parameter. Here $\lambda$ is a model hyperparameter and we set it at $0.1$ for
this study. This objective function can be minimized using any state-of-the-
art continuous optimizer. An excellent characteristic of the augmented
Lagrangian approach is its ability to accurately approximate the solution of a
constrained problem using the solution of unconstrained problems by gradually
increasing the penalty parameter $\rho$ but not to infinity. Hence we have to
gradually increase $\rho$ to minimize the unconstrained augmented Lagrangian
and the value of the $\alpha$ also be updated accordingly. The rule for
updating $\rho$ and $\alpha$ based on the value of $h(W^{i})$ is given by:

$\rho^{i+1}=\left\\{\begin{matrix}(1+\beta)\rho^{i},&\textup{if}\;\;h(W^{i})>\gamma
h(W^{i-1})\\\ \rho^{i},&\textup{otherwise}\end{matrix}\right.$ (9)
$\alpha^{i+1}=\alpha^{i}+\rho^{i}h(W^{i})$ (10)

, where $\beta>0$ and $\gamma<1$ are hyperparameters, and we find that
($\beta=0.1,\gamma=0.25$) work better. The overall process of learning causal
DAG is outlined in Algorithm 1.

Input: Multivariate Time Series Data $X=\\{x^{1},x^{2},x^{3},...,x^{n}\\}$

Output: Adjacency Matrix of the causal graph $W$

1 Instantiate and compile the Causal Neural Network model

2for  _each iteration_ do

3 Train the Neural Network model

4 Compare the acyclicity loss with the previous iteration

5 Update the penalty coefficient of the Lagrangian method ($\rho,\alpha$)

6

Return the adjacency matrix $W$

Algorithm 1 TS-CausalNN Algorithm

##  4 Experimental Setup

We mention the dataset and evaluation criterion used for performance
comparison in this section. Our model is developed on TensorFlow Keras and all
experiments are conducted on Google Colab Runtime with CPU for easy
reproducibility. Fixed seed value is used for randomized data to make the
experimental results reproducible.

###  4.1 Synthetic Datasets

To evaluate the performance of our proposed causal discovery method we have
used synthetic datasets. As we know the ground truth causal graph for the
synthetic datasets, we can measure and compare the learned causal graph
easily. We generate a time series dataset consisting of four non-linear
variables using Gaussian white noise $\varepsilon$ (Dataset-1). The
mathematical description of each variable is given in Equations 10 to 13. The
non-linear characteristic is incorporated in the generation of the synthetic
dataset to mimic the dynamic properties of real world natural system data. The
corresponding true causal graph for the time-series data is given in Figure
4a.

$S1_{t}=2\\{cos(\frac{t}{10})+log(|S1_{t-2}-S1_{t-5}|+1)\\}+0.1\varepsilon 1$
(11)
$S2_{t}=12e^{\frac{S1^{2}_{t-1}}{2}}-4e^{\frac{S1^{2}_{t}}{2}}+\varepsilon 2$
(12) $S3_{t}=-10.5e^{\frac{-S1^{2}_{t-1}}{2}}+\varepsilon 3$ (13)
$S4_{t}=-11.5e^{\frac{-S1^{2}_{t-1}}{2}}+13.5e^{\frac{-S3^{2}_{t-1}}{2}}+1.2e^{\frac{-S4^{2}_{t-1}}{2}}-5e^{\frac{-S3^{2}_{t}}{2}}+\varepsilon
4$ (14)

Figure 4: Causal graph of (a) our synthetic datasets and (b) the real world
Turbulence Kinetic Energy (TKE) dataset.

We have also generated another version of the synthetic dataset (Dataset-2)
from the previous causal graph without the cosine part of the data generation
formulas. For this case, we used the Poisson distribution to incorporate
random noise into the dataset. The goal is to find the performance and
robustness of the causal discovery methods for the time series dataset without
having wave-shape-like features with non-Gaussian noise. The formula used to
generate this dataset is available in Appendix 0.A.

Before applying the proposed causal discovery method, we performed some data
preprocessing to ensure the quality and consistency of the input data.
Different variables of the dataset have different scales of values. All time
series data were normalized within the scale of 0 to 1 to mitigate the impact
of scale differences. Then $l_{max}$ previous timesteps data of each current
timestep is concatenated in front of it to make the dataset size
$(T,n,(l_{max}+1))$. The preprocessed data can be applied to the input
convolution 2D layer of the model directly.

###  4.2 Real World Datasets

Two real world Earth/atmospheric science datasets, namely Turbulence Kinetic
Energy (TKE) and Arctic Sea Ice, were used to evaluate our work. These natural
datasets exhibit high variability, non-stationarity, and complex data
interactions. Owing to these characteristics, we evaluate our model
performance on these datasets to assess how well the proposed model performs
on such complex data.

TKE refers to the mean kinetic energy per unit mass associated with eddies in
turbulent flow [17]. It can be measured by the root-mean-square (RMS) velocity
fluctuations (unit: $m^{2}s^{-2}$). The temporal TKE data used in this study
represents the TKE evolution during a typical cumulus-topped boundary layer
day (local time 05:00 – 18:00) over the DOE Atmospheric Radiation Measurement
(ARM) Southern Great Plains Central Facility. This data file is generated from
an idealized numerical simulation using the Weather Research & Forecasting
Model [42] with modifications from the Large-Eddy Simulation (LES) Symbiotic
Simulation and Observation (LASSO) activity, which is developed through the US
Department of Energy’s ARM facility [14, 7]. Besides TKE, the model also
generated the budget terms determining the temporal change in TKE. The major
budget terms include the TKE vertical shear production term
($SH,m^{2}s^{-3}$), the TKE buoyancy production term ($BU$), and the TKE
turbulent and pressure transport term ($TR$). All these budget terms form the
net temporal change term of TKE ($TEND,m^{2}s^{-3}$). If $TEND$ is positive
(negative), TKE will increase (decrease) in the next timestep. Figure 4b
illustrates how these terms relate to one another through a directed graph.

Arctic sea ice is one of the important components of the world’s climate
system that has a great impact on the increase of extreme weather events [5,
40, 47, 24, 25] and the ice is melting rapidly due to various atmospheric
conditions [37, 39]. Huang et al. [19] conducted a causal discovery analysis
to uncover the links between the arctic sea ice and the atmosphere. We use the
same 11 atmospheric variables with the sea ice extent employed in [19]. This
time series data contains monthly averages from 1980 to 2018 over the Arctic
region of 60N. The causal relationship between the atmospheric variables and
the sea ice extent based on the physics and microphysics literature review
discussed in [19] is represented in Appendix 0.G.

###  4.3 Evaluation Metrics

To evaluate the performance of the time series causal discovery methods we
will use three standard evaluation metrics: Structural Hamming Distance (SHD),
F1 Score and False Discovery Rate (FDR). SHD of the directed graph represents
the smallest number of edge corrections required to transform the predicted
causal graph into the true causal graph. The edges can be corrected by adding
new edges, removing existing edges, and changing the direction of an existing
edge in the graph. A lower SHD is expected for better performance of the
causal discovery method. F1 Score calculate the harmonic mean of precision and
recall, providing a balanced measure of performance. The F1 score ranges from
0 to 1 and a higher value means a better prediction of the true graph. FDR
explains the rate of predicted wrong edges from all predicted edges
considering the direction of each edge. This is the inverse measure of the
precision of the prediction.

##  5 Results

In this section, we present the comparative results of the time series causal
discovery between the proposed method and state-of-the-art methods. Several
typical temporal causal discovery challenges have been considered like time-
lagged and contemporaneous causal relations, different noise distributions,
nonlinearity of data, and autocorrelation.

###  5.1 Baselines

To compare the results of the proposed method with baselines we considered
PCMCI+[32], DYNOTEARS [27], NTS-NOTEARS [44], PCMCI [35], NOTEARS-MLP [50],
and DAG-GNN [48]. The first three methods can directly learn causal graphs for
time series data. Though the other three methods were proposed for non-
temporal data we used these methods due to their popularity and widespread
usage in different domains. We transformed the lagged and instantaneous data
into a long sequence so that we could apply the transformed dataset to the
non-temporal methods to find the lagged and current time causal relationships.
For each method, we tuned hyperparameters to get the best evaluation scores.

###  5.2 Quantitative Results

To evaluate the performance of the selected state-of-the-art methods we
compared the predicted causal graph in both the summary and full temporal
graph setting. The qualitative comparison of the baseline methods is reported
in Table 2 and Table 1. The best results are marked in bold text and
underlined values represent the second best score. To get better prediction
results, we have applied each method multiple times on the same dataset and
reported the best result. The comparative analysis of the summary graph (Table
2) shows that the proposed method achieves the best scores for all three
evaluation criteria for both datasets except for the FDR of DAG-GNN method.
For synthetic dataset-1, DAG-GNN method yields a lower score in terms of FDR
whereas its SHD is close to the proposed model. Since, FDR only relies on the
number of predicted edges, the $(0.00^{*})$ FDR score of DAG-GNN cannot be
considered a perfect score as the few edges predicted by the method align with
the true causal graph but it missed predicting the remaining edges (Figure
4a).

Table 1: Comparison of the predicted temporal full causal graph by different methods for synthetic datasets. Method | Dataset-1 | Dataset-2
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 69 | 0.16 | 0.90 | 22 | 0.38 | 0.74
PCMCI+ [32] | 49 | 0.16 | 0.90 | 11 | 0.62 | 0.55
NOTEARS-MLP [50] | 13 | 0.58 | 0.59 | 20 | 0.44 | 0.70
NTS-NOTEARS [44] | 12 | 0.45 | 0.61 | 8 | 0.33 | 0.33
DAG-GNN [48] | 12 | 0.00 | 1.00 | 18 | 0.28 | 0.80
DYNOTEARS [27] | 13 | 0.51 | 0.61 | 18 | 0.25 | 0.80
Proposed | 8 | 0.60 | 0.45 | 7 | 0.46 | 0.25

To compare the results of the predicted full causal graph, we considered both
the directed edges and the time lag of the edge. From Table 1, we can see that
our proposed method got the best F1 score and the second-lowest SHD and FDR
for dataset-1. For dataset-2, the proposed method yields the best results for
all three quality measures. For dataset-1, NTS-NOTEARS generated the target
causal graph with smaller SHD and lower F1 scores compared to the proposed
method. This can be interpreted as the number of edges generated by NTS-
NOTEARS is lower than the number of edges in the true causal graph (low F1
score) and also the 0.5 FDR value means that half of the predicted edges are
wrong. Similarly, zero F1 score of DAG-GNN method means that it was not able
to find any causal links from the true causal graph. From observing the
evaluation results of the full causal graph, we can find that the proposed
model is capable of generating fewer false causal links compared to other
baseline models.

Table 2: Comparison of the predicted summary causal graph by different methods for synthetic datasets. Method | Dataset-1 | Dataset-2
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 10 | 0.54 | 0.62 | 7 | 0.63 | 0.53
PCMCI+ [32] | 10 | 0.54 | 0.62 | 6 | 0.66 | 0.50
NOTEARS-MLP [50] | 6 | 0.25 | 0.50 | 4 | 0.50 | 0
NTS-NOTEARS [44] | 3 | 0.76 | 0.28 | 5 | 0.44 | 0.33
DAG-GNN [48] | 5 | 0.44 | 0.33 | 8 | 0.50 | 0.60
DYNOTEARS [27] | 4 | 0.74 | 0.40 | 7 | 0.36 | 0.60
Proposed | 3 | 0.76 | 0.28 | 4 | 0.60 | 0.25

###  5.3 Qualitative Results

To gain more insights into the discovered causal relationships, we conducted a
qualitative analysis of the predicted causal graphs from synthetic dataset-1.
Figure 5 illustrates the ground truth and predicted causal graphs of the
proposed and DYNOTEARS methods using direct graphs. The nodes in the graph
model the variables in the dataset and each edge means a causal relationship
between the nodes. The self-loop in the causal graph represents that the
previous timestep of the variable has a lagged effect on its present timestep.
From the predicted causal graphs, we can see that the proposed model predicted
the same number of edges to the ground truth graph but failed to predict edge
$S4\to S4$ and estimated a false edge $S3\to S4$. On the other hand, DYNOTEARS
method failed to distinguish the causal effects as mentioned in the ground
truth graph. Instead, DYNOTEARS method predicts bidirectional effects between
the variables $(S1\leftrightarrow S2)$,$(S1\leftrightarrow
S3)$,$(S2\leftrightarrow S3)$, and $(S3\leftrightarrow S4)$. Whereas $S1$ is
the cause variable for the other three variables and there is no causal link
between variable $S2$ and $S3$.

Figure 5: Comparison of ground truth summary causal graph with predicted
graphs from DYNOTEARS and our proposed model for synthetic dataset 1. Figure
6: Visualization of ground truth full temporal causal graph of synthetic
dataset 1 and the predicted graphs from DYNOTEARS and our proposed model. (The
yellow cell means there is a directed edge from the column index to the row
index.)

The full temporal causal graph of the data generation process of synthetic
dataset 1 is visualized in the 1st plot of Figure 6. Here columns represent
the parent variables with the time lag and rows represent the child/target
variables. The adjacency matrix of the true graph is sparse which helps us to
evaluate the performance of time series causal discovery methods for sparse
causal graphs. From the predicted adjacency matrix of the proposed model, we
can see that the causal links $\\{S1(t-2)\to S1(t),S1(t-1)\to S3(t),S1(t-1)\to
S2(t),S3(t-1)\to S4(t)\\}$ are predicted correctly. However, the proposed
model failed to predict the other three true causal links of the dataset. On
the other hand, DYNOTEARS method predicted all edges from the true causal
graph except ($S1(t-5)\to S1(t)$) and predicted more false edges than the
proposed model. The qualitative and quantitative analysis of the predicted
results gives us the intuition that most state-of-the-art models can work
better on predicting dense causal graphs than sparse causal graphs like the
one used in this evaluation process.

###  5.4 Results on Real World Data

We applied the proposed and baseline models to the TKE and Arctic Sea Ice
datasets to generate causal graphs with time lags 5 and 12, referring to the
very fast time evolution of planetary boundary layer turbulence in TKE data
and annual seasonality in sea ice data, respectively. The evaluation results
of the summary graphs generated by baseline models are summarised in Table 3.
For the TKE dataset, the proposed model performed similarly to PCMCI+ method.
Although NOTEARS-MLP has a lower SHD score, its F1 score is lower than the
proposed model. For the Arctic Sea Ice dataset, our proposed method achieves
the best F1 and FDR scores. Although the SHD scores of a few baseline models
were lower than the proposed model, their F1 and FDR values were not as good.

Table 3: Comparison of the predicted summary causal graph by different methods for real world datasets. Method | TKE | Arctic Sea Ice
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 9 | 0.18 | 0.87 | 62 | 0.31 | 0.68
PCMCI+ [32] | 5 | 0.44 | 0.66 | 50 | 0.32 | 0.57
NOTEARS-MLP [50] | 4 | 0.33 | 0.66 | 71 | 0.38 | 0.68
NTS-NOTEARS [44] | 6 | 0.40 | 0.71 | 53 | 0.10 | 0.76
DAG-GNN [48] | 7 | 0.22 | 0.83 | 66 | 0.21 | 0.76
DYNOTEARS [27] | 8 | 0.00 | 1.00 | 65 | 0.21 | 0.75
Proposed | 5 | 0.44 | 0.66 | 54 | 0.49 | 0.56
Figure 7: Visualization of predicted summary causal graphs of the TKE dataset
using PCMCI+, NOTEARS-MLP, and our proposed model.

Due to space limitation, here we only demonstrate the predicted causal summary
graphs by NTS-NOTEARS and the proposed method from the TKE dataset. For causal
graphs generated by baseline methods on real world datasets, refer to
Appendices 0.F and 0.G. We analyze the predicted causal graphs by comparing
the domain graph represented in Figure 4. From the predicted graphs in Figure
7, we can see that the proposed method recovered the correct edges from $BU\to
TEND$, $TEND\to TKE$, and the connection between each pair of nodes is
unidirectional. On the other hand, PCMCI+ method predicted bidirectional edges
for node pairs $BU\to TEND$, and $SH\to TEND$ and missed the important edge
from $TEND$ to $TKE$. On the other hand, NOTEARS-MLP predicted one correct
edge $TEND\to TKE$ though it has a lower SHD score. The time lag 5 is used for
discovering lagged causal relationships for each variable. The predicted full
causal graph can give us an idea about the lagged effect of different
variables in this dataset which is illustrated in Appendix 0.F. By observing
the predicted full and summary graph we see that the proposed model can better
identify the lagged and contemporaneous causal links.

###  5.5 Ablation Study

To verify the effectiveness of the proposed model for temporal causal
discovery, a comparative study of this model with its variant is shown here.
The quantitative results of the generated summary causal graphs of these
models are illustrated in Table 4. For Causal Conv1D, we used 1D variant of
the proposed custom Causal Conv2D layer. To incorporate this Conv1D layer into
the model, we flattened the latent representation of the earlier convolution
layers of the model. The same training and optimization process was utilized
for both models. The ablation study results show that the Custom Conv2D layer
improves the causal graphs learning performance of the proposed model with a
significant margin for each evaluation score. As we are studying multivariate
time series data convolution 2D layers can better learn the inherent features
of the data generation than convolution 1D layers.

Table 4: Performance Comparison between Custom Causal Conv2D and Conv1D models. Dataset | Causal Conv1D | Causal Conv2D
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
Synthetic Datase-1 | 7 | 0.46 | 0.57 | 3 | 0.77 | 0.28
Synthetic Datase-2 | 8 | 0.33 | 0.66 | 4 | 0.60 | 0.25
TKE | 6 | 0.40 | 0.71 | 5 | 0.44 | 0.66
Arctic Sea Ice | 63 | 0.42 | 0.63 | 54 | 0.49 | 0.56

###  5.6 Robustness Analysis

To investigate the performance stability of the proposed model, we considered
different features of datasets including signal-to-noise ratio (SNR), noise
distribution, non-stationary and stationary variable combinations. We
generated three versions of the synthetic dataset-1 (Synthetic-1.1 with lower
SNR, Synthetic-1.2 with medium SNR, and Synthetic-1.3 with higher SNR
comparatively) to evaluate the response of the proposed model to the different
SNRs and applied all baseline models to them. The evaluation results of the
proposed model for these synthetic datasets are summarized in Table 5 and more
results for baseline models are reported in Appendix 0.C. By observing
evaluation results we can see that the performance of the proposed model
increases with the increase of SNRs. The non-stationarity property analysis of
the TKE and Arctic Sea Ice data revealed that all variables of the TKE data
are non-stationary and there is a good mixture of non-stationary and
stationary variables in the Arctic Sea Ice data (Test scores available in
Appendix 0.E). Our proposed model demonstrated good performance for both of
these real world non-stationarity datasets (Section 5.4). The simulated
datasets qualify the impact of various noise distributions on the proposed
model’s performance. For synthetic datasets 1 and 2, respectively, we employed
Gaussian and Poisson noise distributions. The evaluation findings confirm that
our model can function with both types of noise (Sections 5.3 and 5.2) and is
not restricted to only these two types of noise. The comprehensive experiments
with different cases and the evaluation results confirm the robustness of the
proposed model for causal discovery for time series data with diverse
properties.

Table 5: Evaluation of the predicted causal graphs by the proposed method for different SNRs. Dataset | Summary Graph | Full Graph
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
Synthetic-1.1 | 5 | 0.54 | 0.40 | 12 | 0.33 | 0.72
Synthetic-1.2 | 4 | 0.66 | 0.33 | 11 | 0.42 | 0.66
Synthetic-1.3 | 2 | 0.83 | 0.16 | 9 | 0.47 | 0.60

##  6 Related Works

Causal Convolutional Neural Network (CausalCNN). Recently, some studies have
combined the similarities in hierarchical structures within convolutional
neural networks (CNNs) and causal graphs to develop causal networks using CNNs
to serve a variety of purposes. [6] proposed CexCNN - an explanation technique
for CNNs using Pearl’s theory of counterfactual reasoning to improve
predictability. [23] proposed a Knowledge-oriented Convolutional Neural
Network (K-CNN) for causal relation extraction in texts for natural language
processing tasks. [21] use attribute pairs to pre-train CNN models for
determining pair-wise causal direction in data. [41] proposed a novel
Recurrent Convolutional Network (RCN) where they employ recurrence to 3D CNN
for causal processing of videos. Earlier, [26] proposed an attention-based CNN
approach for time series causal discovery, however, their approach had several
limitations including the stationarity assumption in data. It is important to
note that our proposed approach is distinguishable from these methods in both
functionality and usage.

Causal Discovery on Nonlinear and Non-Stationary Data. Statistical causal
discovery methods, like Granger Causality (GC), cannot handle non-linearity or
non-stationarity in data. While some methods extend the traditional causal
discovery method to handle non-linearity [12, 3], some approaches utilize
neural network-for these extensions [45, 1]. Further, some recent research has
made inroads to propose techniques applicable to non-stationary time-series.
[18] proposed Constraint-based causal Discovery from NOnstationary Data (CD-
NOD). Introducing a non-parametric principled framework, they demonstrated the
efficiency of their approach in forecasting non-stationary data. [4] proposed
a probabilistic deep learning approach, State-Dependent Causal Inference
(SDCI), to perform causal discovery in conditionally stationary time-series
data. [46] introduced a functional causal model (FCM) based approach for
causal learning in non-stationary data with general nonlinear relationships.
[9] proposed a constraint-based causal discovery approach for autocorrelated
and non-stationary time series data (CDANs) to identify lagged and
instantaneous causal relationships in autocorrelated and non-stationary time
series data. Most recently, [10] combined Linear Non-Gaussian Acyclic Model
(LiNGAM) and the Just-In-Time (JIT) framework to propose a new causal
discovery method JIT-LiNGAM to identify causal relations in non-linear and
non-stationary data. Recent work also includes deep learning-based methods
which infer the causal graph directly from observational data using
transformers, generative modeling, and adversarial learning techniques [26,
20, 48], overcoming the non-linearity constraint whereas the challenge of
handling non-stationarity still prevails in them.

Limitations of Related Works. While these existing methodologies have
significantly contributed to the field of time series causal discovery,
several challenges persist. The inability to handle high-dimensional non-
stationary datasets, scalability issues, and difficulty in distinguishing
causation from correlation remain active research areas. Additionally, many
methods struggle with interpreting and providing actionable insights,
especially in complex real world applications such as climate science. By
addressing the challenges posed by traditional and machine learning based
approaches, our method aims to provide a comprehensive solution for effective
time series causal discovery.

##  7 Conclusion

We proposed TS-CausalNN, a score-based causal structure learning method for
non-linear and non-stationary time series data using custom 2D causal
convolution layer. The proposed method utilizes the power of the 2D neural
networks to learn the contemporaneous and time-lagged causal relationships of
all temporal variables simultaneously. The theory and computational procedure
of the method are established from the continuous formulation of the
acyclicity constraint of the causal DAG adapted from Zheng et al. [49]. The
parallel causal convolution layer block of the proposed model keeps the causal
contributor of each target/effect variable segregated from other target
variables, which gives better causal stability to the generated DAG structure.
The proposed model is very simple and user-friendly as any prior knowledge of
variable independence or the data generation model is not required and also
there is no assumption regarding error distribution. We conducted extensive
experiments on synthetic and real world complex time series datasets to
demonstrate the performance of the proposed causal discovery model. Based on
the empirical evaluation results, the proposed model demonstrates superior
causal graph learning capability compared to the baseline methods.

The code is available at http://tinyurl.com/TS-CausalNN.

## Acknowledgment

This work was partially supported by the DOE Office of Science Early Career
Research Program. This work was performed under the auspices of the U.S.
Department of Energy (DOE) by LLNL under contract DE-AC52-07NA27344. LLNL-
CONF-846980. Faruque, Ali and Wang were also partially supported by grants
OAC-1942714 and OAC-2118285 from the U.S. National Science Foundation (NSF).

## References

  * [1] Absar, S., Wu, Y., Zhang, L.: Neural time-invariant causal discovery from time series data. In: 2023 International Joint Conference on Neural Networks (IJCNN). pp. 1–8. IEEE (2023)
  * [2] Appiah, M.O.: Investigating the multivariate granger causality between energy consumption, economic growth and co2 emissions in ghana. Energy Policy 112, 198–208 (2018)
  * [3] Bahadori, M.T., Liu, Y.: On causality inference in time series. In: 2012 AAAI Fall Symposium Series (2012)
  * [4] Carles Balsells Rodas, Ruibo Tu, H. Kjellström: Causal discovery from conditionally stationary time-series. arXiv.org (2021)
  * [5] Cohen, J., Screen, J.A., Furtado, J.C., Barlow, M., Whittleston, D., Coumou, D., Francis, J., Dethloff, K., Entekhabi, D., Overland, J., et al.: Recent arctic amplification and extreme mid-latitude weather. Nature geoscience 7(9), 627–637 (2014)
  * [6] Debbi, H.: Causal Explanation of Convolutional Neural Networks, pp. 633–649. Springer International Publishing (2021)
  * [7] Endo, S., Fridlind, A.M., Lin, W., Vogelmann, A.M., Toto, T., Ackerman, A.S., McFarquhar, G.M., Jackson, R.C., Jonsson, H.H., Liu, Y.: Racoro continental boundary layer cloud investigations: 2. large-eddy simulations of cumulus clouds and evaluation with in situ and ground-based observations. Journal of Geophysical Research: Atmospheres 120(12), 5993–6014 (2015)
  * [8] Entner, D., Hoyer, P.O.: On causal discovery from time series data using fci. Probabilistic graphical models pp. 121–128 (2010)
  * [9] Ferdous, M.H., Hasan, U., Gani, M.O.: Cdans: Temporal Causal Discovery from Autocorrelated and Non-Stationary Time Series Data (2023)
  * [10] Fujiwara, D., Koyama, K., Kiritoshi, K., Okawachi, T., Izumitani, T., Shimizu, S.: Causal discovery for non-stationary non-linear time series data using just-in-time modeling. In: van der Schaar, M., Zhang, C., Janzing, D. (eds.) Proceedings of the Second Conference on Causal Learning and Reasoning. Proceedings of Machine Learning Research, vol. 213, pp. 880–894. PMLR (11–14 Apr 2023), https://proceedings.mlr.press/v213/fujiwara23a.html
  * [11] Gerhardus, A., Runge, J.: High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems 33, 12615–12625 (2020)
  * [12] Gerhardus, A., Runge, J.: High-recall causal discovery for autocorrelated time series with latent confounders. Advances in Neural Information Processing Systems 33, 12615–12625 (2020)
  * [13] Gu, J., Wang, Z., Kuen, J., Ma, L., Shahroudy, A., Shuai, B., Liu, T., Wang, X., Wang, G., Cai, J., et al.: Recent advances in convolutional neural networks. Pattern recognition 77, 354–377 (2018)
  * [14] Gustafson, W.I., Vogelmann, A.M., Li, Z., Cheng, X., Dumas, K.K., Endo, S., Johnson, K.L., Krishna, B., Fairless, T., Xiao, H.: The large-eddy simulation (les) atmospheric radiation measurement (arm) symbiotic simulation and observation (lasso) activity for continental shallow convection. Bulletin of the American Meteorological Society 101(4), E462–E479 (2020)
  * [15] Hasan, U., Hossain, E., Gani, M.O.: A survey on causal discovery methods for iid and time series data. Transactions on Machine Learning Research (2023)
  * [16] Heckerman, D.E., Horvitz, E.J., Nathwani, B.N.: Toward normative expert systems: Part i the pathfinder project. Methods of information in medicine 31(02), 90–105 (1992)
  * [17] Hinze, J.O.: Turbulence. McGraw-Hill (1975)
  * [18] Huang, B., Zhang, K., Zhang, J., Ramsey, J., Sanchez-Romero, R., Glymour, C., Schölkopf, B.: Causal discovery from heterogeneous/nonstationary data. The Journal of Machine Learning Research 21(1), 3482–3534 (2020)
  * [19] Huang, Y., Kleindessner, M., Munishkin, A., Varshney, D., Guo, P., Wang, J.: Benchmarking of data-driven causality discovery approaches in the interactions of arctic sea ice and atmosphere. Frontiers in big Data 4, 642182 (2021)
  * [20] Kalainathan, D., Goudet, O., Guyon, I., Lopez-Paz, D., Sebag, M.: Sam: Structural agnostic model, causal discovery and penalized adversarial learning (2018)
  * [21] Karamjit Singh, Garima Gupta, L. Vig, Gautam M. Shroff, P. Agarwal: Deep Convolutional Neural Networks for Pairwise Causality. arXiv.org (2017)
  * [22] Koller, D., Friedman, N.: Probabilistic graphical models: principles and techniques. MIT press (2009)
  * [23] Li, P., Mao, K.: Knowledge-oriented convolutional neural network for causal relation extraction from natural language texts. Expert Systems with Applications 115, 512–523 (1 2019)
  * [24] Luo, D., Chen, X., Dai, A., Simmonds, I.: Changes in atmospheric blocking circulations linked with winter arctic warming: A new perspective. Journal of Climate 31(18), 7661–7678 (2018)
  * [25] Luo, D., Chen, X., Overland, J., Simmonds, I., Wu, Y., Zhang, P.: Weakened potential vorticity barrier linked to recent winter arctic sea ice loss and midlatitude cold extremes. Journal of Climate 32(14), 4235–4261 (2019)
  * [26] Nauta, M., Bucur, D., Seifert, C.: Causal Discovery with Attention-Based Convolutional Neural Networks. Machine Learning and Knowledge Extraction 1(1), 312–340 (jan 7 2019)
  * [27] Pamfil, R., Sriwattanaworachai, N., Desai, S., Pilgerstorfer, P., Georgatzis, K., Beaumont, P., Aragam, B.: Dynotears: Structure learning from time-series data. In: International Conference on Artificial Intelligence and Statistics. pp. 1595–1605. PMLR (2020)
  * [28] Pearl, J.: Probabilistic reasoning in intelligent systems - networks of plausible inference. In: Morgan Kaufmann series in representation and reasoning (1991), https://api.semanticscholar.org/CorpusID:32583695
  * [29] Pearl, J.: Models, reasoning and inference. Cambridge, UK: CambridgeUniversityPress 19(2), 3 (2000)
  * [30] Peters, J., Janzing, D., Schölkopf, B.: Elements of causal inference: foundations and learning algorithms. The MIT Press, Cambridge, MA, USA (2017)
  * [31] Rajapakse, J.C., Zhou, J.: Learning effective brain connectivity with dynamic bayesian networks. Neuroimage 37(3), 749–760 (2007)
  * [32] Runge, J.: Discovering contemporaneous and lagged causal relations in autocorrelated nonlinear time series datasets. In: Peters, J., Sontag, D. (eds.) Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI). Proceedings of Machine Learning Research, vol. 124, pp. 1388–1397. PMLR (03–06 Aug 2020)
  * [33] Runge, J., Bathiany, S., Bollt, E., Camps-Valls, G., Coumou, D., Deyle, E., Glymour, C., Kretschmer, M., Mahecha, M.D., Muñoz-Marí, J., et al.: Inferring causation from time series in earth system sciences. Nature communications 10(1), 2553 (2019)
  * [34] Runge, J., Bathiany, S., Bollt, E.M., Camps-Valls, G., Coumou, D., Deyle, E.R., Glymour, C., Kretschmer, M., Mahecha, M.D., Muñoz-Marí, J., van Nes, E.H., Peters, J., Quax, R., Reichstein, M., Scheffer, M., Scholkopf, B., Spirtes, P., Sugihara, G., Sun, J., Zhang, K., Zscheischler, J.: Inferring causation from time series in earth system sciences. Nature Communications 10 (2019), https://api.semanticscholar.org/CorpusID:189819550
  * [35] Runge, J., Nowack, P., Kretschmer, M., Flaxman, S., Sejdinovic, D.: Detecting and quantifying causal associations in large nonlinear time series datasets. Science Advances 5(11), eaau4996 (2019). https://doi.org/10.1126/sciadv.aau4996
  * [36] Sanford, A.D., Moosa, I.A.: A bayesian network structure for operational risk modelling in structured finance operations. Journal of the Operational Research Society 63, 431–444 (2012)
  * [37] Serreze, M.C., Stroeve, J.: Arctic sea ice trends, variability and implications for seasonal ice forecasting. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 373(2045), 20140159 (2015)
  * [38] Shah, R.D., Peters, J.: The hardness of conditional independence testing and the generalised covariance measure. The Annals of Statistics 48(3) (Jun 2020). https://doi.org/10.1214/19-aos1857
  * [39] Simmonds, I.: Comparing and contrasting the behaviour of arctic and antarctic sea ice over the 35 year period 1979-2013. Annals of Glaciology 56(69), 18–28 (2015)
  * [40] Simmonds, I., Govekar, P.D.: What are the physical links between arctic sea ice loss and eurasian winter climate? Environmental Research Letters 9(10), 101003 (2014)
  * [41] SingH, G., Cuzzolin, F.: Recurrent Convolutions for Causal 3d CNNs. In: 2019 IEEE/CVF International Conference on Computer Vision Workshop (ICCVW). IEEE (10 2019)
  * [42] Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Liu, Z., Berner, J., Wang, W., Powers, J., Duda, M., Barker, D., et al.: A description of the advanced research wrf version 4. NCAR tech. note ncar/tn-556+ str 145 (2019)
  * [43] Spirtes, P., Glymour, C.N., Scheines, R.: Causation, prediction, and search. MIT press, Cambridge, MA, USA (2000)
  * [44] Sun, X., Schulte, O., Liu, G., Poupart, P.: Nts-notears: Learning nonparametric dbns with prior knowledge. In: International Conference on Artificial Intelligence and Statistics. pp. 1942–1964. PMLR (2023)
  * [45] Tank, A., Covert, I., Foti, N., Shojaie, A., Fox, E.B.: Neural granger causality. IEEE Transactions on Pattern Analysis and Machine Intelligence 44(8), 4267–4279 (2021)
  * [46] Wu, T., Wu, X., Wang, X., Liu, S., Chen, H.: Nonlinear Causal Discovery in Time Series. In: Proceedings of the 31st ACM International Conference on Information & Knowledge Management. ACM (oct 17 2022)
  * [47] Yao, Y., Luo, D., Dai, A., Simmonds, I.: Increased quasi stationarity and persistence of winter ural blocking and eurasian extreme cold events in response to arctic warming. part i: Insights from observational analyses. Journal of Climate 30(10), 3549–3568 (2017)
  * [48] Yu, Y., Chen, J., Gao, T., Yu, M.: Dag-gnn: Dag structure learning with graph neural networks. In: International Conference on Machine Learning. pp. 7154–7163. PMLR (2019)
  * [49] Zheng, X., Aragam, B., Ravikumar, P., Xing, E.P.: DAGs with NO TEARS: Continuous Optimization for Structure Learning. In: Advances in Neural Information Processing Systems (2018)
  * [50] Zheng, X., Dan, C., Aragam, B., Ravikumar, P., Xing, E.P.: Learning sparse nonparametric DAGs. In: International Conference on Artificial Intelligence and Statistics (2020)

APPENDIX

##  Appendix 0.A Synthetic Dataset 2

The synthetic dataset-2 is generated using the following equations and the
noise signals used in this dataset are generated by the Poisson distribution.
Here we used the exponential non-linearity in the same causal relationship of
the synthetic dataset-1 and excluded the seasonal features by removing the
cosine function. Also, the signal-to-noise ratio of this dataset is different
from dataset-1.

$S1_{t}=0.7e^{\frac{-S1^{2}_{t-2}\times S1^{2}_{t-5}}{2}}+\varepsilon 1$ (15)
$S2_{t}=2e^{\frac{S1^{2}_{t-1}}{2}}+0.5e^{\frac{S1^{2}_{t}}{2}}+\varepsilon 2$
(16) $S3_{t}=-5.05e^{\frac{-S1^{2}_{t-1}}{2}}+\varepsilon 3$ (17)
$S4_{t}=-1.15e^{\frac{-S1^{2}_{t-1}}{2}}+2.35e^{\frac{-S3^{2}_{t-1}}{2}}+1.5e^{\frac{-S4^{2}_{t-1}}{2}}+3e^{\frac{-S3^{2}_{t}}{2}}+\varepsilon
4$ (18)

##  Appendix 0.B Method Hyperparameters

To find the best hyperparameters for baseline methods we started using the
parameters suggested by the authors and gradually tuned those values to obtain
better evaluation results. The best results obtained with tuned
hyperparameters are reported in the comparative analysis. The parameters used
for generating evaluation results are given here.

  * •

PCMCI

    * –

Conditional Independence Test = ParCorr

    * –

$tau\\_max$ = Maximum time lag

    * –

$pc\\_alpha$ = None [So the model will use the optimal value from the list
$\\{0.05,0.1,0.2,0.3,0.4,0.5\\}$]

    * –

$alpha\\_level$ = 0.01

  * •

PCMCI+

    * –

Conditional Independence Test = ParCorr

    * –

$tau\\_max$ = Maximum time lag

    * –

$pc\\_alpha$ = None [So the model will use the optimal value from the list
$\\{0.001,0.005,0.01,0.025,0.05\\}$]

  * •

NOTEARS-MLP

    * –

lambda1 = 0.01

    * –

lambda2 = 0.01

    * –

rho = 1.0

    * –

alpha = 0.0

    * –

$w\\_threshold$ = 0.3

  * •

NTS-NOTEARS

    * –

lambda1 = 0.0005

    * –

lambda2 = 0.001

    * –

$w\\_threshold$ = 0.3

    * –

rho = 1.0

    * –

alpha = 0.0

    * –

$number\\_of\\_lags$ = Maximum time lag

  * •

DYNOTEARS

    * –

$tau\\_max$ = Maximum time lag

    * –

$w\\_threshold$ = 0.01

    * –

$lambda\\_w$ = 0.05

    * –

$lambda\\_a$ = 0.05

  * •

Proposed Method

    * –

lambda1 = 0.1

    * –

alpha=0.0

    * –

rho = 1.0

    * –

beta = 0.1

    * –

gamma = 0.25

    * –

$w\\_threshold$ = 0.3

##  Appendix 0.C Robustness tests using data with different Signal-to-Noise
ratios

We have generated three versions of the synthetic dataset-1 to analyze the
performance of the proposed model for different signal-to-noise ratios (SNRs)
of datasets. The SNRs of these synthetic dataset versions are given in Table
6. From the SNR values, we can see that the Synthetic dataset-3 has a higher
signal component concerning the noise. And the synthetic dataset-1 has a
comparatively lower signal component.

Table 6: Signal-to-Noise ratio of synthetic datasets. Dataset | S1 | S2 | S3 | S4 | Average
---|---|---|---|---|---
Synthetic-1.1 | 0.10 | 0.005 | 0.63 | 0.01 | 0.25
Synthetic-1.2 | 0.19 | 0.001 | 1.93 | 0.06 | 0.55
Synthetic-1.3 | 0.52 | 0.02 | 1.37 | 0.56 | 0.62

All the baseline models are applied to these synthetic datasets to quantify
the changes in the evaluation performance of the generated causal graphs due
to the different noise magnitudes. The evaluation results for the summary and
full temporal causal graphs are summarized in Table 7 and Table 8. The best
results are marked in bold text and underlined values represent the second
best score. Some baseline models generated empty graphs, so we placed a ”NE”
for those methods in the comparison table instead of mentioning the zero (0)
F1 score and FDR as 1. From Table 7 we can see that the proposed method
achieved the best SHD for all the datasets. For second and third datasets
yielded the best F1 score with the lowest FDR, which means the edges predicted
by the proposed model are mostly true.

Table 7: Comparison of the predicted summary causal graph by different methods for synthetic datasets with different SNR. Method | Synthtic-1.1 | Synthtic-1.2 | Synthtic-1.3
---|---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 10 | 0.54 | 0.62 | 10 | 0.54 | 0.62 | 10 | 0.54 | 0.62
PCMCI+ [32] | 8 | 0.60 | 0.57 | 10 | 0.54 | 0.62 | 9 | 0.54 | 0.62
NOTEARS-MLP [50] | 6 | 0.25 | 0.50 | 8 | 0.20 | 0.75 | 6 | 0.25 | 0.50
NTS-NOTEARS [44] | ne | ne | ne | ne | ne | ne | 6 | 0.25 | 0.50
DAG-GNN [48] | ne | ne | ne | ne | ne | ne | 3 | 0.66 | 0.00
DYNOTEARS [27] | 7 | 0.45 | 0.57 | 7 | 0.45 | 0.57 | 7 | 0.63 | 0.53
Proposed | 5 | 0.54 | 0.40 | 4 | 0.66 | 0.33 | 2 | 0.83 | 0.16

The evaluation results of the predicted full temporal causal graphs by all
baseline methods for the synthetic datasets are shown in Table 8. Carefully
observing the evaluation result we see that except for the zero (0.00) FDR of
DAG-GNN method, the proposed method achieved the best scores in SHD and FDR
evaluation criteria for all three datasets. For the first synthetic dataset,
the proposed method generated the second best F1 score that is lower than
PCMCI+. Overall from the comparative analysis of the evaluation results, we
can conclude that the proposed model worked better for different SNRs.

Table 8: Comparison of the predicted full causal graph by different methods for synthetic datasets with different SNR. Method | Synthtic-1.1 | Synthtic-1.2 | Synthtic-1.3
---|---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 50 | 0.19 | 0.89 | 54 | 0.20 | 0.88 | 61 | 0.18 | 0.89
PCMCI+ [32] | 46 | 0.20 | 0.88 | 56 | 0.17 | 0.90 | 39 | 0.23 | 0.86
NOTEARS-MLP [50] | 26 | 0.31 | 0.80 | 23 | 0.37 | 0.76 | 27 | 0.30 | 0.81
NTS-NOTEARS [44] | ne | ne | ne | ne | ne | ne | 7 | 0.25 | 0.50
DAG-GNN [48] | ne | ne | ne | ne | ne | ne | 10 | 0.00 | 1.00
DYNOTEARS [27] | 24 | 0.20 | 0.86 | 23 | 0.20 | 0.86 | 23 | 0.34 | 0.78
Proposed | 12 | 0.33 | 0.72 | 11 | 0.42 | 0.66 | 9 | 0.47 | 0.60

##  Appendix 0.D Study with datasets with only lagged causal links

###  0.D.1 Datasets

To evaluate the performance of our proposed causal discovery method on the

we have used synthetic datasets. As we know the ground truth causal graph for
the synthetic datasets, we can measure and compare the learned causal graph
easily. We generate a time series dataset consisting of four non-linear
variables using Gaussian white noise $\varepsilon$ (Dataset-1). The
mathematical description of each variable is given in Equations 10 to 13. The
non-linear characteristic is incorporated in the generation of synthetic
dataset to mimic the dynamic properties of real world natural system data. The
corresponding true causal graph for the time-series data is given in Figure
4a.

$S1_{t}=2\\{cos(\frac{t}{10})+log(|S1_{t-2}-S1_{t-5}|+1)\\}+0.1\varepsilon 1$
(19) $S2_{t}=12e^{\frac{S1^{2}_{t-1}}{2}}+\varepsilon 2$ (20)
$S3_{t}=-10.5e^{\frac{-S1^{2}_{t-1}}{2}}+\varepsilon 3$ (21)
$S4_{t}=-11.5e^{\frac{-S1^{2}_{t-1}}{2}}+13.5e^{\frac{-S3^{2}_{t-1}}{2}}+1.2e^{\frac{-S4^{2}_{t-1}}{2}}+\varepsilon
4$ (22)

We have also generated another version of the synthetic dataset (Dataset-2)
from the previous causal graph without the cosine part of the data generation
formulas. For this case, we used the Poisson distribution to incorporate
random noise into the dataset. The goal is to find the performance and
robustness of the causal discovery methods for the time series dataset without
having wave-shape-like features with non-Gaussian noise. The formula used to
generate this dataset is available in Appendix 0.A.

###  0.D.2 Results

To evaluate the performance of the selected state-of-the-art methods we
compared the predicted causal graph in both the summary and full temporal
graph setting. The qualitative comparison of the baseline methods is reported
in Table 2 and Table 1. The best results are marked in bold text and
underlined values represent the second best score. To get better prediction
results, we have applied each method multiple times on the same dataset and
reported the best result. The comparative analysis of the summary graph (Table
2) shows that the proposed method achieves the best scores for all three
evaluation criteria for both datasets except for the FDR of DAG-GNN method.
For synthetic dataset-1, DAG-GNN method yields a lower score in terms of FDR
whereas its SHD is close to the proposed model. Since, FDR only relies on the
number of predicted edges, the $(0.00^{*})$ FDR score of DAG-GNN cannot be
considered a perfect score as the few edges predicted by the method align with
the true causal graph but it missed predicting the remaining edges (Figure
4a).

Table 9: Comparison of the predicted summary causal graph by different methods for synthetic datasets. Method | Dataset-1 | Dataset-2
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 10 | 0.54 | 0.62 | 9 | 0.57 | 0.60
PCMCI+ [32] | 9 | 0.54 | 0.62 | 7 | 0.63 | 0.53
NOTEARS-MLP [50] | 6 | 0.25 | 0.50 | 6 | 0.25 | 0.50
NTS-NOTEARS [44] | 6 | 0.25 | 0.50 | 7 | 0.36 | 0.60
DAG-GNN [48] | 3 | 0.66 | 0.00∗ | 7 | 0.22 | 0.66
DYNOTEARS [27] | 7 | 0.63 | 0.53 | 7 | 0.53 | 0.55
Proposed | 2 | 0.83 | 0.16 | 3 | 0.72 | 0.20

To compare the results of the predicted full causal graph, we considered both
the directed edges and the time lag of the edge. From Table 1, we can see that
our proposed method got the best F1 score and the second-lowest SHD and FDR
for dataset-1. For dataset-2, the proposed method yields the best results for
all three quality measures. For dataset-1, NTS-NOTEARS generated the target
causal graph with smaller SHD and lower F1 scores compared to the proposed
method. This can be interpreted as the number of edges generated by NTS-
NOTEARS is lower than the number of edges in the true causal graph (low F1
score) and also the 0.5 FDR value means that half of the predicted edges are
wrong. Similarly, zero F1 score of DAG-GNN method means that it was not able
to find any causal links from the true causal graph. From observing the
evaluation results of the full causal graph, we can find that the proposed
model is capable of generating fewer false causal links compared to other
baseline models.

Table 10: Comparison of the predicted temporal full causal graph by different methods for synthetic datasets. Method | Dataset-1 | Dataset-2
---|---|---
| SHD | F1 | FDR | SHD | F1 | FDR
PCMCI [35] | 61 | 0.18 | 0.80 | 32 | 0.30 | 0.82
PCMCI+ [32] | 39 | 0.23 | 0.86 | 23 | 0.37 | 0.76
NOTEARS-MLP [50] | 27 | 0.31 | 0.81 | 24 | 0.33 | 0.77
NTS-NOTEARS [44] | 7 | 0.22 | 0.50 | 9 | 0.31 | 0.66
DAG-GNN [48] | 10 | 0.00 | 1.00 | 10 | 0.00 | 1.00
DYNOTEARS [27] | 23 | 0.34 | 0.78 | 28 | 0.17 | 0.88
Proposed | 9 | 0.47 | 0.60 | 7 | 0.46 | 0.50

##  Appendix 0.E Non-Stationarity Test Results for Real World Datasets

The non-stationarity feature of the real world TKE and Arctic Sea Ice datasets
is evaluated using the Augmented Dickey–Fuller test (ADF) and Kwiatkowski-
Phillips-Schmidt-Shin test (KPSS) statistical test methods for time series
data. The ADF test method assumes a null hypothesis: the time series has a
unit root and is not stationary. Then try to reject the null hypothesis and if
failed to be rejected, it suggests the time series is not stationarity. For
the ADF test, if the p-value of a time series is higher than the 0.05 alpha
level the null hypothesis cannot be rejected. So the time series is not
stationary. The KPSS test works in a somewhat similar manner to the ADF test
but assumes an inverse null hypothesis. The null hypothesis of the KPSS method
is that the time series is stationary. If the p-value is less than 0.05 alpha
level, we can reject the null hypothesis and derive that the time series is
not stationary.

Table 11: Stationarity test results for variables of the TKE dataset. Variables | ADF Test | KSSP Test
---|---|---
| P-value | Stationary | P-value | Stationary
SH | 0.32 | No | 0.01 | No
BU | 0.72 | No | 0.01 | No
TEND | 0.66 | No | 0.01 | No
TKE | 0.58 | No | 0.01 | No

The statistical non-stationarity test results for the TKE dataset are given in
Table 11 and the results for the Arctic Sea Ice data are available in Table
12. The non-stationarity test results revealed that the TKE dataset contains
only non-stationary variables and both test methods have agreement on the test
outcome. In the Arctic Sea Ice dataset, the ADF test found 4 non-stationary
variables and the KSSP method found 3 non-stationary variables. Therefore we
can say that the Arctic Sea Ice data have a mixture of both non-stationary and
stationary variables.

Table 12: Stationarity test results for variables of the Arctic Sea Ice dataset. Variables | ADF Test | KSSP Test
---|---|---
| P-value | Stationary | P-value | Stationary
HFLX | 0.10 | No | 0.07 | Yes
CC | 0.00 | Yes | 0.02 | No
SW | 0.00 | Yes | 0.10 | Yes
U10m | 0.00 | Yes | 0.10 | Yes
SLP | 0.00 | Yes | 0.10 | Yes
PRE | 0.00 | Yes | 0.09 | Yes
ICE | 0.75 | No | 0.01 | No
LW | 0.28 | No | 0.10 | Yes
V10m | 0.00 | Yes | 0.10 | Yes
CW | 0.00 | Yes | 0.10 | Yes
GH | 0.01 | Yes | 0.10 | Yes
RH | 0.11 | No | 0.01 | No
Table 13: Variables in the Arctic Sea Ice Data Variable Abbreviation | Full Name
---|---
HFLX | Heat Flux
CC | Cloud Cover
SW | Net Shortwave Flux
U10m | Zonal wind at 10m
SLP | Sea level pressure
PRE | Total Precipitation
ICE | Sea Ice
LW | Net Longwave Flux
V10m | Meridional Wind at 10m
CW | Total Cloud Wate Path
GH | Geopotential Height
RH | Relative Humidity

##  Appendix 0.F More Results of the TKE Data

For the TKE dataset, the ground truth causal graph, Figure 8a, refers to the
domain knowledge based graph generated with the aid of domain scientists.
Here, the most important relationship is between $TEND$ and $TKE$ which is a
unidirectional causal relation. Observing the predicted causal graphs from
baseline methods, we see that the methods either miss predicting this
important edge or predict a bidirectional relation between $TEND$ and $TKE$,
with an exception of NOTEARS which correctly identifies this relation.
Overall, our method and NOTEARS show promising results with NOTEARS having the
lowest mis-prediction while our method succeeds in capturing the most
important causal relation $TEND\rightarrow TKE$, highlighting its potential to
make meaningful scientific contributions.

Figure 8: Visualization of summary causal graphs generated by different
methods for the TKE dataset.

Figure 9: Visualization of full temporal causal graphs generated by different
methods for the TKE dataset with time lag 5.

##  Appendix 0.G More Results on the Sea Ice Data

Here we present the causal graphs predicted by the proposed and baseline
models on the Arctic Sea Ice dataset with a time lag of 12, referring to the
annual 12-month periods respectively. We use the abbreviated terms in the
causal graphs, whereas, the full form of each of the data variables is given
in Table 13. The true domain graph is given in Figure 10a.

From the predicted graphs, we can see that NTS-NOTEARS predicted the lowest
number of edges (one reason behind its low SHD), however, amongst the
predicted edges, the false detection rate (FDR) is very high, for instance
$ICE\leftarrow RH$. Likewise, NOTEARS-MLP predicted the highest number of
edges with comparable overall performance to our proposed method.

Since the total number of edges in the predicted and domain causal graph is
very high, here, we will focus only on the edges corresponding to the
scientific problem of identifying causes of sea ice melt. In this regard, we
see NOTEARS-MLP, DYNOTEARS, NTS-NOTEARS and PCMCI+, fail to identify any
correct causes (incoming edge) of $ICE$, whereas PCMCI, PCMCI+, DAG-GNN fail
to identify any correct effect (outgoing edge) of $ICE$. Our proposed method
correctly predicts 3 out of 4 outgoing edges from $ICE$, further highlighting
its potential to discover complex real world causal relations.

Figure 10: Visualization of summary causal graphs generated by different
methods for the Arctic Sea Ice dataset with time lag 12.