INTERPRETABLE MODELS FOR GRANGER CAUSALITY USING SELF-EXPLAINING NEURAL NETWORKS
Year: 2021

Abstract:
Exploratory analysis of time series data can yield a better understanding of complex dynamical systems. Granger causality is a practical framework for analysing interactions in sequential data, applied in a wide range of domains. In this paper, we propose a novel framework for inferring multivariate Granger causality under nonlinear dynamics based on an extension of self-explaining neural networks. This framework is more interpretable than other neural-network-based techniques for inferring Granger causality, since in addition to relational inference, it also allows detecting signs of Granger-causal effects and inspecting their variability over time. In comprehensive experiments on simulated data, we show that our framework performs on par with several powerful baseline methods at inferring Granger causality and that it achieves better performance at inferring interaction signs. The results suggest that our framework is a viable and more interpretable alternative to sparse-input neural networks for inferring Granger causality.

heading: 1 INTRODUCTION
text: Granger causality (GC) (Granger, 1969) is a popular practical approach for the analysis of multivariate time series and has become instrumental in exploratory analysis (McCracken, 2016) in various disciplines, such as neuroscience (Roebroeck et al., 2005), economics (Appiah, 2018), and climatology (Charakopoulos et al., 2018). Recently, the focus of the methodological research has been on inferring GC under nonlinear dynamics (Tank et al., 2018; Nauta et al., 2019; Wu et al., 2020; Khanna & Tan, 2020; Löwe et al., 2020), causal structures varying across replicates (Löwe et al., 2020), and unobserved confounding (Nauta et al., 2019; Löwe et al., 2020).
To the best of our knowledge, the latest powerful techniques for inferring GC do not target the effect sign detection (see Section 2.1 for a formal definition) or exploration of effect variability with time and, thus, have limited interpretability. This drawback defeats the purpose of GC analysis as an exploratory statistical tool. In some nonlinear interactions, one variable may have an exclusively positive or negative effect on another if it consistently drives the other variable up or down, respectively. Negative and positive causal relationships are common in many real-world systems, for example, gene regulatory networks feature inhibitory effects (Inoue et al., 2011) or in metabolomics, certain compounds may inhibit or promote synthesis of other metabolites (Rinschen et al., 2019). Differentiating between the two types of interactions would allow inferring and understanding such inhibition and promotion relationships in real-world dynamical systems and would facilitate a more comprehensive and insightful exploratory analysis. Therefore, we see a need for a framework capable of inferring nonlinear GC which is more amenable to interpretation than previously proposed methods (Tank et al., 2018; Nauta et al., 2019; Khanna & Tan, 2020). To this end, we introduce a novel method for detecting nonlinear multivariate Granger causality that is interpretable, in the sense that it allows detecting effect signs and exploring influences among variables throughout time. The main contributions of the paper are as follows:
1. We extend self-explaining neural network models (Alvarez-Melis & Jaakkola, 2018) to time series analysis. The resulting autoregressive model, named generalised vector autore-
gression (GVAR), is interpretable and allows exploring GC relations between variables, signs of Granger-causal effects, and their variability through time.
2. We propose a framework for inferring nonlinear multivariate GC that relies on a GVAR model with sparsity-inducing and time-smoothing penalties. Spurious associations are mitigated by finding relationships that are stable across original and time-reversed (Winkler et al., 2016) time series data.
3. We comprehensively compare the proposed framework and the powerful baseline methods of Tank et al. (2018), Nauta et al. (2019), and Khanna & Tan (2020) on a range of synthetic time series datasets with known Granger-causal relationships. We evaluate the ability of the methods to infer the ground truth GC structure and effect signs.

heading: 2 BACKGROUND AND RELATED WORK

heading: 2.1 GRANGER CAUSALITY
text: Granger-causal relationships are given by a set of directed dependencies within multivariate time series. The classical definition of Granger causality is given, for example, by Lütkepohl (2007). Below we define nonlinear multivariate GC, based on the adaptation by Tank et al. (2018). Consider a time series with p variables: {xt}t∈Z+ = {( x1t x 2 t ... x p t )>} t∈Z+ . Assume that causal relationships between variables are given by the following structural equation model:
xit := gi ( x11:(t−1), ..., x j 1:(t−1), ..., x p 1:(t−1) ) + εit, for 1 ≤ i ≤ p, (1)
where xj1:(t−1) is a shorthand notation for x j 1, x j 2, ..., x j t−1; ε i t are additive innovation terms; and gi(·) are potentially nonlinear functions, specifying how the future values of variable xi depend on the past values of x. We then say that variable xj does not Granger-cause variable xi, denoted as xj 6−→ xi, if and only if gi(·) is constant in xj1:(t−1).
Depending on the form of the functional relationship gi(·), we can also differentiate between positive and negative Granger-causal effects. In this paper, we define the effect sign as follows: if gi(·) is increasing in all xj1:(t−1), then we say that variable x
j has a positive effect on xi, if gi(·) is decreasing in xj1:(t−1), then x
j has a negative effect on xi. Note that an effect may be neither positive nor negative. For example, xj can ‘contribute’ both positively and negatively to the future of xi at different delays, or, for instance, the effect of xj on xi could depend on another variable.
Granger-causal relationships can be summarised by a directed graph G = (V, E), referred to as summary graph (Peters et al., 2017), where V = {1, ..., p} is a set of vertices corresponding to variables, and E = { (i, j) : xi −→ xj } is a set of edges corresponding to Granger-causal relationships. Let A ∈ {0, 1}p×p denote the adjacency matrix of G. The inference problem is then to estimateA from observations {xt}Tt=1, where T is the length of the time series observed. In practice, we usually fit a time series model that explicitly or implicitly infers dependencies between variables. Consequently, a statistical test for GC is performed. A conventional approach (Lütkepohl, 2007) used to test for linear Granger causality is the linear vector autoregression (VAR) (see Appendix A).

