Causal Discovery from Temporal Data: An Overview and New Perspectives

Abstract.

Temporal data, representing chronological observations of complex systems, has always been a typical data structure that can be widely generated by many domains, such as industry, medicine and finance. Analyzing this type of data is extremely valuable for various applications. Thus, different temporal data analysis tasks, e.g., classification, clustering and prediction, have been proposed in the past decades. Among them, causal discovery, learning the causal relations from temporal data, is considered an interesting yet critical task and has attracted much research attention. Existing causal discovery works can be divided into two highly correlated categories according to whether the temporal data is calibrated, i.e., multivariate time series causal discovery, and event sequence causal discovery. However, most previous surveys are only focused on the time series causal discovery and ignore the second category. In this paper, we specify the correlation between the two categories and provide a systematical overview of existing solutions. Furthermore, we provide public datasets, evaluation metrics and new perspectives for temporal data causal discovery.

1. Introduction

Temporal data recording the status changing of complex systems is widely collected by different application domains, such as social networks, bioinformatics, neuroscience and finance, etc.. As one of the most popular data structural, temporal data consists of attribute sequences ordered by time. Owing to the rapid development of sensors and computing devices, research works on temporal data analysis are emerging in recent years. Different approaches have been proposed for different tasks such as classification(Ismail Fawaz et al., 2019; Ratanamahatana and Keogh, 2004), clustering(Aghabozorgi et al., 2015; Liao, 2005), prediction(Weigend, 2018), causal discovery(Edinburgh et al., 2021; Krakovská et al., 2018), etc..

Among these tasks, causal discovery recognizing the causal relations between many temporal components has become a challenging yet critical task for temporal data analysis. The learned causal structures could be beneficial for explaining the data generation process and guiding the design of data analysis methods. According to whether the data is calibrated, the temporal data for causal discovery can be categorized into two groups, i.e., multivariate time series (MTS) and event sequences. Therefore, existing causal discovery methods can also be divided into two groups respectively. In this survey, we aim to provide a thoughtful overview and summarize the frontiers of temporal data causal discovery.

MTS data, describing the calibrated states of multiple variables changing over time, is a general kind of temporal data in many domains. Discovering causal relations from MTS could be beneficial to the explainability and robustness of data analysis models. However, the definitions of causal relations are not unique, leading to different solutions. Accordingly, existing works can be grouped into four categories, i.e., constraint-based methods, score-based methods, functional causal model (FCM)-based methods and Granger causality-based methods. Besides, there also exist some perspectives such as Takens’ causality and differential equations. In this paper, we will specify the main idea and recent advances for each category.

Another task discussed in this survey is the causal discovery from event sequences, which infers causal relationships within irregularly and asynchronously observed time series. Specifically, it takes a sequence of different events as the input and outputs a causal graph representing the causal interactions between different events. This task is of great importance since most real-world events cannot emerge within a fixed time interval. In accordance with the MTS task, we classify the corresponding methods into three main categories: constraint-based, score-based, and Granger causality-based methods. Among these three categories, Granger causality-based methods, especially Granger causality-based Hawkes process models, are well-developed since a natural match-up exists between Granger causality and Hawkes processes. We will further describe these approaches in detail within this review.

Table 1. Highlights of existing reviews on causal discovery.

Reviews	Multivariate Time-series					Event Sequence	Highlights
	Constrain-based	Score-based	FCM-based	Granger	Deep Learning
(Glymour et al., 2019)	No	No	No	Yes	No	No	An overview for causal discovery methods with practical issues and insightful guidelines
(Guo et al., 2021)	No	No	No	No	No	No	Causal discovery methods dealing with big data (high-dimensional, mixed data) are reviewed
(Vowels et al., 2023)	No	No	No	Yes	No	No	A more extensive coverage of continuous optimization approaches compared to other surveys
(Chen et al., 2022b)	No	No	No	No	No	No	A wider concept of deep learning causal discovery methods is introduced
(Moraffah et al., 2021)	Yes	No	Yes	Yes	No	No	The first survey covers the current progress to analyze time series from a causal perspective
(Shojaie and Fox, 2021)	No	No	No	Yes	Yes	No	Recent advances including network-form and more general notions of Granger causality
(Assaad et al., 2022b)	Yes	Yes	Yes	Yes	No	No	A recent and comprehensive review for causal discovery in time series with comparative evaluations
Ours	Yes	Yes	Yes	Yes	Yes	Yes	A systematic review of causal discovery in both MTS and event sequence, with new perspectives

Recently, many surveys (Glymour et al., 2019; Guo et al., 2021; Vowels et al., 2023; Chen et al., 2022b; Moraffah et al., 2021; Shojaie and Fox, 2021; Assaad et al., 2022b; Kitson et al., 2021; Deng et al., 2022b; Heinze-Deml et al., 2018) have been published to summarize the progress of causal discovery. We compared the representative reviews and their highlight points in Table 1. As shown, these surveys fall into two lines. Research works in the first line (Glymour et al., 2019; Vowels et al., 2023; Guo et al., 2021; Chen et al., 2022b) discuss the general causal discovery problem in different perspectives. For example, (Glymour et al., 2019) provide a brief review of the computational causal discovery methods. (Vowels et al., 2023) focus on the flurry developments of continuous optimization approaches. To handle big data, both causal inference and causal discovery methods based on machine learning are introduced in (Guo et al., 2021). Moreover, deep learning causal discovery methods are reviewed in different variable paradigms (Chen et al., 2022b), where the causal relations in data are discussed from a broader perspective. In these papers, temporal data was taken as one special application and many data-specified methods are not included. The surveys in the second line focus on temporal data causal discovery. As illustrated in Table 1, causal discovery methods for bivariate time series are reviewed in (Edinburgh et al., 2021; Krakovská et al., 2018). The approaches for causal inference in time series are recently reviewed in (Moraffah et al., 2021; Shojaie and Fox, 2021). The recent work (Assaad et al., 2022b) discusses and comparatively evaluates the existing solutions of time series causal discovery. Nevertheless, causal discovery methods for event sequences are ignored in these reviews. In this paper, we not only provide a thoughtful overview of causal discovery methods of the two kinds of temporal data but also give an analysis of the connections and differences between them.

Next, we first introduce the background and preliminary of the causal discovery problem in Section 2. The recent progress of causal discovery from MTS and event sequences are specified in Section 3 and Section 4 respectively. After that, we provide an overview of the applications of temporal data causal discovery in Section 5 and summarize the available resources in Section 6. At last, we discuss the limitations and new perspectives of recent temporal data causal discovery methods in Section 7. The whole framework of this survey is shown in Figure 1.

for tree= grow’=east, anchor=west, node options=draw, thick, font=, align=center, , edge=semithick, forked edges, l sep=8mm, s sep=8mm, text width=2.3cm, fork sep = 2mm, , [Temporal Causal Discovery, fill=col1, parent, rotate=90,font=, for tree=s sep=2.0mm, [MTS
Causal Discovery
(§ 3, Table 3), font=, for tree=child, fill=col3, text width = 4.0cm, [Constraint-Based Approaches
(§ 3.1), text width = 3.6cm, [With Causal Sufficiency,fill = col3,text width = 2.8cm [ oCSE (Sun et al., 2015), PCGCE (Assaad et al., 2022a), PCMCI (Runge et al., 2019b; Runge, 2020), text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] [Without Causal Sufficiency,fill = col3,text width = 2.8cm [ ANLTSM (Chu and Glymour, 2008), tsFCI (Entner and Hoyer, 2010), SVAR-FCI (Malinsky and Spirtes, 2018), LPCMCI (Gerhardus and Runge, 2020), text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] ] [Score-Based Approaches
(§ 3.2), text width = 3.6cm, for tree = child, fill = col3,text width = 3.6 cm [Combinatorial Search,fill = col3,text width = 2.8cm [ Structural EM (Friedman et al., 1998), Greedy Hill-climbling Search (Peña et al., 2005), Structural Constraints (de Campos and Ji, 2011), etc., text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] [Continuous Optimization,fill = col3,text width = 2.8cm [ DYNOTEARS (Pamfil et al., 2020), NTS-NOTEARS (Sun et al., 2021), IDYNO (Gao et al., 2022), text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] ] [FCM-Based Approaches
(§ 3.3), for tree = child, fill = col3,text width = 3.6 cm [Independent Component Analysis,fill = col3,text width = 2.8cm [ VAR-LiNGAM (Hyvärinen et al., 2008, 2010a), MCD (Schaechtle et al., 2013), NCDH (Wu et al., 2022b), text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] [Additive Noise Model,fill = col3,text width = 2.8cm [ TiMINo (Peters et al., 2013), NBCB (Assaad et al., 2021), text width=5.5 cm,draw=colline1,line width=1.2pt,fill=col1,] ] ] [Granger Causality
Based Approaches
(§ 3.4), for tree = child, fill = col3,text width = 3.6cm [HSIC-Lasso-GC (Ren et al., 2020), (R)NN-GC (Montalto et al., 2015; Wang et al., 2018), MPIR (Wu et al., 2020), NGC (Tank et al., 2022), eSRU (Khanna and Tan, 2020), SCGL (Xu et al., 2019), GVAR (Marcinkevics and Vogt, 2021), TCDF (Nauta et al., 2019), CR-VAE (Li et al., 2023), InGRA (Chu et al., 2020), ACD (Löwe et al., 2022), etc., text width=9.4 cm,draw=colline1 ,line width=1.0pt,fill=col1, ] ] [Others
(§ 3.5), for tree = child, fill = col3,text width = 3.6cm [Information-theoretic Statistics (Schreiber, 2000; Runge et al., 2012a; Sun and Bollt, 2014), Differential Equation Based Methods (Voortman et al., 2010; Bellot et al., 2022), Nonlinear State-space Methods (Sugihara et al., 2012), Logic-based Methods (Kleinberg and Mishra, 2009), Hybrid Methods (Li et al., 2016), etc., text width=9.4 cm,draw=colline1 ,line width=1.0pt,fill=col1, ] ] ] [Event Sequence Causal Discovery
(§ 4), font=, for tree=child, fill=col2, text width = 4.0cm [Multivariate Point Process
(§ 4.1), for tree = child, fill = col2, text width=3.6cm, [Basics: Intensity Function, Log-likelihood, text width=9.4 cm,draw=colline1 ,line width=1.0pt,fill=col1, ] ] [Granger Causality
Based Approaches
(§ 4.2), for tree = child, fill = col2,text width=3.6cm [GLM Point Process, text width=2.8cm, [ GLM Model (Kim et al., 2011), text width=5.5cm,draw=colline2,line width=1.2pt,fill=col1] ] [Hawkes Process, text width=2.8cm, [MLE-SGLP (Xu et al., 2016), THP (Cai et al., 2021), $L_0$Hawkes (Idé et al., 2021), HGEM (Yu et al., 2020), NPHC (Achab et al., 2017), GC-nsHP (Chen et al., 2022a), MDLH (Jalaldoust et al., 2022), etc., text width=5.5cm,draw=colline2,line width=1.0pt,fill=col1] ] [Wold Process, text width=2.8cm, [Granger-Busca (de Figueiredo et al., 2018), VI-MWP (Etesami et al., 2021), text width=5.5cm,draw=colline2,line width=1.0pt,fill=col1] ] [Neural Point Process, text width=2.8cm, [CAUSE (Zhang et al., 2020), text width=5.5cm,draw=colline2,line width=1.0pt,fill=col1] ] ] [Others
(§ 4.3), for tree = child, fill = col2,text width=3.6cm [Constraint-Based Approaches, text width=2.8cm, [ MMP-LR/NI (Bhattacharjya et al., 2022), CA (Meek, 2014), text width=5.5cm,draw=colline2,line width=1.2pt,fill=col1] ] [Score-Based Approaches, text width=2.8cm, [PGEM (Bhattacharjya et al., 2018), text width=5.5cm,draw=colline2,line width=1.0pt,fill=col1] ] ] ] [Applications
(§ 5, Table 5),font=, for tree=child, fill=col5, text width = 4.0cm, [Scientific Endeavors, for tree = child, fill = col5,text width = 3.6cm, [ Earth Science, Neuroscience, Bioinformatics, etc., text width=9.4cm,draw=colline3,line width=1.0pt,fill=col1] ] [Industrial
Implementations, for tree = child, fill = col5,text width = 3.6cm, [Anomaly Detection, Root Cause Analysis, Business Intelligence in Online Systems, Video Analysis, Urban Data Analysis, Clinical Data Analysis, etc., text width=9.4cm,draw=colline3,line width=1.0pt,fill=col1] ] ] [Discussions &
New Perspectives
(§ 7) , font=, for tree=child, fill=col4, text width = 4.0cm, [Challenges & Practical Considerations
(§ 7.1), for tree = child, fill = col4, text width=3.6cm, [1) Non-stationarity, 2) Heterogeneity, 3) Unobserved Confounders, 4) Subsampling, 5) Expert Knowledge, text width=9.4 cm,draw=colline4 ,line width=1.0pt,fill=col1, ] ] [New Perspectives
(§ 7.2), for tree = child, fill = col4,text width = 3.6cm [1) Amortized Paradigm, 2) Supervised Paradigm, 3) Causal Representation Learning, text width=9.4cm,draw=colline4,line width=1.0pt,fill=col1] ] ] ]

Figure 1. Framework for causal discovery from temporal data.

2. Background & preliminaries

This section begins with the definition of key concepts and assumptions in causal discovery, followed by an overview of three causal graph representations applicable to temporal data. Finally, the problem definitions for causal discovery from MTS and event sequences will be presented.

Table 2. Main notations used in this survey.

Notation	Description
$d,E$	number of time-series variate, and of event types, respectively
$x_i^t$	the $i$-th time series at time $t$ in multivariate time series
$N_e(t)$	the number of the event $e$ occurrences before time $t$
$⟂⟂,\cancel{}⟂⟂$	independent, and not independent
$V,U$	the set of endogenous variables, and of exogenous variables, respectively
$𝒢$	causal graph
$Pa(x_i)$	the parent nodes of $x_i$

2.1. Key concepts and assumptions in causal discovery

Some key concepts serve as the foundation for inferring causal relationships from temporal data. We establish this common ground before discussing research works. Afterward, we present formal definitions for the structural causal model, $d$-separation, causal Markov condition, causal identifiability and causal minimality with notations detailed in Table 2.

Figure 2. Basic DAGs and a simple structural causal model.

Structural Causal Model (SCM). Pearl’s comprehensive theory of causality, as presented in (Pearl, 2009), enables us to draw causal conclusions from observations using causal hierarchy (PCH) (Pearl and Mackenzie, 2018). From that, the structural causal model is defined as a graphical representation of causal relationships that captures how interventions on one or more variables affect the values of other variables in the data generation mechanism. Formally, SCM can be represented in a 4-tuple , where $V,U$ denote the set of endogenous and exogenous variables respectively, $P(U)$ is the distribution of exogenous variables, and $F$ represents the set of the mapping function. Specifically, for $f_i\in F$, the model $x_i:=f_i(Pa(x_i),u_i),i=1,\mathrm{\dots },d$ indicates the assignment of the value $x_i$ to a function of its structural parents $Pa(x_i)$ and exogenous variable $u_i$. For each SCM, we can yield a causal graph DAG $G$ by adding one vertex for each $x_i$ and directing edges from each parent variable in $Pa(x_i)$ (the causes) to child $x_i$ (the effect). The relationship of the SCM and the corresponding DAG is shown in figure 2 (a)(b).

d-separation. $d$-separation is a criterion for determining the absence of causal effects between two sets of variables in a graphical model. Two sets of variables are said to be $d$-separated if every path between them is blocked. In formal, a set of variables $𝐒$ d-separates two variables if $𝐒$ blocks all paths between them. For the given causal graph in figure 2 (d)(e)(f), two vertices $x_1,x_3$ are d-separated by the set of vertices $𝐒$ if $x_2\in 𝐒$. As for the relations in figure 2 (c) (a.k.a., a v-structure or collider), $x_1,x_3$ are also d-separated if $x_2$ and none of the descendants of $x_2$ are in set $𝐒$.

$d$-separation is a fundamental concept in causal discovery because it provides a criterion for determining whether two sets of variables are causally related. If two sets of variables are $d$-separated, then there is no direct or indirect causal effect between them, and they can be considered independent given the observed variables. Conversely, if two sets of variables are not $d$-separated, then there may be a direct or indirect causal effect between them that needs to be accounted for when inferring causal relationships from data. Thus, $d$-separation is an essential tool for identifying causal relationships in graphical models.

Causal Markov Condition. In the causal graph of SCM, each variable is independent of its non-effects given its direct causes (Pearl, 2009). In other words, a variable is conditionally independent of its non-effects (i.e., variables that do not directly cause it) given its parents (i.e., variables that directly cause it). This condition plays an essential role in causal inference. It enables the identification of causal effects from non-experimental data. Formally, the causal Markov condition implies the joint distribution can be factorized according to the following decomposition:

$$P(𝐱)=\prod _i^{d}P(x_i|Pa(x_i))$$

Markov Equivalence Class(MEC). Two graphical models belong to the same MEC if they entail the same set of conditional independence relations among the observed variables, regardless of the specific structure of the graph. For example, the causal diagrams in Figure 2 (d)(e)(f) imply the same d-seperation information $x_1⟂⟂x_3|x_2$ and belong to the same MEC. MEC is important because it enables us to identify the minimal set of conditions necessary for inferring causal relationships from non-experimental data.

Causal Identifiability. A causal effect is identifiable if it can be estimated without making any untestable assumptions or invoking additional information beyond the observed variables. This means that all causal graphs in the same MEC represent equivalent causal structures from an observational viewpoint. In general, causal identifiability requires that the causal graph is acyclic and that all backdoor paths between the treatment and outcome variables are blocked. If these conditions are met, the causal effect can be identified using the $do$-calculus or other causal inference techniques. Thus, the prerequisite of causal discovery is that causal relationships can be identifiable.

Causal Minimality. Consider a DAG $𝒢$ and a probability distribution $P$, $P$ is said to satisfy the causal minimality with respect to $𝒢$ if $P$ is Markovian with respect to $𝒢$ but not to any proper subgraph of $𝒢$. It indicates that all the variables are necessary and sufficient to accurately represent the causal relationships while excluding any variables that do not contribute to the causal mechanism. A distribution is minimal with respect to the causal graph if and only if there is no node that is conditionally independent of any of its parents, given the remaining parents. In other words, all the parents are “active”.

Building on the aforementioned concepts, we introduce three assumptions, causal sufficiency, faithfulness, and temporal priority, which are the untestable foundations of causal discovery.

Causal Sufficiency. A set of variables is causally sufficient if all common causes of all variables are observed (Spirtes et al., 2000). This assumption indicates that the causal graph in SCM can reflect the truth data generation process and there is no hidden confounder. Under the assumption of causal sufficiency, the majority of causal discovery algorithms presume that the causal structure can be depicted as a DAG.

Faithfulness. Faithfulness asserts that all conditional independence relations of $P$ that hold in the observed data are entailed by the causal model $𝒢$, and conversely, all conditional independence relations implied by the causal model are also held in the observed data. Note that faithfulness implies causal minimality. If $P$ is faithful and Markovian with respect to $𝒢$, then the causal minimality is satisfied.

Figure 3. Example of the faithfulness and causal minimality.

Intuitively, faithfulness is not easy to understand. We try to clarify it with an example (Peters et al., 2017). As shown in Figure 3, we assume the generation process of ${𝒢}_1$ as a linear Gaussian SCM:

$X$ $:=N_X$

$Y$ $:=aX+N_Y$

$Z$ $:=bY+cX+N_Z$

The noise variables $N_X\sim 𝒩(0,{\sigma }_x^2)$, $N_Y\sim 𝒩(0,{\sigma }_y^2)$ and $N_Z\sim 𝒩(0,{\sigma }_z^2)$ are jointly independent. Let us consider a special case that $a\cdot b+c=0$. In this setting, the variables $X$ and $Z$ are independent. The direction of $Y\to Z$ would be inverted and the causal model ${𝒢}_1$ is degraded to ${𝒢}_2$. According to the definition, ${𝒢}_1$ and ${𝒢}_2$ satisfy the causal minimality. But the faithfulness is violated in this special case, i.e. ${𝒢}_2$ is not a proper subgraph of ${𝒢}_1$. Thus, the probability of this linear Gaussian model is not faithfulness with respective to ${𝒢}_1$. Although ${𝒢}_2$ is a proper subgraph of $\mathscr{H}$, the distribution does not satisfy causal minimality because the probability is not Markovian with respect to $\mathscr{H}$.

While faithfulness is untestable in practice, it is crucial for deriving valid causal inferences from data because it ensures that the model correctly represents the data-generating mechanisms. If this assumption is violated, the causal relationships are uncertain which is a disaster for causal discovery methods (Spirtes et al., 2000).

