%% ------------------- SUPPLEMENTARY MATERIAL ------------------- %%
% Cardiorespiratory temporal causal links by sport (or lack thereof)

% Matlab code - how Granger-causality-related methods were carried out 
% Marcel M³yńczak, 2018

%% References

% [1] Barnett, L. and Seth, A. K. (2014). The MVGC multivariate Granger causality toolbox: 
%     a new approach to granger-causal inference. Journal of Neuroscience Methods 223, 50?68
% [2] Schiatti, L., Nollo, G., Rossato, G., and Faes, L. (2015). Extended granger causality: 
%     a new tool to identify the structure of physiological networks. Physiological Measurement 36, 827

%% Preparation

% mat1 - The first kind of data, sampled at 25Hz
% mat2 - The second kind of data, samplied at 2.5Hz
% mat3 - The third kind of data, beat-by-beat sequences

% Memory allocation for G-causalities 
% (nan, when statistically insignificant)
param = zeros(1,124)*nan;

%% 1. Method - Pairwise-conditioned Granger Causality based on MVGC framework [1]

% Setting coefficients
ntrials   = 1;              
nobs      = size(mat1,1);   
X         = mat1';
regmode   = 'OLS';  
icregmode = 'LWR'; 
acmaxlags = 5000;
momax     = 2;     
tstat     = '';    
alpha     = 0.05;   
mhtc      = 'FDR';  
nvars     = size(mat1,2);

% Calculating information criteria
[~,~,~,moBIC] = tsdata_to_infocrit(X,momax,icregmode);
morder = moBIC;

% Conversion from input data to VAR
[A,SIG] = tsdata_to_var(X,morder,regmode);
assert(~isbad(A),'VAR estimation failed');

% Conversion from VAR to autocovariance matrix
[G,~] = var_to_autocov(A,SIG,acmaxlags);

if(~isempty(G))
   
% Conversion from autocovariance matrix to G-causality   
F = autocov_to_pwcgc(G);

% P-value calculation
pval = mvgc_pval(F,morder,nobs,ntrials,1,1,nvars-2,tstat); 

% Confronting p-value with alpha
sig = significance(pval,alpha,mhtc);

% Filling the parameters vactor
if(sig(1,2)==1) 
   param(1) = F(1,2); 
end
if(sig(1,3)==1) 
   param(2) = F(1,3); 
end
if(sig(2,3)==1) 
   param(3) = F(2,3); 
end
if(sig(2,1)==1) 
   param(4) = F(2,1); 
end
if(sig(3,1)==1) 
   param(5) = F(3,1); 
end
if(sig(3,2)==1) 
   param(6) = F(3,2); 
end

end   

%% 2. Method - Spectral Granger Causality based on MVGC framework [1]

% Setting coefficients
nobs      = size(mat2,1);   
X         = mat2';
fs        = 2.5;    
fres      = [];     
nvars     = size(mat2,2);

% Calculating information criteria
[~,~,~,moBIC] = tsdata_to_infocrit(X,momax,icregmode);
morder = moBIC;

% Conversion from input data to VAR 
[A,SIG] = tsdata_to_var(X,morder,regmode);
assert(~isbad(A),'VAR estimation failed');

% Conversion from VAR to autocovariance matrix
[G,~] = var_to_autocov(A,SIG,acmaxlags);

if(~isempty(G))

% Conversion from autocovariance matrix to G-causality
f = autocov_to_spwcgc(G,fres);

% P-value calculation
pval = mvgc_pval(f,morder,nobs,ntrials,1,1,nvars-2,tstat); 

% Confronting p-value with alpha
sig  = significance(pval,alpha,mhtc);

% Filling the parameters vactor
[t1,t2] = max(f(1,2,:));
if(sig(1,2,t2)==1)
   param(7) = t1;
   param(8) = t2/size(sig,3)*fs/2;
end
[t1,t2] = max(f(1,3,:));
if(sig(1,3,t2)==1)
   param(9) = t1;
   param(10) = t2/size(sig,3)*fs/2;
end
[t1,t2] = max(f(2,3,:));
if(sig(2,3,t2)==1)
   param(11) = t1;
   param(12) = t2/size(sig,3)*fs/2;
end
[t1,t2] = max(f(2,1,:));
if(sig(2,1,t2)==1)
   param(13) = t1;
   param(14) = t2/size(sig,3)*fs/2;
end
[t1,t2] = max(f(3,1,:));
if(sig(3,1,t2)==1)
   param(15) = t1;
   param(16) = t2/size(sig,3)*fs/2;
end
[t1,t2] = max(f(3,2,:));
if(sig(3,2,t2)==1)
   param(17) = t1;
   param(18) = t2/size(sig,3)*fs/2;
end

end

%% 3. Entended Granger Causality based on the framework from [2]

% Setting coefficients
p = 4;
LRmethod=1;

% Granger causality estimation
[GC1,p_val1,Am1,Su1,SigmaR1,Res]=egc_gcMVAR(mat3',p);

% Filling parameters vector
for ip = 1:16
   if(p_val1(ip)<0.05) 
      param(18+ip) = GC1(ip);
   end
end
clear('ip')

% Extended Granger causality estimation
[EGC1,pval_e1,Bm1,B01,Sw1,SigmaR2,ZeroLag1,LRpruned,LR]=egc_gceMVAR(Y',p,Res,LRmethod);

% Filling parameters vector
for ip = 1:16
   if(pval_e1(ip)<0.05) 
      param(34+ip) = EGC1(ip);
   end
end
clear('ip')