heading: 2.2 RELATED WORK

heading: 2.2.1 TECHNIQUES FOR INFERRING NONLINEAR GRANGER CAUSALITY
text: Relational inference in time series has been studied extensively in statistics and machine learning. Early techniques for inferring undirected relationships include time-varying dynamic Bayesian networks (Song et al., 2009) and time-smoothed, regularised logistic regression with time-varying coefficients (Kolar et al., 2010). Recent approaches to inferring Granger-causal relationships leverage the expressive power of neural networks (Montalto et al., 2015; Wang et al., 2018; Tank et al., 2018; Nauta et al., 2019; Khanna & Tan, 2020; Wu et al., 2020; Löwe et al., 2020) and are often based on regularised autoregressive models, reminiscent of the Lasso Granger method (Arnold et al., 2007).
Tank et al. (2018) propose using sparse-input multilayer perceptron (cMLP) and long short-term memory (cLSTM) to model nonlinear autoregressive relationships within time series. Building on this, Khanna & Tan (2020) introduce a more sample efficient economy statistical recurrent unit (eSRU) architecture with sparse input layer weights. Nauta et al. (2019) propose a temporal causal discovery framework (TCDF) that leverages attention-based convolutional neural networks to test for GC. Appendix B contains further details about these and other relevant methods.
Approaches discussed above (Tank et al., 2018; Nauta et al., 2019; Khanna & Tan, 2020) and in Appendix B (Marinazzo et al., 2008; Ren et al., 2020; Montalto et al., 2015; Wang et al., 2018; Wu et al., 2020; Löwe et al., 2020) focus almost exclusively on relational inference and do not allow easily interpreting signs of GC effects and their variability through time. In this paper, we propose a more interpretable inference framework, building on self explaining-neural networks (Alvarez-Melis & Jaakkola, 2018), that, as shown by experiments, performs on par with the techniques described herein.

heading: 2.2.2 STABILITY-BASED SELECTION PROCEDURES
text: The literature on stability-based model selection is abundant (Ben-Hur et al., 2002; Lange et al., 2003; Meinshausen & Bühlmann, 2010; Sun et al., 2013). For example, Ben-Hur et al. (2002) propose measuring stability of clustering solutions under perturbations to assess structure in the data and select an appropriate number of clusters. Lange et al. (2003) propose a somewhat similar approach. Meinshausen & Bühlmann (2010) introduce the stability selection procedure applicable to a wide range of high-dimensional problems: their method guides the choice of the amount of regularisation based on the error rate control. Sun et al. (2013) investigate a similar procedure in the context of tuning penalised regression models.

heading: 2.2.3 SELF-EXPLAINING NEURAL NETWORKS
text: Alvarez-Melis & Jaakkola (2018) introduce self-explaining neural networks (SENN) – a class of intrinsically interpretable models motivated by explicitness, faithfulness, and stability properties. A SENN with a link function g(·) and interpretable basis concepts h(x) : Rp → Rk follows the form
f(x) = g (θ(x)1h(x)1, ..., θ(x)kh(x)k) , (2) where x ∈ Rp are predictors; and θ(·) is a neural network with k outputs. We refer to θ(x) as generalised coefficients for data point x and use them to ‘explain’ contributions of individual basis concepts to predictions. In the case of g(·) being sum and concepts being raw inputs, Equation 2 simplifies to
f(x) = p∑ j=1 θ(x)jxj . (3)
Appendix C lists additional properties SENNs need to satisfy, as defined by Alvarez-Melis & Jaakkola (2018).
A SENN is trained by minimising the following gradient-regularised loss function, which balances performance with interpretability:
Ly(f(x), y) + λLθ (f(x)) , (4) where Ly(f(x), y) is a loss term for the ground classification or regression task; λ > 0 is a regularisation parameter; and Lθ(f(x)) = ∥∥∇xf(x)− θ(x)>Jhx (x)∥∥2 is the gradient penalty, where Jhx is the Jacobian of h(·) w.r.t. x. This penalty encourages f(·) to be locally linear.