Temporal priority. For two variables, temporal priority means that the cause must have occurred before its effect. It is a foundation assumption of causal discovery from temporal data and creates an asymmetric time relationship in causal processes. The temporal priority helps us to establish the direction of a causal relationship when two variables are causally linked. However, if the sampling frequencies of time series are high, the time difference between events associated with the time series may be indistinguishable. In such cases, two events that occurred at different times could be perceived as instantaneous in the observational time series, leading to contemporaneous causal relationships between causes and effects occurring at different time instants.

2.2. Causal Structure for Temporal Data

Figure 4. Causal graphs for discrete-time observations.

For temporal data, the causal relationship can be intuitively defined by the temporal precedence (Eichler, 2012) indicating the causes precede their effects. It reveals the causality asymmetric in time and can be used to orient a causal relation when two variables are known to be causally related. Based on the temporal precedence, there exist three graphical representations of causal structure, i.e., full-time causal graph, window causal graph, and summary causal graph.

As illustrated in Figure 4 (a), the full-time causal graph represents a complete graph of the dynamic systems. For $d$-variate time series $𝐱$, the measurement at each time point $t$ is a vector $(x_1^t,\mathrm{\dots },x_d^t)$. Vertices in full-time causal graphs consist of the set of component $x_1,\mathrm{\dots },x_d$ at each time point $t$ with lag-specific directed links such as $x_i^{t-k}\to x_j^t$. However, it’s usually difficult to discover full-time causal graphs due to the single observation for each series at each time point.

To remedy this problem, window causal graph is proposed. It assumes a time-homogeneous causal structure such that the dynamics of observation vector $𝐱$ are governed by where the function $f$ determines the following observation based on past and the noise ${𝐮}^t$. As illustrated in Figure 4 (b), the window causal graph is represented in a time window, the size of which amounts to the maximum lag in the full-time causal graph.

As shown in Figure 4 (c), each time series component is collapsed into a node to form the summary causal graph. The summary graph represents causal relations between time series without referring to time lags (Peters et al., 2013). In many applications, it is sufficient to model the relations between temporal variables without precisely knowing the interaction between time instants.

For causal discovery from temporal data, most works aim to find the summary causal graph. Nevertheless, summary causal graphs do not always correspond to an SCM, which means they do not enable interventional predictions that are consistent with the underlying time-resolve SCM (Janzing et al., 2018; Rubenstein et al., 2017).

2.3. Problem Definitions

Figure 5. The problem description: (a) causal discovery from MTS data, (b) causal discovery from event sequence data.

As illustrated in Figure 5, causal discovery from temporal data can be divided into two problems, i.e., causal discovery from MTS and causal discovery from event sequence. Next, we formally define them respectively.

Causal Discovery from MTS. Consider a time series with $d$ variables: ${\{{𝐱}^t\}}_{t\in {\mathbb{Z}}^{+}}={\{{(x_1^t{x}_2^t\mathrm{\dots }{x}_d^t)}^{\top }\}}_{t\in {\mathbb{Z}}^{+}}$. Assume that causal relationships between variables are given by the following structural equation model:

$$x_i^t:=f_i(Pa(x_i^t),u_i^t),i=1,\mathrm{\dots },d,$$

where for any $i\in \{1,\mathrm{\dots },d\}$ at time instance $t$, $Pa(x_i^t)$ is the set of direct parents of $x_i^t$ which can be both in the past and at the same time instance. $u_i^t$ denotes the independent noise and can denote either measurement noise or driving noise (Peters et al., 2022) without losing generality. Causal discovery from MTS aims to find either of the two kinds of outputs, i.e., summary causal graph or window causal graph. As for the summary causal graph, the output is the adjacency matrix $A\in {\mathbb{Z}}^{d\times d}$ which summarizes the causal structure, and the $(i,j)$-th entry of the matrix $A$ is $1$ if past observations of $x_i$ enter the structural equation of $x_j^t$ and $0$ otherwise. We say that ‘$x_i$ causes $x_j$’ if $A_{ij}=1$. As for window causal graph with maximum time lag $K$, the output matrices $W$ and $A^k(k\in \{1,\mathrm{\dots },K\})$ correspond to intra-slice and inter-slice edges, respectively. For example, $W_{ij}=1$ denotes the instantaneous dependency $x_i^t\to x_j^t$, while $A_{ij}^k=1$ denotes a lagged dependency $x_i^{t-k}\to x_j^t$ for $k>0$.

Causal Discovery from Event Sequence. For an event sequence: $\{(t_1,e_1),(t_2,e_2),\mathrm{\dots }\}$, $t_i,i=1,2,\mathrm{\dots }$ indicates the time at which the event occurred, while $e_i,i=1,2,\mathrm{\dots }$ stands for the corresponding event type. We aim to discover the causal relationships between different event types. In general, we can construct a causal graph $G=(g_i),i=1,2,\mathrm{\dots },n$, where each node represents a type of event sequence. Our mission is to discover the edge in the causal graph. For example, if there is a directed edge from node $g_j$ to node $g_i$, we say event-type $g_j$ is a cause of event-type $g_i$.

3. Causal Discovery from Multivariate Time Series

In this section, we review causal discovery methods for multivariate time-series data, including constraint-based approaches, score-based approaches, functional causal model-based approaches, Granger causality, and others. The representative algorithms combined with the characteristics are summarized in table 3.

Table 3. Characteristics of causal discovery algorithm reviewed for time-series data, arranged by category.

Section	Method	Causal Graph	Nonlinear	Instantaneous effects	Hidden confounders	Sufficiency Asm.	Markov Asm.	Faithfulness Asm.	Minimality Asm.
Constraint-based	oCSE (2015) (Sun et al., 2015)	Summary	Yes	No	No	Yes	Yes	Yes
	PCGCE (2022) (Assaad et al., 2022a)	Extended	Yes	Yes	No	Yes	Yes	Yes
	PCMCI (2019) (Runge et al., 2019b)	Window	Yes	No	No	Yes	Yes	Yes
	PCMCI+ (2020) (Runge, 2020)	Window	Yes	Yes	No	Yes	Yes	Yes
	ANLTSM (2008) (Chu and Glymour, 2008)	Window	Yes	Yes	Yes	No	Yes	Yes
	tsFCI (2010) (Entner and Hoyer, 2010)	Window	Yes	No	Yes	No	Yes	Yes
	SVAR-FCI (2018) (Malinsky and Spirtes, 2018)	Window	No	Yes	Yes	No	Yes	Yes
	FCIGCE (2022) (Assaad et al., 2022a)	Extended	Yes	Yes	Yes	No	Yes	Yes
	LPCMCI (2020) (Gerhardus and Runge, 2020)	Window	Yes	Yes	Yes	No	Yes	Yes
Score-based	DYNOTEARS (2020) (Pamfil et al., 2020)	Window	No	Yes	No	Yes	Yes	No	No
	NTS-NOTEARS (2021) (Sun et al., 2021)	Window	Yes	Yes	No	Yes	Yes	No	No
	IDYNO (2022) (Gao et al., 2022)	Window	Yes	Yes	No	Yes	Yes	No	No
FCM-Based	VAR-LiNGAM (2008) (Hyvärinen et al., 2008)	Window	No	Yes	No	Yes	Yes	No	Yes
	NCDH (2022) (Wu et al., 2022b)	Summary	Yes	No	No	Yes	Yes	No	Yes
	TiMINo (2013) (Peters et al., 2013)	Summary	Yes	Yes	No	Yes	Yes	No	Yes
	NBCB (2021) (Assaad et al., 2021)	Summary	Yes	Yes	No	Yes	Yes	Yes	Yes
Granger Causality	HSIC-Lasso-GC (2020) (Ren et al., 2020)	Summary	Yes	No	No	No	No	No	No
	(R)NN-GC (2015,2018) (Montalto et al., 2015; Wang et al., 2018)	Summary	Yes	Yes	No	No	No	No	No
	MPIR (2019) (Wu et al., 2020)	Summary	Yes	No	No	No	No	No	No
	NGC (2022) (Tank et al., 2022)	Summary	Yes	No	No	No	No	No	No
	eSRU (2020) (Khanna and Tan, 2020)	Summary	Yes	No	No	No	No	No	No
	SCGL (2019) (Xu et al., 2019)	Summary	Yes	No	No	No	No	No	No
	GVAR (2021) (Marcinkevics and Vogt, 2021)	Summary	Yes	No	No	No	No	No	No
	TCDF (2019) (Nauta et al., 2019)	Window	Yes	Yes	Yes	No	No	No	No
	CR-VAE (2023) (Li et al., 2023)	Summary	Yes	Yes	No	No	No	No	No
	InGRA (2020) (Chu et al., 2020)	Summary	Yes	No	No	No	No	No	No
	ACD (2022) (Löwe et al., 2022)	Summary	Yes	No	Yes	No	No	No	No
Others	DBCL (2010) (Voortman et al., 2010)	Summary	Yes	Yes	Yes	No	Yes	Yes
	NGM (2022) (Bellot et al., 2022)	Summary	Yes	Yes	No	No	No	No	No
	CCM (2012) (Sugihara et al., 2012)	Summary	Yes	No	No	No	No	No	No
	PCTL(c) (2009,2011) (Kleinberg and Mishra, 2009; Kleinberg, 2011)	Summary	Yes	No	No	No	No	No	No

3.1. Constraint-Based Approaches

As a family of causal discovery algorithms, constraint-based approaches rely on statistical tests of conditional independence and are easy to understand and widely used. We first give the main ideas of constraint-based approaches, including general steps and causal assumptions. The detailed methodologies will be categorized into approaches with and without causal sufficiency assumption, and be introduced respectively.

The general steps are: Firstly, it builds a skeleton between variables based on conditional independence. Secondly, it orients the skeleton according to the orientation criterion in the rules. The goal is to construct Completed Partially Directed Acyclic Graphs (CPDAGs) representing the MEC of the true causal diagram. Central to these approaches to derive MEC from observations are the causal assumptions. These methods are usually under the assumptions of causal Markov property and faithfulness, and some also assume causal sufficiency (no unobserved confounders). In this section, we first review the main algorithms and their extensions to time-series data assuming causal sufficiency, then introduce the approaches for conditions when the causal sufficiency assumption is not guaranteed.

3.1.1. Methods with causal sufficiency

In this part, we review methods with causal sufficiency. To reveal the principles of these approaches, we first give a short introduction to methods in the non-temporal setting. Then several popular constraint-based approaches, which originate from the approaches for non-temporal data, for time series are reviewed on the basis of two types of extensions (transfer entropy and momentary conditional independence tests).

As for extracting causal relations from non-temporal data, the Sprites-Glymour-Scheines (SGS) algorithm (Spirtes et al., 1990) is one of the first constraint-based approaches, being proved to be consistent under independently, identically distributed (i.i.d) observations assuming causal sufficiency. However, it suffers from exhausting the test of independence between all nodes. The very large search problem makes it unsuitable in practice. The Peter-Clark (PC) algorithm (Spirtes et al., 2000), which also assumes causal sufficiency, is introduced to reduce unnecessary conditional independence tests and search procedures. Given $d$ non-temporal variables, the detailed procedure of PC algorithm is defined as follows in 3 steps: (1) Firstly, the algorithm starts with a completed undirected graph $G$. (2) Secondly, the algorithm respectively retrieves whether there exist pairs of variables $i$ and $j$ are conditioned on other $n$ variables when $n=0,1,2,\mathrm{\dots },d-2$. If satisfied, remove undirected edges between $i$ and $j$, and update the conditioned variables to the separation set. It proceeds to the pruned skeleton. (3) Finally, it determines the collider (V-structure) to obtain the CPDAG and determines the remaining undirected edges based on other rules.

Although approaches such as SGS and PC are designed in non-temporal settings, constraint-based approaches for time-series data are usually extended from them. We will review recent four popular constraint-based methods, which also assume causal sufficiency, for time-series data. Among these methods, two extensions (Sun et al., 2015; Assaad et al., 2022a) are based on the causal concept of Transfer Entropy, another two (Runge et al., 2019b; Runge, 2020) of them are extended to time series via momentary conditional independence tests.

Extension to time series based on Transfer Entropy. Traditional constraint-based approaches can be extended to the scenario of time series based on the concept of Transfer Entropy. The Transfer Entropy is a model-free measure of temporal causality, of which the definition and variants will be detailed in subsection 3.5.1. Here we view the Transfer Entropy measure as an off-the-shelf part and review two representative approaches from the perspective of constraint-based methodology.

The Optimal Causation Entropy (oCSE) Principle (Sun et al., 2015) is proposed to guide computational and data-efficient causal discovery algorithms from MTS data. It’s based on the theoretical concept of Causation Entropy, a generalization of Transfer Entropy for measuring pairwise relations to network relations of many variables The oCSE method takes a procedure slightly different from that in PC: instead of limiting as much as possible the size of its conditioning set, it conditions since the start on all potential causes which constitute the past of all available nodes. The algorithm is summarized in Algorithm 1, which consists of aggregative discovery of causal nodes, and progressive removal of non-causal nodes. In detail, given node $j$, two procedures are conducted jointly to infer its direct causal neighbors: (1) Firstly, it discovers a superset $Pa(x_j)$ of $j$’s direct causal neighbors aggregately based on the maximization of Causation Entropy. (2) Secondly, it prunes away non-direct causal neighbors based on the Causation Entropy criterion, for example, $i$ is removed from $Pa(x_j)$ if $\mathrm{CE}(x_i^t\to x_j^{t+1}|Pa(x_j^t)\backslash \{x_i^t\})=0$. It’s a computational and sample-efficient algorithm. However, it assumes that the hidden dynamics follow a stationary first-order Markov process as the Causation Entropy only models causal relations with time lags equal to one. Recently, the PCGCE (Assaad et al., 2022a) is proposed to extract extended summary causal graphs for time-series data based on the PC algorithm and the Greedy Causation Entropy, which is a variant of the Causation Entropy.

Input: Multivariate time series $𝐱$ with $d$ dimensions, a significant threshold $\alpha$

Output: The summary causal graph $G$

1 Initialize an empty graph $G$ with $d$ nodes $V$

2 for $j\in \{1,\mathrm{\dots },d\}$ do

3  # Aggregative Discovery of Causal Nodes

4  $z=\mathrm{\infty }$

5  while $z>0$ and card($Pa(x_j)$) do

6   for $x_i\in V\backslash Pa(x_j)$ do

7    Compute the p-value ($z_p$) corresponding to the test $\mathrm{CE}(x_i^t\to x_j^{t+1}|Pa(x_j^t))>0$

8   if $z_p>\alpha$ then add edge $x_i\to x_j$ to $G$

9  # Progressive Removal of Non-Causal Nodes

10  for $x_i\in Pa(x_j)$ do

11   Compute $z$ corresponding to the test $\mathrm{CE}(x_i^t\to x_j^{t+1}|Pa(x_j^t)\backslash \{x_i^t\})=0$

  if $z>\alpha$ then remove edge $x_i\to x_j$ from $G$

Algorithm 1 oCSE

Extension to time series via Momentary Conditional Independence Tests. The PCMCI algorithm (Runge et al., 2019b) leverages a variant of the PC algorithm that flexibly combines linear or nonlinear conditional independence tests and extracts causal relations from time-series data. The goal of the algorithm is to discover the window causal graph. Different from that of PC algorithm, PCMCI starts by constructing a partially connected graph, where all pairs of nodes $(x_i^{t-k},x_j^t)$ are directed as $x_i^{t-k}\to x_j^t$ if $k>0$. This initialization also caters to temporal priority. The algorithm consists of two stages: (1) As done in PC, PCMCI removes all unnecessary edges based on conditional independence. It furthermore removes homologous edges based on the assumption of consistency through time. (2) Momentary Conditional Independence (MCI) is leveraged to deal with autocorrelation, which may lead to spurious correlation. Here, MCI is a measurement, which conditions on the parents of $x_j^t$ and $x_i^{t-k}$ while testing $X_i^{t-k}\cancel{}⟂⟂X_j^t|Pa(X_j^t)\mathrm{\backslash }\{X_i^{t-k}\},X_i^{t-k}$. It also provides an interpretable notion of causal strength from $x_i^{t-k}$ to $x_j^t$. PCMCI has been shown to be consistent and can be flexibly combined with any kind of conditional independence test (linear or nonlinear), such as partial correlation and mutual information. In recent years, there is also a wealth of machine learning approaches on nonparametric tests that address a wide range of independence and dependence types (Zhang et al., 2011; Runge, 2018).

The PCMCI⁺ algorithm (Runge, 2020) extends PCMCI to include the discovery of instantaneous causal relations. Central to the PCMCI⁺ algorithm are two basic ideas that deviate from the origin PC algorithm: First, it conducts the edge removal process separately for lagged and contemporaneous conditioning sets. Second, it leverages MCI to calibrate CI tests under autocorrelation, which is similar to that in PCMCI. The author in (Runge, 2020) also details the curse and blessing of autocorrelation.

3.1.2. Methods without causal sufficiency

Constraint-based approaches without causal sufficiency will be reviewed in this part. In the beginning, we give a brief introduction to the Fast Causal Inference (FCI) Algorithm (Spirtes et al., 2000) for non-temporal data. Then, methods for MTS data consist of two categories: (1) Fast causal inference through time-series models, which is extended from FCI. (2) The methodology via momentary conditional independence tests.

The FCI algorithm is a generalization of the PC algorithm, which can be used in the presence of latent confounders and proven to be asymptotically correct. It utilizes independence tests on the observed data to extract (partial) information on ancestral relationships between the observed variables, thus the goal of the FCI algorithm is to infer the appropriate PAG. The FCI algorithm starts by constructing a complete graph consisting of undirected edges, similar to the PC algorithm. Then iterative conditional independence tests are conducted for the removal of edges. As a result, the FCI algorithm removes edges that are independent, first when conditioning with Sepsets and the with Possible-Dsep sets. For the remaining undirected edges, ten orientation rules are applied recursively. The detailed FCI algorithm, including theoretical analysis, demonstrates the algorithm is sound and complete and can be found in (Zhang, 2008).

Fast Causal Inference Through Time-series Models. A constraint-based method called additive nonlinear time series model (ANLTSM) (Chu and Glymour, 2008) is proposed under the assumption that the effects of hidden confounders are linear and contemporaneous. To escape the curse of dimensionality for nonparametric conditional independence tests, ANLTSM leverages additive regression model, which can be specified as follows :

$$x_j^t=\sum _{1\le i\le d,i\ne j}{a}_{j,i}{x}_i^t+\sum _{1\le i\le d,1\le l\le \tau }{f}_{j,i,l}(x_i^{t-l})+\sum _{r=1}^h{b}_{j,r}{u}_r^t+e^t.$$

Here, $a_{j,i}$ and $b_{j,r}$ are constant values, and $f_{j,i,l}(\cdot )$ denotes the smooth univariate function. The unobserved effects in the form of multi-dimensional Gaussian white noise can be categorized into two types: $e^t$ reflects the latent direct causes of the observed variables, and ${(u_r^t)}_{1\le r\le h}$ denotes latent common causes. And the latent common causes affect the observed variables at the same instant. For $x_i^t$ and $x_j^t$, $u_r^t$ suffices to be stated as a latent common cause if and only if there exists $1\le r\le h$ such that $b_{j,r}{b}_{i,r}\ne 0$. Based on the aforementioned additive regression model, the FCI algorithm is leveraged to identify lagged and instantaneous causal relations. For detecting the instantaneous relations, the conditional independence between $x_i^t$ and $x_j^t$ is first tested given the set $S$ by estimating the conditional expectation $𝔼(x_i^t|x_j^t\cup S)$, then the significance of prediction relationship between $x_i^t$ and $x_j^t$ is checked using statistical tests such as the F-test or the BIC scores, where the insignificance of the predictor implies the conditional independence between $x_i^t$ and $x_j^t$. The lagged causal relations are identified in a similar way. The remaining edges are oriented based on rules. This method is shown to be consistent if the data generation caters to the additive nonlinear time series models. However, the ANLTSM method restricts contemporaneous interactions to be linear, and latent confounders to be linear and contemporaneous.

Another extension of FCI to time-series data is the tsFCI (Entner and Hoyer, 2010) algorithm, where the FCI algorithm is directly applied via a time window. In detail, by assuming the observed time-series data comes from a system at equilibrium, the original time-series data is transformed into a set of samples of the random vector, via a sliding window of size $\tau$. Then considering every component of the transformed vector as a separate random variable, the original FCI algorithm is directly applicable. As the amount of information derived from standard FCI is quite restricted, the temporal priority and time invariance is further incorporated as background knowledge to make more inferences in the orientation phase. However, the tsFCI ignores selection variables and contemporaneous causal relations. Recently, a constrain-based approach named SVAR-FCI (Malinsky and Spirtes, 2018) is proposed that allows for both instantaneous influences and arbitrary latent confounding in the data-generating process. Similar to tsFCI, it also uses time invariance to infer additional edge removals.

