# Learning latent causal relationships in
multiple time series

Jacek P. Dmochowski
Department of Biomedical Engineering
City College of New York
New York, NY 10031
<EMAIL_ADDRESS>

###### Abstract

Identifying the causal structure of systems with multiple dynamic elements is
critical to several scientific disciplines. The conventional approach is to
conduct statistical tests of causality, for example with Granger Causality,
between observed signals that are selected _a priori_. Here it is posited
that, in many systems, the causal relations are embedded in a _latent_ space
that is expressed in the observed data as a linear mixture. A technique for
blindly identifying the latent sources is presented: the observations are
projected into pairs of components – driving and driven – to maximize the
strength of causality between the pairs. This leads to an optimization problem
with closed form expressions for the objective function and gradient that can
be solved with off-the-shelf techniques. After demonstrating proof-of-concept
on synthetic data with known latent structure, the technique is applied to
recordings from the human brain and historical cryptocurrency prices. In both
cases, the approach recovers multiple strong causal relationships that are not
evident in the observed data. The proposed technique is unsupervised and can
be readily applied to any multiple time series to shed light on the causal
relationships underlying the data.

_K_ eywords Granger Causality $\cdot$ blind source separation $\cdot$
unsupervised learning

##  1 Introduction

Several important problems in the sciences are concerned with the activity of
interacting sources. Two prominent examples are the dynamics of the brain [1],
where information processing manifests as a spatiotemporal pattern of regional
activations constrained by anatomical and functional connectivity, and
financial markets [2], where values of assets evolve in concert with the
decisions of agents. In complex systems such as these, it is critical to infer
the underlying structure governing system evolution, and relatedly, to
forecast future outcomes.

Granger Causality [3] is a popular technique for measuring a form of
dependence rooted in the temporal precedence effect. Time series $x(t)$ is
said to cause, in a Granger sense, time series $y(t)$ if the past of $x$
improves the prediction of the present value of $y$ above that of its own
past. Originating in economics[4], Granger Causality has since found extensive
utilization in neuroscience [5, 6], where it has been applied to recordings of
brain activity captured at various spatial and temporal scales to illuminate
neural circuits [5, 7, 8, 9, 10, 11]. Perhaps driven by the ubiquitous
interest in causal interactions, the technique has been adopted by many
disparate fields, including ecology [12], computational biology [13], and
epidemiology [14, 15]. The utility of Granger Causality has been aided by
several extensions and reformulations of the original technique, most notably
a frequency-domain formulation [16, 17] and a generalization to multivariate
time series [18, 19]. Moreover, several approaches to capturing non-linear
causal interactions between multiple time series have been proposed [20, 21,
22, 23]. Conventionally, these different variants of Granger Causality are
measured between observed signals that are selected _a priori_. In other
words, one must specify the identity of the signals being probed, and the
hypothesized direction of causality. Moreover, this approach implicitly
assumes that the underlying causal relationships exist in the native space
defined by the observations (e.g. the sensors).

The central idea proposed here is that, in many systems, the true causal
relations are embedded in a _latent_ source space, and that these latent
sources enter the observations via an unknown linear mixture. Due to the
mixing process, direct application of tools such as Granger Causality to the
observed data may not optimally reveal the dynamics of the system. Rather, the
approach taken here to identify the latent causal sources is to project the
observations into a component space that maximizes the Granger Causality among
_pairs_ of time series: one signal models the “driving” source, and the other
captures the source being “driven”. It is shown that this can be formulated as
a non-convex optimization problem with closed-form expressions for the
objective function and gradient. Importantly, the optimization does not
require access to the mixing process and thus constitutes blind
identification.

To solve the optimization problem, a simple coordinate descent algorithm that
is implemented with standard numerical packages is presented. By simulating a
vector autoregressive (VAR) system with known structure, it is demonstrated
that the proposed technique indeed identifies the underlying sources, their
connections, and the mixing process. To evaluate the proposed approach on
real-world systems, the technique is then applied to data from the human brain
and the cryptocurrency market. In both cases, it is shown that the proposed
technique recovers multiple pairs of signals whose causal strength is
significantly greater than what is found in the observed data.

The distinctions between Granger and true physical causality have been
previously described [24, 25]. In what follows, the terms “causal” and
“causality” are employed for conciseness with the understanding that the
findings presented here pertain to the Granger form of causality.

## Results

### Motivating example

Consider a simple system with two connected sources, $s_{1}$ and $s_{2}$,
where source 1 “Granger causes” source 2, denoted here by $s_{1}\rightarrow
s_{2}$. In the neuroscience context, $s_{1}$ may represent the mass synaptic
activity at a brain region, and $s_{2}$ the activity of a downstream region to
which $s_{1}$ projects. Due to signal mixing (e.g. volume conduction), the
observed signals are modeled as a linear mixture of the two sources:

$\displaystyle\left(\begin{array}[]{cc}x_{1}(t)\\\
x_{2}(t)\end{array}\right)=\left(\begin{array}[]{cc}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{array}\right)\left(\begin{array}[]{cc}s_{1}(t)\\\
s_{2}(t)\end{array}\right),$ (7)

where the 2-by-2 mixing matrix is assumed to be invertible and where sensor
noise has been omitted for the sake of this illustrative example. Note that
$s_{1}$ and $s_{2}$ are mixed together in the captured signals, potentially
confounding the measure of Granger Causality between $x_{1}$ and $x_{2}$.
Given only the observations, the goal is to identify the “driving” signal
$y(t)\approx s_{1}(t)$ and the “driven” signal $z(t)\approx s_{2}(t)$. Writing
(7) in matrix notation as $\mathbf{x}(t)=\mathbf{A}\mathbf{s}(t)$, $s_{1}$ is
exactly recovered if $y(t)={\mathbf{w}^{*}}^{T}\mathbf{x}(t)$, where
$\mathbf{w}^{\ast}$ is a column vector whose elements are the first row of
$\mathbf{A}^{-1}$ and T is the transpose operation. Similarly, $s_{2}$ is
recovered as $z(t)={\mathbf{v}^{*}}^{T}\mathbf{x}(t)$, where
$\mathbf{v}^{\ast}$ is the second row of $\mathbf{A}^{-1}$. The projection
vectors ${\mathbf{w}^{*}}$ and ${\mathbf{v}^{*}}$ undo the mixing process by
combining the observed signals to form latent components that approximate the
underlying sources. The problem considered here is whether it is possible to
recover $s_{1}$ and $s_{2}$ _without_ access to the mixing process
$\mathbf{A}$. Below, a novel criterion for blind source separation that
maximizes the Granger Causality between pairs of component signals is
proposed.

### Maximizing latent Granger Causality

Given an observable, centered random process ${\bf x}(t)\in\mathbb{R}^{D}$,
the goal is to identify latent variables $y(t)={\bf w}^{T}{\bf x}(t)$ and
$z(t)={\bf v}^{T}{\bf x}(t)$ such that the Granger Causality from $y$ to $z$
is maximized. Namely, it is desired to solve the following optimization
problem:

$\max_{{\bf w},{\bf v}}\mathcal{G}_{y\rightarrow z}$ (8)

where

$\mathcal{G}_{y\rightarrow
z}=1-\frac{E\\{\epsilon_{f}^{2}\\}}{E\\{\epsilon_{r}^{2}\\}}$ (9)

is termed the “strength of causality” [3] from $y$ to $z$, $\epsilon_{f}$ is
the residual of a linear regression predicting $z$ from the history of both
$z$ _and_ $y$ (i.e., the “full” model), and $\epsilon_{r}$ is the residual
when regressing $z$ onto only its past (the “reduced” model).
$\mathcal{G}_{y\rightarrow z}$ is bounded between 0 and 1, with
$\mathcal{G}_{y\rightarrow z}=0$ indicating that $y$ does not aid in the
prediction of $z$, and $\mathcal{G}_{y\rightarrow z}=1$ denoting a zero-error
estimate of $z$ from the past of itself and $y$. The optimization in (8) is
aimed at identifying two projection vectors, $\mathbf{w}\in\mathbb{R}^{D}$ and
$\mathbf{v}\in\mathbb{R}^{D}$, such that the resulting pair of latent
variables maximize the strength of causality (9).

In what follows, a history of $L$ samples is assumed, and the temporal
apertures of $y$ and $z$ are defined as:

$\displaystyle\mathbf{y}_{p}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}y(t-1)&\ldots&y(t-L)\end{array}\right)^{T}$
$\displaystyle\mathbf{z}_{p}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{ccc}z(t-1)&\ldots&z(t-L)\end{array}\right)^{T}.$

To arrive at a form of (9) that can be optimized using gradient-based
techniques, note that the minimum mean squared error (MMSE) corresponding to
the full and reduced models are given by [26]:

$\displaystyle\Phi_{f}=E\\{\epsilon_{f}^{2}\\}=\sigma_{z}^{2}-\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}$
$\displaystyle\Phi_{r}=E\\{\epsilon_{r}^{2}\\}=\sigma_{z}^{2}-\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{q},$
(12)

where $\sigma_{z}^{2}=E\left\\{z^{2}(t)\right\\}$ is the mean power of $z$,
$E\left\\{\cdot\right\\}$ denotes mathematical expectation,

$\displaystyle\mathbf{r}=E\left\\{z(t)\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)\right\\}~{}~{}~{}~{}~{}~{}~{}\mathbf{q}=E\left\\{z(t)\mathbf{z}_{p}(t)\right\\}$

are $2L$ and $L$ dimensional covariance vectors between $z$ and the temporal
apertures of the full and reduced models, respectively, and where

$\displaystyle\mathbf{R}=E\left\\{\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)^{T}\right\\}~{}~{}~{}\mathbf{Q}=E\left\\{\mathbf{z}_{p}(t)\mathbf{z}^{T}_{p}(t)\right\\}$

are $2L$-by-$2L$ and $L$-by-$L$ covariance matrices of the predictors in the
full and reduced models, respectively. Importantly,
$\sigma_{z}^{2}=\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}$, $\mathbf{r}$,
$\mathbf{q}$, $\mathbf{R}$ and $\mathbf{Q}$ can each be expressed in terms of
the projection vectors $\mathbf{w}$ and $\mathbf{v}$ and the spatiotemporal
statistics of the observations (see _Supplementary Note 1_):

$\displaystyle\mathbf{r}=\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{q}=\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right)$
(17)
$\displaystyle\mathbf{R}=\left(\begin{array}[]{c}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\left(\begin{array}[]{cc}\mathbf{I}_{L}\otimes\mathbf{v}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}\end{array}\right)~{}~{}~{}~{}~{}~{}~{}\mathbf{Q}=\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right),$
(22)

where

$\displaystyle\mathbf{\Sigma}_{1:L}=\left(\begin{array}[]{cccc}\mathbf{\Sigma}(1)&\mathbf{0}&\ldots&\mathbf{0}\\\
\mathbf{0}&\mathbf{\Sigma}(2)&\ldots&\mathbf{0}\\\
\vdots&\mathbf{0}&\ddots&\vdots\\\
\mathbf{0}&\ldots&\ldots&\mathbf{\Sigma}(L)\\\ \end{array}\right)$

is an $LD$-by-$LD$ block covariance matrix where
$\mathbf{\Sigma}(\tau)=E\left\\{\mathbf{x}(t)\mathbf{x}^{T}(t-\tau)\right\\}$
is the lagged covariance of the observations,

$\displaystyle\tilde{\mathbf{\Sigma}}=\left(\begin{array}[]{cccc}\mathbf{\Sigma}(0)&\mathbf{\Sigma}(-1)&\ldots&\mathbf{\Sigma}(-L+1)\\\
\mathbf{\Sigma}(1)&\mathbf{\Sigma}(0)&\ldots&\mathbf{\Sigma}(-L+2)\\\
\vdots&\mathbf{\Sigma}(1)&\ddots&\vdots\\\
\mathbf{\Sigma}(L-1)&\ldots&\ldots&\mathbf{\Sigma}(0)\\\ \end{array}\right)$

is an $LD$-by-$LD$ block Toeplitz matrix, $\otimes$ denotes the Kronecker
product, $\mathbf{1}_{K}$ is a column vector of $K$ ones, and $\mathbf{I}_{K}$
is the $K$-by-$K$ identity matrix. Substituting (17) into (Maximizing latent
Granger Causality) and the resulting expressions into (9), one arrives at the
following expression for the strength of causality between latent sources $y$
and $z$:

$\displaystyle\mathcal{G}_{y\rightarrow z}=1-$
$\displaystyle\frac{\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}-\left(\begin{array}[]{@{}c@{}}\mathbf{1}_{L}^{T}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{L}^{T}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}^{T}\right)\left(\mathbf{I}_{2L}\otimes\mathbf{v}\right)\left[\left(\begin{array}[]{@{}c@{}}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\left(\begin{array}[]{@{}c@{}c@{}}\mathbf{I}_{L}\otimes\mathbf{v}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}\end{array}\right)\right]^{-1}\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{@{}c@{}}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)}{\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}-\left(\mathbf{1}_{L}^{T}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}^{T}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)\left[\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)\right]^{-1}\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right)}$
(33)