heading: 3 METHOD
text: We propose an extension of SENNs (Alvarez-Melis & Jaakkola, 2018) to autoregressive time series modelling, which is essentially a vector autoregression (see Equation 11 in Appendix A) with generalised coefficient matrices. We refer to this model as generalised vector autoregression (GVAR). The GVAR model of order K is given by
xt = K∑ k=1 Ψθk (xt−k)xt−k + εt, (5)
where Ψθk : Rp → Rp×p is a neural network parameterised by θk. For brevity, we omit the intercept term here and in following equations. No specific distributional assumptions are made on the additive innovation terms εt. Ψθk (xt−k) is a matrix whose components correspond to the generalised coefficients for lag k at time step t. In particular, the component (i, j) of Ψθk (xt−k) corresponds to the influence of xjt−k on x i t. In our implementation, we use K MLPs for Ψθk(·) with p input units and p2 outputs each, which are then reshaped into an Rp×p matrix. Observe that the model defined in Equation 5 takes on a form of SENN (see Equation 3) with future time series values as the response, past values as basis concepts, and sum as a link function.
Relationships between variables x1, ..., xp and their variability throughout time can be explored by inspecting generalised coefficient matrices. To mitigate spurious inference in multivariate time series, we train GVAR by minimising the following penalised loss function with the mini-batch gradient descent:
1
T −K T∑ t=K+1 ‖xt − x̂t‖22 + λ T −K T∑ t=K+1 R (Ψt) + γ T −K − 1 T−1∑ t=K+1 ‖Ψt+1 −Ψt‖22 , (6) where {xt}Tt=1 is a single observed replicate of a p-variate time series of length T ; x̂t =∑K k=1 Ψθ̂k (xt−k)xt−k is the one-step forecast for the t-th time point by the GVAR model; Ψt is
a shorthand notation for the concatenation of generalised coefficient matrices at the t-th time point:[ Ψθ̂K (xt−K) Ψθ̂K−1 (xt−K+1) ... Ψθ̂1 (xt−1) ] ∈ Rp×Kp; R (·) is a sparsity-inducing penalty term; and λ, γ ≥ 0 are regularisation parameters. The loss function (see Equation 6) consists of three terms: (i) the mean squared error (MSE) loss, (ii) a sparsity-inducing regulariser, and (iii) the smoothing penalty term. Note, that in presence of categorically-valued variables the MSE term can be replaced with e.g. the cross-entropy loss.
The sparsity-inducing term R(·) is an appropriate penalty on the norm of the generalised coefficient matrices. Examples of possible penalties for the linear VAR are provided in Table 4 in Appendix A. These penalties can be easily adapted to the GVAR model. In the current implementation, we employ the elastic-net-style penalty term (Zou & Hastie, 2005; Nicholson et al., 2017) R(Ψt) = α ‖Ψt‖1 + (1− α) ‖Ψt‖ 2 2, with α = 0.5.
The smoothing penalty term, given by 1T−K−1 ∑T−1 t=K+1 ‖Ψt+1 −Ψt‖ 2 2, is the average norm of the difference between generalised coefficient matrices for two consecutive time points. This penalty term encourages smoothness in the evolution of coefficients w.r.t. time and replaces the gradient penalty Lθ (f (x)) from the original formulation of SENN (see Equation 4). Observe that if the term is constrained to be 0, then the GVAR model behaves as a penalised linear VAR on the training data: coefficient matrices are invariant across time steps.
for the loss function in Appendix D.

heading: 3.1 INFERENCE FRAMEWORK
text: Once neural networks Ψθ̂k , k = 1, ...,K, have been trained, we quantify strengths of Grangercausal relationships between variables by aggregating matrices Ψθ̂k (xt) across all time steps into
summary statistics. We aggregate the obtained generalised coefficients into matrix S ∈ Rp×p as follows:
Si,j = max 1≤k≤K
{ medianK+1≤t≤T (∣∣∣∣(Ψθ̂k (xt))i,j ∣∣∣∣)} , for 1 ≤ i, j ≤ p. (7)
Intuitively, Si,j are statistics that quantify the strength of the Granger-causal effect of xi on xj using magnitudes of generalised coefficients. We expect Si,j to be close to 0 for non-causal relationships and Si,j 0 if xi → xj . Note that in practice S is not binary-valued, as opposed to the ground truth adjacency matrix A, which we want to infer, because the outputs of Ψθ̂k(·) are not shrunk to exact zeros. Therefore, we need a procedure deciding for which variable pairs Si,j are significantly different from 0.
To infer a binary matrix of GC relationships, we propose a heuristic stability-based procedure that relies on time-reversed Granger causality (TRGC) (Haufe et al., 2012; Winkler et al., 2016). The intuition behind time reversal is to compare causality scores obtained from original and time-reversed data: we expect relationships to be flipped on time-reversed data (Haufe et al., 2012; Winkler et al., 2016). Winkler et al. (2016) prove the validity of time reversal for linear finite-order autoregressive processes. In our work, time reversal is leveraged for inferring stable dependency structures in nonlinear time series.
Algorithm 1 summarises the proposed stability-based thresholding procedure. During inference, two separate GVAR models are trained: one on the original time series data, and another on timereversed data (lines 3-4 in Algorithm 1). Consequently, we estimate strengths of GC relationships with these two models, as in Equation 7, and choose a threshold for matrixS which yields the highest agreement between thresholded GC strengths estimated on original and time-reversed data (lines 5-9 in Algorithm 1). A sequence ofQ thresholds, given by ξ = (ξ1, ..., ξQ), is considered where the i-th threshold is an ξi-quantile of values in S. The agreement between inferred thresholded structures is measured (line 7 in Algorithm 1) using balanced accuracy score (Brodersen et al., 2010), denoted by BA (·, ·), equal-to the average of sensitivity and specificity, to reflect both sensitivity and specificity of the inference results. Other measures can be used for quantifying the agreement, for example, graph similarity scores (Zager & Verghese, 2008). In this paper, we utilise BA, because considered time series have sparse GC summary graphs and BA weighs positives and negatives equally. In practice, trivial solutions, such as inferring no causal relationships, only self-causal links or all possible causal links, are very stable. The agreement for such solutions is set to 0. Thus, the procedure assumes that the true causal structure is different from these trivial cases. Figure 6 in Appendix E contains an example of stability-based thresholding applied to simulated data.
Algorithm 1: Stability-based thresholding for inferring Granger causality with GVAR.
Input: One replicate of multivariate time series {xt}Tt=1; regularisation parameters λ and γ ≥ 0; model order K ≥ 1; sequence ξ = (ξ1, ..., ξQ), 0 ≤ ξ1 < ξ2 < ... < ξQ ≤ 1.
Output: Estimate Â of the adjacency matrix of the GC summary graph. 1 Let {x̃t}Tt=1 be the time-reversed version of {xt} T t=1, i.e. {x̃1, ..., x̃T } ≡ {xT , ...,x1} . 2 Let τ (X, χ) be the elementwise thresholding operator. For each component ofX , τ (Xi,j , χ) = 1, if |Xi,j | ≥ χ, and τ (Xi,j , χ) = 0, otherwise. 3 Train an order K GVAR with parameters λ and γ by minimising loss in Equation 6 on {xt}Tt=1 and compute S as in Equation 7. 4 Train another GVAR on {x̃t}Tt=1 and compute S̃ as in Equation 7. 5 for i = 1 to Q do 6 Let κi = qξi(S) and κ̃i = qξi(S̃), where qξ(X) denotes the ξ-quantile of X . 7 Evaluate agreement
ςi = 1 2 [ BA ( τ (S, κi) , τ ( S̃>, κ̃i )) + BA ( τ ( S̃>, κ̃i ) , τ (S, κi) )] .
8 end 9 Let i∗ = argmax1≤i≤Q ςi and ξ∗ = ξi∗ .
10 Let Â = τ (S, qξ∗(S)). 11 return Â.
To summarise, this procedure attempts to find a dependency structure that is stable across original and time-reversed data in order to identify significant Granger-causal relationships. In Section 4, we demonstrate the efficacy of this inference framework. In particular, we show that it performs on par with previously proposed approaches mentioned in Section 2.2.

