2.1 Formal setting

Consider the usual Hilbert space, $L_2(\Omega )$ , of finite variance random variables over a probability space $(\Omega, \mathcal{F}, \mathbb{P})$ having inner product $\langle{x},{y}\rangle = \mathbb{E}[xy]$ . We will work with a discrete time and wide-sense stationary (WSS) $N$ -dimensional vector valued process $x(n)$ (with $n \in \mathbb{Z}$ ) where the $N$ elements take values in $L_2$ . We suppose that $x(n)$ has zero mean, $\mathbb{E} x(n) = 0$ , and has absolutely summable matrix valued covariance sequence $R(m) \overset{\Delta }{=} \mathbb{E} x(n)x(n-m)^{\mathsf{T}}$ , and hence an absolutely continuous spectral density.

We will also work frequently with the spaces spanned by the values of such a process.

(1) \begin{equation} \begin{aligned} \mathcal{H}_n^{(x)} &=\mathsf{cl\ } \Bigl \{\sum _{m = 0}^p A(m)^{\mathsf{T}} x(n-m), \ A(m) \in \mathbb{R}^N, p \in \mathbb{N} \Bigr \} \subseteq L_2(\Omega )\\ H_n^{(x)} &= \{a x(n)\ |\ a \in \mathbb{R}\} \subseteq L_2(\Omega ) \end{aligned} \end{equation}

where the closure is in mean-square. We will often omit the superscript $(x)$ , which should be clear from context. These spaces are separable, and as closed subspaces of a Hilbert space they are themselves Hilbert. We will denote the spaces generated in analogous ways by particular components of $x$ as, for example, $\mathcal{H}_n^{(i, j)}$ , $\mathcal{H}_n^{(i)}$ or by all but a particular component as $\mathcal{H}_n^{(-j)}$ .

From the Wold decomposition theorem (Lindquist & Picci, Reference Lindquist and Picci2015), every WSS sequence has the moving average $\mathsf{MA}(\infty )$ representation

(2) \begin{equation} x(n) = c(n) + \sum _{m = 0}^\infty A(m) v(n-m) \end{equation}

where $c(n)$ is a purely deterministic sequence, $v(n)$ is an uncorrelated sequence and $A(0) = I$ . We will assume that $c(n) = 0\ \forall n$ , which makes $x(n)$ a regular non-deterministic sequence. This representation can be inverted (Lindquist & Picci, Reference Lindquist and Picci2015) to yield the Vector Autoregressive $\mathsf{VAR}(\infty )$ form

(3) \begin{equation} x(n) = \sum _{m = 1}^\infty B(m) x(n-m) + v(n) \end{equation}

The Equations (2) and (3) can be represented as $x(n) = \mathsf{A}(z)v(n) = \mathsf{B}(z)x(n) + v(n)$ via the action of the operators

\begin{equation*}\mathsf {A}(z) \;:\!=\; \sum _{m = 0}^\infty A(m)z^{-m}\end{equation*}

and

\begin{equation*}\mathsf {B}(z) \;:\!=\; \sum _{m = 1}^\infty B(m)z^{-m}\end{equation*}

where the operator $z^{-1}$ is the backshift operator acting on $\ell _2^N(\Omega )$ , the space of square summable sequences of $N$ -vectors having elements in $L_2(\Omega )$ that is:

(4) \begin{equation} \mathsf{B}_{ij}(z)x_j(n) \;:\!=\; \sum _{m = 1}^\infty B_{ij}(m)x_j(n - m) \end{equation}

The operator $\mathsf{A}(z)$ contains only non-positive powers of $z$ and is therefore called causal, and if $A(0) = 0$ then it is strictly causal. An operator containing only non-negative powers of $z$ , for example, $\mathsf{C}(z) = c_0 + c_1z + c_2z^2 + \cdots$ is anti-causal. These operators are often referred to as filters in signal processing, terminology which we often adopt as well.

Finally, we have the well known inversion formula

(5) \begin{equation} \mathsf{A}(z) = (I - \mathsf{B}(z))^{-1} = \sum _{k = 0}^\infty \mathsf{B}(z)^k \end{equation}

The above assumptions are quite weak, and are essentially the minimal requirements needed to have a non-degenerate model. The strongest assumption we require is finally that $\Sigma _v$ is a diagonal and positive-definite matrix. Equivalently, the Wold decomposition could be written with $A(0) \ne I$ and $\Sigma _v = I$ , and we would be assuming that $A(0)$ is diagonal. This precludes the possibility of instantaneous causality (Granger, Reference Granger1969) or of processes being only weakly (as opposed to strongly) feedback free, in the terminology of Caines (Reference Caines1976). We formally state our setup as a definition:

2.2 Granger causality

Definition 2.2 (Granger Causality). We will say that component $x_j$ Granger-causes component $x_i$ (with respect to $x$ ) and write $x_j \overset{\text{GC}}{\rightarrow } x_i$ if

(6) \begin{equation}{\xi [{x_i(n)}\ |\ {\mathcal{H}_{n - 1}}]} \lt{\xi [{x_i(n)}\ |\ {\mathcal{H}^{(-j)}_{n - 1}}]} \end{equation}

where $\xi [x \ |\ \mathcal{H}] \;:\!=\; \mathbb{E} (x -{\hat{\mathbb{E}}[x\ |\ \mathcal{H}]})^2$ is the linear mean-squared estimation error and

\begin{equation*}{\hat {\mathbb {E}}[x\ |\ \mathcal {H}]} \;:\!=\; \text {proj}_{\mathcal {H}}(x)\end{equation*}

denotes the (unique) projection onto the Hilbert space $\mathcal{H}$ .

This notion captures the idea that the process $x_j$ provides information about $x_i$ that is not available from elsewhere. The notion is closely related to the information theoretic measure of transfer entropy, indeed, if the distribution of $v(n)$ is known to be Gaussian then they are equivalent (Barnett et al., Reference Barnett, Barrett and Seth2009).

The notion of conditional orthogonality is the essence of Granger causality, and enables us to obtain results for a fairly general class of WSS processes, rather than simply $\mathsf{VAR}(p)$ models.

