Causal Discovery from Temporal Data: An Overview and New Perspectives

Most methods in this subclass are temporal variants of the Peter-Clark (PC) algorithm [171], which assumes causal sufficiency. Endeavors [6, 20, 155, 157] for MTS data extend PC’s ability to model temporal context by considering autocorrelation, non-stationarity, and high dimensionality.

(1) PCMCI and its Variants. Runge et al. [157] propose PCMCI, which starts by constructing a partially connected graph, where all pairs of nodes \((x^{t-k}_i, x^t_j)\) are directed as \(x^{t-k}_i \rightarrow x^t_j\) accounting for temporal priority. The algorithm comprises two stages: (i) PCMCI eliminates redundant edges based on conditional independence and further removes edges by assuming consistency across time. (ii) The application of Momentary Conditional Independence (MCI) [147] addresses auto-correlation and prevents the occurrence of spurious correlations. In this context, MCI serves as a measure that conditions on the parents of \(x^t_j\) and \(x^{t-k}_i\) while assessing \(x^{t-k}_i \not\! \perp \!\!\! \perp x^t_j | Pa(x^t_j) \backslash \lbrace x^{t-k}_i\rbrace , Pa(x^{t-k}_i)\). And \(\not\! \perp \!\!\! \perp\) denotes non-independence.

Recent years have witnessed notable improvement endeavors in PCMCI. In PCMCI+, Runge [155] expands PCMCI’s capability to model contemporaneous relations. In Bagged-PCMCI\(+\), Debeire et al. [43] enhance the stability of results through bootstrap sampling and assigns confidence to edges. In Regime-PCMCI, Saggioro et al. [158] focus on regime-dependent causal relations, and iteratively optimize regimes and causal graphs in each regime. In Mapped-PCMCI, Tibau et al. [182] extract causal relations from MTS at the grid level. In F-PCMCI, Castri et al. [26] employ a feature selection approach to enhance PCMCI’s performance in handling a large number of variables.

(2) oCSE and its Variants. Sun et al. [177] propose oCSE to guide computational and efficient MTS-CD. It builds upon the theoretical framework of Causation Entropy (CE) [176], an extension of Transfer Entropy (TE) [160] that enables the quantification of network relationships among multiple variables rather than just pairwise relations. Slightly different from PC, oCSE encompasses the past of all accessible nodes, rather than minimizing the size of the conditioning set by conditioning on all potential causes.

oCSE is a computational and sample-efficient algorithm. However, it assumes that the hidden dynamics follow a stationary first-order Markov process as the CE only models causal relations with time lags equal to one. Recently, Assaad et al. [6] propose PCGCE to extract extended summary causal graphs based on PC and the Greedy Causation Entropy, which is a variant of CE modeling more intricate historical information. In PC-PMIME, Arsac et al. [5] combine PC with an information-theoretic measure (i.e., PMIME) designed to detect direct coupling in time-series.

(3) CTPC and its Variants. Bregoli et al. [20] propose Continuous-Time PC (CTPC), which is the first constraint-based approach for the structure learning of CTBNs. On the basis of PC, CTPC consists of a proper set of statistics testing conditional independence for time-series data with continuous-time observations. Compared to its score-based counterpart [137], CTPC has better accuracy in terms of structure learning when variables in the CTBN have more than two values. Recently, Markov Blanket-based Continuous-Time PC (MB-CTPC) [185] adapts CTPC and only tests dependencies relevant to the Markov blanket of specified variables, which is helpful in reducing computational time in downstream tasks such as multivariate time-series classification.

Prior methodologies within this subclass involve temporal adaptations of the FCI algorithm [171], which does not assume causal sufficiency. However, recent years have witnessed the emergence of alternative methods that do not rely on the FCI premise. In this subsection, we will begin by providing a concise overview of FCI before exploring its variants for MTS-CD. Subsequently, we will review methods that operate independently of FCI.

FCI offers a broader applicability by accommodating hidden confounders. It’s demonstrated to be accurate, sound, complete, and consistent from a theoretical standpoint [210]. It accomplishes this property through (i) the partial ancestral graphs for representing the presence of hidden confounders, and (ii) the incorporation of supplementary rules [209] for edge orientations.

The variants of FCI extend its ability to model temporal information in MTS.

(1) ANLTSM and its Variants. Chu and Glymour [33] propose ANLTSM, which is a constraint-based additive nonlinear time-series model. It presumes the following data generation process to escape the curse of dimensionality for nonparametric conditional independence tests. The FCI algorithm is then leveraged to identify lagged and contemporaneous causal relations. ANLTSM imposes constraints on contemporaneous interactions, requiring them to be linear, and on hidden confounders, stipulating them to be linear and contemporaneous. Following that, there are subsequent endeavors, building upon a model-defined data generation process. In SVAR-FCI, Malinsky and Spirtes [125] assume that the data generation process for MTS adheres to a structural vector autoregressive (SVAR) model and takes into account contemporaneous influences.

