## Abstract

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.

## 1 Introduction

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 *x _{t}* to

*y*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).

_{t}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 *x _{t}* 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 *Q _{s}* is positive semidefinite.

### 2.2 Steady-State Kalman Filter, Innovations State-Space Models, and the Discrete Algebraic Ricatti Equation

*z*. It is given by Kailath, Sayeed, and Hassibi (2000, theorem 9.2.1), where

_{t}*e*is the innovations sequence of variance , is the Kalman gain sequence, and

_{t}*P*is the state error variance matrix generated from the Ricatti equation, .

_{t}*P*will obey the so-called discrete algebraic Ricatti equation (DARE), where and is the corresponding steady-state Kalman gain. With some clever algebra, (Kailath et al., 2000), the DARE can be rewritten (the Ricatti equation can be similarly rewritten) as where and .

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.

*V*and Kalman gain

*K*. This steady-state filter provides a new state-space representation of the data sequence. We refer to it as an innovations state-space (ISS) model with parameters . We summarize this:

Given the SS model, equation 2.1, with parameters , then provided N,St,De hold:

See appendix A.

**Remarks.**

Henceforth an ISS model with parameters will be required to have

*V*positive definite, detectable, and controllable so that*A*−*KC*is stable.It is well known that any VARMA model can be represented as an ISS model and vice versa (Solo, 1986; Hannan & Deistler, 1988).

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).

**Examples.**

*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*z*is stationary._{t}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*z*can be perfectly predicted from the past (Hannan & Poskitt, 1988; Green, 1988), something that is not realistic in the fMRI application._{t}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.

## 3 Submodels

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 *x _{t}* from the ISS model for

*z*. We have simply where, . We need to calculate the covariance matrix: . We find , and . This leads to:

_{t}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:

## 4 Downsampling

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.

*z*with sampling multiple

_{t}*m*. Let

*t*denote the fine timescale and

*k*the coarse timescale so . The downsampled signal is . To develop the SS model for , we iterate the SS model above to obtain Now set and denote sampled signals, . Then we find where . We now use result

^{16}to find the ISS model corresponding to this SS model.

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 *A ^{m}*. 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 *m*th 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*y*does not weakly Granger-cause (dn-wgc)_{t}*x*if, for all_{t}*t*, Otherwise, we say*y*weakly Granger causes (wgc)_{t}*x*. Because of the elementary identity, the equality of variance matrices in the definition also ensures the equality of predictions: . 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._{t} - •
*Definition: Strong Granger causality*. Under WSS, we say*y*does not strongly Granger cause (dn-sgc)_{t}*x*if, for all_{t}*t*, Otherwise we say*y*strongly Granger-causes (sgc)_{t}*x*. Again, equality of the variance matrices ensures equality of predictions, . This definition agrees with Caines (1976) and Solo (1986)._{t} - •
*Definition: FBI (feedback interconnected).*If*x*Granger-causes_{t}*y*and_{t}*y*Granger-causes_{t}*x*, then we say are feedback interconnected._{t} - •
*Definition: UGC (unidirectionally Granger causes).*If*x*Granger-causes but_{t}*y*dn-gc_{t}*x*we say_{t}*x*unidirectionally Granger-causes ._{t}

### 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:

*y*dn-wgc_{t}*x*iff ._{t}*y*dn-sgc_{t}*x*iff and ._{t}

We can now state a new SS version of this result:

For the stationary ISS model, equations 5.1 and 5.2:

*y*dn-wgc_{t}*x*iff ._{t}*y*dn-sgc_{t}*x*iff and ._{t}

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

*y*to predict

*x*beyond using just the past of

*x*; similarly introduce . Thus, for example, is a directional measure of the influence of on

*x*. Next, define the instantaneous influence measure: . These are then joined in the fundamental decomposition (Geweke, 1982), where, . Geweke (1982) then proceeds to decompose these measures in the frequency domain. Thus, the frequency domain GEM for the dynamic influence of

_{t}*y*on

_{t}*x*is given by Geweke (1982), and is assembled (following Geweke, 1982) as follows.

_{t}^{3}, Clearly, with , .

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 *x _{t}* and

*y*. 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:

_{t}Now pulling all this together with the help of result ^{18}, we have an extension of the results of Geweke (1982) to the state-space/VARMA case.

For the joint ISS model, equation 5.1:

and .

*y*dn-wgc_{t}*x*iff, which holds iff , that is, iff ._{t}*y*dn-sgc_{t}*x*iff and that is, iff and , that is, iff i.e. iff ._{t}

A very nice nested hypothesis testing explanation of the decomposition 5.3 is given by Parzen in the discussion to Geweke (1982).

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: So to test for strong GC, we put these together: Together with similar asymptotics for , 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 *z _{t}* 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.

Result a has also been obtained in the independent work of Seth, Chorley, and Barnett (2013) by different methods.^{1}

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.