heading: 3.1.1 COMPUTATIONAL COMPLEXITY
text: Our inference framework differs from the previously proposed cMLP, cLSTM (Tank et al., 2018), TCDF (Nauta et al., 2019), and eSRU (Khanna & Tan, 2020) w.r.t. computational complexity. Mentioned methods require training p neural networks, one for each variable separately, whereas our inference framework trains 2K neural networks. A clear disadvantage of GVAR is its memory complexity: GVAR has many more parameters, since every MLP it trains has p2 outputs. Appendix F provides a comparison between training times on simulated datasets with p ∈ {4, 15, 20}. In practice, for a moderate order K and a larger p, we observe that training a GVAR model is faster than a cLSTM and eSRU.

heading: 4 EXPERIMENTS
text: The purpose of our experiments is twofold: (i) to compare methods in terms of their ability to infer the underlying GC structure; and (ii) to compare methods in terms of their ability to detect signs of GC effects. We compare GVAR to 5 baseline techniques: VAR with F -tests for Granger causality1 and the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995) for controlling the false discovery rate (FDR) (at q = 0.05); cMLP and cLSTM (Tank et al., 2018)2; TCDF (Nauta et al., 2019)3; and eSRU (Khanna & Tan, 2020)4. We particularly focus on the baselines that, similarly to GVAR, leverage sparsity-inducing penalties, namely cMLP, cLSTM, and eSRU. In addition, we provide a comparison with dynamic Bayesian networks (Murphy & Russell, 2002) in Appendix I. The code is available in the GitHub repository: https://github.com/i6092467/GVAR.
4.1 INFERRING GRANGER CAUSALITY
We first compare methods w.r.t. their ability to infer GC relationships correctly on two synthetic datasets. We evaluate inferred dependencies on each independent replicate/simulation separately against the adjacency matrix of the ground truth GC graph, an example is shown in Figure 2. Each method is trained only on one sequence. Unless otherwise mentioned, we use accuracy (ACC) and balanced accuracy (BA) scores to evaluate thresholded inference results. For cMLP, cLSTM, and eSRU, the relevant weight norms are compared to 0. For TCDF, thresholding is performed within the framework based on the permutation test described by Nauta et al. (2019). For GVAR, thresholded matrices are obtained by applying Algorithm 1. In addition, we look at the continuously-valued inference results: norms of relevant weights, scores, and strengths of GC relationships (see Equation 7). We compare these scores against the true structure using areas under receiver operating characteristic (AUROC) and precision-recall (AUPRC) curves. For all evaluation metrics, we only consider offdiagonal elements of adjacency matrices, ignoring self-causal relationships, which are usually the easiest to infer. Note that our evaluation approach is different from those of Tank et al. (2018) and Khanna & Tan (2020); this partially explains some deviations from their re-
sults. Relevant hyperparameters of all models are tuned to maximise the BA score or AUPRC (if a model fails to shrink any weights to zeros) by performing a grid search (see Appendix H for details about hyperparameter tuning). In Appendix M, we compare the prediction error of all models on held-out data.
1As implemented in the statsmodels library (Seabold & Perktold, 2010). 2https://github.com/iancovert/Neural-GC. 3https://github.com/M-Nauta/TCDF. 4https://github.com/sakhanna/SRU_for_GCI.

heading: 4.1.1 LORENZ 96 MODEL
text: A standard benchmark for the evaluation of GC inference techniques is the Lorenz 96 model (Lorenz, 1995). This continuous time dynamical system in p variables is given by the following nonlinear differential equations:
dxi dt = ( xi+1 − xi−2 ) xi−1 − xi + F, for 1 ≤ i ≤ p, (8)
where x0 := xp, x−1 := xp−1, and xp+1 := x1; and F is a forcing constant that, in combination with p, controls the nonlinearity of the system (Tank et al., 2018; Karimi & Paul, 2010). As can be seen from Equation 8, the true causal structure is quite sparse. Figure 2 shows the adjacency matrix of the summary graph for this dataset (for other datasets, adjacency matrices are visualised in Appendix G). We numerically simulate R = 5 replicates with p = 20 variables and T = 500 observations under F = 10 and F = 40. The setting is similar to the experiments of Tank et al. (2018) and Khanna & Tan (2020), but includes more variables.
Table 1 summarises the performance of the inference techniques on the Lorenz 96 time series under F = 10 and F = 40. For F = 10, all of the methods apart from TCDF are very successful at inferring GC relationships, even linear VAR. On average, GVAR outperforms all baselines, although performance differences are not considerable. For F = 40, the inference problem appears to be more difficult (Appendix J investigates performance of VAR and GVAR across a range of forcing constant values). In this case, TCDF and cLSTM perform surprisingly poorly, whereas cMLP, eSRU, and GVAR achieve somewhat comparable performance levels. GVAR attains the best combination of accuracy and BA scores, whereas cMLP has the highest AUROC and AUPRC. Thus, on Lorenz 96 data, the performance of GVAR is competitive with the other methods.