Methodology Via Momentary Conditional Independence Tests. It’s found that the original FCI algorithm and its temporal variants suffer from low recall in the autocorrelated time-series case due to the low effect size of conditional independence tests (Gerhardus and Runge, 2020). Some researchers aim to extend PCMCI in the presence of unobserved confounding variables to tackle the aforementioned issues. In (Gerhardus and Runge, 2020), the Latent PCMCI (LPCMCI) algorithm is proposed. Central to the LPCMCI algorithm are two ideas that: First, based on the analysis of the effect size in causal discovery, it uses parents of variables as default conditions and non-ancestors are not tested in the condioning sets, which not only avoids inflated false positives but also reduce the sets to be tested. Second, it introduces the notions of middle marks and LPCMCI-PAGs as an unambiguous causal interpretation to facilitate the early orientation of edges. And the LPCMCI algorithm is proven to be order-independent, sound and complete.

3.2. Score-Based Approaches

Another family of causal discovery approaches is based on score function. The main ideas of score-based approaches will first be introduced, including (dynamic) Bayesian Network, characteristics of score-based approaches compared to their constraint-based counterpart, model scoring, and model search. Then, we will review combinatorial search approaches and continuous optimization approaches for MTS, respectively.

3.2.1. Basics of score-based approaches

The score-based approaches are motivated by the idea that graph structures encoding the wrong (conditional) independence will also result in poor model fit. In the score-based approaches, the causal structure is attached to the concept of Bayesian Network (BN) or Dynamic Bayesian Network (DBN) (Dean and Kanazawa, 1989; Murphy, 2002) dealing with temporal data. In light of this, the score-based methods can generate and probabilistically score multiple models, and then output the most probable one. This contrasts with the constrained-based approaches, which derive and output a single model without quantification regarding how likely it’s to be correct. And the faithfulness assumption is diluted in the scored-based approaches by applying a goodness-of-fit measurement instead of a conditional independence test. The problem of learning a BN or DBN from observations can be therefore formulated as: given a set of instances, find the network that best matches them, i.e., optimize the objective functions. It consists of two elements: model scoring and model search.

Model scoring. Common objective functions fall under two categories: the Bayesian scores which focus on goodness-of-fit and allow the incorporation of prior knowledge, and information-theoretic scores which explicitly consider model complexity, aiming to avoid over-fitting, in addition to the goodness-of-fit (Kitson et al., 2021). The family of Bayesian score functions contains Bayesian Dirichlet equivalent (BDe) score (Heckerman et al., 1995), K2 score (Kayaalp and Cooper, 2013), and so on. The most widely used information-theoretic scores include the Bayesian Information Criterion (BIC) (Neath and Cavanaugh, 2012) and the Akaike Information Criterion (AIC) (Burnham and Anderson, 2004).

Model search / Optimization. The score-based approaches cast the problem of searching causal structure $G$ into an optimization program using the aforementioned score functions $S$. The ultimate goal is therefore stated as (Peters et al., 2017):

$$\widehat{G}={\mathrm{argmin}}_{G\mathrm{over}𝐱}S(D,G),$$

where $D$ represents the empirical data for variables $𝐱$. Traditionally, it’s a combinatoric graph-search problem, and the solution is generally sub-optimal as finding a globally optimal network is known to be NP-hard (Chickering, 1995). A line of works, such as Greedy Equivalence Search (GES) (Chickering, 2002) involve local heuristics owing to the large search space of graphs. However, they still suffer from the curse of dimensionality and suboptimal problems. Recently, an algebraic result characterizing the acyclicity constraint is leveraged in structure learning, which turns the combinatoric problems into a continuously optimizing problem (Zheng et al., 2018, 2020), which can be reformulated as:

where $𝐀$ denotes the adjacency matrix, and $h$ is the function used to enforce acyclicity in the inferred structure. The original acyclicity constraint function is implemented as $h(𝐀)=\mathrm{tr}(e^{𝐀\odot 𝐀})-d$ in (Zheng et al., 2018). It relies on the augmented Lagrangian method (ALM) (Yurkiewicz, 1985) to solve the continuous constrained optimization problem. Various works have further adopted the continuous constrained formulation in neural networks to extract nonlinear causal relations (Zheng et al., 2020; Yu et al., 2019; Gao et al., 2021).

In the context of time series, the ultimate goal of score-based approaches is to learn the structure of DBN. A DBN is a probabilistic network where variables are time series, and it can be decomposed into a prior network and a transition network. A prior network provides dependencies between variables in a given time stamp, and a transition network provides dependencies over time. Therefore, a DBN represents contemporaneous and time-delayed effects in the same framework. Based on this extension to time series, we review the score-based methods following a similar paradigm from combinatoric search to continuous constrained optimization.

3.2.2. Combinatorial search approaches

To conduct the combinatorial search based on scoring function from MTS data efficiently, researches have developed various approaches including structural expectation-maximization (Friedman et al., 1998), cross-validation (Peña et al., 2005), and the decomposition of score functions (de Campos and Ji, 2011).

In (Friedman et al., 1998), the author first utilizes Structural Expectation-Maximization (Structural EM) algorithm (Friedman, 1997, 1998), which is originally a standard algorithm for inferring BN, to learn DBN from longitudinal data. The Structural EM algorithm, combining structural and parametric modification with a single EM process, can be shown to find local optima defined by score functions.

In (Peña et al., 2005), the $K$-fold cross-validation (CV) is leveraged as a computationally feasible scoring criterion for learning DBN. Given the observational data $D$, which is randomly split into $K$ folds $D^1,\mathrm{\dots },D^K$ of approximately equal size, the CV value of a model $G$ is formulated as $\frac{1}{T}{\sum }_{k=1}^K\log p(D^k|G,{\widehat{\theta }}^k)$. And a greedy hill-climbing search is used to estimate $E[\log p(D_{T+1}|G,\widehat{\theta })]$. The procedure starts from the empty graph and updates it gradually by applying the highest scoring single edge additional or removal available. Experiments show that the scoring methods based on cross-validation lead to models generalizing better than those based on BIC of BDe for a wide range of sample sizes.

Based on the score functions that are decomposable, the paper (de Campos and Ji, 2011) uses structural constraints to cast the problem of structure learning in DBN into a corresponding augmented BN, and presents a branch-and-bound algorithm to guarantee global optimality. The decomposed form of the optimal goal can be formalized as:

$$(G^{0\ast },G^{\prime \ast })={\mathrm{argmax}}_{G^0,G^{\prime }}(S_{D_0}(G^0)+S_{D_{1:T}}(G^{\prime }))=({\mathrm{argmax}}_{G^0}{S}_{D_0}(G^0)+{\mathrm{argmax}}_{G^{\prime }}{S}_{D_{1:T}}(G^{\prime })),$$

where $G^0$ and $G^{\prime }$ correspond to the prior network and the transition network respectively. Structural constraints, as a way to reduce the search space, specify where arcs may or may not be included. Because of the branch-and-bound properties, the algorithm can be stopped at the best current solution and an upper bound for the global optimum. The proposed method is shown to be able to handle larger data sets than before, benefiting from the branch-and-bound algorithm and structural constraints.

3.2.3. Continuous optimization approaches

Owing to the recent contribution of NOTEARS (Zheng et al., 2018), the score-based learning of DAGs can be reformulated as a continuous constrained optimization problem, which inspires various works (Zheng et al., 2020; Yu et al., 2019; Gao et al., 2021; Ng et al., 2022b, 2019; Lachapelle et al., 2020) in structure learning. At the heart of this line of the method is an algebraic characterization of acyclicity expressed as a constraint function, which is further leveraged to minimize the least square loss while enforcing acyclicity. In the context of time series, some works have also adopted this continuous constrained formulation to support structure learning and causal discovery (Pamfil et al., 2020; Sun et al., 2021; Hsieh et al., 2021; Gao et al., 2022).

DYNOTEARS, introduced in (Pamfil et al., 2020), captures linear relations from time-series data via a continuous optimization approach. It models the data in the following standard SVAR way:

$${𝐱}^t={𝐱}^t𝐖+{𝐱}^{t-1}{𝐀}^1+\mathrm{\dots }+{𝐱}^{t-p}{𝐀}^p+{𝐮}^t,$$

where $p$ is the order of SVAR model, $𝐮$ is a vector of centered error variables. To guarantee the identifiability in SVAR models, the error terms ${𝐞}^t$ are assumed either non-Gaussian or standard Gaussian, i.e., ${𝐮}^t\sim 𝒩(0,I)$, as the identifiability is proven to hold on the two cases (Hyvärinen et al., 2010b; Peters et al., 2017). $𝐖$ and $𝐀$ are weighted adjacency matrices, which correspond to intra-slice edges (contemporaneous relationship) and inter-slice edges (time-lagged relationship), respectively. The SEM can further takes the compact form:${𝐗}^t={𝐗}^t𝐖+{𝐗}^{(t-p):(t-1)}𝐀+𝐔$. The procedure of structure learning revolves around minimizing the least-squares loss subject to an acyclicity constraint, which gives the following optimization problem:

${\min }_{𝐖,𝐀}f(𝐖,𝐀)\mathrm{s}.\mathrm{t}.$ $𝐖\mathrm{is}\mathrm{acyclic},$

$\mathrm{where}f(𝐖,𝐀)={\displaystyle \frac{1}{2n}}||{𝐗}^t-{𝐗}^t𝐖-{𝐗}^{(t-p):(t-1)}𝐀$ $|{|}_F^2+{\lambda }_{𝐖}||𝐖|{|}_1+{\lambda }_{𝐀}||𝐀|{|}_1.$

To sidestep the key difficulty of solving the optimization problem under the acyclicity constraint, DYNOTEARS follow the work in (Zheng et al., 2018), where the trace exponential function $h(𝐖)=\mathrm{tr}(e^{𝐌\odot 𝐌})-d$ is leveraged as an equivalent formulation of acyclicity. The continuous constrained optimization problem is translated via the augmented Lagrangian method into unconstrained problems of the form:

$${\min }_{𝐖,𝐀}F(𝐖,𝐀),\mathrm{where}F(𝐖,𝐀)=f(𝐖,𝐀)+\frac{\rho }{2}h{(𝐖)}^2+\alpha h(𝐖).$$

Towards the optimization of the above smooth augmented objective, two solving approaches are presented separately. The first approach is to use standard solvers such as L-BFGS-B (Zhu et al., 1997). An alternative approach is a two-stage procedure similar to those in (Hyvärinen et al., 2010b), where we can rewrite the equation as $𝐳=\mathrm{𝐳𝐖}+𝐔$ and derive the estimate of $𝐖$ by using static NOTEARS to the error term $𝐳$.

NTS-NOTEARS (Sun et al., 2021) is a recent advance that adopts the continuous constrained formulation. Compared to DYNOTEARS, which is a linear autoregressive model, NTS-NOTEARS is able to extract both linear and non-linear relations among variables. It achieves this by leveraging 1D convolutional neural networks (CNNs), which exploit a sequential topology in the input data and are thus well-suited neural function approximation models for temporal data. $d$ CNNs, each of which the first layer is a 1D convolutional layer with $m$ kernels, are trained jointly where the $j$-th CNN predicts the expectation of targeted variable $x_t^j$ at the specific time $t$ given preceding and contemporaneous input variables. Each CNN can be viewed as a Markov blanket of the target variable. The dependence of child variables on their parents in DBN is given as follows:

$$𝔼[x_j^t|Pa(x_j^t)]={\mathrm{CNN}}_j(\{{𝐱}^{t-k}:1\le k\le K\},{𝐱}_{-j}^t),$$

where parents $Pa(x_j^t)$ are derived from the trained CNNs, and ${𝐱}_{-j}^t$ denotes all variables at time step $t$ except $x_j$. In light of NOTEARS-MLP (Zheng et al., 2020) (a non-linear and NN-based extension of NOTEARS), the dependency strength of an edge in DBN is estimated in the following way:

$$W_{ij}^k={\Vert {\varphi }_{i,j}^k\Vert }_L^2\mathrm{for}k=1,\mathrm{\dots },K+1.$$

In detail, the $x_i^{t-k}$ belongs to the parent set $Pa(x_j^t)$ on the condition that the estimated dependency strength is larger than threshold weight $W_{ij}^k>W_{thres}^k$. The optimization procedure follows a similar way as DYNOTEARS. It’s also worth noticing that NTS-NOTEARS shows prior knowledge of variable dependencies that can be transformed as additional optimization constraints and incorporated into the L-BFGS-B solver.

To handle both observational and interventional data, an algorithm, called IDYNO (Gao et al., 2022), is proposed recently. It first introduces a non-linear objective through neural networks to model complex dynamics, then modifies an objective and general solution approach to handle different distributions on intervention targets.

We can find that it’s a powerful methodology for score-based structure learning to use continuous optimization and avoid the explicit combinatoric traversal of possible causal structures. The past several years have also witnessed numerous applications and extensions of this methodology. However, some boundaries and limitations are further discussed in (Kaiser and Sipos, 2022; Reisach et al., 2021; Ng et al., 2022a), including the influence of data scale and the convergence condition of the augmented Lagrangian method. We recommend you take these issues into consideration for further developments and applications of this family of methods.

3.3. FCM-Based Approaches

The two families of methods above either face the inseparability of the MEC or the need for large samples to confirm causal faithfulness. Causal discovery can also be conducted based on Functional Causal Models (FCM) (Pearl et al., 2000), which is also known as SCM in 2.1 and describes a causal system via a set of equations. Recent years have witnessed the proliferation of FCM-based approaches for both temporal and non-temporal data. In this subsection, we first introduce the main ideas of FCM-based approaches, including the functional causal model and the usage of noise in orienting causal relations. Then two families of FCM-based approaches, i.e., methods using independent component analysis and additive noise model, will be reviewed, respectively.

In FCM, each variable is explained by an equation in terms of its direct causes and some additional noise. For example, the function $x_j=f_j(x_i,u_j)$ explains the causal link $x_i\to x_j$ with some additional noise $u_j$. One basic idea of the FCM-based causal discovery approaches is that statistical noise can be a valuable source of insight, which caters to recent discoveries (Climenhaga et al., 2021) challenging the orthodoxy that the noise should be treated as a nuisance. To be specific, causal relationships can be identified and estimated with the help of noise.

3.3.1. Methods using independent component analysis

In this part, we first introduce the basic idea of this family of methods by reviewing the original algorithm in non-temporal setting (Shimizu et al., 2006). Then, methods for MTS data will be detailed (Hyvärinen et al., 2008, 2010a; Schaechtle et al., 2013; Wu et al., 2022b).

LiNGAM (Shimizu et al., 2006) is a typical FCM-based causal discovery algorithm in non-temporal setting, and has the following assumptions: (1) a linear data generation process, (2) non-Gaussian disturbances, (3) no unobserved confounders. In the LiNGAM model, the relations among observations can be formulated as $𝐱=\mathrm{𝐁𝐱}+𝐮$, where $𝐱,𝐁,𝐮$ respectively denote the vector of variables, the adjacency matrix of the causal graph and the noise vector. The equation can be rewritten as $𝐱=\mathrm{𝐀𝐮}$, where $𝐀={(𝐈-𝐁)}^{-1}$. For the equation, the independent component analysis (ICA) method (Stone, 2004) can be used to estimate $𝐀$, and causal relationships $𝐁$ can be derived. Along this line, DirectLiNGAM (Shimizu et al., 2011) further leverages the regression model to ensure the original models to converge to the correct solution in a controlled number of steps. Extensions of LiNGAM to time series are as follows.

As a temporal extension of LiNGAM, VAR-LiNGAM (Hyvärinen et al., 2008, 2010a) estimates the structural autoregressive (SVAR) models by leveraging non-Gaussianity property. SVAR models reflect both instantaneous and time-delayed causal effects and are among the most prevalent tools in empirical economics to analyze dynamic phenomena (Moneta et al., 2013). In VAR-LiNGAM, a representation of time series is a combination of SVAR and SEM, which is defined as:

(SVAR) $${𝐱}^t=\sum _{k=0}^{\tau }{𝐁}^k{𝐱}^{t-k}+{𝐮}^t$$

where ${𝐁}^k$ is the $n\times n$ matrix of the causal effects between the variables $𝐱$ with time lag $k$. And ${𝐮}^t$ are random processes modeling the external influences or ‘disturbances’, which are assumed to be independent, temporally uncorrelated and non-Gaussian. To estimate the above model, a classic least-squares estimation of the autoregressive (AR) model (time lag $k>0$) is combined, which is formalized as:

(VAR) $${𝐱}^t=\sum _{k=1}^{\tau }{𝐌}^k{𝐱}^{t-k}+{𝐧}^t$$

Based on the SVAR and VAR formalization, the basic idea of VAR-LiNGAM is that we can estimate ${𝐌}^k$ of VAR model in a classic least-square fashion consistently and efficiently. And we can deduce the estimate of instantaneous causal effect through LiNGAM analysis. As for the time-delayed effect, it can be derived from reparametrization. The ensuing method in detail is defined as follows in four steps: (1) Firstly, fit the regressions and denote the least-squares estimates of the AR matrices by ${\widehat{𝐌}}^k$. (2) Secondly, compute the residuals, i.e., ${\widehat{𝐧}}^t={𝐱}^t-{\sum }_{k=1}^{\tau }{\widehat{𝐌}}^k{𝐱}^{t-k}$. (3) Thirdly, perform LiNGAM analysis (Shimizu et al., 2006) based on the equation ${\widehat{𝐧}}^t={𝐁}^0{\widehat{𝐧}}^t+{𝐞}^t$ to derive the estimate of instantaneous causal effect ${\widehat{𝐁}}^0$. (4) Finally, compute the estimates of the time-delayed causal effect ${\widehat{𝐁}}^k(k>0)$ as ${\widehat{𝐁}}^k=(𝐈-{\widehat{𝐁}}^0){\widehat{𝐌}}^k$. The VAR-LiNGAM model degenerates to the LiNGAM model if the order of the autoregressive part is set to zero, i.e., $\tau =0$. And an intensive application of this approach in empirical economics can be found in (Moneta et al., 2013).

The VAR-LiNGAM is extended to the identification and estimate of causal models under time-varying situations (Huang et al., 2015), where Gaussian Process regression is further leveraged to automatically model how the causal model change over time. In (Lanne et al., 2017), the initial VAR-LiNGAM is generalized to the condition where the inferred graphs can contain cycles. And the proposed model is demonstrated theoretically to be identifiable. Another algorithm based on LiNGAM, called the Multi-Dimensional Causal Discovery (MCD), is proposed in (Schaechtle et al., 2013). MCD can efficiently discover causal dependencies in multi-dimensional settings, such as time-series data, by integrating data decomposition and projection.

To get rid of constraints of linear (Hyvärinen et al., 2008, 2010a) or additive assumptions(Peters et al., 2013), an FCM-based algorithm named Nonlinear Causal Discovery via HM-NICA (NCDH) is recently proposed in (Wu et al., 2022b) to extract general nonlinear relations from time series. At the heart of this algorithm, a nonlinear ICA algorithm is leveraged as a measurement of nonlinear relationships. The observations are assumed to be generated by mutually independent latent components:

$$𝐱=𝐟(𝐒)\mathrm{where}𝐟={(f_1,f_2,\mathrm{\dots },f_d)}^T\mathrm{and}𝐒={(S_1,S_2,\mathrm{\dots },S_d)}^T.$$

Similar to that in linear ICA, $𝐒$ contains components that are independent of each other, and the goal of nonlinear ICA is to recover $𝐒$ from $𝐱$. NCDH first leverages the nonlinear ICA combined with HMM (Hälvä and Hyvärinen, 2020) to separate latent noises. As a remedy for the permutation uncertainty of ICA, a series of independence tests are conducted to determine the corresponding relations between the observed variables and the separated noises. A recursive search algorithm is finally taken to extract the causal relations.

3.3.2. Methods using additive noise model

In reality, there’re many non-linear causal relationships that violate the assumption of LiNGAM family methods. Despite recent advances (such as NCDH) extracting causal relations in general nonlinear conditions, their usages are restricted. Another family of FCM-based approaches is based on the additive noise model (ANM) with nonlinear function, which is suitable in more general settings. In this part, the main ideas of methods using ANM will be given firstly. Then we will introduce the detailed methods for MTS data.

