ECONOMY STATISTICAL RECURRENT UNITS FOR INFERRING NONLINEAR GRANGER CAUSALITY
Year: 2020

Abstract:
Granger causality is a widely-used criterion for analyzing interactions in largescale networks. As most physical interactions are inherently nonlinear, we consider the problem of inferring the existence of pairwise Granger causality between nonlinearly interacting stochastic processes from their time series measurements. Our proposed approach relies on modeling the embedded nonlinearities in the measurements using a component-wise time series prediction model based on Statistical Recurrent Units (SRUs). We make a case that the network topology of Granger causal relations is directly inferrable from a structured sparse estimate of the internal parameters of the SRU networks trained to predict the processes’ time series measurements. We propose a variant of SRU, called economy-SRU, which, by design has considerably fewer trainable parameters, and therefore less prone to overfitting. The economy-SRU computes a low-dimensional sketch of its high-dimensional hidden state in the form of random projections to generate the feedback for its recurrent processing. Additionally, the internal weight parameters of the economy-SRU are strategically regularized in a group-wise manner to facilitate the proposed network in extracting meaningful predictive features that are highly time-localized to mimic real-world causal events. Extensive experiments are carried out to demonstrate that the proposed economy-SRU based time series prediction model outperforms the MLP, LSTM and attention-gated CNN-based time series models considered previously for inferring Granger causality.

text: Granger causality is a widely-used criterion for analyzing interactions in largescale networks. As most physical interactions are inherently nonlinear, we consider the problem of inferring the existence of pairwise Granger causality between nonlinearly interacting stochastic processes from their time series measurements. Our proposed approach relies on modeling the embedded nonlinearities in the measurements using a component-wise time series prediction model based on Statistical Recurrent Units (SRUs). We make a case that the network topology of Granger causal relations is directly inferrable from a structured sparse estimate of the internal parameters of the SRU networks trained to predict the processes’ time series measurements. We propose a variant of SRU, called economy-SRU, which, by design has considerably fewer trainable parameters, and therefore less prone to overfitting. The economy-SRU computes a low-dimensional sketch of its high-dimensional hidden state in the form of random projections to generate the feedback for its recurrent processing. Additionally, the internal weight parameters of the economy-SRU are strategically regularized in a group-wise manner to facilitate the proposed network in extracting meaningful predictive features that are highly time-localized to mimic real-world causal events. Extensive experiments are carried out to demonstrate that the proposed economy-SRU based time series prediction model outperforms the MLP, LSTM and attention-gated CNN-based time series models considered previously for inferring Granger causality.

heading: 1 INTRODUCTION
text: The physical mechanisms behind the functioning of any large-scale system can be understood in terms of the networked interactions between the underlying system processes. Granger causality is one widely-accepted criterion used in building network models of interactions between large ensembles of stochastic processes. While Granger causality may not necessarily imply true causality, it has proven effective in qualifying pairwise interactions between stochastic processes in a variety of system identification problems, e.g., gene regulatory network mapping (Fujita et al. (2007)), and the mapping of human brain connectome (Seth et al. (2015)). This perspective has given rise to the canonical problem of inferring pairwise Granger causal relationships between a set of stochastic processes from their time series measurements. At present, the vast majority of Granger causal inference methods adopt a model-based inference approach whereby the measured time series data is modeled using with a suitable parameterized data generative model whose inferred parameters ultimately reveal the true topology of pairwise Granger causal relationships. Such methods typically rely on using linear regression models for inference. However, as illustrated in the classical bivariate example by Baek & Brock (1992), linear model-based Granger causality tests can fail catastrophically in the presence of even mild nonlinearities in the measurements, thus making a strong case for our work which tackles the nonlinearities in the measurements by exploring new generative models of the time series measurements based on recurrent neural networks.

heading: 2 PROBLEM FORMULATION
text: Consider a multivariate dynamical system whose evolution from an initial state is fully characterized by n distinct stochastic processes which can potentially interact nonlinearly among themselves. Our goal here is to unravel the unknown nonlinear system dynamics by mapping the entire network of pairwise interactions between the system-defining stochastic processes, using Granger causality as the qualifier of the individual pairwise interactions.
In order to detect the pairwise Granger causal relations between the stochastic processes, we assume access to their concurrent, uniformly-sampled measurements presented as an n-variate time series x = {xt : t ∈ N} ⊂ Rn. Let xt,i denote the ith component of the n-dimensional vector measurement xt, representing the measured value of process i at time t. Motivated by the framework proposed in Tank et al. (2017), we assume that the measurement samples xt, t ∈ N are generated sequentially according to the following nonlinear, component-wise autoregressive model:
xt,i = fi (xt−p:t−1,1,xt−p:t−1,2, . . . ,xt−p:t−1,n) + et,i, i = 1, 2, . . . n, (1)
where xt−p:t−1,j , {xt−1,j , xt−2,j , . . . , xt−p,j} represents the most recent p measurements of the jth component of x in the immediate past relative to current time t. The scalar-valued component generative function fi captures all of the linear and nonlinear interactions between the n stochastic processes up to time t− 1 that decide the measured value of the ith stochastic process at time t. The residual ei,t encapsulates the combined effect of all instantaneous and exogenous factors influencing the measurement of process i at time t, as well as any imperfections in the presumed model. Equation 1 may be viewed as a generalization of the linear vector autoregressive (VAR) model in the sense that the components of x can be nonlinearly dependent on one another across time. The value p is loosely interpreted to be the order of the above nonlinear autoregressive model.