Definition 2.3 (Conditional Orthogonality). Consider three closed subspaces of a Hilbert space $\mathcal{A}$ , $\mathcal{B}$ , $\mathcal{X}$ . We say that $\mathcal{A}$ is conditionally orthogonal to $\mathcal{B}$ given $\mathcal{X}$ and write $\mathcal{A} \perp \mathcal{B}\ |\ \mathcal{X}$ if

\begin{equation*} \langle {a - {\hat {\mathbb {E}}[a\ |\ \mathcal {X}]}},{b - {\hat {\mathbb {E}}[b\ |\ \mathcal {X}]}}\rangle = 0\ \forall a \in \mathcal {A}, b \in \mathcal {B} \end{equation*}

An equivalent condition is that (see Lindquist & Picci (Reference Lindquist and Picci2015) Proposition 2.4.2)

\begin{equation*} {\hat {\mathbb {E}}[\beta \ |\ {\mathcal {A} \vee \mathcal {X}}]} = {\hat {\mathbb {E}}[\beta \ |\ \mathcal {X}]}\ \forall \beta \in \mathcal {B} \end{equation*}

The following theorem captures a number of well known equivalent definitions or characterizations of Granger causality:

2.3 Granger causality graphs

We establish some graph theoretic notation and terminology, collected formally in definitions for the reader’s convenient reference.

Our principal object of study will be a graph determined by Granger causality relations as follows.

Definition 2.6 (Causality graph). We define the Granger causality graph $\mathcal{G} = ([N], \mathcal{E})$ to be the directed graph formed on $N$ vertices where an edge $(j, i) \in \mathcal{E}$ if and only if $x_j$ Granger-causes $x_i$ (with respect to $x$ ). That is,

\begin{equation*}(j, i) \in \mathcal {E} \iff j \in {pa(i)} \iff x_j \overset {\text {GC}}{\rightarrow } x_i\end{equation*}

The edges of the Granger causality graph $\mathcal{G}$ can be given a general notion of “weight” by associating an edge $(j, i)$ with the strictly causal operator $\mathsf{B}_{ij}(z)$ (see Equation (4)). Hence, the matrix $\mathsf{B}(z)$ is analogous to a weighted adjacency matrix Footnote ¹ for the graph $\mathcal{G}$ . And, in the same way that the $k^{\text{th}}$ power of an adjacency matrix counts the number of paths of length $k$ between nodes, $(\mathsf{B}(z)^k)_{ij}$ is a filter isolating the “action” of $j$ on $i$ at a time lag of $k$ steps, this is exemplified in the inversion formula (5).

From the $\mathsf{VAR}$ representation of $x(n)$ there is clearly a close relationship between each node and its parent nodes, the relationship is quantified through the sparsity pattern of $\mathsf{B}(z)$ . The work of Caines & Chan (Reference Caines and Chan1975) starts by defining causality (including between groups of processes) via the notion of feedback free processes in terms of $\mathsf{A}(z)$ . The following proposition provides a similar interpretation of the sparsity pattern of $A(z)$ in terms of the causality graph $\mathcal{G}$ .

Proposition 2.1 (Ancestor Expansion). The component $x_i(n)$ of $x(n)$ can be represented in terms of its parents in $\mathcal{G}$ :

(7) \begin{equation} x_i(n) = v_i(n) + \mathsf{B}_{ii}(z)x_i(n) + \sum _{k \in{pa(i)}}\mathsf{B}_{ik}(z)x_k(n) \end{equation}

Moreover, $x_i(\cdot )$ can be expanded in terms of its ancestor’s $v(n)$ components only:

(8) \begin{equation} x_i(n) = \mathsf{A}_{ii}(z)v_i(n) + \sum _{k \in{\mathcal{A}({i})} \setminus \{i\}}\mathsf{A}_{ik}(z)v_k(n) \end{equation}

where $\mathsf{A}(z) = \sum _{m = 0}^\infty A(m)z^{-m}$ is the filter from the Wold decomposition representation of $x(n)$ , Equation (2).

Proof. Equation (7) is immediate from the $\mathsf{VAR}(\infty )$ representation of Equation (3) and Theorem 2.1, we are left to demonstrate (8).

From Equation (3), which we are assuming throughout the paper to be invertible, we can write

\begin{equation*} x(n) = (I - \mathsf {B}(z))^{-1} v(n) \end{equation*}

where necessarily $(I - \mathsf{B}(z))^{-1} = \mathsf{A}(z)$ due to the uniqueness of the Wold decomposition. Since $\mathsf{B}(z)$ is stable we have

(9) \begin{equation} (I - \mathsf{B}(z))^{-1} = \sum _{k = 0}^\infty \mathsf{B}(z)^k \end{equation}

Invoking the Cayley-Hamilton theorem allows writing the infinite sum of (9) in terms of finite powers of $\mathsf{B}$ .

Let $S$ be a matrix with elements in $\{0, 1\}$ which represents the sparsity pattern of $\mathsf{B}(z)$ , from Lemma D.1 $S$ is the transpose of the adjacency matrix for $\mathcal{G}$ and hence $(S^k)_{ij}$ is non-zero if and only if $j \in{gp_{k}({i})}$ , and therefore $\mathsf{B}(z)^k_{ij} = 0$ if $j \not \in{gp_{k}({i})}$ . Finally, since ${\mathcal{A}({i})} = \bigcup _{k = 1}^N{gp_{k}({i})}$ we see that $\mathsf{A}_{ij}(z)$ is zero if $j \not \in{\mathcal{A}({i})}$ .

Therefore

\begin{align*} x_i(n) &= [(I - \mathsf{B}(z))^{-1}v(n)]_i\\ &= \sum _{j = 1}^N \mathsf{A}_{ij}(z) v_j(n)\\ &= \mathsf{A}_{ii}(z) v_i(n) + \sum _{\substack{j \in{\mathcal{A}({i})} \\ j \ne i}} \mathsf{A}_{ij}(z) v_j(n) \end{align*}

This completes the proof.