It’s demonstrated in (Hoyer et al., 2008) that the true causal structure can be identified in the ANM with nonlinear functions if the causal minimality condition holds. In ANM, if $x_i\to x_j$, we have $x_j=f(x_i)+u_j$, and the cause $x_i$ and additive noise $u_j$ are independent. If the noise $u$ is subject to non-Gaussian distribution and $f(\cdot )$ is a linear function. In the bivariate case $x_i\to x_j$, we can fit regression models in causal and anti-causal directions, the true orientation can be inferred by testing the independence with residuals. As for the multivariate case, a pairwise strategy can be adopted (Mooij et al., 2009). The correctness of this algorithm is discussed in (Peters et al., 2014).

In (Peters et al., 2013), the Time Series Models with Independent Noise (TiMINo) is proposed, which is a causal discovery method for time series based on ANM. It inputs time-series data and outputs either a summary time graph or remains undecided, which avoids leading to wrong causal conclusions when the model is mis-specified or the data is insufficient. It leverages a similar method as that in non-temporal and multivariate setting (Mooij et al., 2009). In detail, it tries to fit the structural equation models for time series, which can be formulated as follows:

$$x_j^t=f_j(Pa{(x_j^{\tau })}^{t-\tau },\mathrm{\dots },Pa{(x_j^1)}^{t-1},Pa{(x_j^0)}^t,u_j^t),$$

where error terms $u_j^t$ are jointly independent over variable index $j$ and time index $t$. There are several options for fitting methods $f$ such as linear models, generalized additive models, and Gaussian process regression models. For inferring causal relations in the additive noise model, independence tests such as cross-correlations and HSIC (Gretton et al., 2007) can be leveraged.

There are some drawbacks to those functional causal models, such as VAR-LiNGAM and TiMINo. It’s illustrated that those methods are not well scalable across the increase of node numbers (Glymour et al., 2019), and those performances are not promising without a large sample size (Malinsky and Danks, 2018). To overcome those drawbacks, a Noise-Based / Constraint-Based (NBCB) approach is proposed in (Assaad et al., 2021), where the constraint-based approach is further leveraged based on the original additive noise model for time-series data. In detail, the potential causes of each time series are detected by an additive noise model which is similar to that in TiMINo. Unnecessary causal relations are pruned using temporal causal entropy, which is an extension to causation entropy (Sun et al., 2015) measuring the (conditional) dependencies between two-time series.

3.4. Granger Causality Based Approaches

Granger causality is a popular tool for analyzing time-series data in many real-world applications. There exist many causal discovery approaches developed on the basis of Granger causality. In this subsection, we first introduce definitions of Granger causality. Before delving into detailed methods, two categories of Granger causality models for MTS (model-free and model-based) will be given and compared. Due to the superiority of model-based approaches in more general conditions, the rest of this part will focus on two recent advances in model-based approaches: (1) methods based on kernels (3.4.3), and (2) methods based on neural networks (3.4.4).

3.4.1. Basics of Granger causality

Granger causality analysis, which is first proposed in (Granger, 1969), is a powerful method that determines cause and effect based on predictability. A time series $x_i$ Granger-causes $x_j$ if past values of $x_i$ provide unique, statistically significant information about future values of $x_j$. According to this proposition, $x_i$ is defined to be ‘causal’ for $x_j$ if

where denotes the optimal prediction of $x_j^t$ given the history of all relevant information . Here indicates excluding the information of from . The above definition seems general and does not have specific modeling assumptions, whereas there are also various forms of definition for Granger causality based on different model specifications and statistical tools for better representation power and the convenience of inference, such as autoregression model (in Granger’s original paper (Granger, 1969)) and so on. And if all relevant variables are observed and no instantaneous connections exist, Granger causal relations are equivalent to causal relations in underlying DAGs (Peters et al., 2013, 2017).

3.4.2. Early approaches for MTS

Earlier methods for identifying Granger causality were limited to bivariate settings. Specifically, a well-documented (Lütkepohl, 1982) issues for Granger causal analysis in bivariate settings is that the causal findings may be misleading without adjusting for all relevant covariates. On the one hand, it’s necessary to account for more variables to prevent identifying incorrect Granger causal relations (Shojaie and Fox, 2021). On the other hand, MTS widely exist among various fields. Inferring Granger causal relations in MTS, which is also termed graphical Granger causality or network Granger causality in some literature, has become a hot research topic. Various graphical Granger causal analysis models for MTS can be divided into two categories, namely model-free and model-based approaches.

Model-free Methods. The mainstream of model-free approaches for multivariate Granger causality are based on predictability and need to estimate the conditional probability density functions (CPDFs) (Bai et al., 2010). In (Diks and Wolski, 2016), the estimates of the CPDFs are provided, and the the bivariate Diks-Panchenko nonparametric causality test is extended to the multivariate case. By introducing conditional variables into the marginal probability density functions, the copula-based Granger causality model (Hu and Liang, 2014; Kim et al., 2020) can also be extended to multivariate case. Besides, model-free measures such as transfer entropy and directed information (Amblard and Michel, 2011), are able to detect nonlinear dependencies. The definitions and some properties of these model-free estimators will be detailed in 3.5.1. Model-free methods can deal with nonlinear Granger causal relations well. However, these estimators suffer from high variance and require large amounts of data for reliable estimation, and also suffer from curse of dimensionality when the number of variables grows. Thus, in the complex real-word scenarios where is nonlinear and high-dimensional, the utilization of model-free methods to some extend are limited.

Model-based Methods. In contrast to model-free counterparts, model-based methods are computationally efficient and therefore more suitable for inferring Granger causal relations in high-dimensional conditions. The model-based inference approach is adopted by the vast majority of Granger causal whereby the measured time series is modeled by a suited parameterized data generative model. And the inferred parameters ultimately reveal the true topology of Granger causality. Earlier methods along this line are typically using the popular vector autoregressive (VAR) model under the assumption of linear time-series dynamics. For $d$-variate time series $𝐱$, the VAR model is defined as:

$${𝐱}^t=\sum _{k=1}^{\tau }{A}^k{𝐱}^{t-k}+{𝐮}^t,$$

where $A^k$ is a $d\times d$ matrix that specifies how lag $k$ affects the future evolution of the series and ${𝐮}^t$ denotes zero mean noises. In the VAR model, as a straightforward extension form the bivariate case (Granger, 1969), time series $i$ does not Granger-cause time series $j$ if and only if for all time lag $k$, the component $(j,i)$ of $A^k$ equals zero. Thus the Granger causal analysis reduces to determine which entries in $A^k$ are zero over all lags. There are also abundant research works (Arnold et al., 2007; Lozano et al., 2009a; Shojaie and Michailidis, 2010; Basu et al., 2015) reducing the computational complexity via the Lasso penalty and its variants for Granger causal analysis in high-dimensional time series, which are also termed as Lasso Granger causality (Lasso-GC). For these methods, the problem of Granger causal series selection can be generally formulated as follows based on least square loss:

$${\min }_{A^1,\mathrm{\dots },A^{\tau }\in {ℝ}^{d\times d}}\sum _{t=\tau +1}^T{\Vert {𝐱}^t-\sum _{k=1}^{\tau }{A}^k{𝐱}^{t-k}\Vert }_2^2+\lambda R(𝐀),$$

where $R(\cdot )$ is the sparsity-inducing regularizer and has various implementations as shown in table 4. Different penalty terms induce different sparsity patterns in $A^1,\mathrm{\dots },A^{\tau }$, thus inducing different heuristics and constraints in the Granger causal series selection. Except for Lasso-GC, another line of works based on VAR models in multivariate setting worth mentioning is the conditional Granger causality index (CGCI) (Geweke, 1982). For variable $X,Y$ and conditional variables $Z$, by comparing residuals errors of the reduced and full models ${\mathrm{CGCI}}_{X\to Y|Z}=\ln \frac{\mathrm{var}({ϵ}_{Y|Z})}{\mathrm{var}({ϵ}_{Y|XZ})}$, a distinction between direct and indirect causality in multivariate systems can be made based on CGCI. Along this line of works, mBTS-CGCI is proposed in (Siggiridou and Kugiumtzis, 2016) based on a modified backward-in-time selection (mBTS) to limit the order of VAR models, thus can be better applied to high-dimensional scenarios.

Table 4. Common sparsity-inducing penalty terms, described by (Nicholson et al., 2017; Marcinkevics and Vogt, 2021)

Model Structure	Penalty Function
Basic Lasso	${\Vert 𝐀\Vert }_1$
Elastic net	$\alpha {\Vert 𝐀\Vert }_1+(1-\alpha ){\Vert 𝐀\Vert }_2^2,\alpha \in (0,1)$
Lag group Lasso	${\sum }_{k=1}^{\tau }{\Vert {𝐀}^k\Vert }_F$
Component-wise Lasso	${\sum }_{p=1}^d{\sum }_{k=1}^{\tau }{\Vert {({𝐀}^{k:\tau })}_p\Vert }_2$
Element-wise Lasso	${\sum }_{p=1}^d{\sum }_{q=1}^d{\sum }_{k=1}^{\tau }{\Vert {({𝐀}^{k:\tau })}_{p,q}\Vert }_2$
Lag-weighted Lasso	${\sum }_{k=1}^{\tau }{k}^{\alpha }{\Vert {𝐀}^k\Vert }_1,\alpha \in (0,1)$

Although the model-based approaches, compared to the model-free counterparts, take advantage of efficiently processing high-dimensional time series, the fundamental issue of these approaches is the model misspecification. Especially, the notion of multivariate Granger causality based on the vanilla VAR model assumes time series follows linear dynamics, whereas many interactions in real-world applications are inherently nonlinear. Recently, many model-based approaches, which are compatible with nonlinear causal relations, have emerged and can be grouped into two categories: methods based on kernel, and methods based on neural networks. As the generation on Granger causality, fundamental venation and development orientations have been reviewed above and prospected from classic documents. In the following part of this subsection, due to their ability to be leveraged in complex real-world scenarios, we will detail the recent advances of model-based methods in nonlinear and high-dimensional settings, especially new perspectives from neural networks.

3.4.3. Recent advances based on kernels

To extract nonlinear causal relations in a model-based approach, establishing a nonlinear parameter model is a common strategy. A line of works extend Granger causality to kernel methods (Ancona et al., 2004; Marinazzo et al., 2008b, a; Sindhwani et al., 2013; Ren et al., 2020). In (Ancona et al., 2004), Granger causality is extended to bivariate nonlinear cases by means of radial basis functions. Furthermore, a Granger causality analysis model is put forward (Marinazzo et al., 2008b) based on the theory of reproducing kernel Hilbert spaces (RKHS). The key idea is to embed data into a Hilbert space and search for nonlinear relations in that space. This method is then generalized to the multivariate case in (Marinazzo et al., 2008a). In (Sindhwani et al., 2013), a matrix-valued extension of the kernel method is proposed, imposed on a dictionary of vector-valued RKHS. The algorithm is for high-dimensional nonlinear multivariate regression, and can naturally lead to nonlinear generalization of graphical Granger Causality. Recently, an algorithm based on Hilbert-Schmidt independence criterion Lasso Granger causality (HSIC-Lasso-GC) (Ren et al., 2020) is proposed.

3.4.4. Recent advances based on neural networks

Neural networks are able to represent nonlinear, complex, and non-additive interactions between variables. In this part, recent advances of Granger causal methods based on neural networks will be reviewed, including non-uniform embedding (Montalto et al., 2015; Wang et al., 2018), information regularization (Wu et al., 2020), component-wise neural network modeling (Tank et al., 2017, 2022; Khanna and Tan, 2020), low-rank approximation (Xu et al., 2019), self-explaining networks (Marcinkevics and Vogt, 2021), attention mechanisms (Nauta et al., 2019; Schwab et al., 2019), recurrent variational autoencoders (Li et al., 2023), and inductive modeling (Chu et al., 2020; Löwe et al., 2022). Besides, as illustrated in Fig 6, existing NN-based Granger causality approaches can be categorized into four groups: parameter-based (Tank et al., 2022; Khanna and Tan, 2020), attention-based (Nauta et al., 2019; Chu et al., 2020), self explanation-based (Marcinkevics and Vogt, 2021), and relational encoding-based (Löwe et al., 2022).

Figure 6. Four types of NN-based Granger causal discovery methods according to the means to derive causal weights.

DL-extensions with Non-uniform Embedding. A feature selection procedure termed as a non-uniform embedding (NUE) is proposed in NN-GC (Montalto et al., 2015) to identify the significant Granger causes in the MLP model. By greedily adding lagged components of predictor time series as input, an MLP is updated iteratively. A predictor time series is claimed a significant Granger cause of the target time series if at least one of its lagged components is added when the procedure is terminated. In RNN-GC (Wang et al., 2018), the NUE is extended by replacing MLPs with gated RNN models, However, as this technique requires training and comparing many candidate models, it’s costly in high-dimensional settings.

DL-extensions with Information Regularization. For extracting nonlinear dynamics, a method with Minimum Predictive Information Regularization (MPIR) (Wu et al., 2020) is introduced. It leverages learnable corruption for predictor variables and minimizes a mutual information-regularized risk, which combines the benefits of the Granger causality paradigm with deep learning models. In MPIR, the author states that the naive way to combine neural nets with Granger causality suffers from two major drawbacks: instability and inefficiency. The solution is to encourage each ${𝐱}_i^{t-K:t-1}$ to provide as little information to $x_j^t$ as possible while maintaining good prediction via learned corruption, replacing the naive way which predicts $x_j^t$ with one ${𝐱}_i^{t-K:t-1}$ missing at a time. The risk is defined as follows:

$$R_{𝐗,x_j}[f_{\theta },𝐧]=E_{{𝐗}^{t-1},x_j^t,𝐮}[{(x_j^t-f_{\theta }({\widetilde{𝐗}}_{(𝐧)}^{t-K:t-1}))}^2]+\lambda \cdot \sum _{p=1}^{d}I({\tilde{X}}_{i(n)}^{t-K:t-1};X_i^{t-K:t-1}),$$

where ${\widetilde{𝐗}}_{(𝐧)}^{t-K:t-1}:={𝐗}^{t-K:t-1}+𝐧\odot 𝐞$ (or its element-wise representation, ${\tilde{X}}_{i(n)}^{t-K:t-1}:=X_i^{t-K:t-1}+n_i\cdot e_i,i=1,2,\mathrm{\dots },d$) are the noise-corrupted inputs with learnable noise amplitudes $𝐧$ and $u_j\sim N(\mathrm{𝟎},𝐈)$. And $W_{pq}=I({\tilde{X}}_{i(n^{\ast })}^{t-K:t-1};X_i^{t-K:t-1})$ is the minimum predictive information at the minimization of $R_{𝐗,x_j}[f_{\theta },𝐧]$, which contains causal information and measures the predictive strength of variable $i$ for predicting variable $j$, conditioned on all the other observed variables. To be specific, $W_{ij}=0$ if $x_i⟂⟂x_j$. Besides, as it’s inefficient to estimate the mutual information term with a large dimension, an upper bound is derived as an alternative optimization goal. Instead of training many candidate models and suffering from instability and inefficiency, this framework only requires training $d$ models separately.

DL-extensions with Component-wise NN Modeling. Another NN-based approach to measure nonlinear Granger causality is component-wise modeling. A component-wise framework is proposed in (Tank et al., 2017), which can be viewed as a generalization of the linear VAR model. In detail, the generation procedure of each variable can be written as follows:

$${𝐱}_j^t:=g_j({𝐱}_1^{1:(t-1)},\mathrm{\dots },{𝐱}_i^{1:(t-1)},\mathrm{\dots },{𝐱}_d^{1:(t-1)})+u_j^t,\mathrm{for}1\le j\le d,$$

where $g_j(\cdot )$ is a continuous function, based on regularized neural networks implementation, specifying how the past values of $𝐱$ determine the future values of variable ${𝐱}_j$. In this context, the time series ${𝐱}_i$ is Granger non-causal for time series ${𝐱}_j$ (${𝐱}_i\nrightarrow {𝐱}_j$) if and only if $g_j(\cdot )$ is invariant to ${𝐱}_i^{1:(t-1)}$, which can be defined as:

$$g_j({𝐱}_1^{1:(t-1)},\mathrm{\dots },{𝐱}_i^{1:(t-1)},\mathrm{\dots },{𝐱}_d^{1:(t-1)})=g_j({𝐱}_1^{1:(t-1)},\mathrm{\dots },{𝐱}_{i^{\prime }}^{1:(t-1)},\mathrm{\dots },{𝐱}_d^{1:(t-1)}),$$

for all $({𝐱}_1^{1:(t-1)},\mathrm{\dots },{𝐱}_d^{1:(t-1)})$ and all ${𝐱}_i^{1:(t-1)}\ne {𝐱}_{i^{\prime }}^{1:(t-1)}$. We will introduce two methods (Tank et al., 2022; Khanna and Tan, 2020) based on this framework, respectively.

Neural Granger Causality (NGC) is proposed in (Tank et al., 2022) to infer nonlinear Granger causality using structured MLP and LSTM with sparse input layer weights, which are termed as component-wise MLP (cMLP) and component-wise LSTM (cLSTM), respectively. In the cMLP, each nonlinear output $g_j$ is modeled with a separate MLP as to easily disentangle the effects from inputs to outputs. The input matrix of the first layer provides information for penalized selection of Granger causality. To be specific, in the first layer of $g_j(\cdot )$

$$h_1^t=\sigma (\sum _{k=1}^{\tau }{W}_1^k{𝐱}^{t-k}+b_1),$$

if the $i$-th column of weight matrix $W_1^k$ contains zeros for all time lag $k$, then time series $i$ does not Granger-cause series $j$. Analogously to the VAR type methods, the Granger causal series are selected by the following encoding selection(Tank et al., 2017) procedure:

$${\min }_{𝐖}\sum _{t=\tau }^T(x_j^t-g_j(x_{(t-1):(t-\tau )}))+\lambda \sum _{i=1}^{d}R({(W_1)}_{:i}),$$

where sparse inducing penalty $R(\cdot )$ is implemented through group lasso penalty, which extracts causal relations without requiring precise lag specification. As for the cLSTM, it sidesteps the lag selection problem and the Granger causal information can also be easily interpreted in the vanilla LSTM model. The input matrix, which is slightly different from that in MLP, is defined as $W^1={({(W^f)}^{\top },{(W^{in})}^{\top },{(W^o)}^{\top },{(W^c)}^{\top })}^{\top }$, controlling how the past time series affect the forgot gates, input gates, output gates, and cell updates. Granger-causal series can be selected based on a group lasso penalty across columns of $W^1$. In the end, to optimize the non-convex optimization objectives in either cMLP or cLSTM, the proximal gradient descent(Parikh et al., 2014) is used, which leads to the exact zeros in the input matrix. This property in the optimization procedure meets the requirement for interpreting Granger non-causality in the framework. To infer the network topology of Granger causality, $d$ models need to be trained with each variable as a response.

Another sample-efficient architecture economy-SRU (eSRU) is proposed in (Khanna and Tan, 2020). It leverages Statistical Recurrent Units (SRUs) (Oliva et al., 2017) to model the observed time-series data. Here SRUs are a special type of RNNs designed for MTS with time-delayed and nonlinear dependencies and therefore also suited for extracting the network topology of nonlinear Granger causal relations. To be specific, it suffers less from the vanishing and exploding gradient issues owing to an ungated architecture and is able to model both short and long-term temporal dependency among multivariate time series by maintaining multi-time scale summary statistics. Similar to model-based approaches like cLSTM, the measure of Granger causal relationships can be derived from the input-layer weight parameters of the SRUs. However, due to the common issue of data scarcity in the causal inference problem, the original framework suffers from overfitting. Additionally, two modifications are implemented as a remedy for overfitting in eSRU.

DL-extensions with Low-rank Approximation. The scalable causal graph learning (SCGL) framework is proposed in (Xu et al., 2019). The authors first deconstruct data nonlinearity into two types (i.e. univariate-level and multivariate-level nonlinearity), which are modeled separately. The key idea of SCGL is that learning the full size of the adjacency matrix $A\in {ℝ}^{d\times d}$ would be unscalable when the size of variables $d$ is quite large. In practice, the relationship of variables is low-rank in hidden space (Zorzi and Chiuso, 2017; Chiuso and Pillonetto, 2012). Therefore, it’s natural to approximate $A$ via a $k$-rank decomposition, where . The low-rank approximation reduces the noise influence in causal discovery and provides interpretability in downstream time series analysis (Huang et al., 2020c).

DL-extensions with Self-explaining Networks. For better interpretability, the generalized vector autoregression (GVAR) model (Marcinkevics and Vogt, 2021) is proposed. It’s based on an extension of self-explaining neural networks (Alvarez-Melis and Jaakkola, 2018). The self-explaining neural networks are inherently interpretable models motivated by restricted properties, and follow the form:

