MVGC bootstrap demo
Demonstrates confidence interval construction using a nonparametric bootstrap on generated VAR data for a 5-node network with known causal structure (see var5_test). Pairwise-conditional Granger causalities are estimated and confidence intervals constructed using both the theoretical and bootstrap distributions.
Contents
References
[1] L. Barnett and A. K. Seth, The MVGC Multivariate Granger Causality Toolbox: A New Approach to Granger-causal Inference, J. Neurosci. Methods 223, 2014 [ preprint ].
[2] D. A. Freedman, Bootstrapping regression models, Ann. Stats., 9(6), 1981.
Parameters
ntrials = 10; % number of trials nobs = 100; % number of observations per trial nsamps = 100; % number of bootstrap samples regmode = 'OLS'; % VAR model estimation regression mode ('OLS', 'LWR' or empty for default) acmaxlags = 1000; % maximum autocovariance lags (empty for automatic calculation) tstat = 'chi2'; % statistical test for MVGC: 'F' for Granger's F-test, 'chi2' for Geweke's chi2 test or leave empty for default alpha = 0.05; % significance level for all statistical tests seed = 0; % random seed (0 for unseeded)
Generate VAR data
% Seed random number generator. rng_seed(seed); % Get VAR coefficients for 5-node test network. AT = var5_test; nvars = size(AT,1); % number of variables % Residuals covariance matrix. SIGT = eye(nvars); % Generate VAR time series data with normally distributed residuals for % specified coefficients and covariance matrix. ptic('\n*** var_to_tsdata... '); X = var_to_tsdata(AT,SIGT,nobs,ntrials); ptoc;
VAR model estimation and autocovariance calculation
morder = size(AT,3); % actual model order; on real data - i.e. with no generative model % available - use information criteria to estimate (see 'mvgc_demo') % Calculate VAR model ptic('*** tsdata_to_var... '); [A,SIG] = tsdata_to_var(X,morder,regmode); ptoc; % Check for failed regression assert(~isbad(A),'VAR estimation failed'); % Now calculate autocovariance according to the VAR model, to as many lags % as it takes to decay to below the numerical tolerance level, or to acmaxlags % lags if specified (i.e. non-empty). ptic('*** var_to_autocov... '); [G,res] = var_to_autocov(A,SIG,acmaxlags); ptoc; % Report and check for errors. fprintf('\nVAR check:\n'); disp(res); % report results... assert(~res.error,'bad VAR'); % ...and bail out if there are errors
Granger causality estimation
% Time-domain pairwise conditional causalities ptic('*** autocov_to_pwcgc... '); F = autocov_to_pwcgc(G); ptoc; % Theoretical confidence intervals. [FTUP,FTLO] = mvgc_confint(alpha,F,morder,nobs,ntrials,1,1,nvars-2,tstat); % Critical GC value. FTCRIT = mvgc_cval(alpha,morder,nobs,ntrials,1,1,nvars-2,tstat);
Bootstrap
ptic('\n*** bootstrap_tsdata_to_pwcgc\n'); FSAMP = bootstrap_tsdata_to_pwcgc(X,morder,nsamps); ptoc('*** bootstrap_tsdata_to_pwcgc took ',[],1); % (We should really check for failed bootstrap estimates here.) % Bootstrap (empirical) confidence intervals. [FSUP,FSLO] = empirical_confint(alpha,FSAMP); % Note: we haven't calculated a bootstrap critical GC value; for this we would % require an (empirical) null GC distribution, which we could obtain by running % a permutation test (see e.g. 'permtest_tsdata_to_pwcgc').
Plot PWCGC estimates with theoretical and bootstrap confidence intervals
figure(1); clf subplot(2,1,1); plot_confints(F,FTUP,FTLO,FTCRIT); title(sprintf('Theoretical distribution\nconfidence intervals at alpha = %g',alpha)); subplot(2,1,2); plot_confints(F,FSUP,FSLO); title(sprintf('Bootstrap distribution\nconfidence intervals at alpha = %g',alpha));