%Matlab routines to evaluate Granger Causality results on simulated data 

%Routine to simulate data while varying the lags between causally interacting signals. 
load S.mat     %load of the signal S.X, a matrix of 20 rows (20 signals) with 12800 columns (signal values for 12800 time points) and the value of the sampling rate S.fc  
S.X=Rumore_sasq(S.X,1); % add the noise to the signal, see below S.X=Rumore_sasq(S.X,1); 
rit(1)=-9; %-9 is the starting value of the lag of the designed structure
Amp=0.0; % no connection strength factor was added to the connection strength factor of the starting structure
incr=1;
for k=1:40 
rit(k+1)=rit(k)+incr; %the adding value for the lags ranged from -9 to 30 time points with a step of 1
X=S.X(:,1:5120); % each signal of 12800 time pints was cut to 5120 time points
X=Fai_scambi_e_ritardi(X,rit(k),Amp); %see the function below which defines the designed structure for the first and second example
X=X(1:19,:); % the 20th latent signal was not included 
% code to perform Granger causality, see Seth, 2010   
end

%Routine to simulate data while varying the connection strength between causally interacting signals
load S.mat
  S.X=Rumore_sasq(S.X,1); 
rit=0; %no value was added to the starting values of the lags between interacting signals 
Amp(1)=-0.2; %starting value of the connection strength
incr=0.01
for k=1:40
Amp(k+1)=Amp(k)+incr; %a value which ranged from -0.2 to 0.2 with a step of 0.01 was added to each connection strength factor of all the causal links 
X=S.X(:,1:5120);
X=Fai_scambi_e_ritardi(X,rit,Amp(k));
X=X(1:19,:); % the 20th latent signal was not included 
% code to perform Granger causality, see Seth, 2010   
end

%Routine to simulate data while varying the noise
load S.mat
S.X=Rumore_sasq(S.X,0.6); 
noise(1)=0.00; %no noise was added to the starting noise in the signals 
Amp=0.0;
rit=0;
incr=0.02;
for k=1:40
noise(k+1)=noise(k)+incr; % the added noise ranged from 0 to 0.78 time the standard deviation of the amplitude of each signal
X=S.X(:,1:5120);
X=Fai_scambi_e_ritardi(X,rit,Amp);
X=X(1:19,:); 
X=Rumore_sasq(X,noise(k));
% code to perform Granger causality, see Seth, 2010   
end

%Routine to simulate data while varying the number of time points in the temporal window 
load S.mat
 S.X=Rumore_sasq(S.X,1); 
Amp=0.0; 
incr=256;
lun(1)=2*incr; %the lenght of the temporal window increased of 256 time points for each loop
rit=0;
 for k=1:length(S.X)/256-10;
lun(k+1)=lun(k)+incr; % the length of the temporal window ranged from 512 to 10496 time points
X=S.X(:,1:lun(k));
X=Fai_scambi_e_ritardi(X,rit,Amp);
X=X(1:19,:); 
% code to perform Granger causality, see Seth, 2010   
end

function X=Rumore_sasq(X,rum)
    X=X';
    for ch=1:size(X,2)
X(:,ch) = X(:,ch)+rum*std(X(:,ch))*rand(1);
    end
  X=X';    
  
  
  function  X=Ritardi_sas(X,ch_p,ch_a,lag,Amp)
[ri,co]=size(X);
for k=1:co-lag
X(ch_a,k+lag)=X(ch_a,k+lag)+Amp*X(ch_p,k);  
end

  % First example of the simulation
    function X=Fai_scambi_e_ritardi(X,rit,Amp)
X=Ritardi_sas(X,2,4,rit+10,Amp+0.4);  %from 2 to 4       lag=10 connection strength=0.4     
X=Ritardi_sas(X,8,6,rit+14,Amp+0.3);  %from 8 to 6       lag=14 connection strength=0.3
X=Ritardi_sas(X,11,18,rit+14,Amp+0.3);    %from 11 to 18       lag=14 connection strength=0.3   
X=Ritardi_sas(X,11,9,rit+16,Amp+0.35);  %from 11 to 9       lag=16 connection strength=0.35
X=Ritardi_sas(X,4,15,rit+13,Amp+0.45);  %from 4 to 15       lag=13 connection strength=0.45
X=Ritardi_sas(X,4,2,rit+12,Amp+0.5);  %from 4 to 2       lag=12 connection strength=0.5
X=Ritardi_sas(X,17,7,rit+15,Amp+0.55);  %from 17 to 7       lag=15 connection strength=0.25

X=Ritardi_sas(X,10,20,rit+11,0.5);  %from 10 to the latent 20th signal lag=2 connection strength=0.5
X=Ritardi_sas(X,20,13,rit+14,2*Amp+0.6);  %from 20 to 13       lag=3 connection strength=0.6   

  % Second example of the simulation
  function X=Fai_scambi_e_ritardi(X,rit,Amp)
X=Ritardi_sas(X,1,3,rit+14,Amp+0.2);      
X=Ritardi_sas(X,3,6,rit+11,0.05);
X=Ritardi_sas(X,5,8,rit+10,0.2);
X=Ritardi_sas(X,10,9,rit+12,0.25);
X=Ritardi_sas(X,13,5,rit+13,0.05);
X=Ritardi_sas(X,7,19,rit+16,0.1);
X=Ritardi_sas(X,19,3,rit+15,0.15);
X=Ritardi_sas(X,2,4,rit+10,Amp+0.4);  
X=Ritardi_sas(X,8,6,rit+14,Amp+0.3);  
X=Ritardi_sas(X,11,18,rit+14,Amp+0.3);    
X=Ritardi_sas(X,11,9,rit+16,Amp+0.35);  
X=Ritardi_sas(X,4,15,rit+13,Amp+0.45); 
X=Ritardi_sas(X,4,2,rit+12,Amp+0.5);  
X=Ritardi_sas(X,17,7,rit+15,Amp+0.55);  
X=Ritardi_sas(X,10,20,rit+11,0.5);  
X=Ritardi_sas(X,20,13,rit+14,2*Amp+0.4);