The gradient of (33) has a closed-form that is derived in _Supplementary Note
2_. Conventional optimization tools may then be employed to learn projection
vectors $\mathbf{w}^{\ast}$ and $\mathbf{v}^{\ast}$ that maximize the Granger
Causality between resulting latent signals
$y(t)={\mathbf{w}^{\ast}}^{T}\mathbf{x}(t)$ and
$z(t)={\mathbf{v}^{\ast}}^{T}\mathbf{x}(t)$.

### Identifying latent causal structure

The objective in (33) is non-convex, since $\mathcal{G}_{y\rightarrow
z}\left(\mathbf{w},\mathbf{v}\right)=\mathcal{G}_{y\rightarrow
z}\left(a\mathbf{w},b\mathbf{v}\right)$ for arbitrary real scalars $a$ and
$b$. This follows from the fact that the residual error when predicting $z$
from $y$ is equivalent to that when predicting $bz$ from $ay$ – any scaling
factors will be accommodated by the temporal filter that predicts the driven
signal from the driving signal. Thus, the technique presented here is only
able to identify the latent sources up to a scaling factor. As with other
blind source separation techniques such as independent components analysis
(ICA) [27, 28], it is not possible to recover the scale or sign of the
underlying sources.

Another potential ambiguity when optimizing (33) is related to a known
property of multivariate Granger Causality [18, 19]. Namely, the strength of
causality between $y$ and $z$ is invariant to mixtures of $y$ and $z$ in the
driving signal, such that $\mathcal{G}_{y\rightarrow
z}=\mathcal{G}_{ay+bz\rightarrow cz}$. This means that, without appropriate
modifications to the objective function, maximizing (33) will only identify
$z$. To resolve this ambiguity, one can utilize the concept of time-reversed
Granger Causality [29, 30]. Notice that if $y\rightarrow z$ in
$\mathbf{x}(t)$, then $z\rightarrow y$ in $\mathbf{x}(-t)$. Thus, while the
ambiguity in forward time occurs on $y$, it occurs in reversed time on $z$.
One can therefore combine forward and reversed time into a single objective
function according to:

$\max_{{\bf w},{\bf v}}\mathcal{G}_{y\rightarrow z}+\mathcal{G}_{z\rightarrow
y}^{\mathrm{tr}},$ (34)

where $\mathcal{G}_{z\rightarrow y}^{\mathrm{tr}}$ is the strength of
causality between $z(-t)$ and $y(-t)$.

The non-convexity of the objective function means that a local minimizer of
(34) is not guaranteed to represent a global minimum. Many approaches to non-
convex optimization have been developed, including the use of multiple
starting points [31] and stochastic gradient descent [32]. Here, a grouped
coordinate descent algorithm [33] that maximizes $\mathbf{v}$ and $\mathbf{w}$
in an alternating fashion is proposed: instead of combining $\mathbf{v}$ and
$\mathbf{w}$ into a single model parameter and performing a 2$D$ dimensional
optimization, the driving and driven signals are learned in tandem. This
reduces the dimensionality of the problem, partitions the variables in a
natural manner, and is shown empirically to converge to optima that recover
the causal structure underlying the data.

As the cost function is non-convex, there are potentially several pairs of
projection vectors that locally maximize the strength of causality (34) and
thus may yield meaningful latent sources. In order to recover $P$ pairs of
components $\\{y_{i}(t),z_{i}(t)\\},i=1,\ldots,P$, here it is proposed to
repeat the optimization after the first iterate, but not before removing the
contribution of the driving signal $y_{1}(t)$ from the observed data. This
takes the form of a spatiotemporal regression such that any signals that are
correlated with $y_{1}(t)$ or its lagged versions $y_{1}(t-l),l=1,\ldots,L$
are removed. Given that this should include $z_{1}(t)$, the driven signal is
not explicitly removed. This procedure is repeated until the desired number of
component pairs $P$ is obtained. The proposed algorithm is described in
_Supplementary Note 3_.

In what follows, the proposed approach is evaluated on synthetic and real-
world data. The primary criterion employed to assess performance is the
strength of causality (9) among the recovered pairs of components relative to
the strength of causality between observed signals, or those formed by
conventional component analysis techniques. Where possible, the fidelity of
the recovered signals compared to the underlying sources is measured.
Moreover, the recovered components and associated projection vectors are
interpreted based on what is known about the system being investigated (i.e.,
neural dynamics, the cryptocurrency market) to further assess the behavior of
the proposed approach.

### Recovering the causal structure of a three-element system

To test the proposed method’s ability to recover the causal structure embedded
in multiple time series, a series of empirical evaluations was conducted on
synthetic data. Access to the system’s ground-truth structure permitted
measuring the fidelity of the recovered signals with respect to the latent
sources. The data was generated according to a VAR(3) process whose parameters
matched those employed by Stokes and Purdon [34], where $s_{1}\rightarrow
s_{2}$ and $s_{2}\rightarrow s_{3}$. Projection of the three connected sources
to a $D=4$ dimensional observation vector followed as
$\mathbf{x}(t)=\mathbf{A}\mathbf{s}(t)$, where the elements of 4-by-3 mixing
matrix $\mathbf{A}$ were randomly drawn from the uniform distribution
$A_{ij}\sim U[0,1]$. The proposed technique was employed to recover $P=2$
pairs of components.

The latent sources, observed data, and recovered signals of a single
realization are depicted in Fig 2a,b, and c, respectively. The goal of the
proposed approach is to recover the $s_{1}\rightarrow s_{2}$ relationship in
($y_{1},z_{1}$), and the $s_{2}\rightarrow s_{3}$ link in ($y_{2},z_{2}$).
Notice that $s_{2}$ is both a driven signal as well as a driving signal, and
thus the components $z_{1}$ and $y_{2}$ are aiming to capture the _same_
signal. The strength of causality among all pairs of latent sources is
illustrated in Fig 2d, where the order dependence inherent to Granger
Causality is evident in the asymmetry of the matrix (rows correspond to the
driving signals, columns to the driven). The underlying strengths of causality
were measured as: $\mathcal{G}_{s_{1}\rightarrow s_{2}}=0.11\pm 0.0011$ and
$\mathcal{G}_{s_{2}\rightarrow s_{3}}=0.10\pm 0.0010$ (mean $\pm$ sem across
$n=100$ random realizations).

The strength of causality measured among pairs of observed signals was
markedly lower, with a maximum (across all pairs) strength of causality of
$0.066\pm 0.002$, significantly lower than the underlying latent causality
($p=3.9\times 10^{-18}$ comparing to $\mathcal{G}_{s_{1}\rightarrow s_{2}}$,
$p=4.0\times 10^{-18}$ comparing to $\mathcal{G}_{s_{2}\rightarrow s_{3}}$,
Wilcoxon signed rank test, $n=100$). The strengths of causality between
observed signals is depicted for a single realization in Fig 2e, where the
strongest connection was $\mathcal{G}_{x_{1}\rightarrow x_{3}}=0.025$. Notice
that the system structure (i.e., two connected pairs) is no longer apparent,
as the mixing process has both obscured and dampened the underlying causal
relationships.

The causality matrix of the $P=2$ recovered components is shown in Fig 2f
(rows and columns correspond to driving and driven components, respectively).
Two strong connections are readily apparent: $y_{1}\rightarrow z_{1}$ and
$y_{2}\rightarrow z_{2}$. The magnitudes of these causal relationships closely
matched those of the latent sources: $\mathcal{G}_{y_{1}\rightarrow
z_{1}}=0.11\pm 0.0014$ and $\mathcal{G}_{y_{2}\rightarrow z_{2}}=0.10\pm
0.0010$. These values were significantly larger than the maximum causality
among all pairs of observed variables (pair 1: $p=4.7\times 10^{-18}$, pair 2:
$p=4.0\times 10^{-18}$). To determine whether the recovered components
captured the underlying sources, the mixing matrix was estimated by regressing
the driving and driven signals onto the observations. The true mixing matrix
is depicted for a single realization in Fig 2g. The recovered matrices
exhibited a large correlation with the true values ($r^{2}=0.98\pm 0.004$,
shown for a single realization in Fig 2h). Moreover, the time series of
recovered components faithfully tracked the dynamics of the latent sources:
$r^{2}=0.98\pm 0.007$ between $s_{1}$ and $y_{1}$ (Fig 2i), $r^{2}=0.96\pm
0.014$ between $s_{2}$ and $z_{1}$ (Fig 2j), $r^{2}=0.98\pm 0.004$ between
$s_{2}$ and $y_{2}$ (Fig 2k), and $r^{2}=0.99\pm 0.002$ between $s_{3}$ and
$z_{2}$ (Fig 2l). Note that latent source $s_{2}$ was indeed captured by both
$z_{1}$ and $y_{2}$.

### Identifying latent causal connections in the brain

Scalp electroencephalogram (EEG) signals, sometimes referred to as “brain
waves”, arise from the coordinated activity of a large number of neurons in
the cerebral cortex [35]. At any time instant, the set of electric potentials
registered by scalp electrodes is a linear mixture of dipolar generators [36]
(Fig 3a). It was hypothesized that Granger causal relations would be most
strongly observed at the level of these neural generators, estimated by
spatially filtering the EEG [37]. To test this, the proposed technique was
applied to a previously collected EEG data set where $n=12$ healthy
participants viewed television advertisements that originally aired during the
2012 and 2013 _Super Bowl_ football matches [38].

To identify the level of Granger causality among the captured $D=64$ signals,
the strength of causality was measured for all pairs of electrodes (Fig 3b).
The strongest relationship was found between left centroparietal electrode
“CP1” and right centroparietal electrode “CP4”, with $\mathcal{G}=0.073$ (Fig
3b). In order to determine whether conventional spatial filtering approaches
recover stronger causal relationships than those found among electrodes, the
observed data was decomposed with both principal components analysis (PCA) and
independent components analysis (ICA). Surprisingly, the strength of causality
among pairs of PCs and ICs was not larger than that found in the raw
electrodes: a maximum value of $\mathcal{G}=0.067$ was found between principal
components 8 and 3 (Fig 3c), and a maximum of $\mathcal{G}=0.022$ between
independent components 6 and 3 (Fig 3d). Next, the proposed method was
employed to recover $P=3$ pairs of latent components. The strength of
causality among the recovered pairs was substantially larger, with
$\mathcal{G}_{y_{1}\rightarrow z_{1}}=0.32$, $\mathcal{G}_{y_{2}\rightarrow
z_{2}}=0.16$, and $\mathcal{G}_{y_{3}\rightarrow z_{3}}=0.18$, for pairs 1, 2,
and 3, respectively (Fig 3e). The presence of more than two-fold increases in
the strength of causality at multiple component pairs is consistent with the
notion that the underlying causal relationships occur in a latent subspace of
the data.

The coefficients of the spatial filter weights learned by the proposed method
represent the scalp regions expressing the driving and driven signals. For
pair 1, the causing signal $y_{1}$ exhibited peak expression over the right
temporo-parietal region, while the driven signal $z_{1}$ had peak expression
over the left central electrodes (Fig 3f). This indicates that, during this
task, activity over the right temporo-parietal cortex temporally preceded
activity over the left central region. To further interpret the learned
components, the power spectrum of the driving and driven signals were
measured. The power spectrum of scalp EEG is typically segregated into
distinct frequency bands, with a large body of literature documenting
associations between cognitive states and activity in specific bands [39].
Both $y_{1}$ and $z_{1}$ showed high levels of power in the delta band (1-3
Hz), and moderate levels of alpha band (8-13 Hz) power (Fig 3f). The spatial
topographies of the next strongest pair showed peak expression over the left
parieto-occipital ($y_{2}$) and right temporo-parietal regions ($z_{2}$),
indicating inter-hemispheric connectivity (Fig 3g). An interesting pattern
arose in the power spectra of the components: the driving signal was marked by
low delta power and high alpha power, while the driven signal exhibited the
opposite pattern (i.e., high delta power and a notable absence of alpha
power). This result is consistent with previous findings of an inverse
correlation between alpha and delta waves, hypothesized to arise from
thalamocortical inhibition of the brain stem [40, 41]. The topography of
driving signal $y_{3}$ exhibited activation over the left occipital and right
centro-temporal regions, while the corresponding driven signal $z_{3}$ was
concentrated over the left occipital region (Fig 3h). As observed in pair 2,
the driving signal showed a high ratio of alpha-to-delta power, while a low
alpha-to-delta ratio was detected in the driven signal.

To formally test whether the proposed method recovers stronger causal
relations than those found with conventional approaches, a two-way ANOVA
(method $\times$ component) was conducted. For this analysis, the strength of
causality was measured separately for each subject, yielding $n=12$ repeated
measures. A large main effect of method was identified ($F(3)=11.53$,
$p=9.23\times 10^{-7}$; Fig S2). There was no main effect of component
($p=0.59$) and no significant interaction ($p=0.98$). Follow-up tests showed
that the main effect of method was driven by significantly larger strengths of
causality with the proposed method ($\mathcal{G}=0.10\pm 0.031$, $0.088\pm
0.012$, and $0.081\pm 0.014$ for the first three components, means $\pm$ sem
across $n=12$ subjects) relative to the three most connected electrode pairs
($\mathcal{G}=0.049\pm 0.0084$, $0.049\pm 0.0072$, $0.049\pm 0.0079$;
$p=0.034$; $p=4.9\times 10^{-4}$, and $p=0.034$ for components 1, 2, and 3,
respectively; Wilcoxon signed rank test, $n=12$), the three most connected
principal component pairs ($\mathcal{G}=0.050\pm 0.0054$, $0.049\pm 0.0058$,
$0.046\pm 0.0046$; $p=0.034$, $p=0.0049$, and $p=0.0093$), and the three most
connected independent component pairs ($\mathcal{G}=0.033\pm 0.0056$,
$0.030\pm 0.0035$, $0.030\pm 0.0031$; $p=0.016$, $p=4.9\times 10^{-4}$, and
$p=0.0024$). Thus, the proposed technique detected causal relationships whose
magnitude was significantly larger than those measured with conventional
approaches.

