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

% Toy Model 2
% Baccala & Sameshima. Partial directed coherence: a new concept in neural
% structure determination. Biol. Cybern. 84:463-474, 2001.
%                http://dx.doi.org/10.1007/PL00007990
% Example 3 (pag.468)
%                       VAR(3) with feedback

% disp('======================================================================');
% disp('      Model 2: Linear five-dimensional VAR[3] Model')
% disp(' Baccala & Sameshima. Biol.Cybern.84:463-74, 2001. Example 3 (pag.468)')
% disp('             x1-->x2  x1-->x3 x1-->x4 x4-->x5 x5-->x4');
% disp('======================================================================');


if nargin < 2, flgPrint = 0; end;

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

p = 3; % True model order

N=nDiscard+nPoints+p;

randn('state', sum(100*clock))
ei=randn(5,N);
x1=zeros(1,N);
x2=zeros(1,N);
x3=zeros(1,N);
x4=zeros(1,N);
x5=zeros(1,N);
% Variables initialization
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;

chLabels = []; % or
%chLabels = {'x_1';'x_2';'x_3';'x_4';'x_5'};


for t=5:N,
    x1(t) = 0.95*sqrt(2)*x1(t-1) - 0.9025*x1(t-2) + ei(1,t);
    x2(t) = 0.5*x1(t-2) + ei(2,t);
    x3(t) = -0.4*x1(t-3) + ei(3,t);
    x4(t) = -0.5*x1(t-2) + 0.25*sqrt(2)*x4(t-1) + 0.25*sqrt(2)*x5(t-1) + ei(4,t);
    x5(t) = -0.25*sqrt(2)*x4(t-1) + 0.25*sqrt(2)*x5(t-1) + ei(5,t);
end;


y=[x1' x2' x3' x4' x5']; % data must be organized column-wise

u=y(nDiscard+1:nDiscard+nPoints,:);

%==========================================================================
%                     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;

% VAR model as well as performing GCT and ICT.
[Tr_gct, pValue_gct, Tr_igct, pValue_igct] = gct_alg(u,A,pf,gct_signif, ...
                                                     flgPrintResults);

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

% Estimate PDC and asymptotic statistics
c = asymp_pdc(u,A,pf,nFreqs,metric,alpha);

% 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;


%Call MVGC toolbox for cMVGC calculation
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;



%Use MVGC toolbox with true model order.
if morder == p, % If the estimated model order == true order, just copy the result.
    
    f = e;
else
    %Use MVGC toolbox with true model order,which in this is p = 3.
    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(['Using true model order p = ', int2str(IP), ...
              '. Further computation unnecessary.'])
    g = c;
else
    disp(['using true model order p = ', int2str(p), '.'])
    
    fs=1; maxIP = p; % p = 3
    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;

    % VAR model as well as performing GCT and ICT
    [Tr_gct, pValue_gct, Tr_igct, pValue_igct] = gct_alg(u,A,pf,gct_signif, ...
        flgPrintResults);
        
    % Estimate PDC and asymptotic statistics
    g=asymp_pdc(u,A,pf,nFreqs,metric,alpha);
    
    % 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;