(2) tsFCI. Entner and Hoyer [52] propose time-series FCI (tsFCI) for MTS, which directly applies FCI based on a time window. In detail, the initial time-series data is transformed into a collection of samples of the random vector using a sliding window of size K. Subsequently, by treating each element of the transformed vector as an individual random variable, FCI can be directly applied. Nevertheless, tsFCI overlooks selection variables and contemporaneous causal relationships.

In recent years, several methods have emerged that employ alternative approaches to handle hidden confounders, distinct from FCI.

(3) TS-ICD. Rohekar et al. [153] propose time-series iterative causal discovery (TS-ICD) which serves as an efficient approach for MTS-CD while accounting for hidden confounders. The core concept of TS-ICD encompasses two primary elements: (i) Drawing inspiration from [152], TS-ICD adopts an iterative refinement strategy for a graph recovered through preceding iterations with smaller conditioning sets, thereby enhancing statistical power and promoting stability in contrast to FCI. (ii) In the context of temporal data, TS-ICD first learns long-term temporal relations then short term ones, and contemporaneous relations are inferred lastly.

(4) CDANs. In order to address the presence of changing modules (i.e., non-stationarity), Ferdous et al. [54] propose a framework for causal discovery from autocorrelated and non-stationary time-series (CDANs). Their rationale is grounded on the notion that disregarding the presence of changing modules may give rise to erroneous or misleading causal links in MTS-CD. Specifically, CDANs incorporates a proxy variable that captures the variability of changing modules over time. Subsequently, kernel-based conditional independence tests are employed to identify contemporaneous relationships between the proxy variables and the observed variables.

(5) LPCMCI. In order to expand the applicability of PCMCI method to situations involving latent confounders, Gerhardus and Runge [60] propose an extension known as LPCMCI. It introduces novel concepts, namely middle marks and LPCMCI-PAGs, which provide an unambiguous causal interpretation to facilitate the orientation of edges in the presence of hidden confounders. Furthermore, it mitigates the issue of inflated false positives.

The intuition of score-based approaches is to find the best graph structure measuring its fitness to data [106, 159]. The graph structures in these approaches are often modeled as Bayesian Networks or their variants [56, 151, 184] for temporal data, which provide causal explanations under causal edge assumptions [141]. It consists of two elements: (i) model scoring, and (ii) model search.

Model Scoring. Score functions play a crucial role in ranking the inferred graphs based on data. Common objective functions can be categorized into two groups [103], i.e., the Bayesian scores and information-theoretic scores. The Bayesian scores are based on both goodness-of-fit and the incorporation of prior knowledge, including BDe scores, K2 scores, and so on [98]. The information-theoretic scores take the model complexity into account to avoid over-fitting, including BIC and AIC scores [22].

Model Search/Optimization. This process can be additionally divided into two stages, i.e., parent set identification and structure optimization [159]. In the first stage, parent set identification generates a ranked list of suitable candidate parent sets for each variable, with higher-scoring sets listed first. Most approaches concentrate on the second stage, which cast the problem of searching causal structure \(\mathcal {G}\) into an optimization problem using the aforementioned score functions S. As stated in [145], the ultimate goal is: \(\hat{\mathcal {G}} = \mathrm{argmin}_{\mathcal {G}\ \mathrm{over} \ \mathbf {x}} S(\mathcal {D},\mathcal {G})\), where the empirical data \(\mathcal {D}\) represents variables \(\mathbf {x}\). Traditionally, it involves a combinatorial graph search, which is known to be a time-consuming problem. And currently the state-of-the-art on structure optimization is represented by hill climbing and tabu search [37, 162]. Nonetheless, causal discovery algorithms remain vulnerable to the challenges posed by the curse of dimensionality, which can result in suboptimal solutions that fail to fully capture the true underlying causal relationships. Recently, in the context of NOTEARS [213], an algebraic breakthrough has been utilized to characterize the acyclicity constraint in structure learning. This advancement transforms the previously discrete combinatoric search problem into a continuously optimizing problem. Consequently, the convertible relationship \((\iff)\) can be formulated as follows:where \(\mathbf {A}\) denotes the adjacency matrix, and \(h(\cdot)\) function enforces acyclicity in structure learning. The original acyclicity constraint function is implemented as \(h(\mathbf {A}) = \mathrm{tr}(e^{\mathbf {A} \circ \mathbf {A} }) - d\), where \(\circ\) is the Hadamard product. The equation restricts that there exist no loop within any length of steps based on the Taylor expansion [213]. It relies on the augmented Lagrangian method to solve the continuous constrained optimization problem.

In the context of MTS, we review the score-based methods following a similar paradigm from combinatoric search to continuous optimization.

It’s non-trivial to conduct the combinatorial graph search as it’s reported to be an NP-hard problem [31, 174]. Researchers have developed various score-based approaches, including those for DBNs [41, 56] and CTBNs [36, 115, 137, 184].

