Granger causality

Granger causality is the method of inferring directional interaction between time series, which in the context of neuroimaging can be derived from statistical parametric mapping of fMRI/MEG/EEG measurements. Most Granger causality values are calculated using auto-regressive (AR) modeling of the time series. Here are some examples in Matlab.

Codes

The multivariate auto-regressive modeling toolbox (ARFIT) is very useful and has been extensively applied in many studies. Our codes are built on top of this toolbox.

Please download all codes in this page here.
Github repository

The function etc_granger.m calculates the Granger causality values and the associated statistics

Examples

Simulating bi-variate AR model time series
Date: Feb-23-2010
Purpose: Generate bi-variate time series and calculate the Granger causality
Description: This is the example of Granger causality values from bi-variate time series generated from a 1-st order AR model in Roebroeck et al. (NeuroImage 2005) *without* convolving the time series with hemodynamic response function. Mathematically, it is

[x(t+1), y(t+1)]^T=[-0.9, 0; I, -0.9][x(t), y(t)]^T+[ex(t), ey(t)], where [ex(t) ey(t)]^T are bi-variate noise with a covariance matrix S=[1 0; 0 1];

I is the coupling strength and it is set I=0.5 for strong influence. The model clearly suggests a causal modulation from x -> y.

Matlab code: granger_ar_sim.m

Simulating the bi-variate AR model time series
Date: Feb-23-2010
Purpose:Generate dynamic bi-variate time series and calculate the dynamic Granger causality
Description: This is the example of Granger causality values from bi-variate time series generated from a 2nd order AR model. Mathematically, it is

[x(t+2), y(t+2)]^T=[-0.9, 0; I, -0.9][x(t+1), y(t+1)]^T+[-0.3, 0; I*0.5, -0.3][x(t), y(t)]^T+[ex(t), ey(t)], where [ex(t) ey(t)]^T are bi-variate noise with a covariance matrix S=[1 0; 0 1];

I is the coupling strength and it is set I=0.5 for strong influence. Such bi-variate time series is only active between time sample 300 and 600 among total 1000 samples. Other time instants are two white noises. Specifically, the final time series [X(t), Y(t)]^T is [X(t), Y(t)]^T = env(t).*[x(t), n(t)]^T + (1-env(t)).*[xn(t), yn(t)]^T, where env(t) is a time course with value close to 1 between sample 300 and 600 (sigmoidal ramping) and 0 else where. [xn(t) yn(t)]^T are white Gaussian noise with unit variance.

The goal here is to estimate the appropriate window length and AR model order using the SURE metric, suggested by Lin et al (Human Brain Mapp 2009). Subsequently, dynamic Granger causality is estimated using the chosen AR model order and window length.

Matlab code: dgranger_ar_sim.m