function DEMO

osc = input('Choose a set of oscillators (type 1,2,3,4 or 5) \n or 3 uncoupled white noise processes (6), or a set of shifted signals (7): ');

N=10000;

if osc == 1                                                                 % Choose oscillator system
    M=5;
    P=4;
    x=zeros(N,M);
    x(1:P,:)=rand(P,M);
    w=randn(N,M);
    for t=P+1:N
        x(t,:) = [0.6*x(t-1,1)+0.65*x(t-2,2)+w(t,1), ...
            0.5*x(t-1,2)-0.3*x(t-2,2)-0.3*x(t-4,3)+0.6*x(t-1,4)+w(t,2), ...
            0.8*x(t-1,3)-0.7*x(t-2,3)-0.1*x(t-3,5)+w(t,3), ...
            0.5*x(t-1,4)+0.9*x(t-2,3)+0.4*x(t-2,5)+w(t,4), ...
            0.7*x(t-1,5)-0.5*x(t-2,5)-0.2*x(t-1,3)+w(t,5)];
    end

elseif osc == 2 % Baccala-Sameshima, 2001, Fig.2
    M=5;
    P=3;
    x=zeros(N,M);
    x(1:P,:)=rand(P,M);
    w=randn(N,M);
    for t=P+1:N
        x(t,:) = [0.95*sqrt(2)*x(t-1,1)-0.9025*x(t-2,1)+w(t,1), ...
            0.5*x(t-2,1)+w(t,2), ...
            -0.4*x(t-3,1)+w(t,3), ...
            -0.5*x(t-2,1)+0.25*sqrt(2)*x(t-1,4)+0.25*sqrt(2)*x(t-1,5)+w(t,4), ...
            -0.25*sqrt(2)*x(t-1,4)+0.25*sqrt(2)*x(t-1,5)+w(t,5)];
    end
elseif osc == 3 % Baccala-Sameshima, 2001, Fig.4f
    M=5;
    P=2;
    x=zeros(N,M);
    x(1:P,:)=rand(P,M);
    w=randn(N,M);
    for t=P+1:N
        x(t,:) = [0.95*sqrt(2)*x(t-1,1)-0.9025*x(t-2,1)+0.5*x(t-2,5)+w(t,1), ...
            -0.5*x(t-1,1)+w(t,2), ...
            0.4*x(t-2,2)+w(t,3), ...
            -0.5*x(t-1,3)+0.25*sqrt(2)*x(t-1,4)+0.25*sqrt(2)*x(t-1,5)+w(t,4), ...
            -0.25*sqrt(2)*x(t-1,4)+0.25*sqrt(2)*x(t-1,5)+w(t,5)];
    end
elseif osc == 4
    M=3;
    P=4;
    x=zeros(N,M);
    x(1:P,:)=rand(P,M);
    w=randn(N,M);
    for t=P+1:N
        x(t,:) = [0.6*x(t-1,1)+0.65*x(t-2,2)+w(t,1), ...
            0.5*x(t-1,2)-0.3*x(t-2,2)-0.3*x(t-4,3)+w(t,2), ...
            0.8*x(t-1,3)-0.7*x(t-2,3)+w(t,3)];
    end
elseif osc == 5
    M=3;
    P=2;
    x=zeros(N,M);
    x(1:P,:)=rand(P,M);
    w=randn(N,M);
    for t=P+1:N
        x(t,:) = [0.95*sqrt(2)*x(t-1,1)-0.9025*x(t-2,1)+w(t,1), ...
            -0.5*x(t-1,1)+w(t,2), ...
            0.4*x(t-2,2)+w(t,3)];
    end
elseif osc == 6 % Winterhalder et al 2005, Fig. 5
    M=3;
    x=zeros(N,M);
    w=randn(N,M);
    for t=1:N
        x(t,:) = [w(t,1), 500*w(t,2), 500*w(t,3)];
    end