$$f(𝐱)=g(\theta {(𝐱)}_{1}h{(𝐱)}_1,\mathrm{\dots },\theta {(𝐱)}_{k}h{(𝐱)}_k),$$

where $g(\cdot )$ and $𝐡(𝐱)$ denote a link function and the interpretable basis concepts, respectively. Combined with the vector autoregression model, which is often specified in Granger causal inference, the GVAR model is given by

$${𝐱}^t=\sum _{l=1}^{\tau }{\mathrm{\Psi }}_{{\theta }_l}({𝐱}^{t-l}){𝐱}^{t-l}+{𝐮}^t,$$

where ${\mathrm{\Psi }}_{{\theta }_l}:{ℝ}^d\to {ℝ}^{d\times d}$ is a neural network parameterized by ${\theta }_l$, of which the output is the matrix corresponding to the strength of influence. In detail, the strength of influence $x_i^{t-l}\to x_j^t$ is measured by the component $(j,i)$ of ${\mathrm{\Psi }}_{{\theta }_l}({𝐱}^{t-l})$. The loss function consists of three terms: the MSE loss, a sparsity-inducing regularization (can be chosen from 4), and the smooth penalty, which is defined as follows:

$$\frac{1}{T-\tau }\sum _{t=\tau +1}^T{\Vert {𝐱}^t-{\widehat{𝐱}}^t\Vert }_2^2+\frac{\lambda }{T-\tau }\sum _{t=\tau +1}^{T}R({\mathrm{\Psi }}_t)+\frac{\gamma }{T-\tau -1}\sum _{t=\tau +1}^{T-1}{\Vert {\mathrm{\Psi }}_{t+1}-{\mathrm{\Psi }}_t\Vert }_2^2,$$

here ${\{{𝐱}^t\}}_{t=1}^T$ is the observed d-variate time series whereas ${\widehat{𝐱}}^t$ is the one-step forecast made by the GVAR model. Now that the interpreting matrices for each time point $t$ can be derived via ${\mathrm{\Psi }}_{{\widehat{\theta }}_k}({𝐱}^t)$, the signs of Granger causal effects and their variability in time can also be assessed. Furthermore, a procedure of GVAR based on the heuristics of time-reversed Granger causality (Winkler et al., 2016), which expects the relationships to be flipped on time-reversed data, is leveraged to improve the stability of the inferred structures. Compared to the aforementioned methods, such as cMLP, cLSTM, eSRU and MPIR, another key difference is that these methods require training $d$ neural networks, whereas GVAR requires training $2\tau$ networks.

DL-extensions with Attention Mechanisms. The temporal causal discovery framework (TCDF) is introduced in (Nauta et al., 2019), which utilizes attention-based dilated CNN. This framework consists of $d$ independent attention-based CNNs with the same architecture but different target variable $X_j$. For each target variable, a neural network is proposed to derive prediction, attention scores and kernel weights. Intuitively, a high attention score on $X_i$ while forecasting $X_j$ indicates the former contains prediction information towards the latter. A permutation-based procedure is additionally provided for evaluating variable importance and identifying significant causal links. TCDF can discover self-causation and time delays between cause and effect. Besides, by assuming that the bidirectional causal relations can not be instantaneous, it can also detect the presence of hidden confounders with equal delays.

Besides, an interpretable multi-variable LSTM with mixture attention is proposed in IMV-LSTM (Guo et al., 2019, 2018) to extract variable importance knowledge. And it’s widely used as a baseline for causal discovery in multivariate time series. However, the topic on attention and its interpretation is to some extent still a controversial and inconclusive topic (Jain and Wallace, 2019; Wiegreffe and Pinter, 2019; Grimsley et al., 2020). Especially in the context of Granger causal explanation, the naive-trained soft attention mechanisms are noted (Sundararajan et al., 2017; Schwab et al., 2019; Chu et al., 2020) to provide no incentive to yield accurate attributions. In (Schwab et al., 2019), Granger causal attention weights are introduced based on the measures named as the mean Granger-causal error. The decrease in error when adding $i$ can be computed as:$\mathrm{\Delta }{\epsilon }_{X,i}={\epsilon }_{X\mathrm{\backslash }\{i\}}-{\epsilon }_X$, given the auxiliary prediction error ${\epsilon }_X,{\epsilon }_{X\mathrm{\backslash }\{i\}}$ with and without any information from the $i$-th variable. Then the Granger-causal attention factor can be computed as:${\omega }_i(X)=\frac{\mathrm{\Delta }{\epsilon }_{X,i}}{{\sum }_{j=1}^d\mathrm{\Delta }{\epsilon }_{X,j}}$. The attention factor ${\omega }_p(X)$ is able to capture Granger causality, which is zero if $p$-th time series is Granger noncausal for the target series.

DL-extensions with Recurrent Variational Autoencoders. Recently, the causal recurrent variational autoencoder (CR-VAE) (Li et al., 2023) is proposed, where a generative model incorporates Granger causal learning into the data generation process. By preventing encoding future information before decoding, the encoder of CR-VAE obeys the principle of Granger causality To be specific, given time lag $\tau$, a CR-VAE model can be written as:

$${\widehat{𝐱}}^{t-\tau :t}=D_{\theta }({𝐱}^{t-\tau :t-1},E_{\psi }({𝐱}^{t-2\tau -1:t-\tau -1}))+{ϵ}^t,$$

where $E_{\psi },D_{\theta }$ represent encoder and decoder. Another distinct to the classical recurrent VAE is that the CR-VAE leverages a multi-head decoder where the $i$-th head is designed for generating ${𝐱}_i$. Besides, an error-compensation module is leveraged to capture instantaneous effects. The CR-VAE is not only able to extract causal relations, but also conduct the data-generating process in a transparent manner benefiting from the learned causal matrix.

DL-extensions with Inductive Modeling. The problem of methods with inductive modeling is slightly different from the above methods, where MTS data from massive individuals, which entails different causal mechanisms but shares common structures, is collected. The goal is to train a model on samples with heterogenous structures to discover Granger causal relations from each individual. Two approaches with inductive modeling are reviewed here.

An inductive Granger causal modeling (InGRA) is proposed in (Chu et al., 2020), combined with Granger causal attention (Schwab et al., 2019) and prototype learning. As there often exist real-world scenarios where massive multivariate time series data is collected from heterogeneous individuals sharing commonalities. Instead of training one or a set of models for each individual, InGRA trains a global model for individuals potentially having different Granger causal structures, devoid of sample inefficiency and over-fitting issues. Firstly, the Granger causal attention mechanism is leveraged to quantify variable-wise contributions toward prediction. As the Granger causal attention is not robust enough to reconstruct Granger causal topology from limited data of a single individual, InGRA secondly leverages prototype learning, of which the key idea is to solve problems for new inputs based on similarity to prototypical cases, to detect common causal structures. As a result, the Granger causal relations and strengths between the $d-1$ exogenous variables and the target variable are inferred.

A framework termed amortized causal discovery (ACD) is proposed in (Löwe et al., 2022), which aims to train a single model to infer causal relations across samples with different underlying causal graphs but shared dynamics. It’s an encoder-decoder framework, in which the encoder function is defined to infer Granger causal relations of the input sample whereas the decoder function learns to predict the next time-step given the inferred causal relations. In the implementation, a graph neural network is applied to the amortized encoder, and ACD models the functions using variational inference, which is based on the widely used neural relational inference (NRI) model (Kipf et al., 2018). Besides, to derive a causal interpretation of the inferred edges, the proof is provided in ACD to relate the zero-edge function to Granger causality. As a result, the causal relations of previous unseen samples can be inferred without refitting the model.

3.5. Others

The aforementioned four categories of approaches have been the subjects of many endeavors in causal discovery research. For the sake of completeness, we present five types of methods that are distinct from the above approaches in this subsection, including causality based on information-theoretic statistics, causal models based on differential equations, nonlinear state-space methods, logic-based methods, and hybrid methods.

3.5.1. Causality based on information-theoretic statistics

Causal relationships in MTS can be measured based on information-theoretic statistics. As a model-free measure, it’s widely used in constraint-based approaches (3.1.1) and Granger causal models (3.4.2). However, its definitions and characteristics have not been detailed. In this part, we will first introduce Transfer Entropy (Schreiber, 2000), which is the original concept of information-theoretic statistics of causality, and then its variants.

Transfer Entropy (Schreiber, 2000) is a measure of information flow or effective coupling between two processes, regardless of the actual functional relationship. Instead of model-based criterion, which shares the problem that the model might be misspecified, as a model-free measure, it can be combined with a variant of specific structure learning methods. In detail, the Transfer Entropy from $i$ to $j$ (with time lag) can be expressed as:

$$\mathrm{TE}(X_i^t\to X_j^{t+1})=h(X_j^{t+1}|X_j^t)-h(X_j^{t+1}|X_j^t,X_i^t),$$

where $h(\cdot |\cdot )$ denotes the conditional entropy. Here the term $h(X_j^{t+1}|X_j^t)$ measures the uncertainty of $X_j^{t+1}$ given information about $X_j^t$, and $h(X_j^{t+1}|X_j^t,X_i^t)$ measures the uncertainty of $X_j^{t+1}$ given information about both $X_j^t$ and $X_i^t$. Therefore, we can understand Transfer Entropy $\mathrm{TE}(X_i^t\to X_j^{t+1})$ in the causal view of reduction of uncertainty about future dynamics of $X_j$ when the current dynamics of $X_i$ is given addition to that of $X_j$. For Gaussian variables, the equivalence between Transfer Entropy and Granger causality is demonstrated in (Barnett et al., 2009). Furthermore, Transfer Entropy is reformulated into a decomposition form and embedded into the framework of graphical models for multivariate in (Runge et al., 2012b). In (Runge et al., 2012a), the causal coupling strength for multivariate time series is quantified based on a variant of transfer entropy.

Although some utilization in multivariate scenarios, Transfer Entropy suffers from the pairwise limitation. And it’s reported to fail to distinguish between direct and indirect causality in networks (Sun et al., 2015) . As a remedy to pairwise limitation, Causation Entropy (Sun and Bollt, 2014), a model-free information theoretic statistic for inferring causality, is introduced. In detail, the Causation Entropy from the set of nodes $I$ to the set of nodes $J$ conditioning on the set of nodes $C$ is defined as follows:

$$\mathrm{CE}(X_{𝐈}^t\to X_{𝐉}^{t+1}|X_{𝐂}^t)=h(X_{𝐉}^{t+1}|X_{𝐂}^t)-h(X_{𝐉}^{t+1}|X_{𝐂}^t,X_{𝐈}^t),$$

here $I,J,C$ are all subset of nodes $\{1,2,\mathrm{\dots },d\}$. As a type of conditional mutual information, Causation Entropy is a generalization of Transfer Entropy for measuring pairwise relations to network relations of many variables. And similar to the equivalence relations between Transfer entropy and Granger Causality, Causation Entropy also generalizes Granger Causality and Conditional Granger Causality when applied to Gaussian variables. However, according to its definition, this concept assumes that the hidden dynamics follow a stationary first-order Markov process as the Causation Entropy only models causal relations with time lags equal to one. Recently, to measure any lagged or instantaneous relations, an extension of Causation Entropy, named Greedy Causation Entropy, is proposed in (Assaad et al., 2022a).

3.5.2. Causal models based on differential equations

Differential equations are a commonly used modeling tool in many fields, and are especially useful if measurements can be done on the relevant time scale. Compared to the aforementioned causal models, this type of approach is specifically designed to model systems that can be well represented by differential equations (Peters et al., 2022). In this part, we will first review the relationships between differential equations and causal models, for both discrete and continuous time. The first difference-based causal discovery framework will be introduced. Then, we will give the recent advances in this type of method.

There is abundant literature (Peters et al., 2022; Bongers et al., 2018) (Schölkopf, 2019; Mooij et al., 2013; Rubenstein et al., 2018) discussing the relationship between differential equations and structural causal models. For discrete time, a difference-based causal discovery framework is first proposed in (Voortman et al., 2010). The cross-temporal restriction is satisfied, where all causation across time is due to a derivative $\dot{x}$ causing a change in its integral $x$. This characteristic makes the difference-based causal model a restricted form of dynamic SEMs. And difference-based causality learner (DBCL) is leveraged to extract difference-based causal models from data, which is proven to be able to identify the presence or absence of feedback loops. For continuous time, several theoretical endeavors have also been made to derive a causal interpretation of dynamic systems by both ordinary differential equations (ODEs) (Mooij et al., 2013; Rubenstein et al., 2018; Blom et al., 2019; Pfister et al., 2019) and stochastic differential equations (SDEs) (Hansen and Sokol, 2014; Mogensen et al., 2018).

More recently, under a dynamic causal system where the multivariate time series are irregularly-sampled (in infinitesimal interval of time), an algorithm called neural graphical model (NGM) is proposed in (Bellot et al., 2022). In many applications, the underlying causal system of interest can be represented as a dynamic structural model as follows:

$$d𝐱(t)=𝐟(𝐱(t))dt+d𝐰(t),𝐱(0)={𝐱}_0,t\in [0,T],$$

where $𝐰(t)$ is a $d$-dimensional standard Brownian motion, ${𝐱}_0$ is a Gaussian random variable independent of $𝐰(t)$, and the function $𝐟$ describes the causal graph $G$. NGM is a learning algorithm based on penalized Neural Ordinary Differential Equations (neural-ODE). The recovery of causal graph can be cast to penalized optimization problems of the form:

$${\min }_{{𝐟}_{\theta }}\frac{1}{n}\sum _{i=1}^n{\Vert 𝐱(t_i)-\widehat{𝐱}(t_i)\Vert }_2^2,\mathrm{subject}\mathrm{to}{\rho }_{n,T}({𝐟}_{\theta })\mathrm{and}\widehat{𝐱}(t)={𝐟}_{\theta }(\widehat{𝐱}(t_i))dt,$$

where the observation of the systems are at irregular time points .

3.5.3. Nonlinear state-space methods

In this part, we will first introduce the basics of nonlinear state-space methods, including the Takens theorem and Convergent Cross Mapping algorithm. Then variants and recent advances of the original algorithm will be given, to tackle the challenges such as high sensitivity to noise, large sample demands, inconsistent results, and misidentifications.

The state space reconstruction theory proposed by Takens (Takens, 1981) provides a theoretical basis for analyzing the dynamic characteristics of nonlinear systems. Based on this theory, another approach for determining causality, known as Convergent Cross Mapping (CCM), was first proposed in (Sugihara et al., 2012). Developed for coupled time series, this method leverages Takens’ theorem via state space reconstruction. In detail, given two time series $x_1^t$ and $x_2^t$, the attractor manifolds ${ℳ}_{x_1},{ℳ}_{x_2}$ are first reconstructed using $x_1^t$ and $x_2^t$, respectively. Secondly, causality can be detected by measuring the correspondence between ${ℳ}_{x_1}$ and ${ℳ}_{x_2}$, to be specific, by testing whether every local neighborhood defined on one manifold is preserved in the other. Figure 7 gives the illustration of CCM. This methodology has been successfully applied in many fields (Hirata et al., 2016; Ye et al., 2015) where nonlinear systems are dynamically coupled.

However, there exist issues for the original CCM method, such as high sensitivity to observation noise, a requirement for a relatively large number of observations, and inconsistent results under different optimal algorithms. To overcome these challenges, variants of CCM based on time-lagged analysis (Ye et al., 2015), deep Gaussian process (Feng et al., 2019) reservoir computing (Huang et al., 2020b) and neural ODE (Brouwer et al., 2021) were proposed. Besides, most CCM-based approaches have been originally developed for bivariate analysis. Although the same procedures may be used multiple times to ascertain the causal network among multivariate time series, the performance is not guaranteed under high-dimensional conditions (Huang et al., 2020a). Misidentifying indirect causations as direct ones performs one of the key challenges in multivariate settings. Recently, partial cross mapping (PCM), which combines CCM with partial correlation, was proposed (Leng et al., 2020) to eliminate indirect causal influences.

Figure 7. The illustration of the cross convergent mapping procedure.

3.5.4. Logic-based methods

Another type of methodology, used for causal inference and causal discovery in time-series data, is based on logic formulas. The original algorithm of this type of approach will first be introduced and combined with its semeiology and definition of potential causality. Then we will give its variants and recent advances.

In logic-based methods, temporal data can be thought of as observations of the sequence of states the system has occupied and is referred to as traces in model checking. This line of research originates from work in (Kleinberg and Mishra, 2009), where causal relationships are described in terms of temporal logic formulas. To be specific, it first leverages logic, Probabilistic Computation Tree Logic (PCTL), to define prima facie (potential) causality based on temporal priority and the uplift of conditional probability. Given the notation in the original work, the prima facie cause is defined if the following conditions all hold: (1) $F_{>0}^{\le \mathrm{\infty }}c$, (2) $c{\rightsquigarrow }_{\ge p}^{\ge 1,\le \mathrm{\infty }}e$, and (3) , implying that there may exist any number of transitions between $c$ and $e$ and the sum of a set of path probabilities are at least $p$. To separate the underlying prima facie (potential) causes into genuine and spurious causes, the notion of $ϵ$-insignificant cause is introduced by computing the average difference in probabilities for each prima facie cause of an effect in relation to all other prima facie causes of the effect:

$${ϵ}_{avg}(c,e)=\frac{{\sum }_{x\in X\mathrm{\backslash }c}{ϵ}_x(c,e)}{|X|},$$

where ${ϵ}_x(c,e)=P(e|c\wedge x)-P(e|\neg c\wedge x)$. A prima facie cause $c$ is an $ϵ$-insignificant cause of $e$ if . The value of $ϵ$ is chosen based on empirical null hypothesis testing by assuming: (1) data contains two classes, significant and insignificant, (2) the significant class is relatively small to the insignificant class. And false discovery rate control is implemented simultaneously. This methodology has also applications in fields (Kleinberg, 2013).

To expand the methodology to the condition where both discrete and continuous components exist, PCTLc is introduced in (Kleinberg, 2011) to express temporal and probabilistic properties involving discrete and continuous variables, and the significance of relationship in the continuous case is validated via conditional expectation of an effect instead of conditional probability. Besides, a variant (Huang and Kleinberg, 2015) of this logic-based approach was proposed to improve the accuracy of causal discovery and enable faster computation of causal significance, by showing the computational complexity can be reduced under several conditions. Following this line of temporal logic form, a recent work (da Costa and Dasgupta, 2021) combines the idea of decision trees and reconsiders the problem of causal discovery to extract temporal causal sequence relationships from real-time time series.

3.5.5. Hybrid methods: combining score-based and constraint-based approaches

Hybrid approaches are proposed for the benefit of combining the strengths of both constraint-based (3.1) and score-based (3.2) approaches. We cover two parts of hybrid methods, including methods based on max-min hill-climbing heuristics, and methods incorporating the conditional independence tests to improve the local search.

Some researchers develop hybrid approaches based on max-min hill-climbing heuristics (Tsamardinos et al., 2006; Li and Ngom, 2013; Li et al., 2016). As hybrid local learning methods, Max-Min approaches fuse concepts from both constraint-based techniques to limit the space of potential structures and search-and-score Bayesian methods to search for an optimal structure. They are originally leveraged in the structure learning of BN for static data (Tsamardinos et al., 2006). The Max-Min hill-climbing Bayesian network (MMHO-DBN), introduced in (Li et al., 2016), learns the structure of DBN based on an extension of the max-min hill-climbing heuristic and is leveraged in the modeling of real gene expression time-series data.

There are also hybrid approaches that combine conditional-independence tests and local search to improve the criterion score (Ogarrio et al., 2016; Malinsky and Spirtes, 2018). Greedy FCI (GFCI) (Ogarrio et al., 2016) is a hybrid score that combines features of GES with FCI. SVAR-GFCI (Malinsky and Spirtes, 2018) extends this method to causal structure learning from time series. In (Sanchez-Romero et al., 2019), both a variant of the PC-stable algorithm referred to as Fast Adjacency Skewness (FASK), and a hybrid two-step algorithm is proposed for extracting causal relations for time-series data.

4. Causal Discovery from Event Sequences

An important assumption in multivariate time series is that the timestamps are discrete and the time intervals are fixed. However, in the real-world scenario, the vast majority of events will not occur at fixed intervals. In consequence, we need to come up with some methods to deal with these irregular and asynchronous data. We can construct event sequences as $\{(t_1,e_1),(t_2,e_2),\mathrm{\dots }\}$, where the first dimension represents the time at which the corresponding event happens, and the second dimension stands for the corresponding event type. In this section, we will focus on inferring causal relationships in event sequences. First, the multivariate point process is introduced, which is preliminary for causal discovery in event sequences. Then, we review approaches based on the Granger causal model, which are well-developed. Lastly, other approaches including constraint-based and score-based methods are given.

4.1. Multivariate Point Process