heading: 4.1.2 SIMULATED FMRI TIME SERIES
text: Another dataset we consider consists of rich and realistic simulations of blood-oxygen-leveldependent (BOLD) time series (Smith et al., 2011) that were generated using the dynamic causal modelling functional magnetic resonance imaging (fMRI) forward model. In these time series, variables represent ‘activity’ in different spatial regions of interest within the brain. Herein, we consider R = 5 replicates from the simulation no. 3 of the original dataset. These time series contain p = 15 variables and only T = 200 observations. The ground truth causal structure is very sparse (see Appendix G). Details about hyperparameter tuning performed for this dataset can be found in Appendix H.2. This experiment is similar to one presented by Khanna & Tan (2020).
Table 2 provides a comparison of the inference techniques. Surprisingly, TCDF outperforms other methods by a considerable margin (cf. Table 1). It is followed by our method that, on average, outperforms cMLP, cLSTM, and eSRU in terms of both AUROC and AUPRC. GVAR attains a BA score comparable to cLSTM. Importantly, eSRU fails to shrink any weights to exact zeros, thus, hindering the evaluation of accuracy and balanced accuracy scores (marked as ‘NA’ in Table 2). This
experiment demonstrates that the proximal gradient descent (Parikh & Boyd, 2014), as implemented by eSRU (Khanna & Tan, 2020), may fail to shrink any weights to 0 or shrinks all of them, even in relatively simple datasets. cMLP seems to provide little improvement over simple VAR w.r.t. AUROC or AUPRC. In general, this experiment promisingly shows that GVAR performs on par with the techniques proposed by Tank et al. (2018) and Khanna & Tan (2020) in a more realistic and data-scarce scenario than the Lorenz 96 experiment.

heading: 4.2 INFERRING EFFECT SIGN
text: So far, we have only considered inferring GC relationships, but not the signs of Granger-causal effects. Such information can yield a better understanding of relations among variables. To this end, we consider the Lotka–Volterra model with multiple species ( Bacaër (2011) provides a definition of
the original two-species system ) , given by the following differential equations:
dxi
dt = αxi − βxi ∑ j∈Pa(xi) yj − η ( xi )2 , for 1 ≤ i ≤ p, (9)
dyj
dt = δyj ∑ k∈Pa(yj) xk − ρyj , for 1 ≤ j ≤ p, (10)
where xi correspond to population sizes of prey species; yj denote population sizes of predator species; α, β, η, δ, ρ > 0 are fixed parameters controlling strengths of interactions; and Pa(xi), Pa(yj) are sets of Granger-causes of xi and yj , respectively. According to Equations 9 and 10,
Figure 3: Simulated two-species Lotka–Volterra time series (top) and generalised coefficients (bottom). Prey have a positive effect on predators, and vice versa.
the population size of each prey species xi is driven down by ∣∣Pa(xi)∣∣ predator species (negative effects), whereas each
predator species yj is driven up by ∣∣Pa(yj)∣∣ prey populations (positive effects).
We simulate the multi-species Lotka–Volterra system numerically. Appendix K contains details about simulations and the summary graph of the time series. To infer effect directions, we inspect signs of median generalised coefficients for trained GVAR models. For cMLP, cLSTM, TCDF, and eSRU, we inspect signs of averaged weights in relevant layers. For VAR, we examine coefficient signs. For the sake of fair comparison, we restrict all models to a maximum lag of K = 1 (where applicable). In this experiment, we focus on BA scores for positive ( BApos ) and negative ( BAneg ) relationships. Appendix L provides another example of detecting effect signs with GVAR, on a trivial benchmark with linear dynamics.
Table 3 shows the results for this experiment. Linear VAR does not perform well at inferring the GC structure, however, its coefficient signs are strongly associated with true signs of relationships. cMLP provides a considerable improvement in GC inference, and surprisingly its input weights are informative about the signs of GC effects. cLSTM fails to shrink any of the relevant weights to zero; furthermore, the signs of its weights are not associated with the true signs. Although eSRU performs better than VAR at inferring the summary graph, its weights are not associated with effect signs at
all. TCDF performs poorly in this experiment, failing to infer any relationships apart from selfcausation. Our model considerably outperforms all baselines in detecting effect signs, achieving nearly perfect scores: it infers more meaningful and interpretable parameter values than all other models.
These results are not surprising, because the baseline methods, apart from linear VAR, rely on interpreting weights of relevant layers that, in general, do not need to be associated with effect signs and are only informative about the presence or absence of GC interactions. Since the GVAR model follows a form of SENNs (see Equation 2), its generalised coefficients shed more light into how the future of the target variable depends on the past of its predictors. This restricted structure is more intelligible and yet is sufficiently flexible to perform on par with sparse-input neural networks.
In addition to inferring the summary graph, GVAR allows inspecting variability of generalised coefficients. Figure 3 provides an example of generalised coefficients inferred for a two-species Lotka– Volterra system. Although coefficients vary with time, GVAR consistently infers that the predator population is driven up by prey and the prey population is driven down by predators. For the multispecies system used to produce the quantitative results, inferred coefficients behave similarly (see Figure 12 in Appendix K).

heading: 5 CONCLUSION
text: In this paper, we focused on two problems: (i) inferring Granger-causal relationships in multivariate time series under nonlinear dynamics and (ii) inferring signs of Granger-causal relationships. We proposed a novel framework for GC inference based on autoregressive modelling with selfexplaining neural networks and demonstrated that, on simulated data, its performance is promisingly competitive with the related methods of Tank et al. (2018) and Khanna & Tan (2020). Proximal gradient descent employed by cMLP, cLSTM, and eSRU often does not shrink weights to exact zeros and, thus, prevents treating the inference technique as a statistical hypothesis test. Our framework mitigates this problem by performing a stability-based selection of significant relationships, finding a GC structure that is stable on original and time-reversed data. Additionally, proposed GVAR model is more amenable to interpretation, since relationships between variables can be explored by inspecting generalised coefficients, which, as we showed empirically, are more informative than input layer weights. To conclude, the proposed model and inference framework are a viable alternative to previous techniques and are better suited for exploratory analysis of multivariate time series data.
In future research, we plan a thorough investigation of the stability-based thresholding procedure (see Algorithm 1) and of time-reversal for inferring GC. Furthermore, we would like to facilitate a more comprehensive comparison with the baselines on real-world data sets. It would also be interesting to consider better-informed link functions and basis concepts (see Equation 2). Last but not least, we plan to tackle the problem of inferring time-varying GC structures with the introduced framework.