heading: 2.1 GRANGER CAUSALITY IN NONLINEAR DYNAMICAL SYSTEMS
text: We now proceed to interpret Granger causality in the context of the above component-wise time series model. Recalling the standard definition by Granger (1969), a time series v is said to Granger cause another time series u if the past of v contains new information above and beyond the past of u that can improve the predictions of current or future values of u. For x with its n components generated according to equation 1, the concept of Granger causality can be extended as suggested by Tank et al. (2018) as follows. We say that series j does not Granger cause series i if the componentwise generative function fi does not depend on the past measurements in series j, i.e., for all t ≥ 1 and all distinct pairs xt−p:t−1,j and x′t−p:t−1,j ,
fi (xt−p:t−1,1, . . . ,xt−p:t−1,j , . . . ,xt−p:t−1,n) = fi ( xt−p:t−1,1, . . . ,x ′ t−p:t−1,j , . . . ,xt−p:t−1,n ) . (2) From equation 1, it is immediately evident that under the constraint in equation 2, the past of series j does not assert any causal influence on series i, in alignment with the core principle behind Granger causality. Based on the above implication of equation 2, the detection of Granger noncausality between the components of x translates to identifying those components of x whose past is irrelevant to the functional description of each individual fi featured in equation 1.
Note that any reliable inference of pairwise Granger causality between the components of x is feasible only if there are no unobserved confounding factors in the system which could potentially influence x. In this work, we assume that the system of interest is causally sufficient (Spirtes & Zhang (2016)), i.e., none of the n stochastic processes (whose measurements are available) have a common Granger-causing-ancestor that is unobserved.

heading: 2.2 INFERRING GRANGER CAUSALITY USING COMPONENT-WISE RECURRENT MODELS
text: We undertake a model-based inference approach wherein the time series measurements are used as observations to learn an autoregressive model which is anatomically similar to the componentwise generative model described in equation 1 except for the unknown functions fi replaced with their respective parameterized approximations denoted by gi. Let Θi, 1 ≤ i ≤ n denote the complete set of parameters encoding the functional description of the approximating functions {gi}ni=1. Then, the pairwise Granger causality between series i and the components of x is deduced from Θi which is estimated by fitting gi’s output to the ordered measurements in series i. Specifically, if the estimated Θi suggests that gi’s output is independent of the past measurements in series j, then
we declare that series j is Granger noncausal for series i. We aim to design the approximation function gi to be highly expressive and capable of well-approximating any intricate causal coupling between the components of x induced by the component-wise function fi, while simultaneously being easily identifiable from underdetermined measurements.
By virtue of their universal approximation property (Schäfer & Zimmermann (2006)), recurrent neural networks or RNNs are a particularly ideal choice for gi towards inferring the pairwise Granger causal relationships in x. In this work, we investigate the use of a special type of RNN called the statistical recurrent unit (SRU) for inferring pairwise Granger causality between multiple nonlinearly interacting stochastic processes. Introduced by Oliva et al. (2017), an SRU is a highly expressive recurrent neural network designed specifically for modeling multivariate time series data with complex-nonlinear dependencies spanning multiple time lags. Unlike the popular gated RNNs (e.g., long short-term memory (LSTM) (Hochreiter & Schmidhuber (1997)) and gated recurrent unit (GRU)) (Chung et al. (2014)), the SRU’s design is completely devoid of the highly nonlinear sigmoid gating functions and thus less affected by the vanishing/exploding gradient issue during training. Despite its simpler ungated architecture, an SRU can model both short and long-term temporal dependencies in a multivariate time series. It does so by maintaining multi-time scale summary statistics of the time series data in the past, which are preferentially sensitive to different older portions of the time series x. By taking appropriate linear combinations of the summary statistics at different time scales, an SRU is able to construct predictive causal features which can be both highly component-specific and lag-specific at the same time. From the causal inference perspective, this dual-specificity of the SRU’s predictive features is its most desirable feature, as one would argue that causal effects in reality also tend to be highly localized in both space and time.
The main contributions of this paper can be summarized as follows:
1. We propose the use of statistical recurrent units (SRUs) for detecting pairwise Granger causality between the nonlinearly interacting stochastic processes. We show that the entire network of pairwise Granger causal relationships can be inferred directly from the regularized block-sparse estimate of the input-layer weight parameters of the SRUs trained to predict the time series measurements of the individual processes.
2. We propose a modified SRU architecture called economy SRU or eSRU in short. The first of the two proposed modifications is aimed at substantially reducing the number of trainable parameters in the standard SRU model without sacrificing its expressiveness. The second modification entails regularizing the SRU’s internal weight parameters to enhance the interpretability of its learned predictive features. Compared to the standard SRU, the proposed eSRU model is considerably less likely to overfit the time series measurements.
3. We conduct extensive numerical experiments to demonstrate that eSRU is a compelling model for inferring pairwise Granger causality. The proposed model is found to outperform the multi-layer perceptron (MLP), LSTM and attention-gated convolutional neural network (AG-CNN) based models considered in the earlier works.