*B*. It is straightforward to check that this has spectral factorization, where and . Note that and . This factorization corresponds to the model where are independent white noises of variances , respectively. We see that the NMP filtering has introduced spurious dynamics and reversed the GCS.

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 *x _{t}*. The Wiener filter estimate of

*x*is where (Kailath et al., 2000) where is the autocovariance generating function of

_{t}*x*. 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

_{t}*k*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

_{x}*y*, we conclude that deconvolution will distort the GCS since the filters are noncausal.

_{t}## 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:

Forwards causality:

If

*y*dn-sgc_{t}*x*, then dn-sgc ._{t}If

*y*dn-wgc_{t}*x*, then in general wgc ._{t}

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.

We use the partitioned expressions in the discussion leading up to result ^{18}. We also refer to the discussion leading up to theorem ^{8}.

^{18}and the definition of dn-sgc, The other decomposition is , and Next we note from theorem

^{5}that for all . Thus, we deduce We can now write Based on equation 7.4, we now introduce the ISS model for , where is the innovations sequence. Using this, we introduce the estimator of and the estimation error . Below we show We thus rewrite the model for as Now we can construct an ISS model for where is the innovations sequence. In view of equations 7.1, 7.2, and 7.5, and are uncorrelated for all . Thus, we have constructed the joint ISS model: From this we deduce that dn-sgc as required.

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.

### 8.2 Computation

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 . |
---|---|---|---|---|---|---|---|---|---|---|

1.38 | 1.66 | 1.41 | 1.17 | 0.99 | 0.86 | 0.551 | 0.15 | 0.00 | 0.01 | |

0.20 | 0.25 | 0.29 | 0.31 | 0.32 | 0.32 | 0.29 | 0.11 | 0.00 | 0.01 |

m . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 10 . | 20 . | 30 . | 40 . |
---|---|---|---|---|---|---|---|---|---|---|

1.38 | 1.66 | 1.41 | 1.17 | 0.99 | 0.86 | 0.551 | 0.15 | 0.00 | 0.01 | |

0.20 | 0.25 | 0.29 | 0.31 | 0.32 | 0.32 | 0.29 | 0.11 | 0.00 | 0.01 |

#### 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 . |
---|---|---|---|---|---|---|---|---|---|---|

0.93 | 0.88 | 0.77 | 0.683 | 0.62 | 0.57 | 0.42 | 0.13 | 0.00 | 0.01 | |

1.05 | 1.82 | 2.01 | 1.80 | 1.53 | 1.30 | 0.75 | 0.18 | 0.00 | 0.02 |

m . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 10 . | 20 . | 30 . | 40 . |
---|---|---|---|---|---|---|---|---|---|---|

0.93 | 0.88 | 0.77 | 0.683 | 0.62 | 0.57 | 0.42 | 0.13 | 0.00 | 0.01 | |

1.05 | 1.82 | 2.01 | 1.80 | 1.53 | 1.30 | 0.75 | 0.18 | 0.00 | 0.02 |

#### 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 . |
---|---|---|---|---|---|---|---|---|---|---|

0.02 | 0.28 | 0.38 | 0.43 | 0.45 | 0.46 | 0.31 | 0.45 | 0.41 | 0.17 | |

2.93 | 2.38 | 1.62 | 1.24 | 1.02 | 0.86 | 0.37 | 0.53 | 0.45 | 0.18 |

m . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 10 . | 20 . | 30 . | 40 . |
---|---|---|---|---|---|---|---|---|---|---|

0.02 | 0.28 | 0.38 | 0.43 | 0.45 | 0.46 | 0.31 | 0.45 | 0.41 | 0.17 | |

2.93 | 2.38 | 1.62 | 1.24 | 1.02 | 0.86 | 0.37 | 0.53 | 0.45 | 0.18 |

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

*z*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 ._{t}

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 *z _{t}* 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 .

In Figure 1 we show plots of canonical HRFs (with default parameter choices) used in SPM and Glover (1999). They exhibit three main features found in practice:

- Refractory period: Characterized by a time to onset
Overshoot—characterized by a peak time

Undershoot

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).

### 9.2 Deconvolution

As noted in section 1, the current practice is to first apply a deconvolution to the BOLD signals before carrying out a GC analysis. As discussed in section 9.3, this makes sense only for DS-F.

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.

#### 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.

#### F-DS

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.

#### DS-F

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.

## 10 Conclusion

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:

Computations of submodels via the DARE (see theorems

^{1}and^{5})Computation of downsampled models via the DARE (see theorem

^{4})Reliable computation of GEMs via the DARE (see theorems

^{6}and^{7})Effect of filtering on GEMs (see theorem

^{8})

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.

### A.1 Stabilizability

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 *P _{t}* 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.

### A.3 Detectability

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 *all**pass* 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.

## Acknowledgments

This work was partially supported by NIH grant P41EB015896.

## References

## Note

^{1}

The main parts of this theorem were announced in a research workshop at the 2011 Human Brain Mapping meeting