heading: ACKNOWLEDGMENTS
text: We thank Djordje Miladinovic and Mark McMahon for valuable discussions and inputs. We also acknowledge Jonas Rothfuss and Kieran Chin-Cheong for their helpful feedback on the manuscript.

heading: A LINEAR VECTOR AUTOREGRESSION
text: Linear vector autoregression (VAR) (Lütkepohl, 2007) is a time series model conventionally used to test for Granger causality (see Section 2.1). VAR assumes that functions gi(·) in Equation 1 are linear:
xt = ν + K∑ k=1 Ψkxt−k + εt, (11)
where ν ∈ Rp is the intercept vector; Ψk ∈ Rp×p are coefficient matrices; and εt ∼ Np (0,Σε) are Gaussian innovation terms. Parameter K is the order of the VAR model and determines the maximum lag at which Granger-causal interactions occur. In VAR, Granger causality is defined by zero constraints on the coefficients, in particular, xi does not Granger-cause xj if and only if, for all lags k ∈ {1, 2, ...,K}, (Ψk)j,i = 0. These constraints can be tested by performing, for example, F -test or Wald test.
Usually a VAR model is fitted using multivariate least squares. In high-dimensional time series, regularisation can be introduced to avoid inferring spurious associations. Table 4 shows various sparsity-inducing penalties for a linear VAR model of order K (see Equation 11), described by Nicholson et al. (2017). Different penalties induce different sparsity patterns in coefficient matrices Ψ1,Ψ2, ...,ΨK . These penalties can be adapted to the GVAR model for the sparsity-inducing term R(·) in Equation 6.
B INFERRING GRANGER CAUSALITY UNDER NONLINEAR DYNAMICS
Below we provide a more detailed overview of the related work on inferring nonlinear multivariate Granger causality, focusing on the recent machine learning techniques that tackle this problem.
Kernel-based Methods. Kernel-based GC inference techniques provide a natural extension of the VAR model, described in Appendix A, to nonlinear dynamics. Marinazzo et al. (2008) leverage reproducing kernel Hilbert spaces to infer linear Granger causality in an appropriate transformed feature space. Ren et al. (2020) introduce a kernel-based GC inference technique that relies on regularisation – Hilbert–Schmidt independence criterion (HSIC) Lasso GC.
Neural Networks with Non-uniform Embedding. Montalto et al. (2015) propose neural networks with non-uniform embedding (NUE). Significant Granger causes are identified using the NUE, a feature selection procedure. An MLP is ‘grown’ iteratively by greedily adding lagged predictor components as inputs. Once stopping conditions are satisfied, a predictor time series is claimed a significant cause of the target if at least one of its lagged components was added as an input. This technique is prohibitively costly, especially, in a high-dimensional setting, since it requires training and comparing many candidate models. Wang et al. (2018) extend the NUE by replacing MLPs with LSTMs.
Neural Granger Causality. Tank et al. (2018) propose inferring nonlinear Granger causality using structured multilayer perceptron and long short-term memory with sparse input layer weights, cMLP and cLSTM. To infer GC, p models need to be trained with each variable as a response. cMLP and
cLSTM leverage the group Lasso penalty and proximal gradient descent (Parikh & Boyd, 2014) to infer GC relationships from trained input layer weights.
Attention-based Convolutional Neural Networks. Nauta et al. (2019) introduce the temporal causal discovery framework (TCDF) that utilises attention-based convolutional neural networks (CNN). Similarly to cMLP and cLSTM (Tank et al., 2018), the TCDF requires training p neural network models to forecast each variable. Key distinctions of the TCDF are (i) the choice of the temporal convolutional network architecture over MLPs or LSTMs for time series forecasting and (ii) the use of the attention mechanism to perform attribution. In addition to the GC inference, the TCDF can detect time delays at which Granger-causal interactions occur. Furthermore, Nauta et al. (2019) provide a permutation-based procedure for evaluating variable importance and identifying significant causal links.
Economy Statistical Recurrent Units. Khanna & Tan (2020) propose an approach for inferring nonlinear Granger causality similar to cMLP and cLSTM (Tank et al., 2018). Likewise, they penalise norms of weights in some layers to induce sparsity. The key difference from the work of Tank et al. (2018) is the use of statistical recurrent units (SRUs) as a predictive model. Khanna & Tan (2020) propose a new sample-efficient architecture – economy-SRU (eSRU).
Minimum Predictive Information Regularisation. Wu et al. (2020) adopt an informationtheoretic approach to Granger-causal discovery. They introduce learnable corruption, e.g. additive Gaussian noise with learnable variances, for predictor variables and minimise a loss function with minimum predictive information regularisation that encourages the corruption of predictor time series. Similarly to the approaches of Tank et al. (2018); Nauta et al. (2019); Khanna & Tan (2020), this framework requires training p models separately.
Amortised Causal Discovery & Neural Relational Inference. Kipf et al. (2018) introduce the neural relational inference (NRI) model based on graph neural networks and variational autoencoders. The NRI model disentangles the dynamics and the undirected relational structure represented explicitly as a discrete latent graph variable. This allows pooling time series data with shared dynamics, but varying relational structures. Löwe et al. (2020) provide a natural extension of the NRI model to the Granger-causal discovery. They introduce a more general framework of the amortised causal discovery wherein time series replicates have a varying causal structure, but share dynamics. In contrast to the previous methods (Tank et al., 2018; Nauta et al., 2019; Khanna & Tan, 2020; Wu et al., 2020), which in this setting, have to be retrained separately for each replicate, the NRI is trained on the pooled dataset, leveraging shared dynamics.