Proposition 2.1 is ultimately about the sparsity pattern in the Wold decomposition matrices $A(m)$ since $x_i(n) = \sum _{m = 0}^\infty \sum _{j = 1}^N A_{ij}(m)v_j(n - m)$ . The proposition states that if $j \not \in{\mathcal{A}({i})}$ , then $\mathsf{A}_{ij}(z) = 0$ .

2.4 Pairwise granger causality

Recall that Granger causality in general must be understood with respect to a particular universe of observations. If $x_j \overset{\text{GC}}{\rightarrow } x_i$ with respect to $x_{-k}$ , it may not hold with respect to $x$ . For example, $x_k$ may be a common ancestor which when observed, completely explains the connection from $x_j$ to $x_i$ . In this section we study pairwise Granger causality, and seek to understand when knowledge of pairwise relations is sufficient to deduce the fully conditional relations of $\mathcal{G}$ .

This notion is of interest for a variety of reasons. From a purely conceptual standpoint, we will see how the notion can in some sense capture the idea of “flow of information” in the underlying graph, in the sense that if $j \in{\mathcal{A}({i})}$ we expect that $j \overset{\text{PW}}{\rightarrow } i$ . It may also be useful for reasoning about the conditions under which unobserved components of $x(n)$ may or may not interfere with inference in the actually observed components. Finally, motivated from a practical standpoint to analyze causation in large systems, practical estimation procedures based purely on pairwise causality tests are of interest since the computation of such pairwise relations is substantially easier than the full conditional relations.

The following propositions are essentially lemmas used for the proof of the upcoming Proposition 2.4, but remain relevant for providing intuitive insight into the problems at hand.

The previous result can still be strengthened significantly; notice that it is possible to have some $k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}$ where still $j \overset{\text{PW}}{\nrightarrow } i$ , an example is furnished by the three node graph $k \rightarrow i \rightarrow j$ where clearly $k \in{\mathcal{A}({i})}\cap{\mathcal{A}({j})}$ but $j \overset{\text{PW}}{\nrightarrow } i$ . We must introduce the concept of a confounding variable, which effectively eliminates the possibility presented in this example.

A similar proposition as the following also appears as [Datta-Gupta, (Reference Datta-Gupta and Mazumdar2014), Prop. 5.3.2]:

Proof. The proof is by way of contradiction. To this end, suppose that $j$ is a node such that:

1. (a) $j \not \in{\mathcal{A}({i})}$

2. (b) every $k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}$ every $k \rightarrow \cdots \rightarrow j$ path contains $i$ .

Firstly, notice that every $u \in \big ({pa(j)} \setminus \{i\}\big )$ necessarily inherits these same two properties. This follows since if we also had $u \in{\mathcal{A}({i})}$ then $u \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}$ so by our assumption every $u \rightarrow \cdots \rightarrow j$ path must contain $i$ , but $u \in{pa(j)}$ so $u \rightarrow j$ is a path that doesn’t contain $i$ , and therefore $u \not \in{\mathcal{A}({i})}$ ; moreover, if we consider $w \in{\mathcal{A}({i})} \cap{\mathcal{A}({u})}$ then we also have $w \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}$ so the assumption implies that every $w \rightarrow \cdots \rightarrow j$ path must contain $i$ . These properties therefore extend inductively to every $u \in \big ({\mathcal{A}({j})} \setminus \{i\}\big )$ .

In order to deploy a recursive argument, define the following partition of $pa(u)$ , for some node $u$ :

\begin{align*} C_0(u) &= \{k \in{pa(u)}\ |\ i \not \in{\mathcal{A}({k})},{\mathcal{A}({i})} \cap{\mathcal{A}({k})} = \emptyset, k \ne i\}\\ C_1(u) &= \{k \in{pa(u)}\ |\ i \in{\mathcal{A}({k})} \text{ or } k = i\}\\ C_2(u) &= \{k \in{pa(u)}\ |\ i \not \in{\mathcal{A}({k})},{\mathcal{A}({i})} \cap{\mathcal{A}({k})} \ne \emptyset, k \ne i\} \end{align*}

We notice that for any $u$ having the properties $(a), (b)$ above, we must have $C_2(u) = \emptyset$ since if $k \in C_2(u)$ then $\exists w \in{\mathcal{A}({i})} \cap{\mathcal{A}({k})}$ (and $w \in{\mathcal{A}({i})} \cap{\mathcal{A}({u})}$ , since $k \in{pa(u)}$ ) such that $i \not \in{\mathcal{A}({k})}$ and therefore there must be a path ${w \rightarrow \cdots \rightarrow k}\rightarrow u$ which does not contain $i$ , contradicting property $(b)$ . Moreover, for any $u \in \big ({\mathcal{A}({j})} \setminus \{i\}\big )$ and $k \in C_0(u)$ , Proposition 2.3 shows that $H_n^{(i)} \perp \mathcal{H}_{n-1}^{(j)} |\mathcal{H}_{n-1}^{(i)}$ .

In order to establish $j \overset{\text{PW}}{\nrightarrow } i$ , choose an arbitrary element of $\mathcal{H}_{t - 1}^{(j)}$ and represent it via the action of a strictly causal filter $\Phi (z)$ , that is, $\Phi (z) x_j(n) \in \mathcal{H}_{n-1}^{(j)}$ , by Theorem 2.1 it suffices to show that

(10) \begin{equation} \langle{x_i(n)},{\Phi (z)x_j(n) -{\hat{\mathbb{E}}[{\Phi (z)x_j(n)}\ |\ {\mathcal{H}_{t - 1}^{(i)}}]}}\rangle = 0 \end{equation}

Denote $e_j \;:\!=\; e_j^{S(j)}$ , we can write $x_j(n) = e_j^{\mathsf{T}} x_{S(j)}(n)$ , and therefore from Equation (D3) there exist strictly causal filters $\Gamma _s(z)$ and $\Lambda _{sk}(z)$ (defined for ease of notation) such that