For learning DBNs, Friedman et al. [56] propose Structural Expectation-Maximization (Structural-EM). EM is concerned with learning BN models when hidden variables are assumed to be present [55]. And Structural-EM can learn DBN given even partial observations of variables over time. Yu et al. [205] propose an influence score that advances the capability of DBN to identify the \(+/-\) sign, such as activation or repression in biological data, and the relative magnitude between interactive temporal variables. In [151], Robinson and Hartemink introduce a non-stationary dynamic Bayesian network specifically designed for scenarios where the assumption that time-series data is generated by a stationary process is violated. Grzegorczyk and Husmeier [69, 70] propose a series of works learning DBNs under non-stationary scenarios. In B&B, Campos and Ji [41] first cast the problem of structure learning in DBN into a corresponding augmented BN using structural constraints, then use a branch-and-bound approach to guarantee global optimality. The B&B algorithm exhibits two notable properties regarding its computational complexity. First, it is an anytime and exact method, whereby the algorithm can be terminated at any point to return the current best solution along with an upper bound on the globally optimal value. Second, the algorithm can be easily parallelized and can provide a linear speedup in the number of tasks.

For learning CTBNs, Nodelman et al. [137] first propose a score-based approach called Continuous-Time Search and Score (CTSS). CTSS defines a Bayesian score for evaluating different graph candidates, and leverages a greedy hill-climbing search to find the structure with a high score. As there are no acyclicity constraints for CTBNs, the heuristic search process is much more efficient than that for DBNs [36, 137]. For non-stationary data, Vila and Stella [184] propose non-stationary CTBN (nsCTBN) with score functions and the corresponding optimization algorithm covering different non-stationary scenarios in a systematical manner. Recently, Linzner and Koeppl [115] introduce conditional CTBN (cCTBN) for interventional data in an active learning manner.

The majority of approaches in this subclass are temporal adaptations derived from NOTEARS, as discussed in Section 3.2.1. In the context of time-series data, the reviewed methods enhance the abilities of NOTEARS by incorporating temporal context.

(1) DYNOTEARS. Pamfil et al. [140] propose DYNOTEARS, which captures linear causal relations from MTS based on a continuous optimization approach. A structural vector autoregressive model assumption is introduced, which reflects both contemporaneous and time-lagged causal effects:where K is the model order, \(\mathbf {u}\) is a vector of centered error variables. To guarantee the identifiability, the error terms \(\mathbf {u}^t\) are assumed either non-Gaussian or standard Gaussian to hold identifiability [145]. \(\mathbf {W}\) and \(\mathbf {A}\) are weighted adjacency matrices, which correspond to intra-slice edges (contemporaneous relationships) and inter-slice edges (time-lagged relationships), respectively. The procedure of structure learning revolves around minimizing the least-squares loss subject to an acyclicity constraint, and the optimization problem is defined as follows:where n is the number of timesteps, and \(\lambda _{\mathbf {W}}, \lambda _{\mathbf {A}}\) stand for weights of regularization. The acyclicity constraint is imposed on the contemporaneous relationships \(\mathbf {W}\). As a result, the window causal graph can be estimated by optimizing \(\mathbf {W}\) and \(\mathbf {A}\).

(2) NTS-NOTEARS. To extract both linear and nonlinear relations among variables, Sun et al. [178] propose NTS-NOTEARS which leverages 1D convolutional neural networks. The optimization procedure follows a similar way as DYNOTEARS. Besides, prior knowledge of variable dependencies can be transformed as additional optimization constraints.

(3) IDYNO and its Variants. To handle both observational and interventional MTS data, Gao et al. [57] propose an interventional extension of DYNOTEARS (IDYNO), which modifies the score function to handle interventional targets. In Temporal Causal Discovery based on Intervention (TECDI) [110], Li et al. propose scores for imperfect intervention.

Nevertheless, it is important to acknowledge the boundaries and limitations of continuous optimization methods discussed in [96, 136, 149], such as the impact of the data scale and the convergence criteria of the augmented Lagrangian method. It is recommended to consider these aspects for future advancements and applications of this methodology.

The methods in this subclass are temporal adaptations derived from the LiNGAM algorithm [165]. Prior to delving into the specific methods for MTS-CD, it is essential to provide a concise overview of LiNGAM.

Shimizu et al. [165] propose causal discovery based on Linear, Non-Gaussian, Acyclic Models (LiNGAM) in non-temporal setting, which has the following assumptions: (i) a linear data generation process, (ii) non-Gaussian disturbances, and (iii) no hidden confounders. In LiNGAM, the relations among observations can be formulated as \(\mathbf {x} = \mathbf {A}\mathbf {x} + \mathbf {u}\), where \(\mathbf {A}\) denotes the adjacency matrix of the causal graph. The equation can be rewritten as \(\mathbf {x}=\mathbf {B}\mathbf {u}\), where \(\mathbf {B} = (\mathbf {I}-\mathbf {A})^{-1}\). Then, LiNGAM considers \(\mathbf {x}\) as observed signals and \(\mathbf {u}\) as underlying independent components. To uncover hidden factors or sources that contribute to the observed signals, independent component analysis (ICA) [173] can be leveraged to estimate a linear transformation from \(\mathbf {u}\) to \(\mathbf {x}\), i.e. \(\mathbf {B}\). Then causal relationships \(\mathbf {A}\) can be derived based on \(\mathbf {B}\). MTS-CD methods based on ICA are reviewed. They extend LiNGAM’s ability to model temporal context, non-linearity, and time-varying situations.