### Probing latent causality in the cryptocurrency market

Finally, the proposed method was tested on a system without an obvious latent
structure: the cryptocurrency market. Historical prices of $D=19$ popular
cryptocurrencies (Fig 4A, individual traces have been standardized), were
employed for the analysis, which sought to identify the $P=3$ strongest causal
relationships.

Among pairs of individual cryptocurrencies, the strength of causality was
quite modest: $0.028\pm 0.023$ (mean $\pm$ sd across all $n=342$ pairs of
currencies), with a maximum value of
$\mathcal{G}_{\mathrm{ETC}\rightarrow\mathrm{QTUM}}=0.12$ (Fig 4B). In
contrast, the proposed technique identified a primary pair of latent
components with a statistically significant strength of causality
($\mathcal{G}_{y_{1}\rightarrow z_{1}}=0.40$, $p<0.001$, non-parametric
permutation test altering the phase of individual cryptocurrency time series),
representing a more than three-fold increase (Fig 4C). A statistically
significant strength of causality was also found for the second pair of
components ($\mathcal{G}_{y_{2}\rightarrow z_{2}}=0.14$, $p=0.008$; Fig 4C).
Note that, even after removing the contribution from the primary driving
signal $y_{1}$, a latent relationship whose causality exceeded that seen in
the observed data was still recovered. The strength of causality exhibited by
the third pair of latent components ($\mathcal{G}_{y_{3}\rightarrow
z_{3}}=0.080$, $p=0.13$) fell short of significance, but nevertheless exceeded
96% of the individual pair values (compare panels B and C in Fig 4).

The dynamics of the driving and driven components of the first pair are
depicted in Fig 4D, where the temporal precedence of $y_{1}$ relative to
$z_{1}$ is visible in the traces. For example, note that the occurrence of the
three prominent peaks in the spring of 2021 is first observed in $y_{1}$ and
shortly after in $z_{1}$ (see Fig 4D inset). The individual currencies with
the largest expression in the driving signal were BNB (Binance Coin) and ETC
(Ethereum Classic), while the largest contributions to the driven signal were
from QTUM and TRX (Fig 4E, color indicates weight of filter used to construct
$y_{1}$ and $z_{1}$). This result indicates that past fluctuations in the
prices of BNB and ETC predict the current prices of QTUM and TRX. The temporal
precedence of $y_{2}$ relative to $z_{2}$ is also evident in the dynamics of
the second pair of latent components (Fig 4F). For example, a sharp dip in
price occurs near May 2021, first in $y_{2}$ and slightly later in $z_{2}$.
Similar to the first pair of latent components, the currencies best expressed
in $y_{2}$ were ETC and BNB. However, unlike ($y_{1},z_{1}$), the driven
signal here most strongly expressed ADA (Cardano) and ETH (Ethereum) (Fig 4G).
The finding of similar driving signals (but distinct driven signals) in the
first two pairs suggests the presence of multiple “links” emanating from the
latent driver. The currencies best expressed in the driving signal of the
third pair were BNB and XRP (Ripple), while the corresponding driven signal
$z_{3}$ best expressed XLM (Stellar) and ETC (Fig 4I).

## Discussion

The distinction between the proposed technique and conventional univariate and
multivariate Granger Causality can be illuminated by the types of queries that
the different approaches address. In the context of the cryptocurrency market,
univariate Granger Causality addresses questions such as “does the price of
Bitcoin exert a causal influence on the price of Ethereum?” Multivariate
Granger Causality is concerned with questions such as “do the prices of
Bitcoin and Cardano (taken as a group) drive the prices of Ethereum and
Ethereum Classic?” Note that, in both cases, one must specify the elements and
direction of the causal relationship being tested. To identify the full
complement of causal links in the system of interest, such a hypothesis
testing approach will generally require a large number of statistical tests.
In contrast, the proposed method automatically identifies paired groups of
cryptocurrencies, with each group defined such that the strength of causality
from the driving group to the driven group is maximized: the elements and
direction of the causal links are learned directly from the data. This
identification may be performed over several iterations, with each iteration
revealing a generally weaker but distinct causal relationship from the
previous. The weights of the learned filters are interpretable: dimensions
with a large magnitude indicate that the corresponding signal is either
driving activity, or being driven, in a latent subspace of the system.

In applications such as EEG or magnetoencephalography (MEG) where the source
space has a clear physical substrate, the learned filters offer clear insight
into the nature of the latent sources. Namely, the cortical generators of the
scalp topographies in Fig 3f-h may be estimated with source localization [36]
to estimate the spatial origin of the latent sources. Causal relationships
that are obscured at the level of the electrodes may be clarified as genuine
connections between cortical sources. The nature of the latent source space is
less apparent in other problems. In financial systems defined by a set of
evolving prices, the latent sources correspond to a set of linked assets whose
dynamics exhibit a temporal dependence on those of a second set. For example,
the occurrence of an external event (e.g. activity on social media) may
produce a change in the value of a certain group of assets. As a consequence,
the value of a second (disparate) group of assets may also be modulated, and
due to the delay between the price movements, a Granger causal relationship
emerges.

Conventional approaches to blind source separation assume that the underlying
sources are statistically independent, perhaps inspired by the “cocktail party
problem” [42] solved by the auditory system. This assumption is exploited by
Independent Components Analysis (ICA) [27, 28], which projects the observed
signals into components to maximize their statistical independence. On the
other hand, the approach proposed here assumes the existence of Granger Causal
sources, and is thus applicable to systems with temporal dependencies among
the signals of interest. Notice that the criteria optimized by ICA and the
proposed technique to perform source separation are opposing. In the context
of brain signals, ICA is seeking to identify decoupled neural sources, while
the method proposed here aims to recover functionally connected brain regions.
More closely related to the proposed method are approaches that combine
Canonical Correlation Analysis [43] with Granger Causality [44, 45] to test
causal relations between pairs of multivariate time series. These approaches
share a feature of the proposed method by forming components of observed data,
but differ importantly in that the data must already be partitioned into
hypothesized driving and driven signals.

One limitation of the proposed technique is the potential difficulty in
identifying causality in data with very high dimensionality (i.e., the number
of observed signals) or very long temporal dependencies between latent
sources. In either case, the covariance matrices required to identify the
latent causal sources may be poorly estimated, potentially leading to
erroneous estimates of latent Granger Causality. To mitigate this, it is
required to assume some prior information about the structure of the observed
signals. For example, a form of Tikhonov regularization [46] equivalent to
adding uncorrelated noise to the measurements was employed here. More
sophisticated approaches to covariance estimation in high dimensions will
improve the performance of the proposed framework.

A challenge with conventional Granger Causality is the potential presence of
exogenous sources that drive two or more observed variables with different
delays. In this event, spurious relationships between the observed signals may
be inferred. To address this, partial Granger Causality [47] may be employed
to measure the relationship that remains after removing the contribution of
the exogenous source. It is interesting to consider how such confounding
sources may affect the behavior of the proposed technique. If the nature of
the confounding source is known _a priori_ , it should be regressed out of the
data prior to deploying the proposed technique. This was performed in the
cryptocurrency example above, where the global market trend was removed prior
to analysis. In the case of an unknown confounding source, the proposed
approach is expected to provide some shielding from spurious inference. This
follows from the utilization of multiple component pairs to separate the
contributions of distinct latent sources. For example, in the case of a strong
confounding source that enters the observed data, the underlying relationship
may appear in the first pair of latent sources, leaving the genuine causal
relationships in subsequent pairs. The technique proposed here is tasked with
capturing all latent sources that produce Granger Causal links, meaningful or
otherwise. This highlights the importance of interpreting the weights of the
learned projection vectors, which may offer clues as to the origin of the
recovered relationship.

Granger Causality is one of several statistical approaches to measuring
causality. Two popular frameworks that have been successfully applied to
dynamic systems are Dynamic Causal Modeling (DCM) [48] and Structural Equation
Modeling [49]. In DCM, a “forward model” that relates the activity of
underlying sources to the observations is specified, with Bayesian model
selection utilized to estimate the parameters of the underlying sources (i.e.,
connectivity). This allows DCM to take advantage of the known structure of the
system, including nonlinear interactions. The approach proposed here, while
also aiming to identify causal structure, is complementary in nature. The
forward model need not be specified beforehand, and the technique functions
not as a statistical test _per se_ but rather a decomposition of the data,
akin to PCA and ICA. Moreover, the knowledge gleaned from the components
recovered by the decomposition may then be employed in a subsequent hypothesis
testing procedure that has been informed by the method’s findings.

## Materials and Methods

All data and source code are provided at dmochow.github.io/gca. Data analysis
was performed in the MATLAB computing environment (Mathworks, Natick MA).

##### Implementation

To solve the optimization problems at each iteration of the grouped coordinate
descent algorithm (see Algorithm 1 in Supplementary Note 3), we employed the
built-in MATLAB function fmincon with the default interior point algorithm
solver. The maximum number of function evaluations was set to $10^{4}$ and the
maximum number of iterations was set to $4000$. Regularization of the block
covariance matrices $\mathbf{\Sigma}_{1:L}$ and $\tilde{\mathbf{\Sigma}}$ was
implemented by limiting the condition number of each matrix to a value of $c$,
where the value of $c$ was selected based on the dimensionality of the
problem, as specified below. Limiting the condition number was implemented by
adding a small diagonal component $\sigma^{2}\mathbf{I}$ to each covariance
matrix, where the value of
$\sigma^{2}=\frac{(\lambda_{1}-\lambda_{LD}c)}{c-1}$ ensures that the
condition number of the covariance matrix is $c$, where $\lambda_{1}$ and
$\lambda_{LD}$ are the largest and smallest eigenvalues of the block
covariance matrix being regularized [50, 51].

Although the closed-form expression for the gradient of
$\mathcal{G}_{y\rightarrow z}$ (see _Supplementary Note 2_) was verified
empirically, it was more efficient to compute the gradient numerically with
finite differences. The numerous Kronecker products and matrix inverse
operations required to evaluate the gradient expression led to longer run
times compared to the finite differences approximation. Moreover, in order to
guarantee that the optimization identified projections with unit norm, a pair
of nonlinear constraints were added, leading to the following constrained
optimization problem:

$\displaystyle\min_{\mathbf{w},\mathbf{v}}~{}-\left[\mathcal{G}(\mathbf{w},\mathbf{v})+\mathcal{G}^{\mathrm{tr}}(\mathbf{v},\mathbf{w})\right]$
$\displaystyle\mathrm{~{}~{}~{}subject~{}to:~{}}\mathbf{w}^{T}\mathbf{w}=1\mathrm{~{}and~{}}\mathbf{v}^{T}\mathbf{v}=1,$
(35)

where $\mathcal{G}(\mathbf{w},\mathbf{v})$ is the strength of causality (9)
between driving signal $\mathbf{w}^{T}\mathbf{x}(t)$ and driven signal
$\mathbf{v}^{T}\mathbf{x}(t)$, and
$\mathcal{G}^{\mathrm{tr}}(\mathbf{v},\mathbf{w})$ is the strength of
causality between driving signal $\mathbf{v}^{T}\mathbf{x}(-t)$ and driven
signal $\mathbf{w}^{T}\mathbf{x}(-t)$. After each iteration of the grouped
coodinate descent, the driving signal
$y(t)={\mathbf{w}^{\ast}}^{T}\mathbf{x}(t)$ and its lagged versions were
regressed out of the data according to:

$\displaystyle\mathbf{x}(t)$ $\displaystyle=$
$\displaystyle\mathbf{x}(t)-\mathbf{B}^{T}\mathbf{y}_{p}(t)$ (36)