\begin{equation*} x_j(n) = \sum _{s \in S(j)} \Gamma _s(z)v_s(n) + \sum _{\substack {s \in S(j) \\ k \in {pa(s)} \cap S(j)^{\mathsf {c}}}} \Lambda _{sk}(z)x_k(n) \end{equation*}

When we substitute this expression into the left hand side of Equation (10), we may cancel each term involving $v_s$ by Lemma D.2, and each $k \in C_0(s)$ by our earlier argument, leaving us with

\begin{equation*} \sum _{\substack {s \in S(j) \\ k \in C_1(s) \cap S(j)^{\mathsf {c}}}}\langle {x_i(n)},{\Phi (z)\Lambda _{sk}(z)x_k(n) - {\hat {\mathbb {E}}[{\Phi (z)\Lambda _{sk}(z)x_k(n)}\ |\ {\mathcal {H}_{t - 1}^{(i)}}]}}\rangle \end{equation*}

Since each $k \in C_1(s)$ with $k \ne i$ inherits properties $(a)$ and $(b)$ above, we can recursively expand each $x_k$ of the above summation until reaching $k = i$ (which is guaranteed to terminate due to the definition of $C_1(u)$ ) which leaves us with some strictly causal filter $F(z)$ such that the left hand side of Equation (10) is equal to

\begin{equation*} \langle {x_i(n)},{\Phi (z)F(z)x_i(n) - {\hat {\mathbb {E}}[{\Phi (z)F(z)x_i(n)}\ |\ {\mathcal {H}_{t - 1}^{(i)}}]}}\rangle \end{equation*}

and this is $0$ since $\Phi (z)F(z)x_i(n) \in \mathcal{H}_{n-1}^{(i)}$ .

An interesting corollary is the following:

It seems reasonable to expect a converse of Proposition 2.4 to hold, that is, $j \in{\mathcal{A}({i})} \Rightarrow j \overset{\text{PW}}{\rightarrow } i$ . Unfortunately, this is not the case in general, as different paths through $\mathcal{G}$ can lead to cancelation (see Figure 1(a)). In fact, we do not even have $j \in{pa(i)} \Rightarrow j \overset{\text{PW}}{\rightarrow } i$ (see Figure 1(b)).

Figure 1. Examples illustrating the difficulty of obtaining a converse to Proposition 2.4. (a) Path cancelation: $j \in{\mathcal{A}({i})} \nRightarrow j \overset{\text{PW}}{\rightarrow } i$ . (b) Cancellation from time lag: $j \in{pa(i)} \nRightarrow j \overset{\text{PW}}{\rightarrow } i$ .

2.5 Strongly causal graphs

In this section and the next we will seek to understand when converse statements of Proposition 2.4 do hold. One possibility is to restrict the coefficients of the system matrix, for example, by requiring that $B_{ij}(m) \ge 0$ . Instead, we think it more meaningful to focus on the defining feature of time series networks, that is, the topology of $\mathcal{G}$ . To this end, we introduce the following notion of a strongly causal graph (SCG).

Examples of strongly causal graphs include directed trees (or forests), the diamond shaped graph containing $i \rightarrow j \leftarrow k$ and $i \rightarrow j' \leftarrow k$ , and Figure 4 of this paper. Importantly, a maximally connected bipartite graph with $2N$ nodes is also strongly causal, demonstrating both that the number of edges of such a graph can still scale quadratically with the number of nodes. It is evident that the strong causal property is inherited by subgraphs.

