function [c,d,e,f,g] = fguo2008(nPoints)

% Guo, Wu, Ding & Feng. Uncovering interactions in frequency domains.
%         PLoS Computational Biology, 4(5):1-10, February 8, 2008.
%           [http://dx.plos.org/10.1371/journal.pcbi.1000087]
%
% Page 2 Toy Model Example 5-dimensional VAR[3] with large common exogenous input

% disp('======================================================================');
% disp('                   Five dimensional linear VAR[3] Model')
% disp('                    with large common exogenous input')
% disp('                      Guo et al. February 8, 2008.')
% disp('                 x1==>x2  x1-->x3 x1==>x4 x4-->x5 x5-->x4 ');
% disp('             Instantaneous causality in between all variables.');
% disp('======================================================================');

%flgPrintScreen = 'Screen'
%flgPrintScreen = 'Print'
if nargin < 2, flgPrint = 0; end;

nDiscard=5000;  % number of points discarded at beginning of series
%nPoints  =   % number of analyzed samples points

N=nDiscard+nPoints;

p = 3; % True model order

% Variables initialization
ei=randn(7,N);
x1=zeros(1,N);
x2=zeros(1,N);
x3=zeros(1,N);
x4=zeros(1,N);
x5=zeros(1,N);

a1=rand; a2=rand; a3=rand; a4=rand; a5=rand;
b1=2; b2=2; b3=2; b4=2; b5=2;
c1=5; c2=5; c3=5; c4=5; c5=5;

% disp(sprintf('a1=%f; a2=%f; a3=%f; a4=%f; a5=%f;',a1,a2,a3,a4,a5))
% disp('b1 = b2 = b3 = b4 = b5 = 2;')
% disp('c1 = c2 = c3 = c4 = c5 = 5;')
% disp('======================================================================');

for t=1:4,
    x1(t)=randn(1); x2(t)=randn(1); x3(t)=randn(1); x4(t)=randn(1);
    x5(t)=randn(1);
end;

for t=5:N,
    x1(t) = 0.95*sqrt(2)*x1(t-1) - 0.9025*x1(t-2) + ei(1,t) + ...
        a1*ei(6,t) + b1*ei(7,t-1) + c1*ei(7,t-2);
    x2(t) = 0.5*x1(t-2) + ei(2,t) + ...
        a2*ei(6,t) + b2*ei(7,t-1) + c2*ei(7,t-2);
    x3(t) =-0.4*x1(t-3) + ei(3,t) + ...
        a3*ei(6,t) + b3*ei(7,t-1) + c3*ei(7,t-2);
    x4(t) =-0.5*x1(t-2) + 0.25*sqrt(2)*x4(t-1) + 0.25*sqrt(2)*x5(t-1) ...
        + ei(4,t) + a4*ei(6,t) + b4*ei(7,t-1) + c4*ei(7,t-2);
    x5(t) =-0.25*sqrt(2)*x4(t-1) + 0.25*sqrt(2)*x5(t-1) ...
        + ei(5,t) + a5*ei(6,t) + b5*ei(7,t-1) + c5*ei(7,t-2);
end;




y=[x1' x2' x3' x4' x5']; % data must be organized column-wise
u=y(nDiscard+1:nDiscard+nPoints,:);

[nSegLength,nChannels]=size(u);

%==========================================================================
%                     DTFn estimation and analysis parameters
%==========================================================================


%run example_analysis_parameters % Setting default paramaters

%==========================================================================
%                     Estimation and analysis parameters
%==========================================================================

fs=1; maxIP=30;
criterion=1; % AIC - Akaike Information Criteria
alg=1; %1 = Nutall-Strand MVAR estimation algorithm
nFreqs=32;
alpha=0.01;     % Significance level
gct_signif = 0.01; % Granger causality test significance level
igct_signif = 0.01; % Instantaneous GCT significance level
MVARadequacy_signif = 0.05; % VAR model estimation adequacy significance
metric='info'; % euc  = original PDC;
% diag = generalized PDC or gPDC;
% info = information PDC or iPDC.


%run example_pre_processing      % Detrending and/or standardization
%==========================================================================
%                    Detrend and normalization options
%==========================================================================
flgDetrend = 1;
flgStandardize = 0;

[nChannels,nSegLength]=size(u);
if nChannels > nSegLength, u=u.';
    [nChannels,nSegLength]=size(u);
end;

if flgDetrend,
    for i=1:nChannels, u(i,:)=detrend(u(i,:)); end;
    disp('Time series were detrended.');
end;

if flgStandardize,
    for i=1:nChannels, u(i,:)=u(i,:)/std(u(i,:)); end;
    disp('Time series were scale-standardized.');
end;



%run example_mvar_estimation     % Estimating VAR and testing adequacy of

%==========================================================================
%                            MVAR model estimation
%==========================================================================
[IP,pf,A,pb,B,ef,eb,vaic,Vaicv] = mvar(u,maxIP,alg,criterion);

disp(['Number of channels = ' int2str(nChannels) ' with ' ...
    int2str(nSegLength) ' data points; MAR model order = ' int2str(IP) '.']);

%==========================================================================
%    Testing for adequacy of MAR model fitting through Portmanteau test
%==========================================================================
h=20; % testing lag
aValueMVAR = 1 - MVARadequacy_signif;
flgPrintResults = 0;
[Pass,Portmanteau,st,ths]=mvarresidue(ef,nSegLength,IP,aValueMVAR,h,...
    flgPrintResults);

%==========================================================================
%         Granger causality test (GCT) and instantaneous GCT
%==========================================================================
flgPrintResults = 0;
[Tr_gct, pValue_gct, Tr_igct, pValue_igct] = gct_alg(u,A,pf,gct_signif, ...
    flgPrintResults);


% VAR model as well as performing GCT and ICT.

%=======Overriding some default parameters for plotting and analysis=======

%[hxlabel,hylabel] = xplot(c,flgPrinting,fs,w_max,chLabels, ...
%                          flgColor,flgScale,flgMax,flgSignifColor)

flgPrinting = [1 1 1 2 0 0 3]; % overriding default setting
% alpha=0.05;         % Significance level
% gct_signif = 0.05;  % Granger causality test significance level
% igct_signif = 0.05; % Instantaneous GCT significance level
% flgColor = 0;

%run example_dtf_analysis        % Estimate DTF and asymptotic statistics
d = asymp_dtf(u,A,pf,nFreqs,metric,alpha);
% Power spectra and coherence calculation




c=asymp_pdc(u,A,pf,nFreqs,metric,alpha);% Estimate PDC and asymptotic statistics
% Power spectra and coherence calculation
c.SS = ss_alg(A, pf, nFreqs);
c.coh = coh_alg(c.SS);

c.Pass = Pass;
c.Portmanteau = Portmanteau;
c.Tr_gct = Tr_gct;
c.pValue_gct = pValue_gct;
c.Tr_igct = Tr_igct;
c.pValue_igct = pValue_igct;
c.Vaicv = Vaicv;
c.criterion = criterion;
c.alg = alg;
c.fs = fs;
c.p = IP;

%Use MVGC tool box
morder = 'AIC';
[pval,F,sig,pvalun,Fun,sigun,Amvgc,morder] = fmvgc(u,alpha,morder,A);

e.pval = pval;
e.F = F;
e.sig = sig;
e.pvalun = pvalun;
e.Fun = Fun;
e.sigun = sigun;
e.A = Amvgc;
e.p = morder;


if morder == p,
    
    f = e;
else
    %Use MVGC toolbox with true model order,which in this is p =1.
    morder = p;
    [pval,F,sig,pvalun,Fun,sigun,Amvgc,morder] = fmvgc(u,alpha,morder);
    
    f.pval = pval;
    f.F = F;
    f.sig = sig;
    f.pvalun = pvalun;
    f.Fun = Fun;
    f.sigun = sigun;
    f.A = Amvgc;
    f.p = morder;
end;
% =========================================================================

if IP == p,
    disp(['Estimated model order = true model = ', int2str(IP), ...
              '. No further computation.'])
    g = c;
else
    disp(['using true model order p = ', int2str(p), '.'])
    
    fs=1; maxIP = p; % p = 1
    criterion = 5; % Fixed order

    alg=1; %1 = Nutall-Strand MVAR estimation algorithm
    nFreqs=32;
    alpha=0.01;     % Significance level
    gct_signif = 0.01; % Granger causality test significance level
    igct_signif = 0.01; % Instantaneous GCT significance level
    MVARadequacy_signif = 0.05; % VAR model estimation adequacy significance
    metric='info'; % euc  = original PDC;
    % diag = generalized PDC or gPDC;
    % info = information PDC or iPDC.
    
    %run example_mvar_estimation     % Estimating VAR and testing adequacy of
    
    %==========================================================================
    %                            MVAR model estimation
    %==========================================================================
    [IP,pf,A,pb,B,ef,eb,vaic,Vaicv] = mvar(u,maxIP,alg,criterion);
    
    disp(['Number of channels = ' int2str(nChannels) ' with ' ...
        int2str(nSegLength) ' data points; MAR model order = ' int2str(IP) '.']);
    
    %==========================================================================
    %    Testing for adequacy of MAR model fitting through Portmanteau test
    %==========================================================================
    h=20; % testing lag
    aValueMVAR = 1 - MVARadequacy_signif;
    flgPrintResults = 0;
    [Pass,Portmanteau,st,ths]=mvarresidue(ef,nSegLength,IP,aValueMVAR,h,...
        flgPrintResults);
    
    %==========================================================================
    %         Granger causality test (GCT) and instantaneous GCT
    %==========================================================================
    flgPrintResults = 0;
    [Tr_gct, pValue_gct, Tr_igct, pValue_igct] = gct_alg(u,A,pf,gct_signif, ...
        flgPrintResults);
    
    
    % VAR model as well as performing GCT and ICT
    
    % Power spectra and coherence calculation
    
    g=asymp_pdc(u,A,pf,nFreqs,metric,alpha);% Estimate PDC and asymptotic statistics
    % Power spectra and coherence calculation
    g.SS = ss_alg(A, pf, nFreqs);
    g.coh = coh_alg(g.SS);
    
    g.Pass = Pass;
    g.Portmanteau = Portmanteau;
    g.Tr_gct = Tr_gct;
    g.pValue_gct = pValue_gct;
    g.Tr_igct = Tr_igct;
    g.pValue_igct = pValue_igct;
    g.Vaicv = Vaicv;
    g.criterion = criterion;
    g.alg = alg;
    g.fs = fs;
    g.p = IP;
end;