(1) VAR-LiNGAM and its Variants. Hyvärinen et al. [86, 87] propose VAR-LiNGAM for MTS-CD. It assumes the utilization of a SVAR model as SCM, incorporating the non-Gaussian nature of causal process. It’s defined as \(\mathbf {x}^t = \sum _{k=0}^K \mathbf {A}^k \mathbf {x}^{t-k} + \mathbf {u}^t\), where \(\mathbf {A}^k\) is the adjacency matrix of the causality between the variables \(\mathbf {x}\) with time lag k. And \(\mathbf {u}^t\) are independent random variables modeling the exogenous influences, which are assumed to be non-Gaussian. To derive causal adjacency matrix \(\mathbf {A}\), VAR-LiNGAM takes a four-step optimization procedure similar to LiNGAM, i.e., (i) regression phase, (ii) residual computation, (iii) contemporaneous relation inference based on LiNGAM, (iv) delayed relation inference based on matrix transformation. It degenerates to the LiNGAM model if the autoregressive order is zero, i.e., \(K = 0\). The VAR-LiNGAM model is further extended to challenging scenarios, such as time-varying situations [83], cyclic structures [107].

(2) NCDH. Wu et al. [192] propose Nonlinear Causal Discovery via HM-NICA (NCDH), which leverages a nonlinear ICA algorithm as a measurement of nonlinear relationships. NCDH first leverages the combination of nonlinear ICA and hidden Markov model to separate latent noises [72]. As a remedy for the permutation uncertainty of ICA, a sequence of independence tests are conducted to determine the corresponding relations between the observed variables and the separated noises.

Except for ICA-based methods, several works in MTS-CD have emerged that take into account the effects of noise distortion or non-linearity in the causal functions.

(1) TiMINo and its Variants. Peters et al. [144] propose time-series Models with Independent Noise (TiMINo), which fits SCM for time-series data. TiMINo outputs either a summary time graph or remains undecided, which avoids leading to wrong causal conclusions when the model is misspecified or the data is insufficient.

(2) NBCB. Assaad et al. [8] propose Noise-Based/Constraint-Based (NBCB) approach. It first leverages an SCM-based method to infer potential causal relationships. The data generation process is defined as \(x^t_j = f_j(Pa(x^{\tau }_j)^{t - K} , \ldots , Pa(x^{0}_j)^{t}) + u^t_j\), and can be estimated with a Gaussian process. Spurious edges are subsequently eliminated by taking into consideration the set of potential parental relationships. Recently, Bystrova and Assaad et al. [23] propose hybrid approaches, combining SCM-based and constraint-based approaches. They aim at exploiting the best of two groups of endeavors and show robustness to assumption violations.

Statistical learning-based methods encompass two main categories: (i) model-free approaches and (ii) model-based approaches. Model-free approaches involve estimating conditional probability density functions [10] or utilizing transfer entropy [101] to infer Granger causal relationships. These approaches exploit the ability to model nonlinear dependencies, but encounter challenges when the number of variables increases due to the curse of dimensionality. In contrast, model-based methods offer computational efficiency in high-dimensional settings by employing parameterized data generative models, such as VAR models [4] or kernel models [150], to effectively capture the dynamics of the measured time-series.

Due to space limitations, we leave details of statistical learning-based methods, including model-free and model-based methods, in supplementary materials.

Neural networks possess the capability to capture intricate, nonlinear, and non-additive relationships between input and output variables. In order to infer Granger causal relationships within nonlinear and high-dimensional contexts, deep learning-based methodologies have been proposed. The derivation of Granger causality within neural networks is based on the following definitions [180].

(1) NGC and its Variants. Tank et al. [180] propose Neural Granger Causality (NGC), which leverages component-wise MLP (cMLP) or component-wise LSTM (cLSTM) with sparse input layer weights to infer nonlinear Granger causality. To be specific, in cMLP, each nonlinear output \(g_j\) is modeled with a separate MLP to distinguish the effects from inputs to outputs. With penalization, the input matrix of the first layer provides information for the selection of Granger causality. In the first layer of \(g_j(\cdot)\), we have \(h_1^t = \sigma (\sum _{k=1}^{\tau }W^{k}_1 \mathbf {x}^{t-k} + b_1)\). If the ith column of weight matrix \(W^k_1\) contains zeros for all time lag k, then time-series \(\mathbf {x}_i\) does not Granger-cause time-series \(\mathbf {x}_j\). The objective function consists of a data-fitting term and a sparse-inducing term:where sparse inducing penalty \(R(\cdot)\) is implemented through group lasso penalty, which extracts causal relations without requiring precise lag specification. As for cLSTM, it sidesteps the lag selection problem and the Granger causal information is inferred from the input weight of LSTM. In economy Statistical Recurrent Units (eSRU), Khanna and Tan [99] extend cLTSM’s ability in data scarcity scenarios based on statistical recurrent units.