An event sequence records the occurrence of one specific type of event (or ‘event type’ for simplicity). Meanwhile, we can characterize an event sequence through a point process. To discover the relationships between different types of events, we consider its high-dimensional cases, which is to model event sequences through Multivariate Point Processes (MPPs). Therefore, our problem can be defined as inputting a set of point processes, where each point process represents an event sequence, and outputting a causal graph $G$ established by different processes. In the causal graph $G$, each node represents a point process, and each directed edge captures a directed interaction from one point process to another. In this part, we will detail MPPs, including their intensity functions and log-likelihood functions.

Intensity Functions of MPPs. A temporal point process is a stochastic or random process composed of a time series of binary events that occur in continuous time (Daley et al., 2003). MPPs are high-dimensional point processes, implying that they can involve multiple types of events. $ℰ=\{1,\mathrm{\dots },E\}$ is the set of event types. The occurring time of these events $\{t_1,t_2,\mathrm{\dots },t_n|t_i\in [0,T]\}$ are unevenly-distributed. The multivariate point process with $E$ types of events can be represented by $E$ counting processes ${\{N_e\}}_{e=1}^E$, where $N_e=\{N_e(t)|t\in [0,T]\}$. The core of a point process is its conditional intensity function, in which the process’s pattern is captured. A type-$u$ intensity function can be defined as the expected instantaneous rate of type-$e$ event’s occurrence given the history:

$${\lambda }_e(t)=\frac{𝔼[d{N}_e(t)|{\mathscr{H}}_t]}{dt}$$

Here ${\mathscr{H}}_t=\{(t_i,e_i)|t_i<t,e_i\in ℰ\}$ represents all types of events happened before time $t$.

Log-likelihood Functions of MPPs. Next, we show the relationship between the intensity function and the Probability Density Function(PDF) of the joint distribution: $f((t_1,e_1),\mathrm{\dots },(t_n,e_n)|(t_0,e_0))$. Using the chain rule, there is $f((t_1,e_1),\mathrm{\dots },(t_n,e_n)|(t_0,e_0))={\prod }_{j=1}^{n}f((t_j,e_j)|{\mathscr{H}}_{t_j})$. Then, we can set up the likelihood function for estimating the joint distribution:

(1) $$\widetilde{L_0}\triangleq \sum _{j=1}^{n}lnf(t_j|e_j,{\mathscr{H}}_{t_j})+\sum _{j=1}^{n}lnf(e_j|{\mathscr{H}}_{t_j})$$

Since the goal is to infer the causal relationship between different events, here we focus on the first term and omit the second term: $L_0\triangleq {\sum }_{j=1}^{n}lnf(t_j|e_j,{\mathscr{H}}_{t_j})$. The intensity function reflects the expectation of the event happening in $[t,t+dt]$ given the information of ${\mathscr{H}}_{n-1}\triangleq {\mathscr{H}}_{t_{n-1}}$. Similar to the calculation of the force of mortality in survival analysis, there holds,

${\lambda }_e(t|{\mathscr{H}}_{n-1})$ $={\displaystyle \frac{𝔼[N_e(t+dt)-N_e(t)|{\mathscr{H}}_{n-1}]}{dt}}$

$={\displaystyle \frac{{\int }_{t_{n-1}}^{t+dt}f(l|e,{\mathscr{H}}_{n-1})𝑑l-{\int }_{t_{n-1}}^{t}f(l|e,{\mathscr{H}}_{n-1})𝑑l}{1-{\int }_{t_{n-1}}^{t}f(l|e,{\mathscr{H}}_{n-1})𝑑l}}$

$=-{\displaystyle \frac{d}{dt}}ln(1-{\displaystyle {\int }_{t_{n-1}}^t}f(l|e,{\mathscr{H}}_{n-1})𝑑l)$

Integrating the equation above and substituting $f$’s expression into 1, we have,

(2) $$L_0=\sum _{j=1}^n\{ln{\lambda }_{e_j}(t_j|{\mathscr{H}}_{j-1})-{\int }_{t_{j-1}}^t{\lambda }_{e_j}(l|{\mathscr{H}}_{j-1})𝑑l\}$$

We have briefly introduced Multivariate Point Processes and constructed likelihood functions to characterize MPPs in the above. Next, we aim to discover the causal relationships within MPPs using Granger-based, as well as constraint-based and score-based methods.

4.2. Granger Causality Based Approaches

In this subsection, we consider the task to infer Granger causalities in event sequences. Similar to that in MTS, we say $e_j$-type events Granger cause $e_i$ if is useful in forecasting $e_i(t)$. The detailed methods can be categorized according to the following model specifications, i.e., GLM point process, Hawkes process, Wold process, and neural point process.

4.2.1. Methods for GLM point processes

We first introduce causal discovery approaches for event sequences which are modeled via Generalized Linear Model (GLM) of point processes (Truccolo et al., 2005). The GLM assumes that the logarithm of the intensity function has a linear format, i.e., $ln{\lambda }_e(t|{\mathscr{H}}_{n-1})={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2$. Specifically, in our mission, the intensity functions follow,

(3) $$ln{\lambda }_i(t|{\gamma }_i,H_i(t))={\gamma }_{i,0}+\sum _{j=1}^J\sum _{m=1}^{M_i}{\gamma }_{i,j,m}{R}_{j,m}(t)$$

Here ${\gamma }_{i,0}$ can be interpreted as the background intensity of event $e_i$, ${\gamma }_{i,j,m}$ is the intensity on type-$e_i$ events triggered by type-$e_j$ events and $R_{j,m}(t)$ is the number of occurrence of $e_j$-type events happened in $[t-mW,t-(m-1)W]$ ($W$ is a small number which refers to the length of time range). By looking at the sign of ${\sum }_{m=1}^{M_i}{\gamma }_{i,j,m}$, we can distinguish whether type-$e_j$ events have excitatory or inhibitory effects on type-$e_i$ events.

To infer the Granger Causality between type-$e_j$ and type-$e_i$ events, we substitute 3 into the likelihood function 2. Next, we follow a simple thought that we can exclude a certain type of event and then infer the Granger causality by comparing its intensity with the original case. Specifically, we obtain both the likelihood of type-$e_i$’s occurrence with and without type-$e_j$’s effect: $L_i({\gamma }_i)$, $L_i({\gamma }_i^j)$. Then, consider that ${\sum }_{m=1}^{M_i}{\gamma }_{i,j,m}$ is an indicator of the effection type, the Granger causality from type-$e_j$ to type-$e_i$ events can be proposed as (Kim et al., 2011):

${\varphi }_{ij}$ $=-sign({\displaystyle \sum _{m=1}^{M_i}}{\gamma }_{i,j,m}){\mathrm{\Gamma }}_{ij}$

$=-sign({\displaystyle \sum _{m=1}^{M_i}}{\gamma }_{i,j,m})log{\displaystyle \frac{L_i({\gamma }_i^j)}{L_i({\gamma }_i)}}$

Apparently, there exists $L_i({\gamma }_i)\ge L_i({\gamma }_i^j)$, hence, ${\mathrm{\Gamma }}_{ij}=log\frac{L_i({\gamma }_i^j)}{L_i({\gamma }_i)}\le 0$. Only if ’’ is satisfied, type-$e_j$ events will be the Granger cause to type-$e_i$ events. In the next step, Kim et al. (Kim et al., 2011) presented a significance test of these causal interactions by conducting $H_0$ hypothesis: ${\theta }_0={\gamma }_i^j$ and $H_1$ hypothesis: ${\theta }_1={\gamma }_i$. Passing through an FDR significance test, the final causal relationships could be estimated by $\widetilde{{\varphi }_{ij}}$: (1) type-$e_j$ events are an excitatory cause of type-$e_i$ events when $\widetilde{{\varphi }_{ij}}>0$, (2) the cause is inhibitory when , (3) there exists no causal relationship between type-$e_j$ and type-$e_i$ events when $\widetilde{{\varphi }_{ij}}=0$.

4.2.2. Methods for Hawkes processes

In this part, we review methods for Hawkes process. As a particular type of point process, the basics of the Hawkes process are first given. Then we detail approaches based on MLE to infer causal relations, including (1) parameterization strategies and (2) regularization methods. Next, we review other estimation approaches, including (1) graphical event models, (2) generalized method of moments, (3) event sequence separation, and (4) minimum description length. We note that there exists a plethora of literature in this category since a natural match-up between Granger causality and Hawkes processes.

The Hawkes process is a type of point process that has a fixed form of intensity function:

(4) $${\lambda }_{e_i}(t)={\mu }_{e_i}+\sum _{e_j=1}^E{\int }_0^t{\varphi }_{e_i{e}_j}(s)𝑑N_{e_j}(t-s)$$

Here, ${\mu }_{e_i}$ is called the baseline intensity, which can only be affected by exogenous events, hence, is a constant over time. And ${\varphi }_{e_i{e}_j}(s)$, the impact function, measures the decay of the excitement on future type-$e_i$ events triggered by historical type-$e_j$ events. That is to say, it captures the endogenous intensity from $e_j$ to $e_i$. Considering the similarity between definitions of $\varphi$ and the Granger causality, we can directly infer Granger causality by analyzing $\varphi$:

Proposition 1.

(Eichler et al. (Eichler et al., 2017), 2017)

$$e_j\text{does not Granger-cause}{e}_i\iff {\varphi }_{e_i{e}_j}(s)=0,\forall s\in R$$

Therefore, we aim to model ${\varphi }_{e_i{e}_j}(t)$ for each event and all $t\in R$. However, due to the complexity and heterogeneity of event sequences, this mission could be extremely difficult to accomplish. Zhou et al. (Zhou et al., 2013) parameterize ${\varphi }_{e_i{e}_j}(s)$ as $a_{e_i{e}_j}g(s)$. By this means, $\varphi$ is split into events-interaction and time-decaying parts.

MLE Approaches. The Maximum Likelihood Estimation (MLE) can be performed for estimating parameters in 4. We take $\lambda$’s expression into 2, which results in the corresponding likelihood function: $L(A,\mu )$. Here, $A$ is composed of $(a_{e_i{e}_j})$, and $\mu$ is built up by ${\mu }_{e_i}$. Next, consider that in real-world scenarios, most events can only influence a small fraction of other events, and the community structures in the influence networks tend to be low-ranked (Zhou et al., 2013), we should add penalty entries to the MLE loss function. Specifically, the following objective function can be constructed in order to achieve matrix $A$’s low rank and sparsity:

$$\underset{A\ge 0,\mu \ge 0}{min}-L(A,\mu )+{\lambda }_1{\Vert A\Vert }_{\ast }+{\lambda }_2{\Vert A\Vert }_1$$

Here, $\parallel \cdot {\parallel }_{\ast }$ is the nuclear norm, of which performance in reducing the matrix’s rank has been proven. And $\parallel \cdot {\parallel }_1$ is the L1 norm. It can enforce matrix A to gain more sparsity. ${\lambda }_1$, ${\lambda }_2$ are parameters to control the strength of these two penalties. We denoted the object function as $f(A,\mu )$. Apart from this, an EM-based algorithm could be conducted to solve the optimization problem for $A$ and $\mu$. In details, Zhou used the surrogate function $Q(A,\mu ;A^{(m)},{\mu }^{(m)})$ as a tight upper bound of $f(A,\mu )$. By optimizing $Q(A,\mu ;A^{(m)},{\mu }^{(m)})$ iteratively, $f(A,\mu )$ was forced to decrease and, thus, successfully optimized. We summarize this in algorithm LABEL:alg:Algorithm1.

(1) Parameterization Strategies: The parameterization method ${\varphi }_{e_i{e}_j}(s)$ = $a_{e_i{e}_j}g(s)$ mentioned above may suffer from bad performance if the data do not fit its strong assumption. Therefore, to gain robustness in different types of event sequences, Xu et al. (Xu et al., 2016) came up with a strategy to choose a family of basic functions and used their linear combination, ${\sum }_{m=1}^M{a}_{e_i{e}_j}^m{\kappa }_m(s)$, to model the targeted intensity function.

However, in NPHC (Achab et al., 2017), Achab et al. put forward that the basic function strategies would have extraordinary computing complexities when there exist too many types of events (i.e., $E$ is large). Given that our goal is only to infer the Granger causality, there is no need to totally parameterize the Hawkes process, and thus, we only need to estimate the corresponding integral ${\int }_0^{+\mathrm{\infty }}{\varphi }_{e{e}^{\prime }}(s)𝑑s$. Achab denoted the integral as $g_{e{e}^{\prime }}$ while $(g_{e{e}^{\prime }})$ formed matrix G. Then, from Eichler’s proof (Eichler et al., 2017) as well as ${\varphi }_{e{e}^{\prime }}(s)>0,\forall s>0$, it is clear that $(g_{e{e}^{\prime }})=0\iff e^{\prime }\text{does not Granger-cause}e$.

Other works considered the underlying topological relationships within event sequences. In THP (Cai et al., 2021), Cai et al. assumed the existence of a hidden undirected graph structure $G_N$ between events. And the corresponding intensity function is formed as ${\lambda }_{e_i}(n,t)={\mu }_{e_i}+{\sum }_{e_j\in E}{(g_{e_i{e}_j}\ast s_{e_i,e_j,t})}_{G_N}(n)$. Here, $g_{e_i{e}_j}$ is the graph convolution kernel that could capture the effect from the graph neighbors. And $s_{e_i,e_j,t}$ is the time convolution kernel representing the sum of the past impact function . This is based on the assumption that the hidden topological structure will not change during the process.

(2) Regularization Methods: In an aforementioned method we presented $A$’s nuclear norm ${\Vert A\Vert }_{\ast }$ and L1 norm ${\Vert A\Vert }_1$ as regularizers. And in basic function methods, a special sparse-group-lasso regularizer (Simon et al., 2013) is applied to fit their summation parameterization. Specifically, Xu et al. (Xu et al., 2016) conducted a group-lasso penalty as well as a lasso penalty simultaneously in order to enforce $a_{e_i{e}_j}^m=0$ for all $m$, i.e., group sparsity, in addition to a regular sparsity for all the entries $a_{e_i{e}_j}^m$. Nevertheless, in $L_0$ Hawkes (Idé et al., 2021), Ide et al. proved that EM-based MLE algorithms with L1-regularization cannot offer sparse solutions mathematically. Hence, their sparse solutions can appear only as numerical artifacts. Sequentially, Ide presented an L0-regularized EM-MLE algorithm to circumvent this problem. Here, the L0-norm ${\Vert A\Vert }_0$ indicates the number of non-zero entries in matrix $A$.

Similar to the aforementioned topological parameterization strategies, Xu et al. (Xu et al., 2016) considered the underlying topological relationships between event types in constructing our regularizers. Specifically, pairwise similarity ${\sum }_{e_i=1}^E{\sum }_{e_j\in C_e}{\Vert a_{e_i\cdot }-a_{e_j\cdot }\Vert }_F^2+{\Vert a_{\cdot e_i}-a_{\cdot e_j}\Vert }_F^2$ could be presented in order to enforce that similar events could have similar intensity functions. However, we must add that this regularizer requires a predefined cluster structure and thus can be optimized.

Prior domain knowledge could be of great use when discovering the Granger causality. Due to the event sequences’ high dimensionality and heterogeneity, existing algorithms regularly suffer from underfitting and poor interpretability. Hence, it is natural to consider adding domain knowledge from humans to the causal-inferring model. In specific, a bottom-up visualization model with user feedback was established (Jin et al., 2021). Jin et al. set up their based model with the traditional MLE method in MLE-SGLP. During the training process, the user could either confirm or remove a causal relation depending on their domain knowledge from the network. And the model will change its optimization target corresponding to the user’s choice. For example, in accordance with the idea in MLE-SGLP (Xu et al., 2016), Jin (Jin et al., 2021) constructed their intensity function as ${\varphi }_{e_i{e}_j}={\sum }_{m=1}^M{a}_{e_i{e}_j}^m{\kappa }_m(s)$, and set $a_{e_i{e}_j}$ as $[a_{e_i{e}_j}^1,\mathrm{\dots },a_{e_i{e}_j}^n]$. Correspondingly, their objective function could be: $\underset{\mu ,\alpha }{argmin}-L+\alpha {\sum }_{e_i,e_j}{\Vert a_{e_i{e}_j}\Vert }_2$. After the user made their choice to either confirm or delete edges in the causal graph $\widehat{G}$, Jin updated the object function as follows:

(5) $$\begin{array}{c}\hfill \underset{\mu ,\alpha }{argmin}-L+{\alpha }_v\sum _{e_i,e_j}{\Vert a_{e_i{e}_j}(\widehat{G})\Vert }_2\\ \hfill \text{s.t.}{a}_{e_i{e}_j}=0\text{for}(e_j\to e_i)\notin \widehat{G}\end{array}$$

(6) $$a_{e_i{e}_j}(\widehat{G})=\{\begin{array}{cc}0;\text{if}(e_j\to e_i)\text{is confirmed}\hfill & \\ a_{e_i{e}_j};\text{otherwise}\hfill & \end{array}$$

Here the constraints in 5 fit the removal operations, and the updates in 6 represent the user’s confirmations.

Other estimation approaches

(1) Graphical Event Models: The aforementioned methods use Maximum Likelihood Estimation to model the Hawkes processes of event sequences. However, these attempts lack interpretability and require fine-tuning processes for parameters to achieve a good performance. Therefore, entirely data-driven, graph-based, and dependency-captured Graphical Event Models (GEMs) could be presented to infer Granger causalities in the event sequences.

We will elaborate more on GEM’s attributes in 4.3. Here, we only focus on its relationship with the Granger causality. Suppose there is a directed graph $𝒢=(ℰ,𝒜)$, in which the edges represent dependencies between different event types. For each event-type $e$, we assumed that its conditional intensity could only be affected by its parent type, i.e., it follows ${\lambda }_e(t|h_t)={\lambda }_e(t|{[h_t]}_{P_a(e)})$, where $P_a(e)\subseteq ℰ$ is $e$’s parent event in a graph $𝒢$, and ${[h_t]}_{P_a(e)}$ is the history of events which types are listed in the set $P_a(e)$. In accordance with 1, there holds,

Proposition 2.

(Granger Causality in GEMs, Yu et al., 2020 (Yu et al., 2020))
For two event types $e_i$ and $e_j$ in $𝒢=(ℰ,𝒜)$, $e_j$ does not Granger-cause $e_i\iff {\varphi }_{e_i{e}_j}(t)=0,\forall t>0\iff e_j\notin Pa(e_i)$

Hence, one can apply traditional score-based structure learning methods to discover the Granger causality. For example, BIC scores can be presented for learning the optimized graph ${𝒢}^{\ast }$. The optimization approach is consistent. At the same time, Yu conducted a Forward-Backward Search (FBS) to learn the parent types of a certain event type independently (Yu et al., 2020). The Forward-Backward Search with BIC scores is proved to be sound and complete for a family of GEMs (Gunawardana and Meek, 2016).

(2) Generalized Method of Moments: In NPHC, the optimization object is a matrix $G=(g_{e{e}^{\prime }})=({\int }_0^{+\mathrm{\infty }}{\varphi }_{e{e}^{\prime }}(s)𝑑s)$. Therefore, the Generalized Method of Moments (GMM) can be used to address this problem (Hall, 2004). Achab et al. presented a GMM-based NPHC algorithm to model the first, second, and third-order cumulants of matrix $G$ (Achab et al., 2017). Afterward, the Granger causality could be directly attained from $G$. This moment estimation approach is proven consistent and robust to certain observation noise (Trouleau et al., 2021). However, this approach might receive poor results in specific datasets, e.g., datasets with long tails. That is mainly due to GMMs’ general issue: they can only capture the information within a statistical distribution’s moments.

(3) Event Sequences Separation: Another intriguing idea is separating the event sequences into multiple sub-sequences and applying the Hawkes Process model in each sub-sequence correspondingly. In GC-nsHP (Chen et al., 2022a), Chen et al. divided the event sequences ${\mathscr{H}}_n=[(t_1,e_1),\mathrm{\dots },(t_n,e_n)]$ into $K$ different patterns, where $K$ should be predefined according to its applying scenario. ’Events’ in the same pattern are supposed to build up a stationary sub-process of ${\mathscr{H}}_n$. Then, $K$ different Hawkes processes were established specifically for $K$ patterns, and the Granger causality can only be learned inside each pattern. Within each iteration, a Viterbi-path-based pattern reassignment algorithm and an EM-MLE-based parameter-updating algorithm were conducted alternately. In the parameter updating part, consider that $X_{t-1}$ and $X_t$ are more likely to be in the same pattern, Chen added a penalty term to help to put the adjacent sequences into the same pattern.

(4) Minimum Description Length: Following the Minimum Description Length (MDL) principle (Rissanen, 1998; Grünwald and Roos, 2019), Jalaldoust et al. conducted a trade-off between the goodness-of-fit and the model complexity (Jalaldoust et al., 2022). In detail, they partitioned the parameter space $\mathrm{\Theta }$ into $\{{\mathrm{\Theta }}_{\gamma }:\gamma \in \mathrm{\Gamma }\}$, defined a luckiness function $v:\mathrm{\Theta }\to ℝ$, and set the normalized maximum likelihood distribution for each model $\gamma \in \mathrm{\Gamma }$ to be:

(7) $$p_{v|\gamma }^{NML}(x)=\frac{ma{x}_{\theta \in {\mathrm{\Theta }}_{\gamma }}p(x|\theta )v(\theta )}{{\int }_{x\in 𝒳}ma{x}_{\theta \in {\mathrm{\Theta }}_{\gamma }}p(s|\theta )v(\theta )𝑑s}$$

The logarithm of the integral can be seen as the model complexity:

(8) $$COMP(M_{\gamma };v)=log{\int }_{x\in 𝒳}ma{x}_{\theta \in {\mathrm{\Theta }}_{\gamma }}p(s|\theta )v(\theta )𝑑s$$

Jalaldoust picked the optimized model ${\widehat{\gamma }}^{MDL}\in \mathrm{\Gamma }$ using

(9) $${\widehat{\gamma }}^{MDL}=\underset{\gamma \in \mathrm{\Gamma }}{\arg \min }{L}_v(\gamma ;x)=\underset{\gamma \in \mathrm{\Gamma }}{\arg \min }(-log\pi (\gamma )-r_v({\widehat{\theta }}_{v|\gamma }(x);x)+COMP(M_{\gamma };v))$$

where $\pi$ is a uniform distribution and $r_v({\widehat{\theta }}_{v|\gamma }(x);x)$ is the goodness-of-fit relevant to $p$, $v$, and ${\mathrm{\Theta }}_{\gamma }$.

Moreover, consider a one-to-one mapping from $\gamma \in \mathrm{\Gamma }$ to a $p\times p$ adjacent matrix of a causal graph within the set of all binary $p\times p$ matrices. By optimizing the 9, one can choose the most appropriate model from their predefined model family, hence, infer the Granger causal relationships between event types.

4.2.3. Methods for Wold processes

While most of the existing algorithms concerning discovering Granger causality from event sequences are based on Hawkes Processes, we can also model these relationships on another type of process - Wold Processes, which bear less complexity in nature. Suppose we denote ${\delta }_i=t_i-t_{i-1}$ to be the waiting time for $i$-th event from the occurrence of $(i-1)$-th event. Wold Processes are built upon a simple assumption that the current waiting time ${\delta }_i$ is only related to the closest past waiting time ${\delta }_{i-1}$. That is to say, the set $\{{\delta }_i,i\in \mathrm{N}\}$ forms a Markov chain. The inherent Markov property within the Wold processes makes them suitable for modeling the dynamics of certain complex systems. Besides, Figueiredo et al. (de Figueiredo et al., 2018)) have measured the correlation between ${\delta }_i$ and ${\delta }_{i-1}$ on certain datasets. The result shows that in most of their datasets, the median Pearson correlation is above 0.7, which is a sign of the adequacy of the Wold model. Accordingly, following Alve et al. (da Silva Alves et al., 2016) and Figueiredo et al. (de Figueiredo et al., 2018)’s idea, the intensity function can be performed as

$${\lambda }_{e_i}(t)={\mu }_{e_i}+\sum _{e_j\in E}\frac{{\alpha }_{e_i{e}_j}}{{\beta }_{e_j}+{\mathrm{\Delta }}_{e_i{e}_j}(t)}$$

based on the BuSca model. Here, ${\mu }_{e_i}$ is the base intensity as in Hawkes Processes. ${\mathrm{\Delta }}_{e_i{e}_j}(t)$ denote the time interval between the last $e_i$ type occurrence and $e_j$ type occurrence on time t. That is, if we define the closest $e_i$ type event before time $t$ happened at time $t_{e_i}$, and $e_j$ type before time $t_{e_i}$ happened at time $t_{e_j}$, there is ${\mathrm{\Delta }}_{e_i{e}_j}(t)=t_{e_i}-t_{e_j}$. Hence, the cross-type entry ${\sum }_{e_j\in E}\frac{{\alpha }_{e_i{e}_j}}{{\beta }_{e_j}+{\mathrm{\Delta }}_{e_i{e}_j}(t)}$ in our intensity function will be larger if ${\mathrm{\Delta }}_{e_i{e}_j}(t)$ decrease. This perfectly matches the fact that if type-$e_i$ events always happen just before type-$e_j$ events occur, we see a greater probability that $e_j$ has a certain effect on $e_i$. ${\alpha }_{e_i{e}_j}$ is the normalizing entry satisfying ${\sum }_{e_j\in E}{\alpha }_{e_i{e}_j}=1$, while ${\beta }_{e_j}$ is the base rate such that when the time interval ${\mathrm{\Delta }}_{e_i{e}_j}(t)$ between two types are infinitesimal at time $t$, the cross-type entry will converge to $\frac{{\alpha }_{e_i{e}_j}}{{\beta }_{e_j}}$.

The Granger causality can be learned through this Wold-based model by inspecting ${\alpha }_{e_i{e}_j}$. Specifically, if ${\alpha }_{e_i{e}_j}\ne 0$, it is considered that $e_j$ Granger-cause $e_i$. Since the approaches to learning the processes may not have an adequately sparse solution, Figueiredo tested the statistical significance of these possible Granger causal relationships and discard those with low significance. Moreover, the Wold-based model can be learned through MCMC, Expectation Maximization (EM) (de Figueiredo et al., 2018) and Variational Inference (Etesami et al., 2021) approaches. The task is to infer the parameters $\{{\alpha }_{e_i{e}_j},{\beta }_{e_j},{\mu }_{e_i}|\forall e_i,e_j\}$ in the intensity functions, which could reveal all the properties inside event sequences. Here we do not elaborate on the details of these learning methods.

4.2.4. Methods for Neural Point Processes

With the rapid development of neural networks, Neural Point Processes(NPPs) have gradually been utilized to model event sequences and infer causal relations. The core idea of these NPP algorithms is to use neural networks to infer the intensity function ${\lambda }_e(t)$. In specific, they encode an event sequence into the hidden state, during which they capture the feature of the sequence. Then, they use decoders to infer the future intensity function. And there are two major types of NPPs. One is based on the autoregressive(AR) model; its hidden states $h_i$ only update when an event occurs. The other follows the hypothesis that the hidden state $h(t)$ changes continuously in time. The continuous-time models hold advantages that they are natural and more suitable for estimating attributes at any time $t$ because of their continuous traits. Nonetheless, this flexibility comes at a cost. Continuous hidden models could suffer from a slower training speed compared with AR-based discrete models. This is because the evolution, likelihood, and sampling processes might demand numerical approximations. In this part, we will first give the basics of NPPs and then introduce how to learn the Granger causality in NPPs.

Basics of NPPs. The general process of using the AR-based NPP model to infer Granger causality (Zhang et al., 2020) is presented as follows. First, we embed each event into a vector $v_i=[\theta (t_i-t_{i_1});V^T{z}_i$], where $\theta (\cdot )$ is a predefined function, V is the embedding function for events’ type, $[\cdot ;\cdot ]$ could be concatenation and $z_i$ could be the one-hot coding for event type $u_i$. Then, we utilize a sequence encoder(e.g. LSTM or GRU) to encode $\{v_j;j\le i\}$ to $h_i=Enc(h_{i-1},v_i)$. Also, there exists a different encoding method in which the encoding is done independently for each i using, e.g., self-attention. This can better capture the long-range dependencies between events but have heavy computing complexity as well.

Next, we aim to decode the hidden state $h_i$ into the intensity function ${\lambda }_e(t)$. To do this, we need to make some assumptions about $\lambda$. For example, we could predict that the intensity function can be divided into the sum of some interaction-related and time-related functions ${\lambda }_e(t)={\sum }_{m=1}^M{a}_{em}{\kappa }_m(t)$ similar to what Xu did (Xu et al., 2016). Sequentially, we could only infer the $a_{em}$ entries since ${\kappa }_m(s)$ can be chosen from a large function family in which the functions can represent a wide variety of time-varied patterns. Hence,

$$\alpha :{ℝ}^{rank(h_i)}\to {ℝ}_{+}^{K\ast S}$$

is the corresponding decoder for this model. Here, $k$ and $S$ are the numbers of event types and basic functions correspondingly. However, the aforementioned method could not fit the continuous hidden state model. Under this circumstance, since $h(t)$ is continuous, thus carrying a lot more time-varying information than $h_i$ do, we can simply define the intensity as:

(10) $${\lambda }_e(t)=g_e(h(t))$$

Here, $g_e:{ℝ}^{rank(h_i)}\to R_{>0}$ is a non-linear function (e.g. softplus function) which maps $h(t)$ to the corresponding intensity function for event-type $e$ at time $t$.

As for the training process, most of the NPPs nowadays use the Maximum Likelihood Estimation(MLE) method as most of the traditional PP methods do. They take the negative log-likelihood of MLE as the objective function and use the neural network to optimize it. Besides, there exist alternative methods to use for learning the MLE. For example, if we set the objective as ${𝔼}_{X\sim p(X)}[f(X)]$, we can model the point process $P(X)$ using $f(X)$ by variational inference or reinforcement learning.

Inferring Granger Causal Relations from NPPs. When it comes to inferring the Granger causality, attribution methods need to be used (Zhang et al., 2020). That is because, in neural methods, most of the algorithms do not follow the parametrization in the Hawkes processes. On the contrary, their goal is to directly model the processes’ intensity function in order to loosen the strictness of Hawkes process and thus gain more accuracy. Since those intensity functions captured all the characteristics of event sequences, we should take full advantage of them. To do that, Zhang et al. first denoted $x_p=[t_1,e_1,\mathrm{\dots },t_p,e_p,t_{p+1}]$, $\underset{¯}{x_p}=[t_1,0,\mathrm{\dots },t_p,0,t_{p+1}]$ as the baseline input, and $f_k(x_p)={\int }_{t_p}^{t_{p+1}}{\lambda }_e(s)𝑑s$ as the impact function (Zhang et al., 2020). For each event type k, we have:

$$f_k(x_p)-f_k(\underset{¯}{x_p})=\sum _{q=1}^p{A}_q(f_k,x_p,\underset{¯}{x_p})$$

where $A_q(f_k,x_p,\underset{¯}{x_p})$ is the attribution(e.g. Integrated Gradients) for the event type of $z_q$. Hence, $A_q(f_k,x_p,\underset{¯}{x_p})$ can be regarded as the contribution of $z_j$-type events to the prediction of k-type events given the history $x_p$. Next, Zhang conducted a normalization on $A_q(f_k,x_p,\underset{¯}{x_p})$ as

$$Y_{e_i,e_j}=\frac{{\sum }_{s=1}^s{\sum }_{p=1}^{n_s}{\sum }_{q=1}^i𝕀(k_q^s=e_j)A_q(f_{e_i},x_p^s,\underset{¯}{x_p^s})}{{\sum }_{s=1}^s{\sum }_{p=1}^{n_s}𝕀(k_q^s=e_j)}$$

Consecutively, the Granger causality between $u_i$ and $u_j$-type events can be inferred from $Y_{e_i{e}_j}$. This method can measure not only the inhibitive causality but also the magnitude of the causality.

Interestingly, some other neural algorithms just model the intensity function as in Hawkes processes. They set $\mu$ and $\alpha$ as matrices and directly put $H$ and $A$ into the neural networks. Since the input structure is much easier, we could add other hypotheses, like the topological structure between events, and let the neural network(in this case, GCN) optimize $H$ and $A$ iteratively. Then, we could directly infer Granger causality from the matrix $A$.

4.3. Other Inferring Approaches

In this section, we will not directly model intensity functions in point processes. Instead, we focus on discovering the relationships between different processes (i.e., different types of events). To do that, we could utilize the Graphical Event Model mentioned before and loosen the assumption that each node follows the Hawkes Processes. Historically, Didelez et al. and Meek et al. first introduced the Graph Event Models to capture dependencies among events. Based on common graph methods, they assumed that an event type’s intensity function is only related to its parental type. GEMs capture dependencies between various types of events over time, providing a general framework to model the dependency in graph methods. Therefore, similar to the stationary as well as the discrete-time case, constraint-based and score-based approaches can be utilized.

4.3.1. Constraint-based methods

Just like the notion independence between different random variables, we can define process independence for point processes:

Definition 1.

(Didelez, 2008 (Didelez, 2008); Bhattacharjya et al., 2022 (Bhattacharjya et al., 2022))
For processes $X,Y,Z$, s.t. $Y\cap Z=0$, $X$ is a process independent of $Y$ given $Z$ if all events in $X$ have conditional intensities such that if historical information of events in $Z$ is known, then those events in process $Y$ do not provide any further information.

Meek et al. (Meek, 2014) and Bhattacharjya et al. (Bhattacharjya et al., 2022) introduced the notion of ${\delta }^{\ast }$-separation, which is based on $d$-separation but released its restriction of not having self-loops and made each self-loop independent of their own history. Then, they proposed a causal dependence assumption with ${\delta }^{\ast }$-separation analogous to the faithfulness assumption. Based on the causal dependence assumption, several constraint-based methods, such as the PC and max-min parents algorithms, are proposed to learn the causal relationships between different types. There are several Process Independence testers to choose from. For example, we have the NI tester:

(11) $$\text{NI score}=\frac{1}{2}\frac{{\sum }_z{({\lambda }_{x|y,z}-{\lambda }_{x|\widehat{y},z})}^2}{{\sum }_z{\lambda }_{x|y,z}+{\lambda }_{x|\widehat{y},z}){}^2}$$

