function [c,d,e,f,g] = model_1(nPoints)
% Baccala & Sameshima. Overcoming the limitations of correlation analysis
% for many simultaneously processed neural structures, Progress in Brain
% Research, 130:33--47, 2001.
%            http://dx.doi.org/10.1016/S0079-6123(01)30004-3
%
% Example 7-dimensional VAR[2] model with loop and feedback

if nargin < 2, flgPrint = 0; end;

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

N=nDiscard+nPoints;

p = 2; % True model order

disp('======================================================================');
disp('                       Linear VAR[2] Model')
disp('        Baccala & Sameshima. Prog Brain Research, 130:33--47, 2001.')
disp('  x1-->x2  x1-->x3  x2-->x3 x3-->x4 x4==>x5 x5-->x1  x5-->x4 x6==>x7');
disp('======================================================================');

randn('state', sum(100*clock))
ei=randn(7,N);
x1=zeros(1,N);
x2=zeros(1,N);
x3=zeros(1,N);
x4=zeros(1,N);
x5=zeros(1,N);
x6=zeros(1,N);
x7=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); x6(t)=randn(1); x7(t)=randn(1);
end;

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


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

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

%==========================================================================
%                     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 1]; % 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 toolbox
%[pval,F,sig,pvalun,Fun,sigun,A,amo] = fmvgc(u,alpha);

morder = 'AIC';
alpha = 0.01;

[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 necessary.'])
    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;


if flgPrint,
    
    d.SS = ss_alg(A, pf, nFreqs);
    d.coh = coh_alg(d.SS);
    %save baccala2001a_ex3_DTF
    
    %        flgColor - 0: white background;
    %                   1: white [.1 1], light-blue [.01 .1], purple [0 .01]
    %                   2: white [.1 1], light-gray [.01 .1], gray [0 .01]
    
    flgColor = 0;
    %        flgScale - 1: [0 1] / {if max(PDC/DTF) > 1}:[0 max(PDC/DTF)]
    
    %                   2: [0 {if max(PDC/DTF) > 1}:max(PDC/DTF)]/[0 1]/[0 .1]/[0 .01]
    %                                     based on flgMax(PDC/DTF/Threshold/CI/all)
    %                   3: [0 max(PDC/DTF/Thr/CI/all)]
    %                                     based on flgMax(PDC/DTF/Threshold/CI/all)
    flgScale = 3;
    %        flgMax - {'PDC'; 'DTF'; 'Thr'; 'CI'; 'TCI'; 'all'} measure used as
    %                   upper limit for Scale. See flgScale.
    flgMax = 'TCI';
    %         flgSignifColor - 0: black line
    %                          1: black / gray -> significant /not signif PDC/DTF
    %                          2: red  /  gray ->      "
    %                          3: red  / green ->      "
    %                          4: red  / black ->      "
    %                          5: black / green
    flgSignifColor = 3;
    %==========================================================================
    %                        DTFn Matrix Layout Plotting
    %==========================================================================
    w_max=fs/2;
    
    strTitle1 = ['5-dimensional linear VAR[3] Model: '];
    strTitle2 = ['[N=' int2str(nSegLength) 'pts; IP=' int2str(d.p) '; ' ...
        datestr(now) ']'];
    strTitle =[strTitle1 strTitle2];
    
    for kflgColor = flgColor,
        h=figure;
        set(h,'NumberTitle','off','MenuBar','none', ...
            'Name', 'Baccala & Sameshima (2001b) Linear model')
        [hxlabel hylabel] = xplot(d,flgPrinting,fs,w_max,chLabels, ...                                    flgColor,flgScale,flgMax,flgSignifCol
            kflgColor,flgScale,flgMax,flgSignifColor)
        
        %   xplot_title(alpha,metric);
        %    [ax,hT]=suplabel( strTitle, 't' );
        %    set(hT,'FontSize',8)
    end;
end;