elseif osc == 7 % Korzeniewska et al 2003, Fig. 1
    M=8;
    P=6;
    x=zeros(N,M);
    x(1:P,1) = randn(P,1);
    a=0.5;
    w=randn(N,1);
    for t=P+1:N
        x(t,1) = a*x(t-1,1) - a*x(t-2,1) + a*x(t-3,1) - a*x(t-4,1) + a*x(t-5,1) - a*x(t-6,1) + w(t);
    end
    x(:,8) = [rand(1,1); x(1:end-1,1)] + randn(N,1);
    x(:,4) = [rand(1,1); x(1:end-1,1)] + randn(N,1);
    x(:,2) = [rand(1,1); x(1:end-1,4)] + randn(N,1); 
    x(:,5) = [rand(1,1); x(1:end-1,8)] + randn(N,1);
    x(:,3) = [rand(1,1); x(1:end-1,8)] + randn(N,1);
    x(:,6) = [rand(1,1); x(1:end-1,3)] + randn(N,1);
    x(:,7) = [rand(1,1); x(1:end-1,3)] + randn(N,1);
end

mu = sum(x,1)./N;
mu = repmat(mu,N,1);
x=x-mu;                                                                     % remove mean

pmin = 1;                                                                   % check model orders from pmin...
pmax = 20;                                                                  % to pmax.
[pbic, pfpe, logdpbic, logdpfpe] = ModelOrder2(x,pmin,pmax);                % calculate optimal model order according to BIC and FPE
disp(['BIC: The optimal model order is: ', num2str(pbic),' (logdp = ', num2str(logdpbic),')'])
disp(['FPE: The optimal model order is: ', num2str(pfpe),' (logdp = ', num2str(logdpfpe),')'])
disp(' ')

P = pbic;                                                                   % Set the order to be the optimal one
[AR,RC,sigma] = mvar(x, P, 2);                                              % Call mvar to perform a MVAR analysis

%disp('Performing VAR stability test')
z = findroots(AR,M,P);
if all(z <= 1) == 1
    disp('Stability criterion is satisfied')
else
    disp('Stability criterion is violated!');
    disp('End of process!')
    return
end

B = eye(M);                                                                 % B is the identity matrix
A = [eye(M),-AR];                                                           % Bring Ar in the appropriate form...
C = sigma(:,M*P+1:M*(P+1));                                                 % and keep only the last M columns of sigma (noise covariance)
[AF,S,H,PDC,COH,DTF,pCOH,dDTF] = mvfreqz2(B,A,C,0:511,1024);                % Call mvfreqz2 (modified by me) to calculate the desired measures
signif = Significance(x,AF,C,0:511,P,1024);                                 % Calculate the significance level of PDC

figure
z=1;
for r=1:M
    for c=1:M
        pdc=squeeze(PDC(r,c,:));
        sig=squeeze(signif(r,c,:));
        coh=squeeze(abs(COH(r,c,:)));
        dtf=squeeze(DTF(r,c,:));
        ddtf=squeeze(sqrt(dDTF(r,c,:)));

        subplot(M,M,z);
        if r==c
            sp = log10(abs(squeeze(S(r,c,:))));
            plot(sp,'g')                                                    % Spectrum - GREEN
            % hold on
            % plot(coh,'c')
            xlim([0 512]);
            set(gca,'XTick',[0,512],'XTickLabel',[0,0.5])
        else
            plot(pdc)                                                       % PDC - BLUE
            hold on
            plot(sig,'r')                                                   % Significance - RED
            plot(dtf,'k')                                                   % DTF - BLACK
            plot(ddtf,'m')                                                  % dDFT - MAGENTA
            plot(coh,'c')                                                   % Coherence - CYAN
            axis([0 512 0 1]);
            set(gca,'XTick',[0,512],'XTickLabel',[0,0.5])
        end
        z=z+1;
    end
end