Figure 2. Reproduced under the creative commons attribution non-commercial license (http://creativecommons.org/licenses/by-nc/2.5) is a known gene regulatory network for a certain strain of E.Coli. The graph has only a single edge violating strong causality.

Example 2.1. Though examples of SCGs are easy to construct in theory, should practitioners expect SCGs to arise in application? While a positive answer to this question is not necessary for the concept to be useful, it is certainly sufficient. Though the answer is likely to depend upon the particular application area, examples appear to be available in biology, in particular, the authors of Shojaie & Michailidis (Reference Shojaie and Michailidis2010) cite an example of the so called “transcription regulatory network of E.coli” (reproduced in Figure 2), and Ma et al. (Reference Ma, Kemmeren, Gresham and Statnikov2014) study a much larger regulatory network of Saccharomyces cerevisiae. These networks, which we reproduce in Figure 3, appear to have at most a small number of edges which violate the strong causality condition.

Figure 3. Reproduced under the creative commons attribution license (https://creativecommons.org/licenses/by/4.0/) is a much larger gene regulatory network. It exhibits similar qualitative structure (a network of hub nodes) as does the much smaller network of Figure 2.

For later use, and to get a feel for the topological implications of strong causality, we explore a number of properties of such graphs before moving into the main result of this section. The following important property essentially strengthens Proposition 2.4 for the case of strongly causal graphs.

In light of Proposition 2.5, the following provides a partial converse to Proposition 2.4, and supports the intuition of “causal flow” through paths in $\mathcal{G}$ .

Proof. We will show that for some $\psi \in \mathcal{H}_{n-1}^{(j)}$ we have

(11) \begin{equation} \langle{\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{x_i(n) -{\hat{\mathbb{E}}[x_i(n)\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}}\rangle \ne 0 \end{equation}

and therefore that $H_n^{(i)} \not \perp \ \mathcal{H}_{n-1}^{(j)}\ |\ \mathcal{H}_{n-1}^{(i)}$ , which by Theorem (2.1) is enough to establish that $j \overset{\text{PW}}{\rightarrow } i$ .

Firstly, we will establish a representation of $x_i(n)$ that involves $x_j(n)$ . Denote by $a_{r + 1} \rightarrow a_r \rightarrow \cdots \rightarrow a_1 \rightarrow a_0$ with $a_{r + 1} \;:\!=\; j$ and $a_0 \;:\!=\; i$ the unique $j \rightarrow \cdots \rightarrow i$ path in $\mathcal{G}$ , we will expand the representation of Equation (7) backwards along this path:

\begin{align*} x_i(n) &= v_i(n) + \mathsf{B}_{ii}(z)x_i(n) + \sum _{k \in{pa(i)}}\mathsf{B}_{ik}(z) x_k(n)\\ &= \underbrace{v_{a_0}(n) + \mathsf{B}_{a_0a_0}(z)x_i(n) + \sum _{\substack{k \in{pa({a_0})} \\ k \ne a_1}}\mathsf{B}_{a_0 k}(z) x_k(n)}_{\;:\!=\;{\widetilde{\alpha }({a_0},{a_1})}} + \mathsf{B}_{a_0a_1}(z)x_{a_1}(n)\\ &={\widetilde{\alpha }({a_0},{a_1})} + \mathsf{B}_{a_0a_1}(z)\big [{\widetilde{\alpha }({a_1},{a_2})} + \mathsf{B}_{a_1a_2}(z)x_{a_2}(n) \big ]\\ &\overset{(a)}{=} \sum _{\ell = 0}^r \underbrace{\Big (\prod _{m = 0}^{\ell - 1} \mathsf{B}_{A(m) a_{m + 1}}(z) \Big )}_{\;:\!=\; F_\ell (z)}{\widetilde{\alpha }({a_\ell },{a_{\ell + 1}})} + \Big (\prod _{m = 0}^{r}\mathsf{B}_{A(m) a_{m + 1}}(z)\Big )x_{a_{r + 1}}(n)\\ &= \sum _{\ell = 0}^r F_\ell (z){\widetilde{\alpha }({a_\ell },{a_{\ell + 1}})} + F_{r + 1}(z) x_j(n) \end{align*}

where $(a)$ follows by a routine induction argument and where we define $\prod _{m = 0}^{-1} \bullet \;:\!=\; 1$ for notational convenience.

Using this representation to expand Equation (11), we obtain the following cumbersome expression:

\begin{align*} &\langle{\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{F_{r + 1}(z)x_j(n) -{\hat{\mathbb{E}}[{F_{r + 1}(z)x_j(n)}\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}}\rangle \\ &- \langle{\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{{\hat{\mathbb{E}}[{\sum _{\ell = 0}^r F_\ell (z){\widetilde{\alpha }({a_\ell },{a_{\ell + 1}})}}\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}}\rangle \\ &+ \langle{\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{\sum _{\ell = 0}^r F_\ell (z){\widetilde{\alpha }({a_\ell },{a_{\ell + 1}})}}\rangle \end{align*}

Note that by the orthogonality principle, $\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]} \perp \mathcal{H}_{n-1}^{(i)}$ , the middle term above is $0$ . Choosing now the particular value $\psi = F_{r + 1}(z)x_j(n) \in \mathcal{H}_{n-1}^{(j)}$ we arrive at

\begin{align*} &\langle{\psi -{\hat{\mathbb{E}}[{\psi }\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{x_i(n) -{\hat{\mathbb{E}}[{x_i(n)}\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}}\rangle \\ &= \mathbb{E}|F_{r + 1}(z)x_j(n) -{\hat{\mathbb{E}}[{F_{r + 1}(z)x_j(n)}\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}|^2\\ &+ \langle{F_{r + 1}(z)x_j(n) -{\hat{\mathbb{E}}[{F_{r + 1}(z)x_j(n)}\ |\ {\mathcal{H}_{n - 1}^{(i)}}]}},{\sum _{\ell = 0}^r F_\ell (z){\widetilde{\alpha }({a_\ell },{a_{\ell + 1}})}}\rangle \end{align*}

Now since $F_{r + 1}(z) \ne 0$ by Theorem 2.1, and $F_{r + 1}(z)x_j(n) \not \in \mathcal{H}_{n-1}^{(i)}$ , we have by the Cauchy-Schwarz inequality that this expression is equal to $0$ if and only if

\begin{equation*} \sum _{\ell = 0}^r F_\ell (z) {\widetilde {\alpha }({a_\ell }, {a_{\ell + 1}})} \overset {\text {a.s.}}{=} {\hat {\mathbb {E}}[{F_{r + 1}(z)x_j(n)}\ |\ {\mathcal {H}_{n - 1}^{(i)}}]} - F_{r + 1}(z)x_j(n) \end{equation*}

or by rearranging and applying the representation for $x_i(n)$ obtained earlier, if and only if

\begin{equation*} x_i(n) \overset {\text {a.s.}}{=} {\hat {\mathbb {E}}[{F_{r + 1}(z)x_j(n)}\ |\ {\mathcal {H}_{n - 1}^{(i)}}]} \end{equation*}

But, this is impossible since $x_i(n) \not \in\mathcal{H}_{n - 1}^{(i)}$ .

We immediately obtain the corollary, which we remind the reader is, surprisingly, not true in a general graph.

Example 2.2. As a final remark of this subsection we note that a complete converse to Proposition 2.4 is not possible without additional conditions. Consider the “fork” system on $3$ nodes (i.e., $2 \leftarrow 1 \rightarrow 3$ ) defined by

\begin{equation*} x(n) = \left [ \begin {array}{c@{\quad}c@{\quad}c} 0 & 0 & 0\\ a & 0 & 0\\ a & 0 & 0\\ \end {array} \right ] x(n - 1) + v(n) \end{equation*}

In this case, node $1$ is a confounder for nodes $2$ and $3$ , but $x_3(n) = v_3(n) - v_2(n) + x_2(n)$ and $2 \overset{\text{PW}}{\nrightarrow } 3$ (even though $x_2(n)$ and $x_3(n)$ are contemporaneously correlated).

If we were to augment this system by simply adding an autoregressive component (i.e., some “memory”) to $x_1(n)$ , for example, $x_1(n) = v_1(n) + b x_1(n - 1)$ then we would have $2 \overset{\text{PW}}{\rightarrow } 3$ since then $x_3(n) = v_3(n) + av_1(n - 1) - bv_2(n - 1) + bx_2(n - 1)$ . We develop this idea further in the next section.

2.6 Persistent systems

In Section 2.5 we obtained a converse to part $(a)$ of Proposition 2.4 via the notion of a strongly causal graph topology (see Proposition 2.6). In this section, we study conditions under which a converse to part $(b)$ will hold.

Definition 2.10 (Lag Function). Given a causal filter $\mathsf{B}(z) = \sum _{m = 0}^\infty b(m)z^{-m}$ define

(12) \begin{align} m_0(\mathsf{B}) &= \text{inf }\{m \in \mathbb{Z}_+\ |\ b(m) \ne 0\} \end{align}

(13) \begin{align} m_{\infty }(\mathsf{B}) &= \text{sup }\{m \in \mathbb{Z}_+\ |\ b(m) \ne 0\} \end{align}

that is, the “first” and “last” coefficients of the filter $\mathsf{B}(z)$ , where $m_\infty (\mathsf{B}) \;:\!=\; \infty$ if the filter has an infinite length, and $m_0(\mathsf{B}) \;:\!=\; \infty$ if $\mathsf{B}(z) = 0$ .

In order to eliminate the possibility of a particular sort of cancelation, an ad hoc assumption is required. Strictly speaking, the persistence condition is not a necessary or sufficient condition for the following, but cases where the following fails to hold, and persistence does hold, are unavoidable pathologies.

Assumption 2.1 (no-cancellation). Fix $i, j \in [N]$ and let $H_i(z)$ be the strictly causal filter such that

\begin{equation*} H_i(z)x_i(n) = {\hat {\mathbb {E}}[{x_i(n)}\ |\ {\mathcal {H}_{t - 1}^{(i)}}]} \end{equation*}

and similarly for $H_j(z)$ . Then define

(14) \begin{equation} T_{ij}(z)= \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\sigma _k^2\mathsf{A}_{ik}(z^{-1}) (1 - H_i(z^{-1})) (1 - H_j(z))\mathsf{A}_{jk}(z) \end{equation}

where $\sigma _k^2 = \mathbb{E} v_k(n)^2$ .

We will say that Assumption 2.1 is satisfied if for every $i, j \in [N]$ , $T_{ij}(z)$ is either constant over $z$ (i.e., each $z^k$ coefficient for $k \in \mathbb{Z} \setminus \{0\}$ is $0$ ), or is neither causal (i.e., containing only $z^{-k}$ terms, for $k \ge 0$ ) or anti-causal (i.e., containing only $z^k$ terms, for $k \ge 0$ ). Put succinctly, $T_{ij}(z)$ must be two-sided.

Proof. Recalling Theorem 2.1, consider some $\psi \in \mathcal{H}_{n-1}^{(j)}$ and represent it as $\psi (n) = F(z)x_j(n)$ for some strictly causal filter $F(z)$ . Then

\begin{align*} &\langle{\psi (n) -{\hat{\mathbb{E}}[{\psi (n)}\ |\ {\mathcal{H}_{t - 1}^{(i)}}]}},{x_i(n) -{\hat{\mathbb{E}}[{x_i(n)}\ |\ {\mathcal{H}_{t - 1}^{(i)}}]}}\rangle \\ &\overset{(a)}{=} \langle{F(z)x_j(n)},{x_i(n) -{\hat{\mathbb{E}}[{x_i(n)}\ |\ {\mathcal{H}_{t - 1}^{(i)}}]}}\rangle \\ &\overset{(b)}{=} \langle{F(z)\big (\mathsf{A}_{jj}(z)v_j(n) + \sum _{k \in{\mathcal{A}({j})}}\mathsf{A}_{jk}(z)v_k(n)\big )},{(1 - H_i(z))\big (\mathsf{A}_{ii}(z)v_i(n) + \sum _{\ell \in{\mathcal{A}({i})}}\mathsf{A}_{i\ell }(z)v_\ell (n)\big )}\rangle \\ &\overset{(c)}{=} \sum _{k \in{\mathcal{A}({i})}\cap{\mathcal{A}({j})}}\langle{F(z)\mathsf{A}_{jk}(z)v_k(n)},{(1 - H_i(z))\mathsf{A}_{ik}(z)v_k(n)}\rangle, \end{align*}

where $(a)$ applies the orthogonality principle, $(b)$ expands with Equation (8) with $H_i(z)x_i(n) ={\hat{\mathbb{E}}[{x_i(n)}\ |\ {\mathcal{H}_{t - 1}^{(i)}}]}$ , and $(c)$ follows by performing cancelations of $v_k(n) \perp v_\ell (n)$ and noting that by the contra-positive of Proposition 2.5 we cannot have $i \in{\mathcal{A}({j})}$ or $j \in{\mathcal{A}({i})}$ .

Through symmetric calculation, we can obtain the expression relevant to the determination of $i \overset{\text{PW}}{\rightarrow } j$ for $\phi \in \mathcal{H}_{n-1}^{(i)}$ represented by the strictly causal filter $G(z)\;:\; \phi (n) = G(z)x_i(n)$

\begin{align*} &\langle{\phi (n) -{\hat{\mathbb{E}}[{\phi (n)}\ |\ {\mathcal{H}_{t - 1}^{(j)}}]}},{x_j(n) -{\hat{\mathbb{E}}[{x_j(n)}\ |\ {\mathcal{H}_{t - 1}^{(j)}}]}}\rangle \\ &= \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{G(z)\mathsf{A}_{ik}(z)v_k(n)},{(1 - H_j(z))\mathsf{A}_{jk}(z)v_k(n)}\rangle \end{align*}

where $H_j(z)x_j(n) ={\hat{\mathbb{E}}[{x_j(n)}\ |\ {\mathcal{H}_{t - 1}^{(j)}}]}$ .

We have therefore

(15) \begin{align} (j \overset{\text{PW}}{\rightarrow } i)\;:\; \exists F(z) \text{ s.t. } \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{F(z)\mathsf{A}_{jk}(z)v_k(n)},{(1 - H_i(z))\mathsf{A}_{ik}(z)v_k(n)}\rangle \ne 0 \end{align}

(16) \begin{align} (i \overset{\text{PW}}{\rightarrow } j)\;:\; \exists G(z) \text{ s.t. } \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{G(z)\mathsf{A}_{ik}(z)v_k(n)},{(1 - H_j(z))\mathsf{A}_{jk}(z)v_k(n)}\rangle \ne 0 \end{align}

The persistence condition, by Corollary D.1, ensures that for each $k \in{\mathcal{A}({i})}\cap{\mathcal{A}({j})}$ there is some $F(z)$ and some $G(z)$ such that at least one of the above terms constituting the sum over $k$ is non-zero. It remains to eliminate the possibility of cancelation in the sum.

The adjoint of a linear filter $C(z)$ is simply $C(z^{-1})$ , which recall is strictly anti-causal if $C(z)$ is strictly causal. Using this, we can write

\begin{align*} &\sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{F(z)\mathsf{A}_{jk}(z)v_k(n)},{(1 - H_i(z))\mathsf{A}_{ik}(z)v_k(n)}\rangle \\ = &\sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{\mathsf{A}_{ik}(z^{-1})(1 - H_i(z^{-1}))F(z)\mathsf{A}_{jk}(z)v_k(n)},{v_k(n)}\rangle \end{align*}

Moreover, it is sufficient to find some strictly causal $F(z)$ of the form $F(z)(1 - H_j(z))$ (abusing notation) since $1 - H_j(z)$ is causal. Similarly for $G(z)$ , this leads to symmetric expressions for $j \overset{\text{PW}}{\rightarrow } i$ and $i \overset{\text{PW}}{\rightarrow } j$ respectively:

(17) \begin{equation} \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{\mathsf{A}_{ik}(z^{-1})(1 - H_i(z^{-1}))F(z)(1 - H_j(z))\mathsf{A}_{jk}(z)v_k(n)},{v_k(n)}\rangle \end{equation}

(18) \begin{equation} \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}}\langle{\mathsf{A}_{ik}(z^{-1})(1 - H_i(z^{-1}))G(z^{-1})(1 - H_j(z))\mathsf{A}_{jk}(z)v_k(n)},{v_k(n)}\rangle \end{equation}

Recall the filter from Assumption 2.1

(19) \begin{equation} T_{ij}(z) = \sum _{k \in{\mathcal{A}({i})} \cap{\mathcal{A}({j})}} \sigma _k^2\mathsf{A}_{ik}(z^{-1})(1 - H_i(z^{-1}))(1 - H_j(z))\mathsf{A}_{jk}(z) \end{equation}

Since each $v_k(n)$ is uncorrelated through time, $\langle{T_{ij}(z)v_k(n)},{v_k(n)}\rangle = \sigma _k^2T_{ij}(0)$ , and therefore we have $j \overset{\text{PW}}{\rightarrow } i$ if $T_{ij}(z)$ is not causal and $i \overset{\text{PW}}{\rightarrow } j$ if $T_{ij}(z)$ it not anti-causal. Moreover, we have $i \overset{\text{PW}}{\nrightarrow } j$ and $j \overset{\text{PW}}{\rightarrow } i$ if $T_{ij}(z)$ is a constant. Therefore, under Assumption 2.1 $j \overset{\text{PW}}{\rightarrow } i \iff i \overset{\text{PW}}{\rightarrow } j$ .

This follows since if $T_{ij}(z)$ is not causal then $\exists k \gt 0$ such that the $z^k$ coefficient of $T_{ij}(z)$ is non-zero, and we can choose strictly causal $F(z) = z^{-k}$ such that (D11) is non-zero and therefore $j \overset{\text{PW}}{\rightarrow } i$ .

Similarly, if $T_{ij}(z)$ is not anti-causal, then $\exists k \gt 0$ such that the $z^{-k}$ coefficient of $T_{ij}(z)$ is non-zero, and we can choose strictly causal $G(z)$ so that $G(z^{-1}) = z^k$ , and then (D12) is non-zero and therefore $i \overset{\text{PW}}{\rightarrow } j$ .

2.7 Recovering $\mathcal{G}$ via pairwise tests

We arrive at the main conclusion of the theoretical analysis in this paper.

Algorithm 1: Pairwise Granger Causality Algorithm

Example 2.4. The set $W$ of Algorithm 1 collects ancestor relations in $\mathcal{G}$ (see Lemma D.5). In reference to Figure 4, each of the solid black edges, as well as the dotted red edges will be included in $W$ , but not the bidirectional green dash-dotted edges, which we are able to exclude as results of confounding. The groupings $P_0, \ldots, P_3$ are also indicated in Figure 4.

Figure 4. Black arrows indicate true parent-child relations. Red dotted arrows indicate pairwise causality (due to non-parent relations), green dash-dotted arrows indicate bidirectional pairwise causality (due to the confounding node $1$ ). Blue groupings indicate each $P_k$ in Algorithm 1.

The algorithm proceeds first with the parent-less nodes $1, 2$ on the initial iteration where the edge $(1, 3)$ is added to $E$ . On the next iteration, the edges $(3, 4), (2, 4), (3, 5)$ are added, and the false edges $(1, 4), (1, 5)$ are excluded due to the paths $1 \rightarrow 3 \rightarrow 4$ and $1 \rightarrow 3 \rightarrow 5$ already being present. Finally, edge $(4, 6)$ is added, and the false $(1, 6), (3, 6), (2, 6)$ edges are similarly excluded due to the ordering of the inner loop.

That we need to proceed backwards through $P_{k - r}$ as in the inner loop on $r$ can also be seen from this example, where if instead we simply added the set

\begin{equation*} D_k' = \{(i, j) \in \Big (\bigcup _{r = 1}^k P_{k - r}\Big ) \times P_k\ |\ i \overset {\text {PW}}{\rightarrow } j \} \end{equation*}

to $E_k$ then we would infer the false positive edge $1 \rightarrow 4$ . Moreover, the same example shows that simply using the set

\begin{equation*} D_k'' = \{(i, j) \in P_{k - 1} \times P_k\ |\ i \overset {\text {PW}}{\rightarrow } j \} \end{equation*}

causes the edge $1 \rightarrow 3$ to be missed.

3.1 Results

Our first simulation results are summarized in Table 1. It is clear that the superior performance of PWGC in comparison to AdaLASSO is as a result of limiting the false discovery rate. It is unsurprising that PWGC exhibits superior performance when the graph is an SCG, but even in the case of more general DAGs, the PWGC heuristic is still able to more reliably uncover the graph structure for small values of $q$ . We would conjecture that for small $q$ , random graphs are “likely” to be “close” to SCGs in some appropriate sense. As $q$ increases, there are simply not enough edges allowed by the SCG topology for it to be possible to accurately recover $\mathcal{G}$ .

Table 1. Simulation results—PWGC vs AdaLASSO

1. Results of Monte Carlo simulations comparing PWGC and AdaLASSO $(N = 50, p = 5, p_{\text{max}} = 10)$ for small samples and when the SCG assumption doesn’t hold. The superior result is bolded when the difference is statistically significant, as measured by scipy.stats.ttest_rel. 100 iterations are run for each set of parameters.

2. LRE: Log-Relative-Error, that is, the log sum of squared errors at each node relative to the strength of the driving noise $\frac{\mathsf{ln\ } \mathsf{tr } \widehat{\Sigma }_v}{\mathsf{ln\ } \mathsf{tr } \Sigma _v}$ . MCC: Matthew’s Correlation Coefficient.

3. Values of $q$ (edge probability) range between $2/ N, 4/ N, 16/ N$ where $2/ N$ has the property that the random graphs have on average the same number of edges as the SCG.

Figures 6 and 7 provide an empirical analysis of how the algorithm scales with increasing $N$ , and for a small number of training samples $T$ . It is surprising that even for general DAGs (which don’t necessarily have the SCG property) performance degrades extremely slowly with respect to $N$ , and that performance is markedly better than random even when $N \gt T$ . These observations are reminiscent of results for the LASSO algorithm applied to linear regression (Tibshirani et al., Reference Tibshirani, Wainwright and Hastie2015). For example, under certain conditions, it holds that

\begin{align*}||\hat {\beta }_{LASSO} - \beta ^\star ||_2^2 = \mathcal {O}\bigl (\frac {s}{T}\ \mathsf {ln\ } N\bigr ) \end{align*}

where $s$ is the number of non-zero parameters of $\beta ^\star$ in the ground-truth linear regression model (Wainwright, Reference Wainwright2009).

Figure 6. Measures support recovery performance as the number of nodes $N$ increases, and the edge proportion as well as the number of samples $T$ is held fixed. Remarkably, the degradation as $N$ increases is limited, it is primarily the graph topology (SCG or non-SCG) as well as the level of sparsity (measured by $q$ ) which are the determining factors for support recovery performance.

Figure 7. Provides a support recovery comparison for very small values of, $T$ typical for many applications. The system has a fixed number $N=50$ nodes.

Additional experimental results, and a detailed description of the heuristic algorithm based on the theory of Section 2.7 is available in Appendix C.

4.1 Classifier methodology and evaluation

The classifier we apply is an off-the-shelf polynomial kernel logistic regression model: sklearn.linear_model.LogisticRegression (Pedregosa et al., Reference Pedregosa, Varoquaux, Gramfort, Michel, Thirion and Grisel2011), using a Nystroem kernel approximation (with 240 components): sklearn.kernel_approximation.Nystroem(kernel = ‘poly’). Hyperparameters are chosen via 10-fold cross validation with Bayesian Optimization (Kumar & Head, Reference Kumar and Head2017), using $80\%$ of the data. The remaining $20\%$ of the data is held out for a final validation step. A polynomial kernel is natural for this task since the input features are binary (the adjacency matrix); a kernel degree of, say, $3$ , indicates that graphical substructures no more complicated than triangles can be used by the classifier. Similarly to our simulation results, we use the MCC as the metric for evaluation.

As a final verification, a label permutation test (“Test 1” in Ojala & Garriga (Reference Ojala and Garriga2010)) is applied to verify that $A^{(i)}$ and ( $z^{(i)})$ are not independent, and to alleviate the fear that the classifier may simply be picking up on spurious patterns in the estimated causality graphs. This test consists of refit our classifier on randomly relabeled data and comparing the results those of the classifier fit on the correctly labeled data. An inability of the model to accurately discriminate between the labels in an incorrectly labeled data set is further evidence that accurate classification is not simply a result of overfitting or of picking up on spurious variation in the data $A^{(i)}$ .

4.2 Results

The final value of the MCC for the alcoholism classification task, computed on a $10$ -fold cross validation is $MCC = 0.305$ , with the individual folds falling into the interval $[0.229, 0.396]$ . The MCC obtained on the final held-out validation set (i.e., data which was not used for picking any model parameters) was $0.29$ . The p-value estimated from the label permutation test satisfies $p \lt 0.01$ and where the average MCC obtained on the validation data, by the 100 classifiers fit on data with randomly permuted labels, was zero through two decimal places. While the performance of the classifier is not exemplary, these results establish clearly that the performance is statistically significantly better than random.

Figure 8. Is a low-dimensional embedding representing the classification of alcoholic and non-alcoholic subjects based entirely on the Granger causality graphs constructed from EEG signals. Figure 8b further illustrates a classification task using EEG causality graphs. In this case, one particular subject can be almost perfectly separated from another particular subject with a binary classifier. These results indicate the causality graph carries consistent and identifiable cross-subject patterns.

The EEG alcoholism data set has other potential classification labels. One natural example is to classify the subject label. That is, to identify the subject to which the EEG example belongs (each subject had multiple EEG recordings taken). We do not report the details as the task is considerably easier and the methodology is similar to that described in Section 4.1. However, a classifier trained on any two specific subjects is able to discern out-of-sample causality graphs between the two subjects with almost perfect accuracy. Specifically, we obtained an MCC of $0.99$ and the p-value estimated by a permutation test satisfied $p \lt 0.001$ .

An illustration of the alcoholism classifier, as well as a particular example from the subject-subject classifier is provided in Figure 8. The visualization is constructed by projecting onto two basis dimensions by partial least squares (a supervised dimensionality reduction technique). These visualizations are provided entirely for qualitative context, but drawing substantive conclusions from them is not appropriate. In particular, the low-dimensional visualization does not necessarily accurately capture the high-dimensional decision surface utilized by the classifiers.