where $\mathbf{B}=\mathbf{Y}_{p}^{\\#}\mathbf{X}$ is the least-squares
solution to the linear system:

$\displaystyle\mathbf{X}=\mathbf{Y}_{p}\mathbf{B}$ (37)

where $D$-by-$T$ matrix
$\mathbf{X}=\left[\begin{array}[]{ccc}\mathbf{x}(1)&\ldots&\mathbf{x}(T)\end{array}\right]$
and $L$-by-$T$ matrix
$\mathbf{Y}_{p}=\left[\begin{array}[]{ccc}\mathbf{y}_{p}(1)&\ldots&\mathbf{y}_{p}(T)\end{array}\right]$
span the spatiotemporal apertures of the observed and driving signals,
respectively. Convergence was assessed after every iteration, and the search
was stopped when the magnitude of change in both $\mathcal{G}$ and
$\mathcal{G}^{\mathrm{tr}}$ was less than $10^{-6}$.

To measure the strength of causality $\mathcal{G}_{f\rightarrow g}$ between
signals $f$ and $g$, the full and reduced regression models predicting $g(t)$
were explicitly learned, and the residuals then used to obtain
$\mathcal{G}_{f\rightarrow g}$ via Eqn. (9).

##### Synthetic VAR data and analysis

Data was generated by explicitly defining the VAR(3) system analyzed
previously by Stokes and Purdon [34]:

$\displaystyle\left[\begin{array}[]{c}s_{1}(t)\\\ s_{2}(t)\\\
s_{3}(t)\end{array}\right]$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{ccc}-0.9&0&0\\\ -0.356&1.212&0\\\
0&-0.3098&-1.3856\end{array}\right]\left[\begin{array}[]{c}s_{1}(t-1)\\\
s_{2}(t-1)\\\ s_{3}(t-1)\end{array}\right]+$ (64)
$\displaystyle\left[\begin{array}[]{ccc}-0.81&0&0\\\ 0.7136&-0.49&0\\\
0&0.50&-0.64\end{array}\right]\left[\begin{array}[]{c}s_{1}(t-2)\\\
s_{2}(t-2)\\\ s_{3}(t-2)\end{array}\right]+$
$\displaystyle\left[\begin{array}[]{ccc}0&0&0\\\ -0.356&0&0\\\
0&-0.3098&0\end{array}\right]\left[\begin{array}[]{c}s_{1}(t-3)\\\
s_{2}(t-3)\\\
s_{3}(t-3)\end{array}\right]+\left[\begin{array}[]{c}\epsilon_{1}(t)\\\
\epsilon_{2}(t)\\\ \epsilon_{3}(t)\end{array}\right],$

where $\epsilon_{i}$, $i=1,2,3$, are independent and identically distributed
innovation processes with standard deviation $\sigma=1$. $M=100$ realizations,
each with a length of $N=5000$ samples, were generated by passing the vector
innovation process through the impulse response (64). Projection of these
latent sources to a four-dimensional observation vector followed as
$\mathbf{x}(t)=\mathbf{A}\mathbf{s}(t)$, where the elements of 4-by-3 mixing
matrix $\mathbf{A}$ were randomly drawn from the uniform distribution
$A_{ij}\sim U[0,1]$. Notice that measurement noise enters the observed data
via the innovation processes $\epsilon_{i}$. The proposed technique was
employed to recover $P=2$ pairs of causal components:

$\displaystyle y_{i}(t)=\mathbf{w}_{i}^{T}\mathbf{x}(t),~{}~{}~{}~{}i=1,2$
$\displaystyle z_{i}(t)=\mathbf{v}_{i}^{T}\mathbf{x}(t),~{}~{}~{}~{}i=1,2,$

where $\mathbf{w}_{i}$ and $\mathbf{v}_{i}$ were estimated with Algorithm 1.
Convergence was observed in under 20 iterations for pair 1, and under 10 for
the second pair (Figure S1).

The optimization was performed with no regularization of the block covariance
matrices ($c=\infty$) and a maximum lag parameter of $L=3$. $P=2$ pairs were
recovered by the optimization. When comparing the fidelity of the recovered
component pairs with the ground-truth latent sources, the order of the $P=2$
pairs was corrected _post hoc_ if it was evident that the $(y_{1},z_{1})$ pair
matched the $s_{2}\rightarrow s_{3}$ relationship. In practice, the order of
the recovered pairs ($s_{1}\rightarrow s_{2}$, $s_{2}\rightarrow s_{3}$) is
insignificant, as the causal structure reflected by the two pairs is agnostic
to their ordering.

To estimate the mixing matrix from the model’s projection vectors $\mathbf{w}$
and $\mathbf{v}$, the driving signals $y_{1}$ and $y_{2}$, as well as the
driven signal $z_{2}$, were individually regressed onto the observation vector
$\mathbf{x}$. This yielded a $D$-dimensional ‘‘forward model’’ for each of the
three signals, which were then compared to the three columns of the true
mixing matrix. When displaying the estimated and true mixing matrix in Fig
2g,h, the sign and scale (L2 norm) of each estimated forward model was
corrected to match that of the ground-truth mixing matrix column.

When testing for significant differences in the strength of causality between
observed signals and those recovered by the proposed method, the Wilcoxon
signed rank test ($n=100$ independent VAR realizations) was employed. The
maximum value across all pairs of observed signals (i.e.,
$\max_{i,j}\mathcal{G}_{x_{i}\rightarrow x_{j}}$) was compared against the
strength of causality of the first two recovered pairs (i.e.,
$\mathcal{G}_{y_{1}\rightarrow z_{1}}$, $\mathcal{G}_{y_{2}\rightarrow
z_{2}}$). The same procedure was employed to test for significant differences
in the strength of causality between observed and ground-truth latent sources
(i.e., $\mathcal{G}_{s_{1}\rightarrow s_{2}}$, $\mathcal{G}_{s_{2}\rightarrow
s_{3}}$)

##### EEG data and analysis

The neural data employed here to demonstrate the utility of proposed method
has been previously described [38]. Briefly, scalp EEG was collected from
$n=12$ subjects freely viewing a set of 30-60 second advertisements originally
broadcast during the 2012 and 2013 SuperBowl. To demonstrate the utility of
the proposed method, data from a single stimulus was employed here (‘‘Work’’,
Bud Light Platinum). The data was acquired with a 64-channel electrode cap
connected to a BioSemi Active Two amplifier and sampled at rate of 512 Hz. A
set of preprocessing steps comprised of high-pass and notch filtering, removal
of eye motion artifacts by linear regression, and artifact rejection with a
power criterion was applied to denoise the acquired signals. All data samples
identified as artifactual by the preprocessing were linearly interpolated from
neighboring samples. The interpolation allowed the computation of block
covariance matrices in the presence of missing data. Moreover, data was
further downsampled to a sampling frequency of 32 Hz in order to reduce the
dimensionality of the ensuing block covariance matrices. The maximum lag
parameter $L$ was set to 16 samples (500 ms), reflecting a tradeoff between
capturing dependencies occurring on the temporal scale of neural dynamics,
while avoiding excessively large covariance matrices. The number of desired
component pairs was set to $P=3$.

EEG signals were mean centered prior to testing the proposed method. The block
covariance matrices $\mathbf{\Sigma}_{1:L}$ and $\tilde{\mathbf{\Sigma}}$ were
regularized such that the condition number of each matrix was limited to
$K=10^{9}$. The maximum number of iterations in the grouped coordinate descent
was set to 50.

To depict the spatial topographies of the latent components, the ‘‘forward-
model’’ [29] conveying the distribution of the latent source on the scalp
$\mathbf{a}_{w}=\mathbf{\Sigma}(0)\mathbf{w}\left(\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{w}\right)^{-1}$
was computed, where $\mathbf{\Sigma}(0)$ is the lag-zero covariance matrix of
the observations $\mathbf{x}(t)$. Power spectra were estimated with the
Thomson multitaper spectral analysis technique employing a time-bandwidth
product of 64. When comparing the proposed technique with principal components
analysis, the strength of causality was measured between all 90 pairs of the
first 10 principal components (the approximate knee point of the data’s
eigenvalue spectrum). Similarly, the strength of causality was calculated
among all pairs of the 10 independent components formed after performing PCA
on the data. The maximum-kurtosis implementation of ICA was employed [52].

To perform two-way ANOVA with method and component as factors, the spatial
filters learned on the subject-aggregated data were applied to the recordings
of individual subjects, yielding $n=12$ independent measures of the strength
of causality obtained with the proposed method. The three electrode pairs with
the largest (subject-aggregated) strength of causality were selected _post
hoc_. Similarly, the three principal and independent component pairs with the
largest strength of causality were selected. The strength of causality values
at the selected pairs were then measured for all subjects and employed in the
ANOVA procedure. Note that the values of strength of causality yielded by the
proposed method were markedly larger (i.e., $\mathcal{G}=0.32$) when evaluated
on the entire (subject aggregated) data set relative to the values obtained
when applying the spatial filters to individual subjects and averaging across
the cohort (i.e., $\mathcal{G}=0.10\pm 0.0031$).

##### Cryptocurrency data and analysis

Publicly available data was obtained from an online database of historical
cryptocurrency prices as captured on the Binance Exchange
(CryptoDataDownload.com). Data was obtained from the following $D=19$
currencies: ADA, BAT, BNB, BTC, BTT, DASH, EOS, ETC, ETH, LINK, LTC, NEO,
QTUM, TRX, USDC, XLM, XMR, XRP, and ZEC. Prices were obtained at the
resolution of one minute, but subsequently downsampled by a factor of 1800 in
order to capture slower dynamics manifesting across half-day segments. The
opening price in each segment (i.e., as opposed to the high, low, or closing
price) was employed for the analysis.

Due to the fact that the proposed method cannot recover the scale of the
latent sources, each currency’s time series was standardized by removing the
mean and dividing by the standard deviation. Furthermore, in order to capture
genuine causal relationships unaffected by exogenous factors not captured in
these currencies, the mean waveform (‘‘global’’ trend) was linearly regressed
out from the multivariate time series with ordinary least squares.

The proposed algorithm was employed with a maximum lag of $L=4$ (i.e., a two-
day temporal aperture), and the $P=3$ strongest pairs of latent components
were computed. Regularization of the block covariance matrices
$\tilde{\mathbf{\Sigma}}$ and $\mathbf{\Sigma}_{1:L}$ was performed by
limiting the condition number of both matrices to $K=1000$. To test for
statistically significant strengths of causality in the recovered component
pairs, a non-parametric test that employs surrogate data generated by
randomizing the phase spectrum of the original data (while preserving its
power spectrum) was employed [53]. This procedure effectively ‘‘shuffles’’ the
time series of the various cryptocurrency prices such that the genuine
temporal dependencies are removed. The strength of causality measured from the
surrogate records then provides a sample of the null distribution to which the
true values were compared. A total of 1000 surrogate data records were formed,
with the p-value measured as the number of records whose strength of causality
exceeded the true value.

To interpret the constituents of the latent souces learned by the proposed
method, the elements of $\mathbf{w}$ and $\mathbf{v}$ were sorted by
magnitude, and the two elements with the largest absolute value were reported
in the text.

## Acknowledgments

The author would like to thank Amilcar Malave for help with figure
preparation. This research was supported by the Weinbaum - Wallace H. Coulter
Fund.