where $y$ and $\widehat{y}$ indicate the parental state where $Y$ has or has not appeared in its window. We also have the LR tester:

(12) $$\text{LR score}=F_{{\chi }_{2^{|Z|}}^2}(-2[log{L}^{\ast }(X|Y,Z)-log{L}^{\ast }(X|Z)])$$

Here, $F(\cdot )$ is the cumulative distribution function of a chi-squared random variable with $2^{|Z|}$ degrees of freedom. Then, we apply a threshold $\tau$ for each tester, that is, when the score is less than $\tau$, there is no causal relationship between type $X$ and type $Y$.

4.3.2. Score-based methods

Similarly, there are score-based methods that can be applied to GEMs. Bhattacharjya et al. (Bhattacharjya et al., 2018) proposed PGEM - a model that assumed its intensity functions are only influenced by whether or not parent types happened in some recent time window. In addition, they used the BIC criterion on conditional intensities ${\lambda }_{x|u}$ to search for the optimal parent sets for each event type, that is, to infer the graph structure in their PGEM model. The graph structure is a representation of the causal relationships between different types of events.

4.3.3. Transfer Entropy

Recall that Transfer Entropy (TE) can be used to discover causal relationships in discrete-time cases. Here, we can also apply TE to event sequences (i.e., point processes) to identify our continuous-time causal relationships. Specifically, Spinney et al. (Spinney et al., 2017) constructed a continuous-time pairwise Transfer Entropy:

(13)

where $N_X$ is the number of events in the target process and $\tau$ is the length of time when there holds the corresponding intensity function and . The processes are independent when ${𝐓}_{Y\to X}=0$. We can define the conditional TE similarly. There are some existing consistent methods for estimating the continuous-time TE and its conditional form (Shorten et al., 2021).

5. Applications

Temporal causal discovery has been widely used in many areas, such as scientific endeavors (earth science (Runge et al., 2019a), neuroscience (Reid et al., 2019; Weichwald and Peters, 2021; Siddiqi et al., 2022), bioinformatics (Sachs et al., 2005)), industrial implementations (anomaly detection (Qiu et al., 2012), root cause analysis (Vuković and Thalmann, 2022; Liu et al., 2021; Assaad et al., 2023), business intelligence in online systems (Arabzadeh et al., 2018), video analysis (Yi et al., 2020)). Table 5 summarizes the application areas and corresponding studies. For scientific research, the learned causal relations should not usually be considered end results but rather starting points and hypotheses for further studies (Mäkelä et al., 2022). As a facilitator, causal discovery can play a supporting role in a multi-stage approach in an industrial setting (Vuković and Thalmann, 2022). In the rest part of this section, we will review three areas including earth science, anomaly detection and root cause, to explain these main workflows of incorporating temporal causal discovery into both scientific endeavors and industrial implementations, respectively.

Table 5. Major studies in temporal causal discovery applications.

Groups	Application areas	Studies
Scientific endeavors	Earth science	Climate change detection and attribution (e.g., (Lozano et al., 2009b)); Quantifying climate interactions (e.g., (Runge et al., 2014)); Latent driving force detection (e.g., (Trifunov et al., 2019; Shadaydeh et al., 2019)); Causality validation between temperature and greenhouse gases (e.g., (Van Nes et al., 2015)).
	Neuroscience	Dynamic causal models for neural connectivity (e.g., (Penny et al., 2004, 2010; Jafarian et al., 2020)); Granger causal models for neural connectivity (e.g., (Kaminski et al., 2001; Stokes and Purdon, 2017; Sheikhattar et al., 2018; Kim et al., 2011)); Causal inference from noninvasive brain stimulation (e.g., (Bergmann and Hartwigsen, 2021)).
	Bioinformatics	Modeling gene regulatory network (e.g., (Li and Ngom, 2013; Li et al., 2016; Verny et al., 2017; Patil and Vaida, 2022; Wu et al., 2022a)).