heading: C PROPERTIES OF SELF-EXPLAINING NEURAL NETWORKS
text: As defined by Alvarez-Melis & Jaakkola (2018), g(·), θ(·), and h(·) in Equation 2 need to satisfy:
1. g(·) is monotonic and additively separable in its arguments; 2. ∂g∂zi > 0 with zi = θ(x)ih(x)i, for all i; 3. θ(·) is locally difference-bounded by h(·), i.e. for every x0, there exist δ > 0 and L ∈ R s.t. if ‖x− x0‖ < δ, then ‖θ(x)− θ(x0)‖ ≤ L ‖h(x)− h(x0)‖;
4. {h(x)i}ki=1 are interpretable representations of x; 5. k is small.

heading: D ABLATION STUDY OF THE LOSS FUNCTION
text: We inspect hyperparameter tuning results for the GVAR model on Lorenz 96 (see Section 4.1.1) and synthetic fMRI time series (Smith et al., 2011) (see Section 4.1.2) as an ablation study for the loss function proposed (see Equation 6). Figures 4 and 5 show heat maps of BA scores (left) and AUPRCs (right) for different values of parameters λ and γ for Lorenz 96 and fMRI datasets, respectively. For the Lorenz 96 system, sparsity-inducing regularisation appears to be particularly important, nevertheless, there is also an increase in BA and AUPRC from a moderate smoothing penalty. For fMRI, we observe considerable performance gains from introducing both the sparsityinducing and smoothing penalty terms. Given the sparsity of the ground truth GC structure and
the scarce number of observations (T = 200), these gains are not unexpected. During preliminary experiments, we ran grid search across wider ranges of λ and γ values, however, did not observe further improvements from stronger regularisation. In summary, these results empirically motivate the need for two different forms of regularisation leveraged by the GVAR loss function: the sparsityinducing and smoothing penalty terms.

heading: E STABILITY-BASED THRESHOLDING: EXAMPLE
text: Figure 6 shows an example of agreement between dependency structures inferred on original and time-reversed synthetic sequences across a range of thresholds (see Algorithm 1). In addition, we plot the BA score for resulting thresholded matrices evaluated against the true adjacency matrix. As can be seen, the peak of stability agrees with the highest BA achieved. In both cases, the procedure described by Algorithm 1 chooses the optimal threshold, which results in the highest agreement with the true dependency structure (unknown at the time of inference).

heading: F COMPARISON OF TRAINING & INFERENCE TIME
text: To compare the considered methods in terms of their computational complexity, we measure training and inference time across three simulated datasets with p ∈ {4, 15, 20} variables and varying time series lengths. This experiment was performed on an Intel Core i7-7500U CPU (2.70 GHz × 4) with a GeForce GTX 950M GPU. All models were trained for 1000 epochs with a mini-batch size of 64. In each dataset, the same numbers of hidden layers and hidden units were used across all models. When applicable, models were restricted to the same order (K). Table 5 contains average training and inference time in seconds with standard deviations. Observe that for the fMRI and Lorenz 96 datasets, GVAR is substantially faster than cLSTM and eSRU.

heading: G GC SUMMARY GRAPHS OF SIMULATED TIME SERIES

heading: H HYPERPARAMETER TUNING
text: In our experiments (see Section 4), for all of the inference techniques compared, we searched across a grid of hyperparameters that control the sparsity of inferred GC structures. Other hyperparameters were fine-tuned manually. Final results reported in the paper correspond to the best hyperparameter configurations. With this testing setup, our goal was to fairly compare best achievable inferential performance of the techniques.
Tables 6, 7, and 8 provide ranges for hyperparameter values considered in each experiment. For cMLP and cLSTM (Tank et al., 2018), parameter λ is the weight of the group Lasso penalty; for TCDF (Nauta et al., 2019), significance parameter α is used to decide which potential GC relationships are significant; eSRU (Khanna & Tan, 2020) has three different penalties weighted by λ1:3. For the stability-based thresholding (see Algorithm 1) in GVAR, we used Q = 20 equally spaced values in [0, 1] as sequence ξ5. For Lorenz 96 and fMRI experiments, grid search results are plotted in Figures 4, 8, and 5. Figure 9 contains GVAR grid search results for the Lotka–Volterra experiment.
5We did not observe high sensitivity of performance w.r.t. ξ, as long as sufficiently many evenly spaced sparsity levels are considered.
H.1 LORENZ 96
H.2 FMRI

heading: I COMPARISON WITH DYNAMIC BAYESIAN NETWORKS
text: We provide a comparison between GVAR and linear Gaussian dynamic Bayesian networks (DBN). DBNs are a classical approach to temporal structure learning (Murphy & Russell, 2002). We use R (R Core Team, 2020) package dbnR (Quesada, 2020) to fit DBNs on all datasets considered in Section 4. We use two structure learning algorithms: the max-min hill-climbing (MMHC) (Tsamardinos et al., 2006) and the particle swarm optimisation (Xing-Chen et al., 2007). Table 9 contains average balanced accuracies achieved by DBNs and GVAR for inferring the GC structure. Not surprisingly, DBNs outperform GVAR on the time series with linear dynamics, but fail to infer the true structure on Lorenz 96, fMRI, and Lotka–Volterra datasets.