## References

  * [1] Olaf Sporns.  Networks of the Brain.  MIT press, 2010.
  * [2] Ruey S Tsay.  Analysis of financial time series.  John wiley & sons, 2005.
  * [3] Clive WJ Granger.  Investigating causal relations by econometric models and cross-spectral methods.  Econometrica: journal of the Econometric Society, pages 424--438, 1969.
  * [4] Clive WJ Granger.  Essays in econometrics: collected papers of Clive WJ Granger, volume 32.  Cambridge University Press, 2001.
  * [5] Mingzhou Ding, Yonghong Chen, and Steven L Bressler.  17 granger causality: basic theory and application to neuroscience.  Handbook of time series analysis: recent theoretical developments and applications, 437, 2006.
  * [6] Anil K Seth, Adam B Barrett, and Lionel Barnett.  Granger causality analysis in neuroscience and neuroimaging.  Journal of Neuroscience, 35(8):3293--3297, 2015.
  * [7] Corrado Bernasconi and Peter KoÈnig.  On the directionality of cortical interactions studied by structural analysis of electrophysiological recordings.  Biological cybernetics, 81(3):199--210, 1999.
  * [8] Maciej Kamiński, Mingzhou Ding, Wilson A Truccolo, and Steven L Bressler.  Evaluating causal relations in neural systems: Granger causality, directed transfer function and statistical assessment of significance.  Biological cybernetics, 85(2):145--157, 2001.
  * [9] Rainer Goebel, Alard Roebroeck, Dae-Shik Kim, and Elia Formisano.  Investigating directed cortical interactions in time-resolved fmri data using vector autoregressive modeling and granger causality mapping.  Magnetic resonance imaging, 21(10):1251--1261, 2003.
  * [10] Alireza Sheikhattar, Sina Miran, Ji Liu, Jonathan B Fritz, Shihab A Shamma, Patrick O Kanold, and Behtash Babadi.  Extracting neuronal functional network dynamics via adaptive granger causality analysis.  Proceedings of the National Academy of Sciences, 115(17):E3869--E3878, 2018.
  * [11] Raul Vicente, Michael Wibral, Michael Lindner, and Gordon Pipa.  Transfer entropy—a model-free measure of effective connectivity for the neurosciences.  Journal of computational neuroscience, 30(1):45--67, 2011.
  * [12] George Sugihara, Robert May, Hao Ye, Chih-hao Hsieh, Ethan Deyle, Michael Fogarty, and Stephan Munch.  Detecting causality in complex ecosystems.  science, 338(6106):496--500, 2012.
  * [13] Justin D Finkle, Jia J Wu, and Neda Bagheri.  Windowed granger causal inference strategy improves discovery of gene regulatory networks.  Proceedings of the National Academy of Sciences, 115(9):2252--2257, 2018.
  * [14] Michael Eichler and Vanessa Didelez.  On granger causality and the effect of interventions in time series.  Lifetime data analysis, 16(1):3--32, 2010.
  * [15] Samantha Kleinberg and George Hripcsak.  A review of causal inference for biomedical informatics.  Journal of biomedical informatics, 44(6):1102--1112, 2011.
  * [16] John Geweke.  Measurement of linear dependence and feedback between multiple time series.  Journal of the American statistical association, 77(378):304--313, 1982.
  * [17] John F Geweke.  Measures of conditional linear dependence and feedback between time series.  Journal of the American Statistical Association, 79(388):907--915, 1984.
  * [18] Adam B Barrett, Lionel Barnett, and Anil K Seth.  Multivariate granger causality and generalized variance.  Physical Review E, 81(4):041907, 2010.
  * [19] Lionel Barnett and Anil K Seth.  The mvgc multivariate granger causality toolbox: a new approach to granger-causal inference.  Journal of neuroscience methods, 223:50--68, 2014.
  * [20] Craig Hiemstra and Jonathan D Jones.  Testing for linear and nonlinear granger causality in the stock price-volume relation.  The Journal of Finance, 49(5):1639--1664, 1994.
  * [21] Nicola Ancona, Daniele Marinazzo, and Sebastiano Stramaglia.  Radial basis function approach to nonlinear granger causality of time series.  Physical Review E, 70(5):056221, 2004.
  * [22] Daniele Marinazzo, Mario Pellicoro, and Sebastiano Stramaglia.  Kernel method for nonlinear granger causality.  Physical review letters, 100(14):144103, 2008.
  * [23] Alex Tank, Ian Covert, Nicholas Foti, Ali Shojaie, and Emily Fox.  Neural granger causality.  arXiv preprint arXiv:1802.05842, 2018.
  * [24] Mariusz Maziarz.  A review of the granger-causality fallacy.  The journal of philosophical economics: Reflections on economic and social issues, 8(2):86--105, 2015.
  * [25] Greta Grassmann.  New considerations on the validity of the wiener-granger causality test.  Heliyon, 6(10):e05208, 2020.
  * [26] Norbert Wiener et al.  Extrapolation, interpolation, and smoothing of stationary time series: with engineering applications, volume 8.  MIT press Cambridge, MA, 1964.
  * [27] Pierre Comon.  Independent component analysis, a new concept?  Signal processing, 36(3):287--314, 1994.
  * [28] Aapo Hyvärinen and Erkki Oja.  Independent component analysis: algorithms and applications.  Neural networks, 13(4-5):411--430, 2000.
  * [29] Stefan Haufe, Vadim V Nikulin, Klaus-Robert Müller, and Guido Nolte.  A critical assessment of connectivity measures for eeg data: a simulation study.  Neuroimage, 64:120--133, 2013.
  * [30] Irene Winkler, Danny Panknin, Daniel Bartz, Klaus-Robert Müller, and Stefan Haufe.  Validity of time reversal for testing granger causality.  IEEE Transactions on Signal Processing, 64(11):2746--2760, 2016\.
  * [31] Zsolt Ugray, Leon Lasdon, John Plummer, Fred Glover, James Kelly, and Rafael Martí.  Scatter search and local nlp solvers: A multistart framework for global optimization.  INFORMS Journal on computing, 19(3):328--340, 2007.
  * [32] Léon Bottou, Frank E Curtis, and Jorge Nocedal.  Optimization methods for large-scale machine learning.  Siam Review, 60(2):223--311, 2018.
  * [33] JC Bezdek, RJ Hathaway, RE Howard, CA Wilson, and MP Windham.  Local convergence analysis of a grouped variable version of coordinate descent.  Journal of Optimization Theory and Applications, 54(3):471--477, 1987.
  * [34] Patrick A Stokes and Patrick L Purdon.  A study of problems encountered in granger causality analysis from a neuroscience perspective.  Proceedings of the national academy of sciences, 114(34):E7063--E7072, 2017.
  * [35] György Buzsáki, Costas A Anastassiou, and Christof Koch.  The origin of extracellular fields and currents—eeg, ecog, lfp and spikes.  Nature reviews neuroscience, 13(6):407--420, 2012.
  * [36] Sylvain Baillet, John C Mosher, and Richard M Leahy.  Electromagnetic brain mapping.  IEEE Signal processing magazine, 18(6):14--30, 2001.
  * [37] Lucas C Parra, Clay D Spence, Adam D Gerson, and Paul Sajda.  Recipes for the linear analysis of eeg.  Neuroimage, 28(2):326--341, 2005.
  * [38] Jacek P Dmochowski, Matthew A Bezdek, Brian P Abelson, John S Johnson, Eric H Schumacher, and Lucas C Parra.  Audience preferences are predicted by temporal reliability of neural processing.  Nature communications, 5(1):1--9, 2014.
  * [39] Wolfgang Klimesch.  Eeg alpha and theta oscillations reflect cognitive and memory performance: a review and analysis.  Brain research reviews, 29(2-3):169--195, 1999.
  * [40] David L Robinson.  The technical, neurological and psychological significance of ‘alpha’,‘delta’and ‘theta’waves confounded in eeg evoked potentials: a study of peak latencies.  Clinical Neurophysiology, 110(8):1427--1434, 1999.
  * [41] David L Robinson.  How brain arousal systems determine different temperament types and the major dimensions of personality.  Personality and Individual Differences, 31(8):1233--1259, 2001.
  * [42] Josh H McDermott.  The cocktail party problem.  Current Biology, 19(22):R1024--R1027, 2009.
  * [43] Harold Hotelling.  Relations between two sets of variates.  In Breakthroughs in statistics, pages 162--190. Springer, 1992.
  * [44] João R Sato, André Fujita, Elisson F Cardoso, Carlos E Thomaz, Michael J Brammer, and Edson Amaro Jr.  Analyzing the connectivity between regions of interest: an approach based on cluster granger causality for fmri data analysis.  Neuroimage, 52(4):1444--1455, 2010.
  * [45] Guorong Wu, Xujun Duan, Wei Liao, Qing Gao, and Huafu Chen.  Kernel canonical-correlation granger causality for multiple time series.  Physical Review E, 83(4):041921, 2011.
  * [46] Gene H Golub, Per Christian Hansen, and Dianne P O’Leary.  Tikhonov regularization and total least squares.  SIAM journal on matrix analysis and applications, 21(1):185--194, 1999.
  * [47] Shuixia Guo, Anil K Seth, Keith M Kendrick, Cong Zhou, and Jianfeng Feng.  Partial granger causality—eliminating exogenous inputs and latent variables.  Journal of neuroscience methods, 172(1):79--93, 2008.
  * [48] Karl J Friston, Lee Harrison, and Will Penny.  Dynamic causal modelling.  Neuroimage, 19(4):1273--1302, 2003.
  * [49] AR McIntosh and F Gonzalez-Lima.  Structural modeling of functional neural pathways mapped with 2-deoxyglucose: effects of acoustic startle habituation on the auditory system.  Brain research, 547(2):295--302, 1991.
  * [50] Arthur E Hoerl and Robert W Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics, 12(1):55--67, 1970.
  * [51] Jemima M Tabeart, Sarah L Dance, Amos S Lawless, Nancy K Nichols, and Joanne A Waller.  Improving the condition number of estimated covariance matrices.  Tellus A: Dynamic Meteorology and Oceanography, 72(1):1--19, 2020\.
  * [52] Mark Girolami and Colin Fyfe.  Negentropy and kurtosis as projection pursuit indices provide generalised ica algorithms.  In Advances in Neural Information Processing Systems Workshop, volume 9, 1996.
  * [53] James Theiler, Stephen Eubank, André Longtin, Bryan Galdrikian, and J Doyne Farmer.  Testing for nonlinearity in time series: the method of surrogate data.  Physica D: Nonlinear Phenomena, 58(1-4):77--94, 1992.
  * [54] Jan R Magnus and Heinz Neudecker.  Matrix differential calculus with applications in statistics and econometrics.  John Wiley & Sons, 2019.

## Supplementary Note 1: Closed-form expressions for latent covariance

In order to obtain a closed-form expression for the strength of causality
$\mathcal{G}_{y\rightarrow z}$ between driving signal $y(t)$ and driven signal
$z(t)$, expressions for the minimum mean squared error (MMSE) are required for
both the reduced and full linear regression models:

$\displaystyle z(t)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{L}h_{l}z(t-l)+\epsilon_{r}(t)~{}~{}~{}\mathrm{(reduced~{}model)}$
(65) $\displaystyle z(t)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{L}g_{1l}z(t-l)+\sum_{l=1}^{L}g_{2l}y(t-l)+\epsilon_{f}(t)~{}~{}~{}\mathrm{(full~{}model)}$
(66)

where $h_{l}$ are the coefficients of the temporal filter predicting the
driven signal from its past in the reduced model, $g_{1l}$ are the
coefficients of the filter predicting the driven signal from its own past in
the full model, and $g_{2l}$ are the temporal filter weights of the filter
predicting the driven signal from the past of the driving signal. The
regression models can be more compactly written in vector notation as:

$\displaystyle z(t)$ $\displaystyle=$
$\displaystyle\mathbf{h}^{T}\mathbf{z}_{p}+\epsilon_{r}(t)$ (67)
$\displaystyle z(t)$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\mathbf{g}_{1}&\mathbf{g}_{2}\end{array}\right]\left[\begin{array}[]{cc}\mathbf{z}_{p}\\\
\mathbf{y}_{p}\end{array}\right]+\epsilon_{f}(t).$ (71)

The coefficients of the filter $\mathbf{h}$ that minimizes the residual in the
reduced model are given by [26]:

$\displaystyle\mathbf{h}=\mathbf{Q}^{-1}\mathbf{q},$ (72)

where $\mathbf{q}=E\left\\{z(t)\mathbf{z}_{p}(t)\right\\}$ is the covariance
vector between the desired signal $z$ and its own past, whose $l$th element is
given by $E\left\\{z(t)z(t-l)\right\\}$, and where
$\mathbf{Q}=E\left\\{\mathbf{z}_{p}(t)\mathbf{z}^{T}_{p}(t)\right\\}$ is the
covariance matrix of $\mathbf{z}_{p}(t)$, where the element at row $i$ and
column $j$ is given by
$E\left\\{z(t-i)z(t-j)\right\\}=E\left\\{z(t)z(t-j+i)\right\\}$ under the
assumption of wide-sense stationary observation data. By substituting (72)
into (67) and solving for the residual power, the corresponding MMSE follows
as:

$\displaystyle\Phi_{r}=E\\{\epsilon_{r}^{2}\\}=\sigma_{z}^{2}-\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{q}.$
(73)

It is required to express $\sigma_{z}$, $\mathbf{q}$, and $\mathbf{Q}$ in
terms of the projection vector $\mathbf{v}$ and the statistics of the observed
data $\mathbf{x}(t)$. The power of the desired signal $z$ in the regression
models is given by:

$\displaystyle\sigma_{z}^{2}=E\\{z^{2}(t)\\}=\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v},$
(74)

where
$\mathbf{\Sigma}(\tau)=E\left\\{\mathbf{x}(t)\mathbf{x}^{T}(t-\tau)\right\\}$
is the lagged covariance matrix of the observed data. Substitution of
$z(t)=\mathbf{v}^{T}\mathbf{x}(t)$ into the expression for $\mathbf{q}$ leads
to:

$\displaystyle\mathbf{q}$ $\displaystyle=$ $\displaystyle
E\left\\{z(t)\mathbf{z}_{p}(t)\right\\}$ (79) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}E\left\\{z(t)z(t-1)\right\\}\\\
E\left\\{z(t)z(t-2)\right\\}\\\ \vdots\\\
E\left\\{z(t)z(t-L)\right\\}\end{array}\right)$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{v}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(2)\mathbf{v}\\\ \vdots\\\
\mathbf{v}^{T}\mathbf{\Sigma}(L)\mathbf{v}\end{array}\right)$ (84)
$\displaystyle=$ $\displaystyle\left(\begin{array}[]{cccc}\mathbf{v}^{T}&&&\\\
&\mathbf{v}^{T}&&\\\ &&\ddots&\\\
&&&\mathbf{v}^{T}\end{array}\right)\left(\begin{array}[]{cccc}\mathbf{\Sigma}(1)&&&\\\
&\mathbf{\Sigma}(2)&&\\\ &&\ddots&\\\
&&&\mathbf{\Sigma}(L)\end{array}\right)\left(\begin{array}[]{c}\mathbf{v}\\\
\mathbf{v}\\\ \vdots\\\ \mathbf{v}\end{array}\right).$ (97)

By utilizing the Kronecker product $\otimes$, one can write (97) as the
following matrix product:

$\displaystyle\mathbf{q}$ $\displaystyle=$
$\displaystyle\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right),$
(98)

where

$\displaystyle\mathbf{\Sigma}_{1:L}=\left(\begin{array}[]{cccc}\mathbf{\Sigma}(1)&\mathbf{0}&\ldots&\mathbf{0}\\\
\mathbf{0}&\mathbf{\Sigma}(2)&\ldots&\mathbf{0}\\\
\vdots&\mathbf{0}&\ddots&\vdots\\\
\mathbf{0}&\ldots&\ldots&\mathbf{\Sigma}(L)\\\ \end{array}\right)$