heading: 3 PROPOSED GRANGER CAUSAL INFERENCE FRAMEWORK
text: In the proposed scheme, each of the unknown generative functions fi, 1 ≤ i ≤ n in the presumed component-wise model of x in (1) is individually approximated by a distinct SRU network. The ith SRU network sequentially processes the time series measurements x and outputs a next-step prediction sequence x̂+i = {x̂i,2, x̂i,3, . . . , x̂i,t+1, . . .} ⊂ R, where x̂i,t+1 denotes the predicted value of component series i at time t + 1. The prediction x̂i,t+1 is computed in a recurrent fashion by combining the current input sample xt at time t with the summary statistics of past samples of x up to and including time t− 1 as illustrated in Figure 1. The following update equations describe the sequential processing of the input time series x within the ith SRU network in order to generate a prediction of xi,t+1.
Feedback: ri,t = h ( W (i)r ui,t−1 + b (i) r ) ∈ Rdr . (3a)
Recurrent statistics: φi,t = h ( W (i) in xt + W (i) f ri,t−1 + b (i) in ) ∈ Rdφ . (3b)
Multi-scale summary statistics: ui,t = [( uα1i,t )T ( uα2i,t )T . . . ( uαmi,t )T ]T ∈ Rmdφ , αj ∈ A,∀j. (3c)
Single-scale summary statistics: uαji,t = (1− αj)u αj i,t−1 + αjφi,t, ∈ R dφ , αj ∈ [0, 1]. (3d) Output features: oi,t = h ( W (i)o ui,t + b (i) o ) ∈ Rdo . (3e)
Output prediction: x̂i,t+1 = ( w(i)y )T oi,t + b (i) y ∈ R. (3f)
The function h in the above updates is the elementwise Rectified Linear Unit (ReLU) operator, h(·) := max(·, 0), which serves as the nonlinear activation in the three dedicated single layer neural networks that generate the recurrent statistics φi,t, the feedback ri,t and the output features oi,t in the ith SRU network. In order to generate the next-step prediction of series i at time t, the ith SRU network first prepares the feedback ri,t by nonlinearly transforming its last hidden state ui,t−1. As stated in equation 3a, a single layer ReLU network parameterized by weight matrix W (i)r and bias b (i) r maps the hidden state ui,t−1 to the feedback ri,t. Another single layer ReLU network parameterized by weight matrices W (i)in ,W (i) f and bias b (i) in takes the input xt and the feedback ri,t and tranforms them into the recurrent statistics φi,t as described in equation 3b. Equation 3d describes how the network’s multi-timescale hidden states uαi,t for α ∈ A = {α1, α2, . . . , αm} ⊂ [0, 1] are updated in parallel by taking exponentially weighted moving averages of the recurrent statistics φi,t corresponding to m different scales in A. A third single layer ReLU network parameterized by W (i)o and b (i) o transforms the concatenated multi-
timescale summary statistics ui,t = [ (uα1i,t ) T (uα2i,t ) T . . . (uαmi,m) T ]T
to generate the nonlinear causal features oi,t which, according to Oliva et al. (2017), are arguably highly sensitive to the input time series measurements at specific lags. Finally, the network generates the next-step prediction of series i as x̂i,t+1 by linearly combining the nonlinear output features in oi,t, as depicted in equation 3f.
For values of scale α ≈ 1, the single-scale summary statistic uαi,t in equation 3d is more sensitive to the recent past measurements in x. On the other hand, α ≈ 0 yields a summary statistic that is more representative of the older portions of the input time series. Oliva et al. (2017) elaborates on how the SRU is able to generate output features (oi,t, 1 ≤ i ≤ n) that are preferentially sensitive to the measurements from specific past segments of x by taking appropriate linear combinations of the summary statistics corresponding to different values of α in A.