Industrial implementations	Anomaly detection	Causal structure as detection reference (e.g., (Qiu et al., 2012; Behzadi et al., 2017; Apte et al., 2021; Yang et al., 2022a)); Detection from imbalanced data (e.g., (Huang et al., 2020c)).
	Root cause analysis	Oscillation propagation tracing in the control loop (e.g., (Landman et al., 2014; Landman and Jämsä-Jounela, 2016; Chen et al., 2017; Lindner et al., 2018)); Alarm flood reduction (e.g., (Wang et al., 2015; Rodrigo et al., 2016; Wunderlich and Niggemann, 2017)); Industrial knowledge combined analysis (e.g., (Landman and Jämsä-Jounela, 2016; Cao et al., 2022; Thambirajah et al., 2009; Winchester et al., 2022)).
	Business intelligence in online systems	User interest prediction (e.g., (Arabzadeh et al., 2018; Hauffa et al., 2019)); Social media analysis (e.g., (Chang et al., 2013; Tsapeli et al., 2017; Kuzma et al., 2021; Chen et al., 2020); Online advertising (e.g., (Nuara et al., 2019; Yao et al., 2022b; Chu et al., 2020)); User-item interaction in recommendation (e.g., (Shang and Sun, 2020)); User activity modeling (e.g., (Li et al., 2017; Yao, 2022)).
	Video analysis	Video analysis and reasoning (e.g., (Yi et al., 2020; Li et al., 2020a)); Interpretable Gait Recognition (e.g., (Balazia et al., 2022)).
	Urban data analysis	Trajectory pattern mining (e.g., (Chu et al., 2016; Yang et al., 2022b)); Traffic flow prediction (e.g., (Li et al., 2015)); Visual urban and causal analytics (e.g., (Deng et al., 2022a)).
	Clinical data analysis	Causal chain discovery (e.g., (Wei et al., 2022)); Hypothesis testing (e.g., (Pandey, 2021)); Stable causal structure learning (e.g., (Rahmadi, 2019)).
	Signal processing	Blind source separation (e.g., (Testi et al., 2020; Testi and Giorgetti, 2021)); Compressed sensing (e.g., (Kathpalia and Nagaraj, 2022)).
	Financial analysis	Causal discovery for financial news (e.g., (Tetereva, 2018; Rambaldi et al., 2015)).
	Military	Battlefield sequential events analysis (e.g., (Li et al., 2022)).
	Robotics and dynamic control systems	Identifying causal structure (e.g., (Baumann et al., 2020)); Causal generalization (e.g., (Sheikhlar et al., 2021)).

Earth science and climate change research: Temporal causal discovery approaches have been widely used in the community of earth science and climate change research (Lozano et al., 2009b; Ebert-Uphoff and Deng, 2012; Runge et al., 2014; Van Nes et al., 2015; Hannart et al., 2016; Runge et al., 2019a; Trifunov et al., 2019). Climate is a complex and chaotic system, incorporating spatio-temporal information. Traditional climate models based on forward simulations have inherent limitation in describing such system due to uncertainties, simplifications, and discrepancies from observed data (Lozano et al., 2009b). Whereas, commonly used data centric methods such as lagged cross-correlation and regression analysis, aiming at deriving insights into interaction mechanisms between climate process, may lead to ambiguous conclusions in the field (Runge et al., 2014). To overcome the aforementioned issues, it’s reasonable to meaningfully characterize causal relationships among parameters of interest and make assertions. Specifically, spatio-temporal Granger modeling via group elastic net is proposed in (Lozano et al., 2009b) to conduct climate change detection and attribution, where the extreme-value theory to model and attribute extreme events in climate, such as severe heatwaves and floods. In (Runge et al., 2014), a graphical Granger model followed by a causal interaction strength measure is proposed to quantify the strength and delay of climate interactions and overcome the possible artifacts from vanilla correlation or regression methods. Another challenge is the existence of unobserved confounders, which may either lead to incorrect attribution or perform as a nonnegligible driving factor. A line of work (Trifunov et al., 2019; Shadaydeh et al., 2019) detect the latent driving force of abnormal event in climate by estimating the causal link intensity between confounded variables. Besides, in climate system some parameters of interest show strong coupling, thus impose difficulties for identification of causal orientation. The convergent cross mapping (CCM) technique, which is designed for strong coupling dynamic systems, is used in (Van Nes et al., 2015) to identify the causality between temperature and greenhouse gases, between which the statistical association is well documented while the causality is different to extract from the observed data. A recent overview of time series causal discovery in the earth system is also provided in (Runge et al., 2019a), where avenues for future work in both method developments and scientific endeavors are depicted.

Industrial temporal anomaly detection: In industrial systems, detecting anomalies in massive temporal data, which is derived from sensors, logs, physical measurements, system settings, etc, is meaningful while challenging. The anomalies can be roughly categorized into univariate anomaly, which has been extensively studied, and dependency anomaly, which is much more challenging to detect but common in real-world applications. As the challenges mainly come from high dimensions and complex dependency in data, methods (Qiu et al., 2012; Behzadi et al., 2017; Apte et al., 2021; Yang et al., 2022a; Huang et al., 2020c) based on temporal causal discovery have played a nonnegligible role in the dependency anomaly detection by providing efficient, robust and interpretable results. Causal discovery can facilitate the detection of the generative mechanisms of an underlying system. The key idea of this family of work is first to construct causal graphs from multivariate time series, and then detect anomalies according to the extracted causal relations. To be specific, in (Qiu et al., 2012; Behzadi et al., 2017), Granger graphical models are built on a reference data set and a testing data set respectively, the distribution differences (such as KL-divergence and Jensen-Shannon divergence) between the two learned models are computed as anomalous measures. In (Apte et al., 2021), the inferred relation based on Granger causality is termed causally anomalous if it violates the domain knowledge or the frequently observed forms. Recently, a causal perspective is also taken in (Yang et al., 2022a) to detect multivariate time series anomalies and leveraged in AIOps applications. In this work, the computation cost is reduced because instead of modeling joint distribution directly, it models factorized distribution modules from learned causal structures, where each corresponds to a local causal mechanism. Besides, as for the imbalanced flight data where the anomalous data points are rare, a time series classification method is proposed in (Huang et al., 2020c) based on nonlinear Granger causality learning.

Root cause analysis in manufacturing process: The root cause analysis is a vital task to ensure process safety and productivity in the industrial context, where the manufacturing processes are temporal and complex scenarios usually composed of multiple process units and a large number of feedback control loops. However, the acceptance of powerful ML methods in this field is hindered due to increasing requirements of fairness, accountability, and transparency (a.k.a., FAT principle (Shin and Park, 2019)), especially in sensitive-use cases (Vuković and Thalmann, 2022). To alleviate this issue, extracting knowledge such as causal relationships is paramount in this field. The last decade has witnessed the proliferation of the causal discovery methods for root cause analysis (Landman et al., 2014; Wang et al., 2015; Rashidi et al., 2018; Liu et al., 2021; Vuković and Thalmann, 2022). For instance, temporal causal discovery approaches such as Granger causality, transfer entropy, and their variants are leveraged to trace the oscillation propagation in the control loop (Landman et al., 2014; Landman and Jämsä-Jounela, 2016; Chen et al., 2017; Lindner et al., 2018). The reduction of alarm flood, which has been recognized as a major cause of industrial incidents, is another aspect of industrial root cause analysis. Among three typical nuisance alarms (i.e., repetitive alarms, standing alarms and consequence alarms (Henningsen and Kemmerer, 1995)), it’s challenging to suppress the consequence alarms and to provide a proper on the condition that the abnormality occurs and propagates. To identify all causal relations between alarms is of help (Hollender and Beuthel, 2007), and a line of work (Wang et al., 2015; Rodrigo et al., 2016; Wunderlich and Niggemann, 2017) leverages causal discovery approaches in this task. Besides, profound industrial knowledge, such as information flow and energy flow, can be combined with causal discovery to eliminate spurious relations (Landman and Jämsä-Jounela, 2016; Cao et al., 2022; Thambirajah et al., 2009).

6. Performance Evaluation

In this section, we give an overview of the benchmark datasets and evaluation metrics used in temporal causal discovery.

6.1. Datasets

We will briefly introduce some of the datasets used in temporal causal discovery, including MTS datasets and event-sequence datasets.

Datasets for MTS causal discovery range from health data to financial data. We discuss some of the commonly used datasets, which are publicly available and with the ground truth of causal graphs.

• Lorenz-96 simulated data: It’s a nonlinear model formulated in (Lorenz, 1996) to simulate climate dynamics. The continuous dynamics in a $d$-dimensional Lorenz model are given by $\frac{\partial {𝐱}_i^t}{\partial t}=-{𝐱}_{i-1}^t({𝐱}_{i-2}^t-{𝐱}_{i+1}^t)-{𝐱}_i^t+F,i\le i\le d$. The system dynamics become increasingly chaotic for higher values of forcing constant $F$. As a standard benchmark, it’s used by (Tank et al., 2022; Khanna and Tan, 2020; Marcinkevics and Vogt, 2021; Chu et al., 2020; Li et al., 2023).

• Linear VAR simulated data: Time series measurements are generated according to the linear VAR model. In (Tank et al., 2022; Khanna and Tan, 2020), it’s used to analyze methods’ performance when the true underlying dynamics are linear.

• CMU Human motion capture (CMU MoCap) data: It’s a data set from CMU MoCap database, containing data about joint angles, body position. Causal discovery methods can be leveraged to extract nonlinear dependencies between different regions of the body (Tank et al., 2022).

• DREAM-3 in Silico Network Inference Challenge: In DREAM-3 IN Silico Network Challenge (Prill et al., 2010), time-series data is simulated using continuous gene expression and regulation dynamics. Five gene regulation networks are to be inferred from gene expression level trajectories recorded. This dataset has been used to evaluate causal discovery algorithms in (Tank et al., 2022; Khanna and Tan, 2020).

• Blood-oxygenation-level dependent (BOLD) imaging data: In this dataset (Smith et al., 2011), time-ordered samples of the BOLD signals measure different brain regions of interest in human subjects. It’s generated using the dynamic causal modeling functional magnetic resonance imaging (fMRI) forward model. In (Khanna and Tan, 2020; Nauta et al., 2019), causal discovery methods are applied to estimate the connections in the human brain based on BOLD imaging data.

• Simulated financial time series: The dataset (Kleinberg, 2013) is created using factor model to describe portfolio’s return depending on three factors and a portfolio-specific error term. Thus the true relationships are known. It’s used by (Nauta et al., 2019).

As for event sequences, datasets range from online behavior to electricity. However, the true information on causal relationships is not accessible under all scenarios.

• MemeTracker: It is a dataset that captures online articles’ website, publication time, and all the hyperlinks within. This data set originally represents how a meme flow on different websites. The domain of the website and the publication time are considered an event type and its occurring time. And the hyperlinks between different websites can be seen as the ground truth of causal relationships. It’s used by (Achab et al., 2017; de Figueiredo et al., 2018; Zhang et al., 2020).

• IPTV viewing records: This dataset(Luo et al., 2014) records the user’s viewing behavior, i.e., what program and when they watch in the IPTV systems. The type of program and the time of watching the program can be deemed as an event type and its occurring time, respectively. It’s used by (Xu et al., 2016; Chen et al., 2022a; Zhang et al., 2020). However, ground-truth causal relationships are not included in this dataset.

• Power grid failure event data: This dataset includes abrupt changes in the voltage or current signals within Phasor Measurement Units (PMUs) as well as each PMU’s ID. The mission of the causal diagnosis task with this dataset is to infer the causalities within the grid (Idé et al., 2021). Since the network topology is not given out of privacy concerns, this is a non-ground-truth task.

• G-7 bonds: This dataset(Demirer et al., 2018) includes the daily return volatility of sovereign bonds of countries in the Group of Seven. The goal of dealing with this dataset is to discover the causal network underneath sovereign bonds (Jalaldoust et al., 2022). Expert knowledge from the domain can be deemed as ground truth.

6.2. Evaluation Metrics

In this part, we will explain different metrics used in the literature. Given the inferred probability of an edge $p(A_{ij})$ thresholded by $thre\in (0,1)$, the set of ground truth edges in causal graph $E_{GT}=\{(i,j):A_{ij}^{\ast }=1\}$, and the set of ground truth missing edges in causal graph $E_{MS}=\{(i,j):A_{ij}^{\ast }=0\}$, the definition and description of commonly used metrics is provided as follows:

• True Positive Rate (TPR): As a ratio of common edges found in the causal discovery results and the ground truth adjacencies over the total number of ground truth edges, the TPR metric is defined as $TPR=\frac{|\{(i,j):p(A_{ij})\ge thre\}\cap E_{GT}|}{|E_{GT}|}$.

• False Positive Rate (FPR): Similar to that in TPR, FPR refers to the ratio of common edges found in the causal discovery results and the ground truth missing adjacencies over the number of ground truth missing edges, which is defined as $FPR=\frac{|\{(i,j):p(A_{ij})\ge thre\}\cap E_{MS}|}{|E_{MS}|}$.

• Area Under the Receiver Operator Curve (AUROC): The Receiver Operator Curve (ROC) is defined as the ratio of TPR and FPR given the threshold $thre$ varies between $0$ and $1$. The area under the ROC (AUROC) is then widely used to assess the performance of causal discovery algorithms.

• Structural Hamming Distance (SHD): SHD is a metric describing the number of edge edition that need to be made to turn the discovered graph to its ground truth counterpart, which sums the number of missing edges, extra edges, and incorrect edges.

7. Discussion and New Perspectives

In this section, we first discuss challenges and practical considerations, including non-stationarity, heterogeneity, unobserved confounders, subsampling, and expert knowledge. Then, two new perspectives of temporal causal discovery are provided, which in our opinion will be a promising avenue for future research.

7.1. Challenges and Practical Considerations

Non-stationarity of data: We are often faced with non-stationarity in practical scenarios, where the probability distributions of temporal variables conditional on their causes or even the causal relations may change across time, especially for temporal data. In this condition, causal discovery approaches presuming a fixed causal model may give misleading results. Whereas, several types of research have shown that non-stationarity contains information for causal discovery (Tian and Pearl, 2001; Peters et al., 2016; Zhang et al., 2017; Huang et al., 2019). Thus, it’s important to properly tackle the non-stationarity in applications. Non-stationarity may result from the change of underlying systems and can be seen as a soft intervention (Korb et al., 2004) done by nature. Following this idea, a line of work (Zhang et al., 2017; Huang et al., 2020d) leverages a surrogate such as time and domain index to account for nonstationarity where the causal relations are changed, and the CD-NOD framework is proposed. Instead of leveraging informative non-stationarity to causal structure learning, another set of research focuses on modeling time-varying relationships (Gao and Yang, 2022). Besides, the approach for slowly varying non-stationary process, such as evolutionary spectral and locally stationary processes, is proposed in (Du and Xiang, 2020).

Heterogeneity of data: In causal discovery for practical applications, the heterogeneity of data lies in two levels: (1) The interacting temporal processes are heterogeneous (having different distributions), for instance, causally related meteorological observations from different stations are influenced by several major weather systems separately (Behzadi et al., 2019). (2) The underlying generating process changes across data sets or different domains (Glymour et al., 2019), for instance stock prices from different markets (Huang et al., 2020d) or individual behaviours in different paradigms (Chu et al., 2020). For the first condition where the heterogeneity exists among temporal variables, the inferred relations of the traditional causal discovery approaches, which have been designed for specific homogeneous data types, may be inaccurate. As a remedy, several variants of Granger causality, based on methods such as generalized linear models and minimum message length, are proposed in (Behzadi et al., 2017, 2019; Hlavácková-Schindler and Plant, 2020). For the second condition, a line of work (Zhang et al., 2017; Huang et al., 2020d) leverages the distribution shift from heterogeneity as a soft intervention to assist causal structure learning, which is similar to that in non-stationary data. Whereas, another line of causal discovery approaches (Chu et al., 2020; Löwe et al., 2022) in the second condition focuses on inductively modeling typical structure in heterogeneous data within an end-to-end framework.

Unobserved confounders: In practice, we are often met with cases where causal sufficiency is violated, i.e., there exist unobserved confounders. This challenging setting may lead to incorrect causal relations (Geiger et al., 2015). As summarized in Table 3, most temporal causal discovery approaches cannot handle unobserved confounders in a straightforward way. Several constraint-based approaches are designed without causal sufficiency and approaches Besides, unobserved confounders are modeled by applying a structural bias in (Löwe et al., 2022). Several recent studies termed as causal representation learning take a new perspective on unobserved confounders. It will be detailed in subsection (7.2.2).

Subsampling: In real-world applications, temporal data, especially time series, may be sampled at a rate lower than the rate of the underlying causal process due to the difficulties in data collection. An ordinary causal discovery algorithm for sub-sampled time series may lead to spurious causal relations and missed ones. Several remarks and approaches (Danks and Plis, 2013; Gong et al., 2015; Plis et al., 2015; Gong et al., 2017; Hyttinen et al., 2016; Tank et al., 2019) are proposed for this issue.

Expert knowledge: Expert knowledge can help the causal discovery process in practice. The approaches of fusing expert knowledge can be categorized into three types (Kitson et al., 2021): (1) Soft constraints: the learning process can be influenced by the knowledge (O’Donnell et al., 2006). (2) Hard constraints: the learnt structure must conform to the enforced requirements (i.e., conditions given with a probability $p=0$ or $p=1$). In (Asvatourian et al., 2020), hard constraints are leveraged in structure learning with a time dependant exposure. Studies in (Sun et al., 2021) add prior knowledge forbidding the existence of intra-slice dependencies, which is helpful to recover edges that are not explicitly encoded by the prior knowledge. (3) Interactive learning: the human input is leveraged in the learning process (Messaoud et al., 2009; Melkas et al., 2021; Zhu et al., 2022; Jin et al., 2021).

7.2. New perspectives

7.2.1. Extension in amortized and supervised paradigms

In the traditional paradigms, causal discovery methods mostly either treat observational data separately or train a distinct model for each individual. These methods do not make full use of the common structure across different samples or supervised information from the datasets whose causal structures are clearly explored, thus suffering from several issues such as the small sample challenge and lack the inductive capability. Recently, causal discovery is conducted in new paradigms to solve this problem. We can roughly categorize them into two groups: methods based on amortized modeling (Chu et al., 2020; Löwe et al., 2022), and methods based on supervised learning (Benozzo et al., 2017; Wang and Kording, 2022). We introduce them in this subsection, which we believe are a promising avenue for future research.

In amortized modeling, a global causal discovery framework is trained for individuals with different causal structures. As for scenarios with temporal data, these approaches have been detailed in 3.4.4 as the deep learning extension of Granger causality with inductive modeling. InGRA (Chu et al., 2020) leverages prototype learning to extract common causal structure while ACD (Löwe et al., 2022) proposes an encoder-decoder framework to conduct amortized causal discovery. These methods make full use of information from massive samples and are able to infer causal relations for newly arrived individuals, which are useful in real-world applications such as e-commerce, social network, and neuroimages.

Another line of work has predominately focused on treating the inference process as a black box and learning the mapping from sample data to causal graph structures via supervised learning. Here the label information is causal structure and can be easily accessed in synthetic datasets. Earlier work (Lopez-Paz et al., 2015; Ton et al., 2021) on learning causal relations by supervised learning is restricted to learning pairwise causal direction where the problem is cast into a classification task to distinguish between $X\to Y$ and $Y\to X$ by using observed samples. It’s later extended to discovery graph structure in (Li et al., 2020b; Petersen et al., 2022). As the labeled information for training is often originated from synthetic data or real-world datasets which have been explored, the requirement of a supervised approach, in which the distributions of training and test data match or highly overlap, is not guaranteed. In (Ma et al., 2022; Ke et al., 2022), methods such as vicinal graph and meta-learning are leveraged in supervised causal discovery to tackle this ‘domain shift’ issue. For the temporal setting, a supervised estimation of Granger causality between time series is proposed in (Benozzo et al., 2017). As a recent advance, a method for learning causal discovery is proposed in (Wang and Kording, 2022) where the learned from large datasets with known causal relations outperform the algorithm in the traditional paradigm when testing on temporal datasets such as fMRI.

7.2.2. Extension in causal representation learning

Extracting the causes of particular phenomena whether explicitly or implicitly from a deep learning black box can be beneficial to the downstream tasks. The aforementioned causal discovery methods focus on inferring relations between observed variables, or start from the premise that the causal variables are given before hand. Although some approaches learn causal relations under unobserved variables. There exist real-world observations (e.g., sensor measurements, image pixels in video) which are not well structured to causal variables to begin with. As a generalization of causal discovery from observed variables, there has recently been a growing interest in causal representation learning (Locatello et al., 2020; Schölkopf et al., 2021; Yang et al., 2021), which aims at learning representation of causal factors in an underlying system. It estimates latent causal variable graphs from observations.

A line of works in causal representation learning identifies independent factors of variations based on disentanglement and Independent Component Analysis (ICA). At the heart of this methodology is the postulation of mutually independent latent factors. It’s hard to identify true latent variables, especially in general nonlinear cases. As a remedy, recent approaches (Locatello et al., 2020; Khemakhem et al., 2020; Hyvärinen and Morioka, 2017, 2016) leverage additional information in multiple views, auxiliary variables, or temporal structure, combined with deep learning methods like VAEs and contrastive learning. A connection between ICA and causality has been recently drawn in (Gresele et al., 2021; Monti et al., 2019). In the context of temporal data, the identifiability of causal variables from temporal sequences is discussed in latent temporal causal process estimation (LEAP) (Yao et al., 2022c). It first provides causal identifiability conditions in a nonparametric, nonstationary setting, and a parametric setting. Then it proposes a learning framework to extract latent causal relations, which extends VAE with a learned causal process network by enforcing the assumed conditions. The non-stationary noise, modeled by flow-based estimators, can be viewed as a soft intervention to aid identification. In line with LEAP, subsequent works (Yao et al., 2022a) extend the identification theory to a more general case.

Another line of work leverage intervention and data augmentation to help to identify latent causal relations. Under data augmentation, it’s demonstrated in (von Kügelgen et al., 2021) that common contrastive learning methods can block-identify causal variables that remain unchanged. For the temporal setting, CITRIS (Lippe et al., 2022b) is proposed. It’s a VAE framework learning causal representation where latent causal factors have possibly been interved on. By using intervention target information for identification, CITRIS is devoid of suffering from functional or distributional form constraints. Besides, causal factors in CITRIS are considered as either scalars or potentially multidimensional vectors, which is more practical in complex scenarios. Along this line of work, instantaneous causal relations are extracted in iCITRIS (Lippe et al., 2022a).

8. Conclusion

Causal discovery in temporal data is fundamental to understanding the dynamics and estimating the causal effects of interest. This article reviews two categories of temporal causal discovery: multivariate time series causal discovery, and event sequence causal discovery. Multivariate time series causal discovery can be categorized into four groups, including constraint-based, score-based, FCM-based, and Granger causal model. Main ideas and recent advances for each type are reviewed. For causal discovery in event sequence, we can classify these algorithms into constraint-based, score-based, and Granger causal models, which are in accordance with multivariate time series causal discovery. We note that Granger causal models are especially well-developed for event sequence due to a natural match-up between Granger causality and Hawkes processes. To bridge the gap between abundant temporal causal discovery algorithms with real-world impacts, we introduce several major studies including scientific endeavors and industrial implementations. We also provide an extensive list of resources, including datasets and metrics, which can be used as a guideline for future research in this field. Whilst many algorithms are offered with theoretical or empirical guarantees, the quality of the inferred relations is dependent on many issues, including non-stationarity, heterogeneity, unobserved confounders, subsampling and expert knowledge. We discuss these challenges and practical considerations. Lastly, we introduce new perspectives of causal discovery, where avenues for future work in amortized modeling, supervised learning, and causal representation learning are depicted.