is an $LD$-by-$LD$ block diagonal covariance matrix, $\mathbf{1}_{L}$ is a
vector of all ones, and $\mathbf{I}_{L}$ is the $L$-by-$L$ identity matrix.

Similarly, the covariance matrix $\mathbf{Q}$ can be written as:

$\displaystyle\mathbf{Q}$ $\displaystyle=$ $\displaystyle
E\left\\{\mathbf{z}_{p}(t)\mathbf{z}^{T}_{p}(t)\right\\}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cccc}E\left\\{z^{2}(t-1)\right\\}&E\left\\{z(t-1)z(t-2)\right\\}&\ldots&E\left\\{z(t-1)z(t-L)\right\\}\\\
E\left\\{z(t-2)z(t-1)\right\\}&E\left\\{z^{2}(t-2)\right\\}&\ldots&E\left\\{z(t-2)z(t-L)\right\\}\\\
\vdots&&\ddots&\vdots\\\
E\left\\{z(t-L)z(t-1)\right\\}&&\ldots&E\left\\{z^{2}(t-L)\right\\}\end{array}\right]$
$\displaystyle=$
$\displaystyle\left[\begin{array}[]{cccc}\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{v}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-1)\mathbf{v}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(-1)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-2)\mathbf{v}]\\\
\vdots&&\ddots&\vdots\\\
\mathbf{v}^{T}\mathbf{\Sigma}(-L+1)\mathbf{v}&&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}\end{array}\right]$
$\displaystyle=$
$\displaystyle\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right),$

where

$\displaystyle\tilde{\mathbf{\Sigma}}=\left(\begin{array}[]{cccc}\mathbf{\Sigma}(0)&\mathbf{\Sigma}(-1)&\ldots&\mathbf{\Sigma}(-L+1)\\\
\mathbf{\Sigma}(1)&\mathbf{\Sigma}(0)&\ldots&\mathbf{\Sigma}(-L+2)\\\
\vdots&\mathbf{\Sigma}(1)&\ddots&\vdots\\\
\mathbf{\Sigma}(L-1)&\ldots&\ldots&\mathbf{\Sigma}(0)\\\ \end{array}\right)$

is an $LD$-by-$LD$ block Toeplitz matrix. The elements of matrix $\mathbf{Q}$
can now be differentiated with respect to the elements of the projection
vector $\mathbf{v}$. As described in Supplementary Note 2, this necessitates
the employment of matrix differentials [54].

A similar development can be performed to derive at closed-form expressions
for the covariance vector $\mathbf{r}$ and covariance matrix $\mathbf{R}$. The
Wiener filter of the full regression model (67) is given by:

$\displaystyle\mathbf{g}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\mathbf{g}_{1}\\\
\mathbf{g}_{2}\end{array}\right)$ (105) $\displaystyle=$
$\displaystyle\mathbf{R}^{-1}\mathbf{r},$ (106)

and the corresponding MMSE is given by:

$\displaystyle\Phi_{f}=E\\{\epsilon_{f}^{2}\\}=\sigma_{z}^{2}-\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}.$
(107)

Substituting $z(t)=\mathbf{v}^{T}\mathbf{x}(t)$ and
$y(t)=\mathbf{w}^{T}\mathbf{x}(t)$ into the definition of $\mathbf{r}$ yields:

$\displaystyle\mathbf{r}$ $\displaystyle=$ $\displaystyle
E\left\\{z(t)\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)\right\\}$ (110) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}E\left\\{z(t)z(t-1)\right\\}\\\
E\left\\{z(t)z(t-2)\right\\}\\\ \vdots\\\ E\left\\{z(t)z(t-L)\right\\}\\\
E\left\\{z(t)y(t-1)\right\\}\\\ E\left\\{z(t)y(t-2)\right\\}\\\ \vdots\\\
E\left\\{z(t)y(t-L)\right\\}\end{array}\right)$ (119) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{v}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(2)\mathbf{v}\\\ \vdots\\\
\mathbf{v}^{T}\mathbf{\Sigma}(L)\mathbf{v}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{w}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(2)\mathbf{w}\\\ \vdots\\\
\mathbf{v}^{T}\mathbf{\Sigma}(L)\mathbf{w}\end{array}\right).$ (128)

It is straightforward to verify that (128) can be factored according to:

$\displaystyle\mathbf{r}=\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right).$ (131)

Finally, the explicit expression for covariance matrix $\mathbf{R}$ is given
by:

$\displaystyle\mathbf{R}$ $\displaystyle=$ $\displaystyle
E\left\\{\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)\left(\begin{array}[]{c}\mathbf{z}_{p}(t)\\\
\mathbf{y}_{p}(t)\end{array}\right)^{T}\right\\}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccccc}\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{v}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-1)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{w}&\mathbf{v}^{T}\mathbf{\Sigma}(1)\mathbf{w}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-1)\mathbf{w}\\\
\mathbf{v}^{T}\mathbf{\Sigma}(-1)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-2)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(-1)\mathbf{w}&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{w}&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(L-2)\mathbf{w}\\\
\vdots&&\ddots&\vdots&&&\\\
\mathbf{v}^{T}\mathbf{\Sigma}(-L+1)\mathbf{v}&&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\mathbf{v}^{T}\mathbf{\Sigma}(-L+1)\mathbf{w}&&\ldots&\mathbf{v}^{T}\mathbf{\Sigma}(0)\mathbf{w}\\\
\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\mathbf{w}^{T}\mathbf{\Sigma}(1)\mathbf{v}&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(L-1)\mathbf{v}&\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{w}&\mathbf{w}^{T}\mathbf{\Sigma}(1)\mathbf{w}&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(L-1)\mathbf{w}\\\
\mathbf{w}^{T}\mathbf{\Sigma}(-1)\mathbf{w}&\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(L-2)\mathbf{v}&\mathbf{w}^{T}\mathbf{\Sigma}(-1)\mathbf{w}&\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{w}&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(L-2)\mathbf{w}\\\
\vdots&&\ddots&\vdots&&&\\\
\mathbf{w}^{T}\mathbf{\Sigma}(-L+1)\mathbf{v}&&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{v}&\mathbf{w}^{T}\mathbf{\Sigma}(-L+1)\mathbf{w}&&\ldots&\mathbf{w}^{T}\mathbf{\Sigma}(0)\mathbf{w}\end{array}\right),$

which can be factored according to:

$\displaystyle\mathbf{R}=\left(\begin{array}[]{c}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\left(\begin{array}[]{cc}\mathbf{I}_{L}\otimes\mathbf{v}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}\end{array}\right).$ (138)

## Supplementary Note 2: Gradient of objective function

The goal is to differentiate the objective function:

$\displaystyle\mathcal{G}_{y\rightarrow z}$ $\displaystyle=$ $\displaystyle
1-\frac{\Phi_{f}}{\Phi_{r}}$ (139) $\displaystyle=$ $\displaystyle
1-\frac{\sigma_{z}^{2}-\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}}{\sigma_{z}^{2}-\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{q}}$

with respect to the projection vectors $\mathbf{v}$ and $\mathbf{w}$. The
derivation below relies on the chain rule, and involves the differentiation of
matrices $\mathbf{R}$ and $\mathbf{Q}$ with respect to vectors $\mathbf{v}$
and $\mathbf{w}$. The reader is referred to Magnus and Neudecker [54] for an
excellent treatment of matrix differentials, which is utilized here.

Taking differentials of both sides of (139) yields:

$\displaystyle d\mathcal{G}_{y\rightarrow z}$
$\displaystyle=\frac{\Phi_{f}d\Phi_{r}-\Phi_{r}d\Phi_{f}}{\Phi_{r}^{2}}$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-d\left(\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{q}\right)\right)-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-d\left(\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}\right)\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-d\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{q}-\mathbf{q}^{T}d\mathbf{Q}^{-1}\mathbf{q}-\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{q}\right)-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-d\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{r}-\mathbf{r}^{T}d\mathbf{R}^{-1}\mathbf{r}-\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{r}\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-2\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{q}+\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{Q}\mathbf{Q}^{-1}\mathbf{q}\right)-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-2\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{r}+\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{R}\mathbf{R}^{-1}\mathbf{r}\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-2\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{q}+\mathrm{tr}~{}\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{Q}\mathbf{Q}^{-1}\mathbf{q}\right)-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-2\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{r}+\mathrm{tr}~{}\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{R}\mathbf{R}^{-1}\mathbf{r}\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-2\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{q}+\mathrm{tr}~{}\mathbf{Q}^{-1}\mathbf{q}\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{Q}\right)-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-2\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{r}+\mathrm{tr}~{}\mathbf{R}^{-1}\mathbf{r}\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{R}\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(d\sigma_{z}^{2}-2\mathbf{q}^{T}\mathbf{Q}^{-1}d\mathbf{q}+\mathrm{vec}~{}\left(\mathbf{Q}^{-1}\mathbf{q}\mathbf{q}^{T}\mathbf{Q}^{-1}\right)^{T}\mathrm{vec}~{}d\mathbf{Q}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{\Phi_{r}}\left(d\sigma_{z}^{2}-2\mathbf{r}^{T}\mathbf{R}^{-1}d\mathbf{r}+\mathrm{vec}~{}\left(\mathbf{R}^{-1}\mathbf{r}\mathbf{r}^{T}\mathbf{R}^{-1}\right)^{T}\mathrm{vec}~{}d\mathbf{R}\right)$
$\displaystyle=\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(2\mathbf{v}^{T}\mathbf{\Sigma}(0)d\mathbf{v}-2\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{J}_{q}d\mathbf{v}+\mathrm{vec}\left(\mathbf{Q}^{-1}\mathbf{q}\mathbf{q}^{T}\mathbf{Q}^{-1}\right)^{T}\mathbf{J}_{Q}d\mathbf{v}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\frac{1}{\Phi_{r}}\left(2\mathbf{v}^{T}\mathbf{\Sigma}(0)d\mathbf{v}-2\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{J}_{r}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right)+\mathrm{vec}\left(\mathbf{R}^{-1}\mathbf{r}\mathbf{r}^{T}\mathbf{R}^{-1}\right)^{T}\mathbf{J}_{R}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right)\right),$ (144)

where $\mathrm{tr}$ is the matrix trace operator, $\mathrm{vec}$ is an
operator that transforms a matrix into a column vector by stacking the
columns, and where the following Jacobian matrices have been defined:

$\displaystyle d\mathbf{q}$ $\displaystyle=$
$\displaystyle\mathbf{J}_{q}d\mathbf{v}$
$\displaystyle\mathrm{vec}~{}d\mathbf{Q}$ $\displaystyle=$
$\displaystyle\mathbf{J}_{Q}d\mathbf{v}$ $\displaystyle d\mathbf{r}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{r}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right)$ $\displaystyle\mathrm{vec}~{}d\mathbf{R}$
$\displaystyle=$
$\displaystyle\mathbf{J}_{R}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right).$

Once closed-form expressions for these Jacobians are derived, it remains to
substitute those expressions into (144).

#### Differential of $\mathbf{q}$

The covariance vector $\mathbf{q}$ is defined by:

$\displaystyle\mathbf{q}=\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right).$
(147)

The goal is to arrive at an expression of the form:

$\displaystyle d\mathbf{q}=\mathbf{J}_{q}d\mathbf{v},$ (148)

where $\mathbf{J}_{q}$ is the Jacobian matrix that one seeks to identify. The
following identity (often termed the ‘‘vec’’ rule) will prove useful
throughout:

$\displaystyle\mathrm{vec}(\mathbf{A}\mathbf{B}\mathbf{C})=\left(\mathbf{C}^{T}\otimes\mathbf{A}\right)\mathrm{vec}\left(\mathbf{B}\right),$

where matrices $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ are defined such
that the conventional matrix product $\mathbf{A}\mathbf{B}\mathbf{C}$ is
valid. Taking differentials of both sides of (147) and applying the vec rule,
one obtains:

$\displaystyle d\mathbf{q}$ $\displaystyle=$
$\displaystyle\left(\left(\mathbf{1}_{L}^{T}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}^{T}\otimes\mathbf{I}_{L}\right)d~{}\mathrm{vec}\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)+\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\mathbf{\Sigma}_{1:L}~{}d~{}\mathrm{vec}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right).$
(149)

The differentials on the right-hand side may be written as [54]:

$\displaystyle
d~{}\mathrm{vec}\left(\mathbf{1}_{L}\otimes\mathbf{v}\right)=\left(\mathbf{1}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{v},$
(150)

and

$\displaystyle
d~{}\mathrm{vec}\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)=\left(\mathbf{I}_{L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{v},$
(151)

where $\mathbf{K}_{DL}$ is a commutation matrix satisfying:

$\displaystyle\mathbf{K}_{mn}\mathrm{vec}\left(\mathbf{D}\right)=\mathrm{vec}\left(\mathbf{D}^{T}\right)$

for $m$-by-$n$ matrix $\mathbf{D}$. In deriving (150) and (151), the property
on p. 206 of Magnus and Neudecker [54] has been invoked to move the
differential outside of the vec operator. Substituting (150) and (151) into
(149), one obtains the required result:

$\displaystyle
d\mathbf{q}=\left[\left(\left(\mathbf{1}_{L}^{T}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}^{T}\otimes\mathbf{I}_{L}\right)\left(\mathbf{I}_{L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)+\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{I}_{D}\right)\right]d\mathbf{v},$

where the Jacobian is identified as:

$\displaystyle\mathbf{J}_{q}=\left(\left(\mathbf{1}_{L}^{T}\otimes\mathbf{v}^{T}\right)\mathbf{\Sigma}_{1:L}^{T}\otimes\mathbf{I}_{L}\right)\left(\mathbf{I}_{L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)+\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\mathbf{\Sigma}_{1:L}\left(\mathbf{1}_{L}\otimes\mathbf{I}_{D}\right).$

#### Differential of $\mathbf{Q}$

The $L$-by-$L$ covariance matrix of predictors in the reduced regression model
is given by:

$\displaystyle\mathbf{Q}=\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right).$

The differential follows as:

$\displaystyle d\mathbf{Q}=\left(\mathbf{I}_{L}\otimes
d\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)+\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)^{T}\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes
d\mathbf{v}\right).$ (152)

By vectorizing both sides of (152) and applying the vec rule to both terms on
the right-hand side, one obtains:

$\displaystyle\mathrm{vec}~{}d\mathbf{Q}=\left[\left(\tilde{\mathbf{\Sigma}}\left(\mathbf{I}_{L}\otimes\mathbf{v}\right)\right)^{T}\otimes\mathbf{I}_{L}\right]\mathrm{vec}\left(\mathbf{I}_{L}\otimes
d\mathbf{v}^{T}\right)+\left(\mathbf{I}_{L}\otimes\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\tilde{\mathbf{\Sigma}}\right)\mathrm{vec}\left(\mathbf{I}_{L}\otimes
d\mathbf{v}\right).$

Again using the property at the bottom of page 206 in Magnus and Neudecker
[54], one obtains the following expression:

$\displaystyle\mathrm{vec}~{}d\mathbf{Q}$ $\displaystyle=$
$\displaystyle\left[\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\tilde{\mathbf{\Sigma}}^{T}\otimes\mathbf{I}_{L}\right]\left(\mathbf{I}_{L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{v}+$
$\displaystyle\left(\mathbf{I}_{L}\otimes\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\tilde{\mathbf{\Sigma}}\right)\left(\mathbf{I}_{L}\otimes\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{v},$

from which the Jacobian is identified as:

$\displaystyle\mathbf{J}_{Q}$ $\displaystyle=$
$\displaystyle\left[\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\tilde{\mathbf{\Sigma}}^{T}\otimes\mathbf{I}_{L}\right]\left(\mathbf{I}_{L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)+\left(\mathbf{I}_{L}\otimes\left(\mathbf{I}_{L}\otimes\mathbf{v}^{T}\right)\tilde{\mathbf{\Sigma}}\right)\left(\mathbf{I}_{L}\otimes\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right).$

#### Differential of r

The $2L$-dimensional covariance vector in the full regression model is given
by:

$\displaystyle\mathbf{r}=\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right).$ (155)

Taking differentials of (155), one obtains:

$\displaystyle d\mathbf{r}=\left(\mathbf{I}_{2L}\otimes
d\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)+\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes
d\mathbf{v}\\\ \mathbf{1}_{L}\otimes d\mathbf{w}\end{array}\right).$ (160)

Applying the vec operator to both sides of (160) yields:

$\displaystyle\mathrm{vec}~{}d\mathbf{r}=\mathrm{vec}\left(\left(\mathbf{I}_{2L}\otimes
d\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)\right)+\mathrm{vec}\left(\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes
d\mathbf{v}\\\ \mathbf{1}_{L}\otimes d\mathbf{w}\end{array}\right)\right).$
(165)

The vec rule can now be applied to the right hand side of (165):

$\displaystyle
d\mathbf{r}=\left(\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)^{T}\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)^{T}\otimes\mathbf{I}_{2L}\right)\mathrm{vec}\left(\mathbf{I}_{2L}\otimes\mathbf{dv}^{T}\right)+\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\mathrm{vec}~{}\left(\mathbf{1}_{L}\otimes
d\mathbf{v}\right)\\\ \mathrm{vec}~{}\left(\mathbf{1}_{L}\otimes
d\mathbf{w}\right)\end{array}\right).$

Applying the property on page 206 of Magnus and Neudecker [54], the resulting
expression may be written as:

$\displaystyle d\mathbf{r}$ $\displaystyle=$
$\displaystyle\left(\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)^{T}\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)^{T}\otimes\mathbf{I}_{2L}\right)\left(\mathbf{I}_{2L}\otimes\mathbf{K}_{D,2L}\right)\left(\mathrm{vec}~{}\mathbf{I}_{2L}\otimes\mathbf{I}_{D}\right)\mathbf{dv}+$
$\displaystyle\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{c}\left(\mathrm{vec}~{}\mathbf{1}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{v}\\\
\left(\mathrm{vec}~{}\mathbf{1}_{L}\otimes\mathbf{I}_{D}\right)d\mathbf{w}\end{array}\right),$

which can then be expressed as:

$\displaystyle d\mathbf{r}$ $\displaystyle=$
$\displaystyle\left(\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)^{T}\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)^{T}\otimes\mathbf{I}_{2L}\right)\left(\mathbf{I}_{2L}\otimes\mathbf{K}_{D,2L}\right)\left(\mathrm{vec}~{}\mathbf{I}_{2L}\otimes\mathbf{I}_{D}\right)\mathbf{dv}+$
$\displaystyle\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{cc}\mathbf{1}_{L}\otimes\mathbf{I}_{D}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}_{L}\otimes\mathbf{I}_{D}\end{array}\right)\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right).$

The Jacobian of $\mathbf{r}$ with respect to $\mathbf{v}$ and $\mathbf{w}$ can
now be identified:

$\displaystyle\mathbf{J}_{r}$ $\displaystyle=$
$\displaystyle\left[\begin{array}[]{cc}\left(\left(\begin{array}[]{c}\mathbf{1}_{L}\otimes\mathbf{v}\\\
\mathbf{1}_{L}\otimes\mathbf{w}\end{array}\right)^{T}\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)^{T}\otimes\mathbf{I}_{2L}\right)\left(\mathbf{I}_{2L}\otimes\mathbf{K}_{D,2L}\right)\left(\mathrm{vec}~{}\mathbf{I}_{2L}\otimes\mathbf{I}_{D}\right)&\mathbf{0}\end{array}\right]+$
$\displaystyle\left(\mathbf{I}_{2L}\otimes\mathbf{v}^{T}\right)\left(\mathbf{I}_{2}\otimes\mathbf{\Sigma}_{1:L}\right)\left(\begin{array}[]{cc}\mathbf{1}_{L}\otimes\mathbf{I}_{D}&\mathbf{0}\\\
\mathbf{0}&\mathbf{1}_{L}\otimes\mathbf{I}_{D}\end{array}\right).$

#### Differential of $\mathbf{R}$

The covariance matrix of predictors in the full regression model is given by
the $2L$-by-$2L$ matrix:

$\displaystyle\mathbf{R}=\left(\begin{array}[]{c}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\left(\begin{array}[]{cc}\mathbf{I}_{L}\otimes\mathbf{v}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}\end{array}\right).$

By following the same development as above, the differential of the elements
of $\mathbf{R}$ is given by:

$\displaystyle\mathrm{vec}~{}d\mathbf{R}$
$\displaystyle=\left[\left(\begin{array}[]{cc}\mathbf{I}_{L}\otimes\mathbf{v}^{T}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}^{T}\right)\otimes\mathbf{I}_{2L}\right]\left(\mathbf{K}_{2,2LD}\otimes\mathbf{I}_{L}\right)^{-1}\left\\{\mathbf{I}_{2}\otimes\left[\left(\mathbf{I}_{2L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\left(\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\right)\otimes\mathbf{I}_{D}\right)\right]\right\\}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right)$ (179)
$\displaystyle+\left[\mathbf{I}_{2L}\otimes\left(\begin{array}[]{c}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\right]\left(\mathbf{I}_{2}\otimes\mathbf{K}_{2,L}\otimes\mathbf{I}_{LD}\right)^{-1}\left[\mathbf{I}_{4,2}\otimes\left(\mathbf{I}_{L}\otimes\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\right]\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right),$ (184)

where $\mathbf{I}_{4,2}$ is a $4$ by $2$ matrix with ones at row 1, column 1
and at row 4, column 2, and zeros elsewhere. Thus, the Jacobian matrix
$\mathbf{J}_{R}$ is identified as:

$\displaystyle\mathbf{J}_{R}$
$\displaystyle=\left[\left(\begin{array}[]{cc}\mathbf{I}_{L}\otimes\mathbf{v}^{T}&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}^{T}\right)\otimes\mathbf{I}_{2L}\right]\left(\mathbf{K}_{2,2LD}\otimes\mathbf{I}_{L}\right)^{-1}\left\\{\mathbf{I}_{2}\otimes\left[\left(\mathbf{I}_{2L}\otimes\mathbf{K}_{DL}\right)\left(\mathrm{vec}~{}\left(\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\right)\otimes\mathbf{I}_{D}\right)\right]\right\\}$
(187)
$\displaystyle+\left[\mathbf{I}_{2L}\otimes\left(\begin{array}[]{c}\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{v}^{T}\\\
\mathbf{1}_{2}^{T}\otimes\mathbf{I}_{L}\otimes\mathbf{w}^{T}\end{array}\right)\left(\mathbf{I}_{2}\otimes\tilde{\mathbf{\Sigma}}\right)\right]\left(\mathbf{I}_{2}\otimes\mathbf{K}_{2,L}\otimes\mathbf{I}_{LD}\right)^{-1}\left[\mathbf{I}_{4,2}\otimes\left(\mathbf{I}_{L}\otimes\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\left(\mathrm{vec}~{}\mathbf{I}_{L}\otimes\mathbf{I}_{D}\right)\right].$
(190)

#### Gradient of objective function

Having identified the Jacobians $\mathbf{J}_{q}$, $\mathbf{J}_{Q}$,
$\mathbf{J}_{r}$, and $\mathbf{J}_{R}$, the final expression for the gradient
of the objective function may now assembled.

$\displaystyle d\mathcal{G}_{y\rightarrow z}$ $\displaystyle=$
$\displaystyle\frac{2}{\Phi_{r}}\left(\frac{\Phi_{f}}{\Phi_{r}}-1\right)\mathbf{v}^{T}\mathbf{\Sigma}(0)d\mathbf{v}$
(196)
$\displaystyle-\frac{2\Phi_{f}}{\Phi_{r}^{2}}\mathbf{q}^{T}\mathbf{Q}^{-1}\mathbf{J}_{q}d\mathbf{v}$
$\displaystyle+\frac{\Phi_{f}}{\Phi_{r}^{2}}\mathrm{vec}\left(\mathbf{Q}^{-1}\mathbf{q}\mathbf{q}^{T}\mathbf{Q}^{-1}\right)^{T}\mathbf{J}_{Q}d\mathbf{v}$
$\displaystyle\frac{2}{\Phi_{r}}\mathbf{r}^{T}\mathbf{R}^{-1}\mathbf{J}_{r}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right)$
$\displaystyle-\frac{1}{\Phi_{r}}\mathrm{vec}\left(\mathbf{R}^{-1}\mathbf{r}\mathbf{r}^{T}\mathbf{R}^{-1}\right)^{T}\mathbf{J}_{R}\left(\begin{array}[]{c}d\mathbf{v}\\\
d\mathbf{w}\end{array}\right),$

from which one identifies the gradient of the objective function as:

$\displaystyle\mathbf{\nabla}\mathcal{G}=-\frac{2\mathcal{G}}{\Phi_{r}}\left(\begin{array}[]{c}\mathbf{\Sigma}(0)\mathbf{v}\\\
\mathbf{0}\end{array}\right)-\frac{2\Phi_{f}}{\Phi_{r}^{2}}\left(\begin{array}[]{c}\mathbf{J}_{q}^{T}\mathbf{Q}^{-1}\mathbf{q}\\\
\mathbf{0}\end{array}\right)+\frac{\Phi_{f}}{\Phi_{r}^{2}}\left(\begin{array}[]{c}\mathbf{J}_{Q}^{T}\mathrm{vec}\left(\mathbf{Q}^{-1}\mathbf{q}\mathbf{q}^{T}\mathbf{Q}^{-1}\right)\\\
\mathbf{0}\end{array}\right)$ (203)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\frac{2}{\Phi_{r}}\mathbf{J}_{r}^{T}\mathbf{R}^{-1}\mathbf{r}-\frac{1}{\Phi_{r}}\mathbf{J}_{R}^{T}\mathrm{vec}\left(\mathbf{R}^{-1}\mathbf{r}\mathbf{r}^{T}\mathbf{R}^{-1}\right).$

## Supplementary Note 3: Grouped coordinate descent algorithm