(2) NN-GC and its variants. Montalto et al. [133] propose neural networks Granger causality (NN-GC), which is a feature selection procedure to identify the significant Granger causes in the MLP model. By greedily adding lagged components of predictor variable as input, an MLP is updated iteratively. A predictor variable is claimed to be a significant Granger cause of the target variable if at least one of its lagged components is considered as the input of neural networks when the procedure is terminated. In RNN-GC [190], the feature selection procedure is extended by replacing MLPs in NN-GC with gated RNN models, However, as this technique requires training and comparing many candidate models, it’s costly in high-dimensional settings.

(3) GVAR. Marcinkevics and Vogt [126] propose the generalized vector autoregression (GVAR) model for better interpretability in Granger causal discovery. It’s based on a variant of self-explaining neural networks [3]. The self-explaining neural networks are inherently interpretable models motivated by restricted properties, and follow the form: \(f(\mathbf {x}) = g(\theta (\mathbf {x})_1h(\mathbf {x})_1, \ldots , \theta (\mathbf {x})_kh(\mathbf {x})_k),\) where \(g(\cdot)\) and \(\mathbf {h}(\mathbf {x})\) denote a link function and the interpretable basis concepts, respectively. Combined with the VAR model for MTS, the GVAR model is given by \(\mathbf {x}^t = \sum _{l=1}^{\tau } \Psi _{\theta _l}(\mathbf {x}^{t-l})\mathbf {x}^{t-l} + \mathbf {u}^t,\) where \(\mathbf {u}^t\) denotes the additive noise term and \(\Psi _{\theta _l}: \mathbb {R}^d \rightarrow \mathbb {R}^{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^{t-l}_i \rightarrow x^t_j\) is measured by the \((i,j)\)-component of \(\Psi _{\theta _l}(\mathbf {x}^{t-l})\). The loss function consists of three terms: the MSE loss, a sparsity-inducing regularization, and the smooth penalty. And variability of signs \((+/-)\) over time can also be captured.

(4) SCGL. To infer Granger causality in high-dimensional settings, Xu et al. [193] propose the scalable causal graph learning (SCGL) framework based on matrix factorization and low-rank approximation. SCGL deconstructs 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 \mathbb {R}^{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 as in [32]. Therefore, SCGL approximates A via a k-rank decomposition, where \(k\lt d\). The low-rank approximation reduces the noise influence in causal discovery and provides interpretability in downstream time-series analysis [85].

(5) TCDF. Nauta et al. [135] propose the temporal causal discovery framework (TCDF), which utilizes attention-based dilated CNN to infer Granger causality. TCDF contains d independent attention-based CNNs for different target variables \((X_1,\ldots ,X_d)\). 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 toward the latter. TCDF evaluates variable importance and identifies significant causal links based on permutation. It 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.

(6) InGRA. In inductive Granger causal modeling (InGRA), Chu et al. [34] combine the Granger causal attention in [161] with prototype learning to infer Granger causality in heterogeneous MTS data. To be specific, the Granger causal attention mechanism is first leveraged to compute contributions of each variable toward prediction. Due to data limitation, learned Granger causal attention is not robust enough. InGRA then leverages prototype learning, of which the key idea is to compute the similarities of new samples to prototypical cases, to detect common causal structures.

(7) ACD. Löwe et al. [122] propose amortized causal discovery (ACD), which aims at training a model to infer Granger causal structures across different samples with various causal structures. Its implicit assumption is that different samples share similar dynamics. ACD is an encoder-decoder framework. The encoder function is defined to infer Granger causal relations of the input sample, and the decoder function forecasts the next time-step given the learned causal relations. As a result, the causal relations of previous unseen samples can be inferred without refitting the model.

(8) CR-VAE. Li et al. [109] propose causal recurrent variational autoencoder (CR-VAE) 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 An error-compensation module is leveraged to capture instantaneous effects. CR-VAE is not only able to infer Granger causality but also conduct the data-generating process in a transparent manner benefiting from the learned causal matrix.

Causality from time-series can be quantified with transfer entropy [160], a metric that measures the flow of information between two dynamic processes. For example, the transfer entropy from i to j (with time lag) can be expressed aswhere \(h(.|.)\) is the conditional entropy. Transfer entropy is capable of capturing nonlinear relationships without the need for specific models. Although the original definition of transfer entropy pertains to bivariate settings, the extension of this concept to multivariate scenarios is an area of interest, such as [6, 156, 176].

For Gaussian variables, Granger causality and transfer entropy are demonstrated to be equivalent [11]. To a certain extent, this family of methods, characterized by its model-free nature, shares commonalities with some constraint-based [7, 177] and Granger causality-based [101] approaches.

This particular family of methods is tailored for MTS systems that can be effectively represented by differential equations. The connection between differential equations and SCMs has been extensively explored in previous literature [143, 154]. In the context of discrete time, Voortman et al. [186] propose difference-based causality learner (DBCL), which can be regarded as a constrained version of dynamic SCMs. In the case of continuous time, theoretical efforts have been made to establish a causal interpretation of dynamic systems using both Ordinary Differential Equations (ODEs) [154] and Stochastic Differential Equations (SDEs) [74, 132].

More recently, neural graphical model (NGM) [13] has been proposed to model dynamic causal systems where the MTS data is irregularly sampled based on Neural Ordinary Differential Equations (neural-ODE). To be specific, the underlying causal system of interest can be represented as a dynamic structural model as follows:where \(\mathbf {w}(t)\) is a d-dimensional standard Brownian motion, \(\mathbf {x}_0\) is a Gaussian random variable independent of \(\mathbf {w}(t)\), and the function \(\mathbf {f}\) describes the causal graph \(\mathcal {G}\). NGM is a learning algorithm based on penalized neural-ODE. It’s suited to irregularly sampled MTS.

This family of methods is well-suited for situations characterized by the dynamic coupling of nonlinear systems. Sugihara et al. [175] propose convergent cross mapping (CCM), which infers causality based on Takens’ state-space theory [179]. In detail, given two time-series \(\mathbf {x}_1^t\) and \(\mathbf {x}_2^t\), the attractor manifolds \(\mathcal {M}_{x_1}, \mathcal {M}_{x_2}\) are first reconstructed. Secondly, causality can be detected by measuring the correspondence between \(\mathcal {M}_{x_1}\) and \(\mathcal {M}_{x_2}\), by testing whether one manifold preserves every local neighborhood defined in the other. Recent advancements [21] have been made to enhance CCM in challenging scenarios.

As a distinct type of point process, we commence by presenting the fundamental aspects of the Hawkes process, encompassing the formulation of the intensity function and the inferential procedures employed to ascertain Granger causality based on the intensity function. Subsequently, we delve into a comprehensive exploration of the intricate methodologies associated with this domain.

The Hawkes process [76], classified as a point process, exhibits temporal dynamics characterized by either intensities that increase abruptly or decay gradually. Notably, this process adheres to a predetermined structure of the intensity function, denoted asIn this context, the term \(\mu _{e_i}\) refers to the baseline intensity, which remains unaffected by endogenous events and therefore remains constant throughout the temporal domain. Conversely, the function \(\phi _{e_ie_j}(s)\) serves as the impact function, quantifying the attenuation of the influence imparted by preceding type-\(e_j\) events on subsequent occurrences of type-\(e_i\) events. Essentially, the impact function encapsulates the internal intensity transfer from \(e_j\) to \(e_i\). Given the resemblance observed in the definitions of \(\phi\) and Granger causality, it is plausible to deduce Granger causality directly through the examination of \(\phi\). To be specific, the following proposition [51] holds:where \(\mathbb {R}\) indicates the set of all real numbers. Consequently, it becomes possible to formulate \(\phi _{e_ie_j}(t)\) or its modified versions as a means to deduce Granger causality from event sequence data modeled by the Multivariate Hawkes process. Notably, numerous studies [2, 25, 27, 81, 88, 91, 93, 111, 194, 206, 214] have emerged within this domain in recent years.

(1) MLE-SGLP. Xu et al. [194] propose MLE-based algorithm with Sparse-Group-Lasso and Pairwise (MLE-SGLP) similarity constraints, which learns Granger causality from Multivariate Hawkes processes with a maximum likelihood estimator and a sparse-group-lasso regularizer. Nonetheless, the MLE-SGLP may encounter limitations due to its elevated computational complexity and limited modeling capability. Parameterization strategies on intensity functions and regularization methods on likelihood functions can be utilized to augment the modeling capacity. While certain techniques such as Generalized Method of Moments and event sequence separation can be leveraged to reduce the computational complexity.

(2) NPHC. Achab et al. [2] argue that causality in the Multivariate Hawkes process can be inferred without estimating the shape of the impact function and they propose Non Parametric Hawkes Cumulant (NPHC). Instead of the detailed form of the intensity function, NPHC infers Granger causality based on its integrals. NPHC has a competitive computational complexity compared with methods such as MLE-SGLP, This is attributed to NPHC’s innovative estimation technique, which focuses on the integrals of intensity functions rather than their explicit forms, thereby optimizing computational efficiency.

(3) \(L_0\)Hawkes. Idé et al. [88] claim that EM-based MLE algorithms with \(L_1\)-regularization cannot offer sparse solutions for Hawkes process mathematically. Hence, they propose \(L_0\)Hawkes, which is an \(L_0\)-regularized EM-based MLE algorithm to circumvent this problem. To be specific, the \(L_0\)-norm \(\Vert A\Vert _0\) indicates counts of non-zero entries in matrix \(\mathbf {A}\). This model exhibits a substantial enhancement in its modeling capacity. Detailed theoretical analysis can be found in [88].

(4) GC-nsHP. Chen et al. [27] propose Granger Causality for non-stationary Hawkes Process (GC-nsHP). It first assumes that a non-stationary Hawkes process can be approximated by utilizing a mixture of non-overlapping, and stationary processes. Subsequently, a dynamic programming-based algorithm is employed to partition the non-stationary process into multiple stationary sub-processes. At last, GC-nsHP leverages an EM-based algorithm to learn Granger causality in each sub-process. It demonstrates a substantial improvement in its computational complexity

(5) THP. Cai et al. [25] propose Topological Hawkes Process (THP) to address the Granger causal discovery problem by incorporating prior knowledge of the underlying topological relationships. For example, in telecommunication networks, alarms as events are propagated across different devices with an underlying topological structure.¹THP perceive the conventional Hawkes process as a temporal convolution and subsequently expand it into the time-graph domain by incorporating graph convolutional networks. Subsequently, an optimization method combining an EM-based algorithm with hill-climbing is employed. THP could significantly improve the model capacity when information about the underlying topological structure is available. However, it also entails a relatively high computational complexity.

(6) MDLH. Jalaldoust et al. [91] introduce the method MDLH, which employs Minimum Description Length (MDL) as the principle to frame the Granger causal problem as a model selection task. MDL, as a practical tool based on information compression, facilitates model selection. In the context of Granger causal discovery, MDLH leverages a Monte-Carlo method for inference. It is especially suitable for estimating Granger causalities from short event sequences.

(7) MMLH. Hlaváčková-Schindler et al. [81] propose the method MMLH, in which the Granger causality among dimensions of multivariate Hawkes processes are seen as a variable selection problem and is approached by minimum message length principle. Similarly as MDLH, it is suitable for estimating Granger causalities from short event sequences.

Wold process is a class of point process defined in terms of a Markovian joint distribution of inter-event times, and is used to model causality in event sequences with reduced complexity [42, 53]. In this part, we will begin by providing an overview of the Wold process, highlighting its differentiation from the Hawkes process. Subsequently, we will delve into a comprehensive examination of the methodologies employed for Granger causal discovery.

Contrasting with the Hawkes process, where the intensity function relies on the entire history of past events, the Wold process exhibits a Markov chain of finite memory for its inter-event time intervals. To be specific, we denote \(\delta _i=t_i-t_{i-1}\) as the waiting time for the ith event from the occurrence of the \((i-1)\)-th event. Wold Processes are constructed based on the fundamental assumption that the current waiting time \(\delta _i\) is solely dependent on the closest preceding waiting time \(\delta _{i-1}\). A commonly employed representation (Busca model [38]) of the intensity function for the Multivariate Wold process can be expressed aswhere \(\mu _{e_i}\) denotes the Poisson rate that characterizes the exogenous events, and \(\sum _{e_j\in E}\frac{\alpha _{e_ie_j}}{\beta _{e_j}+\Delta _{e_ie_j}(t)}\) represents the endogenous Wold rate, which is regulated by the time interval and the strength of interaction. These distinctive attributes contribute to a lower time complexity in estimating the Wold process. Moreover, it’s more suitable in contexts such as communication dynamics [53].

In recent years, studies [42, 53] have focused on inferring Granger causality from the Multivariate Wold process.

(1) Granger-Busca. Figueiredo et al. [42] introduce Granger-Busca, the pioneering approach to learning Granger causality from event sequences modeled by Multivariate Wold processes. Granger-Busca utilizes a Markov Chain Monte Carlo (MCMC) sampling algorithm for parameter optimization, and Granger causality between \(e_i\) and \(e_j\) is established when \(\alpha _{e_ie_j} \ne 0\) in Equation (18). Leveraging the inherent characteristics of the Wold process, Granger-Busca demonstrates significantly improved computational efficiency compared to Hawkes process-based methods. However, Granger-Busca imposes several constraints on the structure. For instance, it necessitates that each node possesses at least one outgoing edge.

(2) Var-Wold. In order to alleviate the assumptions made by Granger-Busca, Etesami et al. [53] introduce Var-Wold, a novel framework that utilizes a Bayesian approach (variational inference) to optimize parameters in multivariate Wold processes. Furthermore, owing to the Markovian nature of the intensity function, the Bayesian approach overcomes the issue of long memory present in the Hawkes process.

With the rapid advancement of neural networks, Neural Point Processes (NPPs) have gained traction in modeling event sequences and inferring causal relationships. In contrast to statistical point processes such as Hawkes processes and Wold processes, NPPs harness the learning capabilities of neural networks to effectively model sequences, often surpassing statistical point processes in terms of predictive performance [24]. At the core of these NPP algorithms lies the concept of employing neural networks to estimate the intensity function \(\lambda _e(t)\). Specifically, they encode the event sequence into a hidden state, capturing its underlying features, and subsequently employ decoders to infer the future intensity function.

Several studies [24, 212] have emerged within this domain in recent years.

(1) CAUSE. Zhang et al. [212] propose CAUSE, which infers Granger causality from NPPs. As NPPs are less interpretable and unable to directly reveal Granger causality, attribution methods are leveraged. This method can also measure the inhibitive causality, in which an event has the ability to inhibit or suppress the occurrence of subsequent events, and the magnitude of the causality.

(2) TNPAR. Cai et al. [24] introduce Topological Neural Poisson Autoregressive (TNPAR) model, which integrates prior knowledge of the underlying topological relationships to learn Granger causality, resembling the approach employed by THP. TNPAR comprises two distinct stages: (i) the generation stage, and (ii) the inference stage. During the generation stage, to model dynamics entailed in event sequences, TNPAR leverages a modified version of the neural Poisson process. Then, an amortized inference procedure is employed to estimate causality.

Table 3 is a comparison between all Granger causality discovery methods on event sequences. Different approaches’ ideas, model capacities, complexities, data volumes, and the scenarios where they can be most effectively utilized are being compared.

Furthermore, there have been studies [167, 168] focusing on the estimation of influences between events within NPPs, akin to the notion of Granger causality within the context of NPPs.

Constraint-based methods rely on the assessment of independence between processes, akin to the examination of conditional independence and d-separation in MTS-CD. Corresponding notions [45] and procedures can be identified in constraint-based methods [1, 17] applied to event sequences.

(1) \(\boldsymbol {\mu }\)-PC. Absar and Zhang [1] propose a constraint-based approach called \(\mu\)-PC, which leverages \(\mu\)-separation [131] and extends PC. \(\mu\)-separation offers a formal graphical depiction of statistical relationships within temporal data. To handle event sequence data, a novel technique for testing conditional independence is introduced, leveraging the Recurrent Marked Temporal Point Process (RMTPP) model [48].

(2) CB-PGEM. Bhattacharjya et al. [17] propose Constraint-Based Proximal Graphical Event Model (CB-PGEM), which introduces process independence tests for causal proximal graphical event model [18] and uses constraint-based structure discovery algorithms (e.g., max-min parents algorithm) for Bayesian networks to deploy these tests.

Score-based ES-CD methods revolve around the inference of causality through the optimization of score functions, a domain that has witnessed several studies [15, 16, 18].

(1) SB-PGEM and its Variants. Bhattacharjya et al. [18] propose Score-Based Proximal Graphical Event Model (SB-PGEM) which employs a score-based optimization approach to estimate the model parameters. The authors utilize BIC as the score to search for optimal time windows and parent sets for each event type. In a more recent work, Bhattacharjya et al. [16] extend SB-PGEM while introducing the Bayesian Gamma score to be the score function to address the challenges of limited data by incorporating background knowledge into the modeling process.

(2) CIR. Disregarding the conditional causal relationship, Bhattacharjya et al. [15] introduce a collection of scores, namely conditional intensity-based scores, along with associated algorithms to estimate the cause-effect associations between event pairs derived from extensive event datasets.

Given realistic demands and the reviewed methods, we assert the necessity and feasibility of a unified causal discovery framework for both MTS and ES data. In real-world scenarios, heterogeneous temporal data from diverse sources are encountered, such as MTS data from sensors and ES data extracted from text posts. Relying solely on either MTS-CD or ES-CD may introduce confounding factors. Existing MTS-CD and ES-CD methods discussed in Section 3 and Section 4 share a similar taxonomy and exhibit fundamental consistency, paving the way for a unified framework in future temporal causal discovery research.

The concept of temporal causal discovery, in its various forms, can empower numerous facets of machine learning and data mining.

We emphasize three new problems of temporal causal discovery that have not been well studied. (i) Online temporal causal discovery. Different from previous methods which train distinct models for individual samples, a comprehensive framework for global causal discovery is trained to accommodate individuals exhibiting diverse causal structures [34, 122]. They effectively utilize information from extensive datasets, enabling the inference of causal relationships for newly encountered individuals. (ii) Supervised temporal causal discovery. Distinct from the conventional frameworks, the inference process is treated as a black box to learn the mapping from sample data to causal graph structures [14, 146], which aligns with the paradigm with most artificial intelligence tasks. (iii) Causal representation learning for temporal data. The aforementioned causal discovery methods focus on inferring relations between observed variables, or start from the premise that the causal variables are given beforehand. As a generalization of temporal causal discovery from observed variables, there has recently been a growing interest in temporal causal representation learning [116, 202], which aims at the learning representation of causal factors in an underlying dynamic system. It estimates latent causal variable graphs from observations.

Recently, we have witnessed great breakthroughs of Large Language Models (LLMs) in many domains[113, 118]. It reveals the synergistic potentials between causal discovery and LLMs in two aspects: (i) How do LLMs affect causal discovery research? (ii) How can causal methods benefit the advances of LLMs?

Owing to LLMs, extant causal discovery methods can be augmented by prior knowledge extracted from a large-scale human corpus dataset. In [100], Kıcıman et al. find that GPT-3.5/4 can get competitive results in the causal discovery task by merely leveraging label information. However, researchers [208] argue that current LLMs relying on human knowledge may not be enough to offer new insights that are beyond limitations of human intuition. Therefore, recent researches [120] reveal that LLMs open a door for the collaborative improvement between data-driven insights from causal discovery algorithms and causal knowledge from human corpus dataset.

Owing to causal research, opportunities are presented to enhance current LLMs in two key aspects. Firstly, it fosters the trustworthiness of LLMs [119, 208] by eliminating bias and generating interpretable results. Secondly, by integrating causal-inspired approaches, LLMs’ capacity for causal reasoning [94, 208] can be further advanced.