The recent interest in the dynamics of networks and the advent, across a range of applications, of measuring modalities that operate on different temporal scales have put the spotlight on some significant gaps in the theory of multivariate time series. Fundamental to the description of network dynamics is the direction of interaction between nodes, accompanied by a measure of the strength of such interactions. Granger causality and its associated frequency domain strength measures (GEMs) (due to Geweke) provide a framework for the formulation and analysis of these issues. In pursuing this setup, three significant unresolved issues emerge. First, computing GEMs involves computing submodels of vector time series models, for which reliable methods do not exist. Second, the impact of filtering on GEMs has never been definitively established. Third, the impact of downsampling on GEMs has never been established. In this work, using state-space methods, we resolve all these issues and illustrate the results with some simulations. Our analysis is motivated by some problems in (fMRI) brain imaging, to which we apply it, but it is of general applicability.
Following the operational development of the notion of (temporal) causality by Granger (1969) and Sims (1972), Granger causality (GC) analysis has become an important part of time series and econometric testing and inference (Hamilton, 1994). It has also been applied in the biosciences (Kaminski & Blinowska, 1991; Bernasconi & Konig, 1999; Ding, Bressler, Yang, & Liang, 2000), climatology (global warming; Suna & Wang, 1996; Kaufmann & Stern, 1997; Triacca, 2005), and, most recently, functional magnetic resonance imaging (fMRI). This last application has stimulated this letter.
Since its introduction into fMRI (Goebel, 2003; Valdes-Sosa, 2004; Roebroeck, Formisano, & Goebel, 2005; Yamashita, Sadato, Okada, & Ozaki, 2005) GC has become the subject of an intense debate: (see Valdes-Sosa, Roebroeck, Daunizeau, & Friston, 2011; Roebroeck, Formisano, & Goebel, 2011) and associated commentary. There are two main issues in that debate that occur more widely in dynamic networks. First is the impact of downsampling on GC. In the fMRI neuroimaging application, causal processes may operate on a timescale on the order of tens of milliseconds, whereas the recorded signals are typically available only on a 1 second timescale, although sampling up to a 100 ms rate is now becoming possible (Feinberg & Yacoub, 2012). So it is natural to wonder if GC analysis on a slower timescale can reveal GC structure (GCS) on a faster timescale. Second is the impact of filtering on GC due to the hemodynamic response function, which relates neural activity to the recorded fMRI signal. Since intuitively GC will be sensitive to time delay, the variability of the hemodynamic response function (HRF), particularly spatially varying delay (called time to onset in the fMRI literature) and time to peak (confusingly sometimes called delay in the fMRI literature) has been suggested as a potential source of problems (Deshpande, Sathian, & Hu, 2010; Handwerker, Gonzalez-Castillo, D’Esposito, & Bandettini, 2012). A referee has, however, pointed out that following David et al. (2008), deconvolution (Glover, 1999) has gained favor (Havlicek, Friston, Jan, Brazdil, & Calhoun, 2011) as a potential means of dealing with this.
An important advance in GC theory and tools was made by Geweke (1982), who provided measures of the strength of causality or influence of one time series on another (henceforth called GEM, for Geweke causality measure), including frequency domain decompositions of them. Subsequently it was pointed out that the GEMs are measures of directional mutual information (Rissannen & Wax, 1987). The GEMs were extended to conditional causality in Geweke (1984). However, GEMs have not found as wide an application as they should have, partly because of some technical difficulties in calculating them, which are discussed further below. But GEMs (and their frequency domain versions) are precisely the tool needed to pursue both GC downsampling and filtering questions.
In the econometric literature, it was appreciated early that downsampling, especially in the presence of aggregation, could cause problems. This was implicit in the work of Sims (1971) and mentioned also in the work of Christiano and Eichenbaum (1987), who gave an example of contradictory causal analysis based on monthly versus quarterly data, and also discussed in Marcet (1991). But precise general conditions under which problems do and do not arise have never been given. We do so in this letter.
Some of the econometric discussion has been framed in terms of sampling of continuous time models (Sims, 1971; Marcet, 1991; Christiano & Eichenbaum, 1987). And authors such as Sims (1971) have suggested that models are best formulated initially in continuous time. While this is a view that I have long shared (and related to the modeling approach advocated by Middleton and Goodwin (1990), we deal with only discrete time models here. To cast our development in terms of continuous time models would require a considerable development of its own without changing our basic message.
The issue at stake, in its simplest form, is the following. Suppose that a pair of (possibly vector) processes exhibit a unidirectional GC relation (a influences b, but b does not influence a), but suppose measurement time series are available only at a slower timescale or as filtered time series or both. Then two questions arise. The first, which we call the forward question, is this: When is the unidirectional GC relation preserved? The second, which we call the reverse question, is harder. Suppose the measurement time series exhibit a unidirectional GC relation. Does that mean the underlying unfiltered faster timescale processes do? The latter question is the more important and so far has received no theoretical attention. The more general form of this issue is to what extent the GC structure (GCS) between two time series is distorted by downsampling and or filtering or both.
In order to resolve these issues, we need to develop some theory and some computational and modeling tools. First, to compute GEMs, one needs to be able to find submodels from a larger model (i.e., one having more time series). Thus, to compute the GEMs between time series , Geweke (1982, 1984) attempted to avoid this by fitting submodels separately to xt to yt and then also fitting a joint model to . Unfortunately this can generate negative values for some of the frequency domain GEMs (Chen, Bressler, & Ding, 2006). Properly computing submodels will resolve this problem, and previous work has not accomplished this (we discuss the attempts in Dufour & Taamouti, 2010, and Chen et al., 2006, below).
Second, one needs to be able to compute how models transform when they are downsampled. This has been done only in special cases (Pandit & Wu, 1983) or by methods that are not computationally realistic. We provide computationally reliable, state-space-based methods for doing this here.
Third, we need to study the effect of downsampling and filtering on GEMs. We will show how to use the newly developed tools to do this.
To sum up, we can say that previous discussions, including those above, as well as Geweke (1978), Telser (1967), Sims (1971), and Marcet (1991), fail to provide general algorithms for finding submodels or models induced by downsampling. Indeed both these problems have remained open problems in multivariate time series in their own right for several decades, and we resolve them here. Further, there does not seem to have been any theoretical discussion of the effect of filtering on GEMs, and we resolve that here also. To do that, it turns out that state-space models provide the proper framework.
Throughout this work, we deal with the dynamic interaction between two vector time series. It is well known in econometrics that if there is a third vector time series involved in the dynamics but not accounted for, then spurious causality can occur for reasons that have nothing to do with downsampling, a situation discussed by Hsiao (1982) and Geweke (1984). Other causes of spurious causality such as observation noise are also not discussed. Of course, the impact of downsampling in the presence of a third (vector) variable is also of interest but will be pursued elsewhere.
Finally, our whole discussion is carried out in the framework of linear time series models. It is of great interest to pursue nonlinear versions of these issues, but that would be a major separate task.
The remainder of the letter, is organized as follows. In section 2, we review and modify some state-space results important for system identification or model fitting and needed in the following sections. In section 3, we develop state-space methods for computing submodels of innovations state-space models. In section 4, we develop methods for transforming state-space models under downsampling. In section 5, we review GC and GEMs and extend them to a state space-setting. In section 6, we study the effect of filtering on GC via frequency domain GEMs. In section 7, we offer a theory to explain when GC is preserved under downsampling. In section 8, we discuss the reverse GC problem showing how spurious causality can be induced by downsampling. In section 9, we set up the framework for the application of these results to fMRI. This includes discussion of the phase properties of HRFs as well as the impact of deconvolution on the assessment of GC. Conclusions are in section 10. There are three appendixes.
1.1 Acronyms and Notation
GC is Granger causality or Granger causes; GCS is Granger causal structure; dn-gc does not Granger cause; GEM is Geweke causality measure. HRF is hemodynamic response function; SS is state space or state-space model; ISS is innovations state-space model; UC is unit circle; VAR is vector autoregression; VARMA is vector autoregressive moving average process; wp1 is with probability 1. LHS denotes left-hand side. The lag or backshift operator is denoted L or ; thus, if xt is a signal, then . L or are causal operators since they operate on the past. or z are noncausal operators since they operate on the future (so ). The causal filter is sometimes denoted h, and its noncausal reflection is .
denotes the values , so . For stationary processes, we have . If are positive semidefinite matrices, then means M−N is positive semidefinite. A stable square matrix is one whose eigenvalues all have modulus .
2 State Space
The computational methods we develop rely on state-space techniques and spectral factorization, the latter being intimately related to the steady-state Kalman filter. In this section, we review and modify some basic results in state-space theory, Kalman filtering, and spectral factorization. Our discussion deals with two vector time series, collected as, .
2.1 State-Space Models
It is common with SS models to take , but for equivalence between the class of VARMA models and the class of state-space models, it is necessary to allow .
Now by matrix partitioning, . So we introduce
N: Noise Condition. R is positive definite.
N implies that Qs is positive semidefinite.
2.2 Steady-State Kalman Filter, Innovations State-Space Models, and the Discrete Algebraic Ricatti Equation
We now introduce two assumptions.
St: Stabilizability condition: is stabilizable (see appendix A)
De: Detectability condition: is detectable.
In appendix A, it is shown De is equivalent to being detectable. Also it holds automatically if A is stable.
Given the SS model, equation 2.1, with parameters , then provided N,St,De hold:
See appendix A.
Henceforth an ISS model with parameters will be required to have V positive definite, detectable, and controllable so that A−KC is stable.
Note that the ISS model with parameters can also be written as the SS model with parameters .
The DARE is a quadratic matrix equation, but it can be computed using the (numerically reliable) DARE command in Matlab as follows. Compute: and then, .
Note that stationarity is not required for this result.
2.3 Stationarity and Spectral Factorization
Given an ISS model with parameters , we now introduce
Assumption EV: Eigenvalue stability condition: A has all eigenvalues with modulus , that is, A is a stability matrix.
With this assumption we can obtain an infinite vector moving average (VMA) representation, an infinite vector autoregressive (VAR) representation, and a spectral factorization.
Before stating these results, we need a definition;
Definition. Minimum phase transfer function. A matrix of causal transfer functions is said to be minimum phase (MP) if its inverse exists and is causal and stable. If is not minimum phase, it is called nonminimum phase (NMP).
Causal and stable. has inverse , which follows from summing the geometric series. Then is clearly causal and is stable since its pole has modulus (i.e., it lies inside the unit circle). Another way to state this is to say that is MP because its zero has modulus (i.e., it lies inside the unit circle, UC).
Causal but not stable. is not minimum phase since is not stable since its pole has modulus . If we rewrite , then . This is stable but not causal since the expansion is in terms of the future. Again we can rephrase this example to say is NMP because its zero is outside the unit circle.
Not causal but stable. is nonminimum phase since is not causal. Note that is stable. We can rephrase this as follows. has two zeros: one is inside the unit circle, and the other is at (i.e., makes vanish and so is outside the unit circle).
The following result is based on Kailath et al. (2000, theorem 8.3.2) and the surrounding discussion.
Write equation 2.3, in operator form. The series is convergent wp1 and in mean square since A is stable.
Rewrite equation as . Then write this in operator form. The series is convergent wp1 and in mean square since A−KC is stable and zt is stationary.
This follows from standard formulas for spectra of filtered stationary time series applied to (a).
From (a),(b) and by (b) is causal and stable, and the result follows.
Result 17 is a special case of a general result that given a full-rank multivariate spectrum , there exists a unique causal stable minimum phase spectral factor with and positive-definite innovations variance matrix such that equation 2.6 holds (Hannan & Deistler, 1988; Green, 1988). In general may have some roots on the unit circle (Hannan & Poskitt, 1988; Green, 1988), but the assumptions in result 17 rule this case out. Such roots mean that some linear combinations of zt can be perfectly predicted from the past (Hannan & Poskitt, 1988; Green, 1988), something that is not realistic in the fMRI application.
Result 17 is also crucial from a system identification or model-fitting point of view. From that point of view, all we can know (from second-order statistics) is the spectrum, and so if, naturally, we want a unique model, the only model we can obtain is the causal stable minimum phase model: the ISS model. The standard approach to SS model fitting is the so-called state-space subspace method (Deistler, Peternell, & Scherer, 1995; Bauer, 2005), and indeed it delivers an ISS model. The alternative approach of fitting a VARMA model (Hannan & Deistler, 1988; Lutkepohl, 1993) is equivalent to getting an ISS model.
We need result 16, however, since when we form submodels, we do not immediately get an ISS model; rather, we must compute it.
Our computation of causality measures requires that we compute induced submodels. In this section, we show how to obtain an ISS submodel from the ISS joint model.
We partition into subvector signals of interest and partition the state-space model correspondingly: and . We first read out an SS submodel for xt from the ISS model for zt. We have simply where, . We need to calculate the covariance matrix: . We find , and . This leads to:
First, we note by partitioning where so that and are both positive definite. Now we need only check conditions N,St,De of result 16. We need to show that is positive definite, is detectable, and is stabilizable; in fact, we show it is controllable.
The first is already established. The second follows trivially since A is stable. We use the PBH test (see appendix A) to check the third.
For implementation in Matlab, positive definiteness in constructing Q can be an issue. A simple resolution is to use a Cholesky factorization of and form and then form .
The matrix from DARE, .
Dufour and Taamouti (2010) discuss a method for obtaining submodels, but it is flawed. First, it requires computing the inverse of the VAR operator. While this might be feasible (analytically) on a toy example, there is no known numerically reliable way to do this in general (computation of determinants is notoriously ill conditioned). Second, it requires the solution of simultaneous quadratic autocovariance equations to determine VMA parameters for which no algorithm is given. In fact, these are precisely the equations required for a spectral factorization of a VMA process. There are reliable algorithms for doing this, but given the flaw already revealed, we need not discuss this approach any further.
Next we state an important corollary:
Any submodel is in general a VARMA model, not a VAR. To put it another way, the class of VARMA models is closed under the forming of submodels, whereas the class of VAR models is not.
This means that VAR models are not generic and is a strong argument against their use. Any vector time series can be regarded as a submodel of a larger-dimensional time series and thus must in general obey a VARMA model. This result (which is well known in time series folklore) is significant for econometrics where VAR models are in widespread use.
For the next section we need:
There are two approaches to the problem of finding the model obeyed by a downsampled process: frequency domain and time domain. While the general formula for the spectrum of a downsampled process has long been known, it is not straightforward to use and has not yielded any general computational approach to finding submodels of parameterized spectra.
The spectral formula, however, has made it very clear that downsampling leads to aliasing, so that higher-frequency information is lost. This strongly suggests that GCS will be affected by downsampling.
The most complete (time domain) work on downsampling seems to be that of Pandit and Wu (1983), although they treat only the first- and second-order scalar cases. There is work in the engineering literature for systems with observed inputs, but that is also limited and in any case not helpful here. We follow an SS route.
Using result 16, we need to show the following. R is positive definite, is detectable, and is stabilizable. In fact, we show controllability. The first holds trivially; the second does as well since A is stable and thus so is Am. For the third, we use the PBH test.
Suppose controllability fails. Then there is a left eigenvector q (possibly complex) with and . Since is positive definite, this delivers for .
Using this, we now find . Iterating this yields . Thus, if is an mth root of , then . Since , we thus conclude that is not controllable. But this is a contradiction since is stable.
In Matlab, we would compute, , yielding and .
More specifically () obeys , where ; .
5 Granger Causality
In this section, we review and extend some basic results in Granger causality. In particular, we extend GEMs to the state-space setting and show how to compute them reliably.
Since the development of Granger causality, it has become clear (Dufour & Renault, 1998; Dufour & Taamouti, 2010) that in general, one cannot address the causality issue with only one-step-ahead measures as commonly used. One needs to look at causality over all forecast horizons. However one-step measures are sufficient when one is considering only two vector time series as we are (Dufour & Renault, 1998, proposition 2.3).
5.1 Granger Causality Definitions
Our definitions of one-step Granger causality naturally draw on Granger (1963, 1969), Sims (1972), and Solo (1986) but are also influenced by Caines (1976), who, drawing on the work of Pierce and Haugh (1977), distinguished between weak and strong GC, or what Caines calls weak and strong feedback free processes. We introduce:
Condition WSS: The vector time series are jointly second-order stationary.
- •Definition: Weak Granger causality. Under WSS, we say yt does not weakly Granger-cause (dn-wgc) xt if, for all t, . This definition agrees with Granger (1969); Caines and Chan (1975), who do not use the designator weak; and Caines (1976) and Solo (1986) who do.
- •Definition: Strong Granger causality. Under WSS, we say yt does not strongly Granger cause (dn-sgc) xt if, for all t, . This definition agrees with Caines (1976) and Solo (1986).
Definition: FBI (feedback interconnected). If xt Granger-causes yt and yt Granger-causes xt, then we say are feedback interconnected.
Definition: UGC (unidirectionally Granger causes). If xt Granger-causes but yt dn-gc xt we say xt unidirectionally Granger-causes .
5.2 Granger Causality for Stationary State-Space Models
If obeys a Wold model of the form where is a one-sided square summable moving average polynomial with , which is partitioned as in equation 5.2, then:
yt dn-wgc xt iff .
yt dn-sgc xt iff and .
We can now state a new SS version of this result:
For the stationary ISS model, equations 5.1 and 5.2:
yt dn-wgc xt iff .
yt dn-sgc xt iff and .
The proof follows immediately from result 18 since .
By the Cayley-Hamilton theorem, we can replace part a of theorem 5 with .
Collecting these equations together gives , which says that the pair is not controllable. Also, we have , which says that the pair is not observable. Thus, the representation of is not minimal.
From a data analysis point of view, we need to embed this result in a well-behaved hypothesis test. The results of Geweke (1982), suitably modified, allow us to do this.
5.3 Geweke Causality Measures for SS Models
Introduce the normalized cross-covariance-based matrix, . Then, using a well-known partitioned matrix determinant formula (Magnus & Neudecker, 1999), we find . This means that the instantaneous causality measure depends only on the canonical correlations (which are the eigenvalues of ) between (Seber, 1984; Kshirsagar, 1972).
To implement these formulas, we need expressions for . To get them, Geweke (1982) fits separate models to each of xt and yt. But this causes positivity problems with (Chen et al., 2006). Instead we obtain the required quantities from the correct submodel obtained in the previous section. We have:
For the joint ISS model, equation 5.1:
yt dn-wgc xt iff, which holds iff , that is, iff .
yt dn-sgc xt iff and that is, iff and , that is, iff i.e. iff .
It is straightforward to see that the GEMs are unaffected by scaling of the variables. This is a problem for other GC measures (Edin, 2010).
- For completeness, we state extensions of the inferential results in Geweke (1982) without proof. Suppose we fit an SS model to data using, for example, so-called state-space subspace methods (Deistler et al., 1995; Bauer, 2005) or VARMA methods in, for example, Lutkepohl (1993). Let be the corresponding GEM estimators. If we denote true values with a superscript 0, we find under some regularity conditions: , we see that the fundamental decomposition, equation 5.3, has a sample version involving a decomposition of a chi-squared into sums of smaller chi-squared statistics.
Chen et al. (2006) also attempt to derive without fitting separate models to . However the proposed procedure to compute involves a two-sided filter and is thus in error. The only way to get is by spectral factorization (which produces one-sided or causal filters), as we have done.
Other kinds of causality measures have emerged in the literature (Kaminski & Blinowska, 1991; Baccala & Sameshima, 2001; Takahashi, Baccala, & Sameshima, 2008; Pascual-Marqui et al., 2014), but it is not known whether they obey the properties in theorems 6 and 7. However these properties are crucial to our subsequent analysis.
6 Effect of Filtering on Granger Causality Measures
Now the import of the frequency domain GEM becomes apparent since it allows us to determine the effect of one-sided (or causal) filtering on GC.
We need to be clear on the situation envisaged here. The unfiltered time series are the underlying series of interest, but we have access only to the filtered time series, so we can only find the GEMs from the spectrum of the filtered time series. What we need to know is when those filtered GEMs are the same as the underlying GEMs. We have:
Suppose we filter zt with a stable, full-rank, one-sided filter . Then:
If is minimum phase, the GEMs (and so GC) are unaffected by filtering.
If has the form where is a scalar all-pass filter and is stable, minimum phase, then the GEMs (and so GC) are unaffected by filtering.
If is nonminimum phase and case b does not hold, then the GEMs (and so GC) are changed by filtering.
We give three examples:
Differential delay. Suppose , where are independent zero mean white noises with variances , respectively, while . So the two series are white noises that exhibit an instantaneous GC. The filtering delays one series relative to the other. Then we have , and we see that while is stable, causal, and invertible; indeed, . Thus, we see that the differential delay has introduced a spurious dynamic GC relation and the original purely instantaneous GC is lost.
Noncausal filter. We use the same setup as in example 6 A except that now where and . We rewrite this as where . Although is noncausal and so NMP, it is the same for both filters and so by theorem 6b will not affect the GEMS. So we can replace with . Both entries in are causal but unstable, and so NMP and the GCS will be altered. In fact, we now show this explicitly.
Wiener deconvolution. When filtered versions of are observed in noise, one can estimate them by using the Wiener fitler (Kailath et al., 2000). The optimal filter uses the joint spectrum, but fMRI practice (Glover, 1999) is to filter each signal separately. This is suboptimal but will still reduce the effect of the noise. Suppose we observe where is a stable filter and is a white noise of variance independent of xt. The Wiener filter estimate of xt is where (Kailath et al., 2000) where is the autocovariance generating function of xt. In fMRI practice (Glover, 1999), is set to 1. This is suboptimal and may lead to performance worse than no filtering. But continuing, we carry out a spectral factorization giving where is a causal minimum phase filter with . The filter thus becomes . Since kx is minimum phase and causal, then is stable and causal and will not affect the GEMS. We are left with the noncausal filter — for example, in the simple case, , the spectral factorization reduces to solving a quadratic equation and one obtains . Thus, . With a similar deconvolution (where ) applied to a noisy filtered version of yt, we conclude that deconvolution will distort the GCS since the filters are noncausal.
7 Downsampling and Forwards Granger Causality
We now consider to what extent GC is preserved under downsampling.
Using the sampled notation of our discussion above and defining , we have the following result:
If yt dn-sgc xt, then dn-sgc .
If yt dn-wgc xt, then in general wgc .
Part a is new, although technically a special case of a result of the authors established in a non-SS framework.
We might consider taking part b as a formalization of longstanding folklore in econometrics (Christiano & Eichenbaum, 1987; Marcet, 1991) that downsampling can destroy unidirectional Granger causality. However, that same folklore is flawed because it failed to recognize the possibility of part a. The folklore is further flawed because it failed to recognize the more serious reverse problem discussed below.
Consider then . The second term vanishes for . The first term vanishes for since is uncorrelated with the past and hence ; for , it vanishes since is uncorrelated with the past. For , it vanishes by the definition of (Kailath et al., 2000).
A perusal of the proof of theorem 12, part a, shows that we cannot construct the block lower triangular joint ISS model. In general, we obtain a full block ISS model.
8 Downsampling and Reverse Granger Causality
We now come to the more serious issue of how downsampling might distort GC. To establish distortion, we simply have to exhibit a numerical example, but that is not as simple as one might hope.
8.1 Simulation Design
Designing a procedure to generate a wide class of examples of spurious causality is not as simple as one might hope. We develop such a procedure for a bivariate vector autoregression of order one: a bivariate VAR(1). On the one hand, this is about the simplest example one can consider; on the other hand, it is general enough to generate important behaviors.
We note that this model can be written as an ISS model with parameters, . Hence, all the computations described above are easily carried out.
But the real issue is how to select the parameters. By a straightforward scaling argument, it is easy to see that we may set without loss of generality. Thus, we need to choose only .
Some reflection shows that there are two issues. First, we must ensure the process is stationary; for the eigenvalues of A, we must have . Second, to design a simulation, we need to choose , but these quantities depend on the parameters in a highly nonlinear way, so it is not obvious how to do this. And five parameters are already too many to pursue this by trial and error.
For the first issue, we have and . Our approach is to select and then find to satisfy . This requires solving a quadratic equation. If we denote the solutions as , then we get two cases: and . This leaves us to select .
In appendix B we show that where . Similarly, where . But we also show that and , so we select , thereby setting a lower bounds to . This seems to be the best one can do, and as we see below, it works quite well. So given , compute and . This gives four cases and multiplied by the two cases above yields eight cases.
This is not quite the end of the story since the values need to be consistent with the values. Specifically the quadratic equation to be solved for must have real roots. Thus, the discriminant must be . So . There are four cases: two with real roots and two with complex roots.
If are real, then we require . This always holds if . If , we have a binding constraint that restricts the sizes of .
If are complex conjugates, then is negative. If , then the condition never holds. If , then there is a binding constraint that restricts the sizes of . In particular, note that if , then one cannot have complex roots for A. We now use this design procedure to illustrate reverse causality.
We describe the steps used to generate the results below. We assume the state-space model for comes in ISS form. Since standard state-space subspace model-fitting algorithms (Larimore, 1983; van Overschee & de Moor, 1996; Bauer, 2005) generate ISS models, this is a reasonable assumption. Otherwise we use result 16 to generate the corresponding ISS model.
8.3 Scenario Studies
We now illustrate the various results above with some bivariate simulations.
Example 1: GEMs Decline Gracefully
Here, for the underlying process, y pushes x much harder than x pushes y. This pattern is roughly preserved with slower sampling, but the relative strengths change as seen in Table 1.
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
Example 2: GEMs Reverse
In this case, the underlying processes push each other with roughly equal strength. But subsampling yields a false picture with x pushing y much harder than the reverse (see Table 2).
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
Example 3: Near Equal Strength Dynamics Becomes Nearly Unidirectional
In this case, the underlying relation is one of near equal strength feedback interconnection. But as Table 3 shows, almost immediately a very unequal relation appears under subsampling which soon decays to a near unidirectional relation.
Example 4: Near Unidirectional Dynamics Becomes Near Equal Strength
In this case a near unidirectional dynamic relation immediately becomes one of significant but unequal strengths and then one of near equal strength (see Table 4).
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
|m .||1 .||2 .||3 .||4 .||5 .||6 .||10 .||20 .||30 .||40 .|
There is nothing pathological about these examples, and using the design procedure developed above, it is easy to generate other similar kinds of examples. They make it emphatically clear that GC cannot be reliably discerned from downsampled data.
In independent work, Seth et al. (2013) have shown simulations where downsampling distorts GC. But their work is purely empirical, and they do not have the theoretical framework developed here to design specific types of distortion as we have done.
9 Application to fMRI
In this section we apply the above results to fMRI, but to do so, we need to describe fMRI data generation schemes (dgs’s); we give three such schemes The following feature is common to all of them.
We suppose there is an underlying neural activity (NA) 2-node circuit operating on a fast timescale of order 50 ms and that it is described by a state-space model. Specifically, the (possibly vector) NA signals at each node are , which we collect together into a joint NA vector signal , which is described by a state-space model such as equation 2.1 or 2.3.
Filter then downsample (F-DS). The BOLD signal recorded at a particular voxel is generated in two stages. The first stage is the same as F. In the second stage, the scanner samples this signal delivering a downsampled (or slow) BOLD signal — for example, if the NA timescale is 50 ms and the downsampled time scale is 200 ms, then .
Downsample then filter (DS-F). Here we suppose the NA zt is downsampled to a slow timescale giving a slow NA . This is then filtered by (a slow spatially varying) HRF to generate the slow BOLD signal .
Note that F is the limiting case of F-DS where fMRI is recorded on the same timescale as NA. Clearly F-DS makes physiological sense, whereas DS-F does not. Indeed F-DS is advocated by Deshpande et al. (2010), who provide further physiological references to back it up. However, below we will see that DS-F seems to be behind some of the processing methods used in practice.
The aim now is to try to solve the reverse GC problem; to determine the GCS of fast NA circuit zt from the GCS of either the fast BOLD circuit or the slow BOLD circuit . To do this, we see from theorem 8, that we need to investigate the phase properties of the HRF.
9.1 HRFs Are Nonminimum Phase
We denote an HRF by and its transfer function as . We assume is bounded input, bounded output stable (Kailath et al., 2000), for which the necessary and sufficient condition is .
Overshoot—characterized by a peak time
Also there are persistent reports in the literature of HRFs with an initial small dip immediately following the refractory period (e.g., Miezin, Maccotta, Ollinger, Petersen, & Buckner, 2000).
A number of results in the control engineering literature show how nonminimum phase (NMP) zeros of the transfer function imply overshoot and undershoot of an impulse response (see Damm & Muhirwa, 2014). But none of these works provide converse results, showing what impulse response shape implies about transfer function zeros.
We give a preliminary general result and then two converse results.
Any HRF with a refractory period is NMP.
We can write the transfer function as where is a causal transfer function. Even if is MP, then is NMP because is not a causal operator.
Certainly the time to onset can be estimated, and then the BOLD signal can be shifted in time to align the start time with the stimulus start time. This effectively corrects the HRF for the refractory period. But this makes sense only if it is known that the refractory period is due to information transmission delay. If it is due to the GCS, then such alignment, being noncausal, will destroy the GCS.
To continue, we now assume the HRF has been corrected for the refractory period so that the corrected transfer function is .
Dip: if and , then the HRF is NMP.
Note that −area above the dip + area under the overshoot − area above the undershoot. The fact that the overshoot area (far) exceeds the sum of the other two areas (so that ) seems to be a universal property of actual HRFs.
We prove that has a real zero outside the UC. We can set , which is real. Then consider that . For , the second term is bounded by and as . Since , then for R large enough, . Since , this means there must exist with as required.
No dip: if and , then the HRF is NMP.
Note that and is satisfied by the default canonical HRFs. Further, because HRFs are very smooth, the successive differences will be very small. can be split into four pieces: a positive sum up to the peak time, a negative sum from there to the zero crossing, a positive sum to the trough of the undershoot, and a negative sum from there back up to the the time axis. We thus expect to be small.
We prove has a complex zero outside the UC. Put . Then consider that . For , the second term is bounded by , which as and so since , then for R large enough, . Since , there must thus exist for which and the result is established.
Since the condition and more generally the possibility of complex zeros outside the UC is so crucial, more flexible modeling of the HRF is called for. Recent approaches include Khalidov, Fadili, Lazeyras, Ville, and Unser (2011), Zafer, Blu, and Ville (2014), Wu et al. (2013), Sreenivasan, Havlicek, and Deshpande (2015) and our approach using Laguerre-polynomials (Cassidy, Long, Rae, & Solo, 2012).
We now see from the deconvolution example 6 C that even if Wiener deconvolution is applied to noise-free data , it will distort the GCS because it necessarily involves noncausal filtering. This is discussed further in section 9.3.
9.3 Reverse GC Problem
We discuss approaches under F, F-DS, and DS-F.
There is only one approach possible here. If the HRF is minimum phase, then the fast filtering will not affect the GCS and the first stage would consist of a GC analysis of the fast BOLD data. That is, fit an SS model to ; find submodels for using DARE (see theorem 1); and compute the GEMS using the expressions in theorem 6. This would deliver the GCS of the fast BOLD circuit and, hence, of the NA circuit. But in practice, the BOLD data will be noisy, and an attempt to reduce the impact of the noise by deconvolution will distort the GCS. If the HRF is nonminimum phase, then this GCS would be distorted and nothing can be done to overcome that.
Again there is only one approach possible here. If the HRF is minimum phase, then the fast filtering will not affect the GCS and the first stage would consist of a GC analysis of the slow BOLD data—that is, fit an SS model to ; find submodels for using DARE (see theorem 1); and compute the GEMS using the expressions in theorem 6. This would deliver the GCS of the slow BOLD circuit. The second stage is to try to get the GCS of the fast BOLD circuit. But now we see that it must fail since the results of the reverse GC study in section 8 demonstrate that we cannot use the GCS of the slow BOLD circuit to discern the GCS of the fast BOLD circuit. So we cannot get at the GCS of the fast NA circuit.
Notice that in F-DS there is no role for deconvolution.
If the HRF, is NMP then the GCS will be distorted and there is nothing we can do to overcome this. If we apply deconvolution, then the analysis above shows it must fail since the Wiener filter is NMP.
There are two approaches here. The direct approach has two stages. If the HRF is minimum phase, then the slow filtering will not affect the GCS and the first stage would consist of a GC analysis of the slow BOLD data—that is, fit an SS model to ; find submodels for using DARE (see theorem 1); and compute the GEMS using the expressions in theorem 6. This would deliver the GCS of the slow NA circuit since by theorem 8, the GEMS would be the same as for the slow BOLD data. The second stage is to try to get the GCS of the fast NA circuit. But now we see that it must fail, since the results of the reverse GC study in section 8 demonstrate that we cannot use the GCS of the slow NA circuit to discern the GCS of the fast NA circuit.
The indirect approach has three stages. If the slow HRF is minimum phase, then in the zeroth stage, we could inverse filter (deconvolve) the slow BOLD data to obtain the slow NA. In the first stage, we obtain its GCS using the procedure described in the first stage of the direct approach. But for the second stage, we are in the same situation as in the direct approach and must fail to discern the GCS of the fast NA.
In practice, slow timescale deconvolution has gained favor recently. It is done using a Wiener filter (Glover, 1999) since this deals with possible noise (which we have ignored in our discussion). But as discussed above, this must fail since the Wiener filter is NMP.
We see here that only in DS-F is there a role for slow timescale deconvolution.
In sum we conclude that under any data-generating scheme, GC analysis must fail.
This letter has given a theoretical and computational analysis of the use of Granger causality (GC) and applied the results to fMRI. There were two main issues: the effect of downsampling and the effect of hemodynamic filtering on Granger causal structure (GCS). To deal with these issues, a number of novel results in multivariate time series and Granger causality were developed via state-space methods as follows:
Theorem 8 showed that nonminimum phase filters distort GCS. Using the results in the first three methods, we were able to develop, in Section 8, a framework for generating downsampling-induced spurious Granger causality on demand and provided a number of illustrations.
In section 9 these results were applied to fMRI. To do this we described three data generating processes for fMRI:
F where fMRI time series are filtered neural activity (NA)
F-DS where NA is filtered, then downsampled to give fMRI time series
DS-F where NA is downsampled, then fitered to give fMRI time series
While F-DS has strong physiological justification, DS-F seems to be behind a number of approaches in the literature. We showed that HRFs are nonminimum phase due to a refractory period. If that period can be corrected for, the phase character of the HRF has to be checked on a case-by-case basis. Certainly default canonical HRFs are nonminimum phase, and so more flexible HRF modeling will be needed.
Then we discussed the reverse GC problem under each data-generating process. In each case, we concluded that GC analysis could not work.
All this leads to the conclusion that linear Granger causality analysis of fMRI data cannot be used to discern neuronal level–driving relationships. Not only is the timescale too slow, but even with faster sampling, the nonminimum phase aspect of the HRF and the noise will still compromise the method.
Future work would naturally include an extension of the Granger causality results to handle the presence of a third vector time series and extensions to deal with time-varying Granger causality. Nongaussian versions could mitigate the nonminimum phase problem to some extent, but there does not seem to be any evidence for the nongaussianity of fMRI data. Extensions to nonlinear Granger causality are currently of great interest but need considerable development.
Appendix A: Stabilizability, Detectability, and DARE
In this appendix we restate and modify for our purposes some standard state-space results. We rely mostly on Kailath et al. (2000, appendixes E and C).
We denote an eigenvalue of a matrix by and a corresponding eigenvector by q. We say is a stable eigenvalue if ; otherwise, is an unstable eigenvalue.
The pair is controllable if there exists a matrix G so that A−BG is stable (i.e., all eigenvalues of A−BG are stable). is controllable iff any of the following conditons hold,
Controllability matrix: has rank n.
Rank test: for all eigenvalues of A.
PBH test: There is no left eigenvector of A that is orthogonal to B, that is, if , then .
The pair is stabilizable if for all unstable eigenvalues of A. Three useful tests for stabilizability are:
PBH test: is stabilizable iff there is no left eigenvector of A corresponding to an unstable eigenvalue that is orthogonal to B, that is, if and , then .
is stabilizable if is controllable.
is stabilizable if A is stable.
A.2 Theorem DARE
See Kailath et al. (2000, theorem E6.1, lemma 14.2.1, section 14.7).
Under conditions N,St,De, the DARE has a unique positive semidefinite solution P, which is stabilizing (i.e., is a stable matrix). Further if we initialize , then Pt is nondecreasing and as .
“ stable” means that A−KC is stable (see below), and this implies that is controllable (see below).
Since , then is positive definite.
We first note (taking limits in) (Kailath et al., 2000, equation 9.5.12) . We then have . So is stable iff A−KC is stable. But then is controllable.
The pair is detectable if is stabilizable.
Suppose is detectable but is not. Then by the PBH test, there is a right eigenvector p of A corresponding to an unstable eigenvalue of A with . But then while which contradicts the detectbility of . The reverse argument is much the same.
Part a follows from the discussion leading to theorem DARE. Part b follows from the remarks after theorem DARE.
Appendix B: GEMs for Bivariate VAR(1)
Appendix C: Spectral Factorization
Suppose is a white noise sequence with . Let be a stable causal, possibly nonminimum phase filter. Then has spectrum where . We can then find a unique causal, stable minimum phase spectral factorization . Let have Cholesky factorizations and set . Then . Since is minimum phase, we can introduce the causal filter . Such a filter is called an allpass filter (Hannan & Deistler, 1988; Green, 1988). Now, or (i.e., a decomposition of a nonminimum phase (matrix) filter into a product of a minimum phase filter and an all-pass filter. We can also write this as showing how the nonminimum phase filter is transformed to yield a spectral factor.
This work was partially supported by NIH grant P41EB015896.
The main parts of this theorem were announced in a research workshop at the 2011 Human Brain Mapping meeting