% MVGC for the V1-V4 interaction
function [f,pval,sig,F,alpha]=granger(X,fmax,fs,momax,fres)
%ntrials = 10; % number of trials
%nobs = 80000; % number of observations per trial
%nvars = 2; % number of variables
F=0;pval=0;sig=0;alpha=0.05;
nobs=length(X(1,:,1));
ntrials=length(X(1,1,:));
nvars=2;
regmode = 'OLS'; % VAR model estimation regression mode ('OLS', 'LWR' or empty for default)
icregmode = 'LWR'; % information criteria regression mode ('OLS', 'LWR' or empty for default)
morder = 'BIC'; % model order to use ('actual', 'AIC', 'BIC' or supplied numerical value)
%momax = 60; % maximum model order for model order estimation (before, it was 20)
acmaxlags = []; % maximum autocovariance lags (empty for automatic calculation) (before, it was 1000)
tstat = ''; % statistical test for MVGC: 'F' for Granger's F-test (default) or 'chi2' for Geweke's chi2 test
alpha = 0.05; % significance level for significance test
mhtc = 'FDR'; % multiple hypothesis test correction (see routine 'significance')
%fs = 10000; % sample rate (Hz) (before, it was 200)
%fres = []; % frequency resolution (empty for automatic calculation)
%seed = 0; % random seed (0 for unseeded)
%--------------------------------------
%now X would be your variable data, with size X(5,1000,10) --> 5 variables, 1000 observations per trial, and 10 trials in total
%X=zeros(nvars,nobs,ntrials);
%----------------------------------------
% Calculate information criteria up to specified maximum model order.
ptic('\n*** tsdata_to_infocrit\n');
[AIC,BIC,moAIC,moBIC] = tsdata_to_infocrit(X,momax,icregmode);
ptoc('*** tsdata_to_infocrit took ');
% Plot information criteria.
% figure(88); clf;
% plot_tsdata([AIC BIC]',{'AIC','BIC'},1/fs);
% title('Model order estimation');
%amo = 10; % actual model order, we don't have it here
fprintf('\nbest model order (AIC) = %d\n',moAIC);
fprintf('best model order (BIC) = %d\n',moBIC);
%fprintf('actual model order = %d\n',amo);
% Select model order.
%if strcmpi(morder,'actual')
% morder = amo;
% fprintf('\nusing actual model order = %d\n',morder);
if strcmpi(morder,'AIC')
morder = moAIC;
fprintf('\nusing AIC best model order = %d\n',morder);
elseif strcmpi(morder,'BIC')
morder = moBIC;
fprintf('\nusing BIC best model order = %d\n',morder);
else
fprintf('\nusing specified model order = %d\n',morder);
end
%morder=25;
fprintf('\nin spite of the above, we use model order = %d\n',morder);
%----------------------------------------------
% Estimate VAR model of selected order from data.
ptic('\n*** tsdata_to_var... ');
[A,SIG] = tsdata_to_var(X,morder,regmode);
ptoc;
% Check for failed regression
assert(~isbad(A),'VAR estimation failed');
% NOTE: at this point we have a model and are finished with the data! - all
% subsequent calculations work from the estimated VAR parameters A and SIG.
%----------------------------------------------------
% The autocovariance sequence drives many Granger causality calculations (see
% next section). Now we calculate the autocovariance sequence G 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,info] = var_to_autocov(A,SIG,acmaxlags);
ptoc;
% The above routine does a LOT of error checking and issues useful diagnostics.
% If there are problems with your data (e.g. non-stationarity, colinearity,
% etc.) there's a good chance it'll show up at this point - and the diagnostics
% may supply useful information as to what went wrong. It is thus essential to
% report and check for errors here.
var_info(info,true); % report results (and bail out on error)
%-------------------------------------------------------
ptic('*** autocov_to_pwcgc... ');
F = autocov_to_pwcgc(G);
ptoc;
% Check for failed GC calculation
assert(~isbad(F,false),'GC calculation failed');
% Significance test using theoretical null distribution, adjusting for multiple
% hypotheses.
pval = mvgc_pval(F,morder,nobs,ntrials,1,1,nvars-2,tstat); % take careful note of arguments!
sig = significance(pval,alpha,mhtc);
%-------------------------------------------------------
% Calculate spectral pairwise-conditional causalities at given frequency
% resolution - again, this only requires the autocovariance sequence.
ptic('\n*** autocov_to_spwcgc... ');
f= autocov_to_spwcgc(G,fres);
ptoc;
% Check for failed spectral GC calculation
assert(~isbad(f,false),'spectral GC calculation failed');
% Plot spectral causal graph.
%figure(23); clf;
%plot_spw(f,fs,[0,fmax]);