To solve the optimization problem in Eqn. (34), a grouped coordinate descent
algorithm was proposed. The procedure is described in Algorithm 1, where
$\mathbf{0}$ is a vector of zeros, $\sigma$ is a small positive number,
$\mathcal{N}(\mathbf{m},\mathbf{C})$ is the multivariate normal distribution
with mean vector $\mathbf{m}$ and covariance matrix $\mathbf{C}$, $\mathbf{X}$
is a $D$-by-$T$ matrix storing the observed data, $\mathbf{Y}_{p}$ is an
$L$-by-$T$ convolution matrix allowing the regression of $\mathbf{y}_{p}(t)$
onto $\mathbf{X}$, $\mathrm{convmatrix}$ is a routine that produces a
convolution matrix, # denotes the Moore-Penrose pseudoinverse, and
$\mathbf{X}_{r}$ is the observed data after regressing out the contribution of
the driving signal(s). The function $\mathcal{G}(\mathbf{w},\mathbf{v})$
evaluates the strength of causality (33) between latent sources
$y(t)=\mathbf{w}^{T}\mathbf{x}(t)$ and $z(t)=\mathbf{v}^{T}\mathbf{x}(t)$
(forward time), and $\mathcal{G}^{\mathrm{tr}}(\mathbf{v},\mathbf{w})$
evaluates the strength of causality between
$y(t)=\mathbf{v}^{T}\mathbf{x}(-t)$ and $z(t)=\mathbf{w}^{T}\mathbf{x}(-t)$
(reversed time).

Algorithm 1 Grouped coordinate descent for maximizing strength of causality
among latent sources.

$p\leftarrow 1$

while $p\leq P$ do

$\mathbf{w}_{t}^{(p)}\leftarrow\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$

$\mathbf{v}_{t}^{(p)}\leftarrow\mathcal{N}(\mathbf{0},\sigma^{2}\mathbf{I})$

repeat

$\mathbf{v}^{\ast}\leftarrow\arg\max_{\mathbf{v}}$
$\mathcal{G}(\mathbf{w}_{t}^{(p)},\mathbf{v})+\mathcal{G}^{\mathrm{tr}}(\mathbf{v},\mathbf{w}_{t}^{(p)})$

$\mathbf{v}_{t}^{(p)}\leftarrow\mathbf{v}^{\ast}$

$\mathbf{w}^{\ast}\leftarrow\arg\max_{\mathbf{w}}\mathcal{G}(\mathbf{w},\mathbf{v}_{t}^{(p)})+\mathcal{G}^{\mathrm{tr}}(\mathbf{v}_{t}^{(p)},\mathbf{w})$

$\mathbf{w}_{t}^{(p)}\leftarrow\mathbf{w}^{\ast}$

until converged

$\mathbf{y}_{p}\leftarrow{\mathbf{w}_{t}^{(p)}}^{T}\mathbf{X}$

$\mathbf{Y}_{p}\leftarrow\mathrm{convmatrix}\left(\mathbf{y}_{p}\right)$

$\mathbf{X}_{r}\leftarrow\mathbf{X}\left(\mathbf{I}-\mathbf{Y}_{p}^{\\#}\mathbf{Y}_{p}\right)$

$\mathbf{X}\leftarrow\mathbf{X}_{r}$

$p\leftarrow p+1$

end while

A MATLAB implementation of Algorithm 1 is provided at dmochow.github.io/gca.

Figure 1:  Extracting latent causal components from multiple time series. (a)
In the framework proposed here, Granger causal relations occur between pairs
of latent sources whose activity governs the system dynamics. Due to a linear
mixing process, the observed data manifests as a superposition of the latent
sources. Consequently, direct application of Granger Causality to the observed
data may not correctly identify the underlying causal structure. Instead, the
approach proposed here consists of projecting the data into _pairs_ of
components that aim to recover the directed causal connections present in the
latent space. (b) The connection weights between the observed and recovered
data are determined by solving an optimization problem that maximizes the
Granger Causality between driving signal components $y$ and driven signal
components $z$. Note that although the latent components are formed
instantaneously once the filters $\mathbf{w}$ and $\mathbf{v}$ are available,
learning of the optimal filter weights requires observations collected over a
temporal aperture of sufficient extent to estimate the data covariance.
Figure 2:  Identifying the latent causal structure of a three-element system.
(a) A simple system with three latent sources was simulated, with
$s_{1}\rightarrow s_{2}$ and $s_{2}\rightarrow s_{3}$. Data from a single
realization is shown. (b) Projection of the source signals to a four-
dimensional observation vector was modeled by a mixing matrix whose values
were drawn from the uniform distribution. (c) Two pairs of latent sources were
recovered by the proposed method. (d) The strength of causality between all
pairs of latent sources, where it is evident that the two connections are
$\mathcal{G}_{s_{1}\rightarrow s_{2}}=0.11\pm 0.0011$ and
$\mathcal{G}_{s_{2}\rightarrow s_{3}}=0.10\pm 0.001$ (mean $\pm$ sem across
100 realizations). (e) When computed between pairs of observed signals, the
strength of causality is visibly dampened, with a maximum value of $0.066\pm
0.002$, occurring at $x_{1}\rightarrow x_{3}$ for the depicted realization.
Notice also that the system’s causal structure is not readily apparent in the
observed data. (f) The strength of causality between recovered components
showed two strong connections: $\mathcal{G}_{y_{1}\rightarrow z_{1}}=0.11\pm
0.0014$ ($p=4.7\times 10^{-18}$, $n=100$, Wilcoxon signed rank test against
the maximum strength of causality in the observed data) and
$\mathcal{G}_{y_{2}\rightarrow z_{2}}=0.10\pm 0.0010$ ($p=4.0\times
10^{-18}$). (g) The mixing matrix relating the latent sources to the observed
signals, shown for a single realization. (h) The estimate of the mixing matrix
as derived by regressing the recovered components onto the observed data,
where $r^{2}=0.98\pm 0.004$. (i) Scatter between $s_{1}$ and $y_{1}$, where
the recovered driving signal of the first pair captured $98\pm 0.7$% of the
variability in latent source 1. (j) Same as (i) but now shown between $s_{2}$
and $z_{1}$, where more than $96\%$ of the variability in latent source 2 was
captured. (k) The driving signal of the second pair explained over 98% of the
variability in the second latent source. Note that both $z_{1}$ and $y_{2}$
captured $s_{2}$. (l) Same as (i) but now shown between $s_{3}$ and $z_{2}$,
where over $99\%$ of the variability in latent source 3 was captured.  Figure
3:  Learning causal relationships in the human brain. A data set consisting of
the scalp EEG of $n=12$ subjects viewing television advertisements was
employed to evaluate the causal relationships recovered by the proposed method
with those from conventional approaches. (a) The scalp EEG manifests as a
linear mixture of interacting electrical sources in the cerebral cortex. (b)
The strength of causality between all pairs of electrode signals, where the
largest value $0.073$ was observed between electrodes CP1 (left centro-
parietal) and CP4 (right centro-parietal). (c) The matrix of Granger Causality
among the first 10 principal components, where a maximum of $0.067$ was
detected between PCs 8 and 3. (d) Same as (c) but now shown for the first 10
independent components, where a maximum of $0.022$ was found between
components 6 and 3. (e) The proposed approach recovered three pairs of
components whose strength of causality was $0.32$, $0.16$, and $0.18$,
respectively, depicted in the diagonal elements of the matrix. (f) The scalp
topographies (top row) and power spectra (bottom row) of the driving and
driven signal in the first pair of recovered components. The driving signal
exhibited a right temporo-parietal distribution and large delta (1-3 Hz)
power, while the driven signal showed peak expression over the left central
electrodes. (g) Same as (f) but now for the second pair of recovered
components, where the driving signal originated over the left parieto-
occipital region and the driven signal over the right temporo-parietal region.
Notice also that the driving and driven signals showed opposing pattterns in
their power spectrum: $y_{2}$ was marked by low delta and high alpha (8-13 Hz)
power, while $z_{2}$ showed the opposite (high delta and low alpha). (h) The
scalp topographies and power spectra of the third pair of recovered
components. The driving signal showed bilateral expression in the parieto-
occipital region, while the driven signal was concentrated over the left
occipital region. Once again, the driving signal showed a high level of alpha
band power, and high delta power was observed in the power spectrum of the
driven signal.

Figure 4:  Probing causality in the cryptocurrency market. (a) Historical
price data was obtained for 19 cryptocurrencies spanning a time period of
approximately 15 months. The individual time series have been standardized to
accommodate a large dynamic range between currencies. (b) The conventional
Granger Causality measure between pairs of individual currencies was modest
(mean $\pm$ sd of $0.028\pm 0.023$), peaking at a value of 0.12 at ETC
$\rightarrow$ LTC. (c) On the other hand, the proposed method recovered two
pairs of latent variables with strong and statistically significant strengths
of causality: $\mathcal{G}_{y_{1}\rightarrow z_{1}}=0.40$ ($p<0.001$, non-
parametric permutation test with 1000 phase-randomized surrogate data records)
and $\mathcal{G}_{y_{2}\rightarrow z_{2}}=0.14$ ($p=0.008$), and a third pair
whose strength of causality fell short of significance but still exceeded
$95\%$ of individual currency pair values ($\mathcal{G}_{y_{3}\rightarrow
z_{3}}=0.080$, $p=0.13$). (d) The dynamics of driving signal $y_{1}$ clearly
precede those of the driven signal $z_{1}$, with the time series peaks
consistently occurring a few samples earlier for $y_{1}$ (see inset). (d) The
weights of the filter employed to construct the driving signal ($\mathbf{w}$,
top row) and those used to construct the driven signal ($\mathbf{v}$, bottom
row) show that the largest contribution to $y_{1}$ was exhibited by the
currencies BNB and ETC, while QTUM and TRX were most strongly expressed by
$z_{1}$ (note that the magnitudes of the depicted weights represent the
contribution of the individual currencies to the latent variables). (f) Past
values of driving signal $y_{2}$ can be seen to predict present values of
driven signal $z_{2}$, with the largest inflection of $y_{2}$ occurring a
short time before the corresponding dip in $z_{2}$. (g) Similar to the first
latent pair, the currencies ETC and BNB contributed most strongly to $y_{2}$,
but the driven signal here best expressed ADA and ETH. (h) The dynamics of the
third pair of latent variables exhibited a dynamic distinct from those of the
first two pairs, with the temporal precedence of $y_{3}$ relative to $z_{3}$
once again visible in the middle region of the curves. (i) The
cryptocurrencies BNB and XRP most strongly contributed to $y_{3}$, while
$z_{3}$ was best expressed by XLM and ETC.  Figure S1:  Convergence of grouped
coordinate descent on VAR(3) process data. Vertical axes depict the objective
function of the proposed approach, namely the strength of causality between
the driving and driven signal. Horizontal axes depict the iteration number in
logarithmic spacing. (a) The strength of causality between $y_{1}$ and
$z_{1}$, measured when holding the driving signal fixed and optimizing the
driven signal. Convergence occurs near iteration 10. (b) Same as (a) but now
measured when holding the driven signal fixed and optimizing the driving
signal. Approximately 20 iterations were required for convergence. (c) The
strength of causality at each iteration when employing the individually
optimized driving and driven signals. (d) Same as (a) but now for the second
pair of components. Convergence occurred after only 3 iterations. (e) Rapid
convergence was also found when optimizing the driving signal of the second
pair. (f) The strength of causality in the second pair of components at each
iteration. In all curves, shading represents the sem across $n=100$
realizations.  Figure S2:  Proposed technique recovers stronger causal
relationships than conventional techniques. EEG data from $n=12$ subjects was
employed to test the hypothesis that the component pairs recovered by the
proposed approach reflect stronger causal links than yielded by existing
approaches. The strength of causality was measured for the three most causally
connected pairs of electrodes, principal components, and independent
components (all selected _post hoc_). The resulting values were compared to
the strength of causality among the three pairs of recovered components
obtained with the proposed approach. A two-way ANOVA with method and component
as factors was conducted, yielding a significant main effect of method
($F(3)=11.53$, $p=9.23\times 10^{-7}$). There was no main effect of component
($p=0.59$), and no interaction ($p=0.98$). Follow-up tests showed that the
main effect of method was driven by significantly larger strengths of
causality with the proposed method ($\mathcal{G}=0.10\pm 0.031$, $0.088\pm
0.012$, and $0.081\pm 0.014$ for the first three components, means $\pm$ sem
across $n=12$ subjects) relative to the three most connected electrode pairs
($\mathcal{G}=0.049\pm 0.0084$, $0.049\pm 0.0072$, $0.049\pm 0.0079$;
$p=0.034$; $p=4.9\times 10^{-4}$, and $p=0.034$ for components 1, 2, and 3,
respectively; Wilcoxon signed rank test, $n=12$), the three most connected
principal component pairs ($\mathcal{G}=0.050\pm 0.0054$, $0.049\pm 0.0058$,
$0.046\pm 0.0046$; $p=0.034$, $p=0.0049$, and $p=0.0093$), and the three most
connected independent component pairs ($\mathcal{G}=0.033\pm 0.0056$,
$0.030\pm 0.0035$, $0.030\pm 0.0031$; $p=0.016$, $p=4.9\times 10^{-4}$, and
$p=0.0024$). This provides evidence that the proposed approach is able to
detect stronger causal relationships than what is found in the observed data
and conventional components analyses.