heading: 3.1 INFERRING PAIRWISE GRANGER CAUSALITY USING SRUS
text: Let Θ(i)SRU , {W (i) f ,W (i) in , b (i) in ,W (i) r , b (i) r ,W (i) o , b (i) o ,w (i) y , b (i) y } denote the complete set of parameters of the ith SRU network approximating fi in the presumed component-wise model of x. From equation 3b, we observe that the weight matrix W (i)in regulates the influence of the individual components of the input time series x on the generation of the recurrent statistics φi,t, and ultimately the next-step prediction of series i. In real-world dynamical systems, the networked interactions are typically sparse which implies that very few dimensions of the input time series x actually play a role in the generation of its individual components. Bearing this property of the networked interactions in mind, we are interested in learning the parameters Θ(i)SRU such that the i th SRU’s sequential
output closely matches with series i’s measurements, while simultaneously seeking a column-sparse estimate of the weight matrix W (i)in .
We propose to learn the parameters Θ(i)SRU of the i th SRU network by minimizing the penalized mean squared prediction error loss as shown below.
Θ̂ (i) SRU := arg min
Θ (i) SRU
1
T − 1 T−1∑ t=1 (x̂i,t − xi,t+1)2 + λ1 n∑ j=1 ‖W (i)in (:, j)‖2. (4)
In the above, the network output x̂i,t depends nonlinearly on W (i) in according to the composite relation described by the updates (3a)-(3f) and W (i)in (:, j) denotes the j th column in the weight matrix W (i)in . The `1-group norm penalty in the objective is known to promote column sparsity in the estimated W (i)in (Simon et al. (2013)). From equation 3b, a straightforward implication of the column vector W (i)in (:, j) being estimated as the all-zeros vector is that the past measurements in series j do not influence the predicted future value of series i. In this case, we declare that series j does not Granger-cause series i. Moreover, the index set supporting the non-zero columns in the estimated weight matrix Ŵ (i)in enumerates the components of x which are likely to Granger-cause series i. Likewise, the entire network of pairwise Granger causal relationships in x can be deduced from the non-zero column support of the estimated weight matrices W (i)in , 1 ≤ i ≤ n in the n SRU networks trained to predict the components of x.
The component-wise SRU optimization problem in equation 4 is nonconvex and potentially has multiple local minima. To solve for Θ̂(i)SRU, we use first-order gradient-based methods such as stochastic gradient descent which have been found to be consistently successful in finding good solutions of nonconvex deep neural network optimization problems (Allen-Zhu et al. (2019)). Since our approach of detecting Granger noncausality hinges upon correctly identifying the all-zero columns of W (i)in , it is important that the first-order gradient based parameter updates used for minimizing the penalized SRU loss ensure that majority of the coefficients in W (i)in iterates become exactly zero after a certain number of iterations. Seeking exact column sparsity in the converged solution of W (i)in , we follow the same approach as Tank et al. (2018) and resort to a first-order proximal gradient descent algorithm to find a regularized solution of the SRU optimization. The gradients needed for executing the gradient descent updates of the SRU network parameters are computed efficiently using the backpropagation through time (BPTT) procedure (Jaeger (2002)).

heading: 4 ECONOMY SRU: A REMEDY FOR OVERFITTING
text: By computing the summary statistics of past measurements at sufficiently granular time scales, an SRU can learn predictive causal features which are highly localized in time. While a higher granularity of α in A translates to a more general SRU model that fits better to the time series measurements, it also entails substantial increase in the number of trainable parameters. Since measurement scarcity is typical in causal inference problems, the proposed component-wise SRU based time series prediction model is usually overparameterized and thus susceptible to overfitting. The typical high dimensionality of the recurrent statistic φt accentuates this issue.
To alleviate the overfitting concerns, we propose two modifications to the standard SRU (Oliva et al. (2017)) aimed primarily at reducing its likelihood of overfitting the time series measurements. The modifications are relevant regardless of the current Granger causal inference context, and henceforth we refer to the modified SRU as Economy-SRU (eSRU).

heading: 4.1 MODIFICATION-I: GENERATING FEEDBACK FROM A LOW-DIMENSIONAL SKETCH OF SUMMARY STATISTICS
text: We propose to reduce the number of trainable parameters in the ith SRU network by substituting the feedback ReLU network parameterized by W (i)r and b (i) r with the two stage network shown in Fig. 2. The first stage implements the linear matrix-vector multiplication operation D (i) r ui,t to generate the output vi,t ∈ Rd ′ r , where D(i)r ∈ Rd ′ r×mdφ is a fixed, full row-rank
matrix with d′r mdφ. The d′r-dimensional output of the first stage can be viewed as a lowdimensional, stable embedding of the multi-timescale summary statistics ui,t. The entries of the constant matrix D(i)r are drawn independently from a zero mean Gaussian distribution with variance 1d′r . The stage-1 processing is based on the premise that for most real-world systems and
Linear ReLU vi,t−1ui,t−1 ri,t
D (i) r W ′(i) r , b ′(i) r
Stage-1 (Encoder)
Stage-2 (Decoder)
tics vi,t to the feedback vector ri,t. The second stage ReLU network is parameterized by weight matrix W ′,(i)r ∈ Rdr×d ′ r and bias b′,(i)r ∈ Rdr . Compared to the standard SRU’s feedback whose generation is controlled by mdφdr + dr trainable parameters, the proposed feedback network has only d′rdr + dr trainable parameters, which is substantially fewer when d ′ r mdφ. Consequently, the modified SRU is less susceptible to overfitting.

heading: 4.2 MODIFICATION-II: GROUPED-SPARSE MIXING OF MULTI-SCALE SUMMARY STATISTICS FOR LEARNING TIME-LOCALIZED PREDICTIVE FEATURES
text: In the standard SRU proposed by Oliva et al. (2017), there are no restrictions on the weight matrix W (i)o parameterizing the ReLU network that maps the summary statistics ui,t to the final predictive features in oi,t. Noting that the number of parameters in themdφ×do sized weight matrix W (i)o usually dominates the overall number of trainable parameters in the SRU, any meaningful effort towards addressing the model overfitting concerns must consider regularizing the weights in W (i)o . In this spirit, we propose the following penalized optimization problem to estimate the parameters Θ
(i) eSRU = (Θ (i) SRU\{W (i) r }) ∪ {W ′(i)r } of the eSRU model equipped with the two-stage feedback
proposed in Section 4.1:
Θ̂ (i) eSRU = arg min
Θ (i) eSRU
1
T−1 T−1∑ t=1 (x̂i,t − xi,t+1)2 +λ1 n∑ j=1 ‖W (i)in (:, j)‖2 +λ2 do∑ j=1 dφ∑ k=1 ‖W (i)o (j,Gj,k)‖2.
(5) Here λ1 and λ2 are positive constants that bias the group sparse penalizations against the eSRU’s fit to the measurements in the ith component series. The term W io(j,Gj,k) (1 ≤ j ≤ do, 1 ≤ k ≤ dφ) denotes the subvector of the jth row in W (i)o obtained by extracting the weight coefficients indexed by set Gj,k. As shown via an example in Fig. 3, the index set Gj,k enumerates the m weight coefficients in the row vector W (i)o (j, :) which are multiplied to the exponentially weighted running averages of kth recurrent statistic φi,t(k) corresponding to the m timescales in A prior to being transformed by the neural unit generating the jth predictive feature in oi,t. Compared to equation 4, the second penalty term in equation 5 promotes a group-sparse solution for W (i)o to the effect that each predictive feature in oi,t depends on only a few components of the recurrent statistic φi,t via their linearly mixed multi-scale exponentially weighted averages. We opine that the learned linear mixtures, represented by the intermediate products W (i)o (j,Gj,k)ui,t(Gk,j), are highly sensitive to certain past segments of the input time series x. Consequently, the output features in oi,t are both time-localized and component-specific, a common trait of real-world causal effects.
1Gaussian random matrices of appropriate dimensions are approximately isometries with overwhelming probability (Johnson & Lindenstrauss (1984)). However, instead of using n independent instantiations of a Gaussian random matrix for initializing D(i)r , 1 ≤ i ≤ n, we recommend initializing them with the same random matrix, as the latter strategy reduces the probability that any one of them is spurious encoder by n-fold.
The above group-sparse regularization of the weight coefficients in W (i)o , combined with the column-sparsity of W (i)in , is pivotal to enforcing that the occurrence of any future pattern in a time series can be attributed to the past occurrences of a few highly time-localized patterns in the ancestral time series. The results of our numerical experiments further confirm that by choosing λ1 and λ2 appropriately, the proposed group-wise sparsity inducing regularization of W (i) o ameliorates overfitting, and the optimization in equation 5 yields an estimate of W (i)in whose column support closely reflects the true pairwise Granger causal relationships between the components of x.

heading: 5 EXPERIMENTS
text: We evaluate the performance of the proposed SRU- and eSRU-based component-wise time series models in inferring pairwise Granger causal relationships in a multivariate time series. The proposed models are compared to the existing MLP- and LSTM-based models in Tank et al. (2018) and the attention-gated CNN-based model (referred hereafter as Temporal Causal Discovery Framework (TCDF)) in Nauta et al. (2019). To ensure parity between the competing models, the maximum size of all the input/hidden/output layers in the different NN/RNN time series models is fixed to 10, unless specified otherwise. The complete list of tuned hyperparameters of the considered models used for different datasets is provided in Appendix G. The performance of each method is qualified in terms of its AUROC (Area Under the Receiver Operating Characteristic curve). Here, the ROC curve illustrates the trade off between the true-positive rate (TPR) and the false-positive rate (FPR) achieved by the methods towards the detection of n2 pairwise Granger causal relationships between the nmeasured processes in the experiment. The ROC curves of SRU and eSRU models are obtained by sweeping through different values of the regularization parameter λ1 in equation 4 and equation 5, respectively. Likewise, the ROCs of component-wise MLP and LSTM models are obtained by varying λ1’s counterpart in Tank et al. (2018). For TCDF, the ROC curve is obtained by varying the threshold that is applied to attention scores of the trained AG-CNN model in Nauta et al. (2019).

heading: 5.1 LORENZ-96 SIMULATIONS
text: In the first set of experiments, the time series measurements x intended for Granger causal inference are generated according to the Lorenz-96 model which has been extensively used in climate science for modeling and prediction purposes (Schneider et al. (2017)). In the Lorenz-96 model of an nvariable system, the individual state trajectories of the n variables are governed by the following set of odinary differential equations:
∂xt,i ∂t = −xt,i−1 (xt,i−2 − xt,i+1)− xt,i + F, 1 ≤ i ≤ n. (6)
where the first and the second terms on the RHS represent the advection and the diffusion in the system, respectively, and the third term F is the magnitude of the external forcing. The system dynamics becomes increasingly chaotic for higher values of F (Karimi & Paul (2010)). We evaluate and compare the accuracy of the proposed methods in inferring pairwise Granger causal relationships between n = 10 variables with Lorenz-96 dynamics. We consider two settings: F = 10 and F = 40 in order to simulate two different strengths of nonlinearity in the causal interactions between
the variables. Here, the ground truth is straightforward i.e., for any 1 ≤ i ≤ n, the ith component of time series x is Granger caused by its components with time indices from i− 2 to i+ 1. In the case of weak nonlinear interactions (F = 10), from Table 1a, we observe that eSRU achieves the highest AUROC among all competing models. The gap in performance is more pronounced when fewer time series measurements (T = 250) are available. In case of stronger nonlinear interactions (F = 40), we observe that both SRU and eSRU are the only models that are able to perfectly recover the true Granger causal network (Table 1b). Surprisingly, the SRU and eSRU models perform poorer when F is small. This could be attributed to the proposed models not sufficiently regularized when fitted to weakly-interacting time series measurements that are less nonlinear.

heading: 5.2 VAR SIMULATIONS
text: In the second set of simulations, we consider the time series measurements x to be generated according to a 3rd order linear VAR model:
xt = A (1)xt−1 + A (2)xt−2 + A (3)xt−3 + wt, t ≥ 1, (7)
where the matrices A(i), i = 1, 2, 3 contain the regression coefficients which model the linear interactions between its n = 10 components. The noise term wt is Gaussian distributed with zero mean and covariance 0.01I . We consider a sparse network of Granger causal interactions with only 30% of the regression coefficients in Ai selected uniformly being non-zero and the regression matrices Ai being collectively joint sparse (same setup as in Bolstad et al. (2011)). All non-zero regression coefficients are set equal to 0.0994 which guarantees the stability of the simulated VAR process.
From Table 2, we observe that all time series models generally achieve a higher AUROC as the number of measurements available increases. For T = 500, the component-wise MLP and the proposed eSRU are statistically tied when comparing their average AUROCs. For T = 1000, eSRU significantly outperforms the rest of the time series models and is able to recover the true Granger causal network almost perfectly.
Table 2: Averaged AUROC for 5 independently generated VAR datasets
MODEL AVERAGE AUROC
T = 500 T = 1000
MLP 0.94± 0.03 0.93± 0.02 LSTM 0.79± 0.12 0.8± 0.09 TCDF 0.77± 0.07 0.78± 0.04 SRU 0.82± 0.06 0.91± 0.04 eSRU 0.93± 0.05 0.98± 0.01
Table 3: Average AUROC corresponding to inferred brain connectivity for 5 human subjects
MODEL AVERAGE AUROC
T = 200
MLP 0.81± 0.04 LSTM 0.70± 0.03 TCDF 0.75± 0.04 SRU 0.78± 0.02 eSRU 0.84± 0.03

heading: 5.3 IN SILICO ESTIMATION OF BRAIN CONNECTIVITY USING BOLD SIGNALS
text: In the third set of experiments, we apply the different learning methods to estimate the connections in the human brain from simulated blood oxygenation level dependent (BOLD) imaging data. Here, the individual components of x comprise T = 200 time-ordered samples of
the BOLD signals simulated for n = 15 different brain regions of interest (ROIs) in a human subject. To conduct the experiments, we use simulated BOLD time series measurements corresponding to the five different human subjects (labelled as 2 to 6) in the Sim-3.mat file shared at https://www.fmrib.ox.ac.uk/datasets/netsim/index.html. The generation of the Sim3 dataset is described in Smith et al. (2011). The goal here is to detect the directed connectivity between different brain ROIs in the form of pairwise Granger causal relationships between the components of x.
From Table 3, it is evident that eSRU is more robust to overfitting compared to the standard SRU and detects the true Granger causal relationships more reliably. Interestingly, a single-layer cMLP model is found to outperform more complex cLSTM and attention gated-CNN (TCDF) models; however we expect the latter models to perform better when more time series measurements are available.

heading: 5.4 DREAM-3 IN SILICO NETWORK INFERENCE CHALLENGE
text: In the final set of experiments, we evaluate the performance of the different time series models in inferring gene regulation networks synthesized for the DREAM-3 In Silico Network Challenge (Prill et al. (2010); Marbach et al. (2009)). Here, the time series x represents the in silico measurements of the gene expression levels of n = 100 genes, available for estimating the gene regulatory networks of E.coli and yeast. A total of five gene regulation networks are to be inferred (two for E.coli and three for yeast) from the networks’ gene expression level trajectories recorded while they recover from 46 different perturbations (each trajectory has 21 time points). All NN/RNN models are implemented with 10 neurons per layer, except for the componentwise MLP model which has 5 neurons per layer. From Table 4, we can observe that the proposed SRU and eSRU models are generally more accurate
compared to the MLP, LSTM, and attention-gated CNN (TCDF) models in inferring the true genegene interactions. For four out of the five gene regulatory networks, either SRU or eSRU was the best performing model among the competing ones.

heading: 6 CONCLUSION
text: In this work, we addressed the problem of inferring pairwise Granger causal relationships between stochastic processes that interact nonlinearly. We showed that the such causality between the processes can be robustly inferred from the regularized internal parameters of the proposed eSRU-based recurrent models trained to predict the time series measurements of the individal processes. Future work includes:
i Investigating the use of other loss functions besides the mean-square error loss which can capture the exogenous and instantaneous causal effects in a more realistic way.
ii Incorporating unobserved confounding variables/processes in recurrent models.
iii Inferring Granger causality from multi-rate time series measurements.

heading: ACKNOWLEDGEMENTS
text: This work is supported by a Singapore Ministry of Education (MOE) Tier 2 Grant (R-263-000-C83112).

heading: A REVIEW OF RELATED WORK
text: Initial efforts in testing for nonlinear Granger causality focused mostly on the nonparameteric approach. Baek & Brock (1992) proposed a general statistical test to detect nonlinear Granger causality between two variables under the assumption that the linear VAR modeling of their time series measurements results in i.i.d residual errors. The test involves computing correlation integral estimators of the conditional probabilities concerning distances between carefully selected lead-lag subsequences of the input bivariate time series. Successive works by Hiemstra & Jones (1994); Diks & Panchenko (2006); Bai et al. (2010); Diks & Wolski (2016) proposed their improved variants of
the Baek-Brock test. Hiemstra & Jones (1994) modified and extended the original Baek-Brock test to allow for weakly dependent residual errors in linear VAR modeling of the time series measurements. Diks & Panchenko (2006) proposed a rectified version of the Heimstra-Jones’ test statistic, making it unbiased while fixing the issue of overrejection of the null hypothesis. Later works by Bai et al. (2010); Diks & Wolski (2016) extended the bi-variate test in Diks & Panchenko (2006) to the multivariate setting. The biggest common drawback of these nonparameteric tests is the large sample sizes required to robustly estimate the conditional probabilities that constitute the test statistic. Furthermore, the prevalent strategy in these methods of testing each one of the variable-pairs individually to detect pairwise Granger causality is unappealing from a computational standpoint, especially when a very large number of variables are involved.
In the model driven approach, the Granger causal relationships are inferred directly from the parameters of a data generative model fitted to the time series measurements. Compared to the nonparameteric approach, the model-based inference approach is considerably more sample efficient, however the scope of inferrable causal dependencies is dictated by the choice of data generative model. Nonlinear kernel based regression models have been found to be reasonably effective in testing of nonlinear Granger causality. Kernel methods rely on linearization of the causal interactions in a kernel-induced high dimensional feature space; the linearized interactions are subsequently modeled using a linear VAR model in the feature space. Based on this idea, Marinazzo et al. (2008) proposes a kernel Granger causality index to detect pairwise nonlinear Granger causality in the multivariate case. In Sindhwani et al. (2013); Lim et al. (2014), the nonlinear dependencies in the time series measurements are modeled using nonlinear functions expressible as sums of vector valued functions in the induced reproducing kernel Hilbert space (RKHS) of a matrix-valued kernel. In Lim et al. (2014), additional smoothness and structured sparsity constraints are imposed on the kernel parameters to promote consistency of the time series fitted nonlinear model. Shen et al. (2016) proposes a nonlinear kernel-based structural VAR model to capture instantaneous nonlinear interactions. The existing kernel based regression models are restrictive as they consider only additive linear combinations of the RKHS functions to approximate the nonlinear dependencies in the time series. Furthermore, deciding the optimal order of kernel based regression models is difficult as it requires prior knowledge of the mimimum time delay beyond which the causal influences are negligible.
By virtue of their universal approximation ability, RNNs offer a pragmatic way forward in modeling of complex nonlinear dependencies in the time series measurements for the purpose of inferring Granger causality. Mutiple recent works by Wang et al. (2018); Duggento et al. (2019); Abbasvandi & Nasrabadi (2019) have investigated the use of different types of RNNs for inferring nonlinear Granger causal relationships. However, they all adopt the same naı̈ve strategy whereby each pairwise causal relationship is tested individually by estimating its causal connection strength. The strength of the causal connection from series j to series i is determined by the ratio of mean-squared prediction errors incurred by unrestricted and restricted RNN models towards predicting series i using the past measurement sequences of all n component including and excluding the jth component alone, respectively. The pairwise testing strategy however does not scale well computationally as the number of component series becomes very large. This strategy also fails to exploit the typical sparse connectivity of networked interactions between the processes which has unlocked significant performance gains in the existing linear methods (Bahadori & Liu (2013); Bolstad et al. (2011)).
In a recent work by Tank et al. (2018), the pairwise Granger causal relationships are inferred directly from the weight parameters of component-wise MLP or LSTM networks fitted to the time series measurements. By enforcing column-sparsity of the input-layer weight matrices in the fitted MLP/LSTM models, their proposed approach returns a sparsely connected estimate of the underlying Granger causal network. Due to its feedforward architecture, a traditional MLP network is not well-suited for modeling ordered data such as a time series. Tank et al. (2018) demonstrated that the MLP network can learn short range temporal dependencies spanning a few time delays by letting the network’s input stage process multi-lag time series data over sliding windows. However, modeling long-range temporal dependencies using the same approach requires a larger sliding window size which entails an inconvenient increase in the number of trainable parameters. The simulation results in Tank et al. (2018) indicate that MLP models are generally outperformed by LSTM models in extracting the true topology of pairwise Granger causality, especially when the processes interact in a highly nonlinear and intricate manner. While purposefully designed for modeling short and long term temporal dependencies in a time series, the LSTM (Hochreiter & Schmidhuber (1997)) is very
general and often too much overparameterized and thus prone to overfitting. While using overparameterized models for inference is preferable when there is abundant training data available to leverage upon, there are several applications where the data available for causal inference is extremely scarce. It is our opinion that using a simpler RNN model combined with meaningful regularization of the model parameters is the best way forward in inferring Granger causal relationships from underdetermined time series measurements. Building on the ideas put forth by Tank et al. (2018), this paper investigates the use of Statistical Recurrent Units (SRUs) towards inferring Granger causality.

heading: B PROXIMAL GRADIENT DESCENT UPDATES FOR ESTIMATING THE
text: REGULARIZED WEIGHT PARAMETERS IN THE SRU AND eSRU MODELS
Noting that the proximal operator corresponding to mixed `1-`2 norm (group-norm) is the groupwise soft-thresholding operator, we use the following proximal gradient-descent updates to minimize the ith SRU’s regularized loss in equation 4:
W (i),t+1 in (:, j) = Sλ1η
( W
(i),t in (:, j)− η∇W (i)in (:,j)li(Θ (i),t SRU )
) , ∀j ∈ [n]. (8)
Here, li(Θti) , 1 T−1 ∑T−1 t=1 (x̂i,t − xi,t+1)
2 is the unregularized SRU loss function, η is the gradientdescent stepsize and Sλ1η is the elementwise soft-thresholding operator defined below.
Sλ1η(w) , { w − λ1η w‖w‖2 , ‖w‖2 > λ1η 0, ‖w‖2 ≤ λ1η , ∀w ∈ Rn. (9)
The columns of weight matrix W (i)i in the i th eSRU model are also updated in exactly the same fashion as above.
Likewise, the jth row of the group-norm regularized weight matrix W (i)o in the eSRU optimization in equation 5 is updated as shown below.
W (i),t+1o (j,Gj,k) = Sλ2η ( W (i),to (j,Gj,k)− η∇W (i)o (j,Gj,k)li(Θ (i),t eSRU) ) , ∀j = 1, 2, . . . do. (10)
The gradient of the unregularized loss function li, 1 ≤ i ≤ n associated with the SRU and eSRU models used in the above updates is evaluated via the backpropagation through time (BPTT) procedure (Jaeger (2002)).

heading: C ABLATION STUDY
text: C.1 CHOICE OF ENCODER FOR eSRU FEEDBACK’S STAGE-1: FIXED OR DATA DEPENDENT
As a possible further enhancement of the proposed eSRU time series model, one may consider learning the encoding map, D(i)r , in the feedback path, as trainable parameters of the ith eSRU. In Table 5, we compare the Granger causality detection performance of this particular eSRU variant and the proposed design wherein D(i)r is taken to be a random matrix with i.i.d. Gaussian entries. The experimental setup is kept the same as in Section 5, and the entries of D(i)r in the eSRU variant are `2-norm penalized during training.
We observe that the performance of these two models is statistically tied, which indicates that the randomly constructed D(i)r is able to distill the necessary information from the high-dimensional summary statistics ui,t−1 required for generating the feedback ri,t. Based on these results, we recommend using the proposed eSRU design with its randomly constructed encoding map D(i)r , because of its simpler design and reduced training complexity.
C.2 IMPACT OF GROUP-SPARSE REGULARIZATION OF W (i)o
In order to highlight the importance of learning time-localized predictive features in detecting Granger causality, we compare the following two time series models:
Once again, we use the same experimental settings as mentioned in Section 5. From Table 6, we observe that barring the Lorenz-96(T=250/500,F=40) datasets, for which nearly perfect recovery of the Granger causal network is achieved, the average AUROC improves consistently for the other datasets by switching from unstructured ridge regularization to the proposed group-sparse regularization of the output weight matrix W (i)o .
D IMPLEMENTATION DETAILS
• Activation function for SRU and eSRU models While the standard SRU proposed by Oliva et al. (2017) uses ReLU neurons, we found in our numerical experiments that using the Exponential Linear Unit (ELU) activation resulted in better performance. The ELU activation function is defined as
ELU(x) =
{ x x > 0
α(ex − 1) x ≤ 0 , α > 0. (11)
In our simulations, the constant α is set equal to one.
• Number of neural layers in SRU model To approximate the generative functions fi in equation 1, we consider the simplest architecture for the SRU networks, whereby the constituent ReLU networks generating the recurrent features, output features and feedback have a single layer feedforward design with equal number of neurons.
• Number of neural layers in Economy-SRU model The ReLU networks used for generating the recurrent and output features in the proposed
eSRU model have a single-layer feedforward design. However, the second stage of eSRU’s modified feedback can be either single or multi-layered feedforward network. Provided that d′r mdφ, a multi-layer implementation of the second stage of eSRU’s feedback can still have fewer trainable parameters overall compared to the SRU’s single layer feedback network. The simulation results in Section 5 are obtained using a two-layer ReLU network in the second stage of eSRU’s feedback for the DREAM-3 experiments, and while using a three-layer design for the Lorenz-96, VAR and NetSim experiments.
• Self-interactions in Dream-3 gene networks The in-silico gene networks synthesized for the DREAM-3 challenge have no selfconnections. Noting that none of the Granger causal inference methods evaluated in our experiments intentionally suppress the self-interactions, the reported AUROC values are computed by ignoring any self-connections in the inferred Granger causal networks.

heading: E PYTHON CODES
text: cMLP & cLSTM models Pytorch implementation of the componentwise MLP and LSTM models are taken from https: //github.com/icc2115/Neural-GC.
Temporal Causal Discovery Framework (TCDF) Pytorch implementation of the attention-gated CNN based Temporal Causal Discovery Framework (TCDF) is taken from https://github.com/M-Nauta/TCDF.
Proposed SRU and Economy-SRU models Pytorch implementations of the proposed componentwise SRU and eSRU models are shared at https://github.com/sakhanna/SRU_for_GCI.

heading: F ROC PLOTS
text: The receiver operating characteristics (ROC) of different Granger causal inference methods are compared in Figures 4-7. Here, an ROC curve represents the trade-off between the true-positive rate (TPR) and the false-positive rate (FPR) achieved by a given method while inferring the underlying pairwise Granger causal relationships.

heading: G TUNED HYPERPARAMETERS
text: Tables 7 to 11 summarize the chosen hyperparameters and configurations of the different NN/RNN models used for generating the results reported in Section 5. For the Dream-3 experiments, the model configurations used for inferring the E.coli-1 gene regulatory network Prill et al. (2010) have been provided in the tables.

Review summaries:
In this paper the authors propose using Statistical Recurrent Units to predict the network for Granger causality. They motivate this choice by the high representation power of SRUs for multivariate time series, the good performance they usually enjoy and as a way to alleviate the vanishing gradient problem. More importantly the particular form of the SRUs gives a very simple predictor and therefore explanability for Granger causality: the authors propose to simply mark serie $i$ as Granger caused by $j$ if the $j$th column of the input mixing matrix of the $g_i$ is non-zero.

the paper attempts to infer Granger causality between nonlinearly interacting stochastic processes from their time series measurements. instead of using MLP/LSTM etc to to model time series measurement, the paper proposed to use component-wise time series prediction model with Statistical Recurrent Units to model the measurements. they consider a low-dimensional version of SRU, which they call economy-SRU. in particular, they use group-wise regularizing to accompany the particular structure of the model to aid interpretability. they compared the performance with existing models with MLP/LSTM and show some gains in a few examples (but not all.) the proposal is interesting, but the experiment section might need further strengthening. currently, the experimental results do not immediately pop out as showing eSRU particularly useful.