heading: J THE LORENZ 96 SYSTEM: FURTHER EXPERIMENTS
text: Figure 10: Inferential performance of GVAR across a range of forcing constant values.
In addition to the experiments in Section 4.1.1, we examine the performance of VAR and GVAR models across a range of forcing constant values F = 0, 5, 10, 25, 50 for the Lorenz 96 system with p = 20 variables. Figure 10 shows average AUPRCs with bands corresponding to the 95% CI for the mean. It appears that for both models, inference is more challenging for lower (< 10) and higher values of F (> 20). This observation is in agreement with the results in Section 4.1.1, where all inference techniques performed worse under F = 40 than under F = 10. Note that herein same GVAR hyperparameters were used across all values of F . It is possible that better inferential performance could be achieved with GVAR after comprehensive hyperparameter tuning.

heading: K THE LOTKA–VOLTERRA SYSTEM
text: The original Lotka–Volterra system (Bacaër, 2011) includes only one predator and one prey species, population sizes of which are denoted by x and y, respectively. Population dynamics are given by the following coupled differential equations:
dx dt = αx− βxy, (12) dy dt = δyx− ρy, (13)
where α, β, δ, ρ > 0 are fixed parameters determining strengths of interactions.
In this paper, we consider a multiple species version of the system, given by Equations 9 and 10 in Section 4.2. We simulate the system under α = ρ = 1.1, β = δ = 0.2, η = 2.75 × 10−5,∣∣Pa(xi)∣∣ = ∣∣Pa(yj)∣∣ = 2, p = 10, i.e. 2p = 20 variables in total, with T = 2000 observations. Figure 11 depicts signs of GC effects between variables in a multi-species Lotka–Volterra with 2p = 20 species and 2 parents per variable. We simulate this system numerically by using the Runge-Kutta method6. We make a few adjustments to the state transition equations, in particular: we introduce normally-distributed innovation terms to make simulated data noisy; during state transitions, we clip all population sizes below 0. Figure 12 shows traces of generalised coefficients inferred by GVAR: magnitudes and signs of coefficients reflect the true dependency structure.
6Simulations are based on the implementation available at https://github.com/smkalami/ lotka-volterra-in-python.

heading: L EFFECT SIGN DETECTION IN A LINEAR VAR
text: Herein we provide results for the evaluation of GVAR and our inference framework on a very simple synthetic time series dataset. We simulate time series with p = 4 variables and linear interaction dynamics given by the following equations:
xt = a1xt−1 + εxt , wt = a2wt−1 + a3xt−1 + εwt ,
yt = a4yt−1 + a5wt−1 + ε y t , zt = a6zt−1 + a7wt−1 + a8yt−1 + ε z t,
(14)
where coefficients ai ∼ U ([−0.8,−0.2] ∪ [0.2, 0.8]) are sampled independently in each simulation; and ε·t ∼ N (0, 0.16) are additive innovation terms. This is an adapted version of one of artificial datasets described by Peters et al. (2013), but without instantaneous effects.
The GC summary graph of the system is visualised in Figure 13. It is considerably denser than for the Lorenz 96, fMRI, and Lotka–Volterra time series investigated in Section 4.
Similarly to the experiment described in Section 4.2, we infer GC relationships with the proposed framework and evaluate inference results against the true dependency structure and effect signs. Table 10 contains average performance across 10 simulations achieved by GVAR with hyperparameter values K = 1, λ = 0.2, and γ = 0.5. In addition, we provide results for some of the baselines (no systematic hyperparameter tuning was performed for this experiment).
GVAR attains perfect AUROC and AUPRC in all 10 simulations. In some cases, stability-based thresholding fails to recover a completely correct GC structure, nevertheless, average accuracy and balanced accuracy scores are satisfactory. Signs of inferred generalised coefficients mostly agree with the ground truth effect signs, as given by coefficients a1:8 in Equation 14.
Not surprisingly, linear VAR performs the best on this dataset w.r.t. all evaluation metrics. Both cMLP and eSRU successfully infer GC
relationships, achieving results comparable to GVAR. However, neither infers effect signs as well as GVAR. Thus, similarly to the experiment in Section 4.2, we conclude that generalised coefficients are more interpretable than neural network weights leveraged by cMLP, TCDF, and eSRU.
To summarise, this simple experiment serves as a sanity check and shows that our GC inference framework performs reasonably in low-dimensional time series with linear dynamics and a relatively dense GC summary graph (cf. Figure 7). Generally, the method successfully infers both the dependency structure and interaction signs.

heading: M PREDICTION ERROR
text: Herein we evaluate the prediction error of models on held-out data. Last 20% of time series points were held out to perform prediction on Lorenz 96, fMRI, and Lotka–Volterra datasets. Root-meansquare error (RMSE) was computed for predictions across R = 5 independent replicates:
RMSE = 1
p p∑ j=1 √√√√∑Tt=1 (x̂jt − xjt) T , (15)
where x̂jt is the one-step forecast made by a model for the t-th point of the j-th variable, and T is the length of the held-out time series segment. Table 11 contains average RMSEs for all models across the considered datasets. In general, RMSEs are not associated with the inferential performance of the models (cf. tables 1, 2, and 3). For example, while TCDF achieves the best inferential performance on fMRI (see Table 2), its prediction error is higher than for cMLP. This ‘misalignment’ between the prediction error and the consistency of variable selection is not surprising and has been discussed before, e.g. by Meinshausen & Bühlmann (2010).

Review summaries:
This paper primarily deals with learning Granger-causal relationships in multivariate time series in the nonlinear dynamics setting. The core method uses vector autoregressive modeling with sparsity inducing regularizers (elastic net and smoothness based fused lasso) along with the recently proposed with self-explaining neural networks (for interpretability). The authors also augment the framework by learning Granger-causal structures that are stable on original and time-reversed data. Exhaustive empirical analysis is done with recent GC baselines.  Some of my concerns with the paper are the following