function xcorrel(method,NetDef,NN,W1,W2,par1,par2,par3)
%  xcorrel
%  ------- 
%          Show six different cross-correlation functions determined from
%          the predictions of a neural network input-output model of a 
%          dynamic system.
%          I.e. a network model which has been generated by nnarx, nnrarx,
%          nnarmax1+2, nnrarmax1+2, or nnoe.
%
%          The following cross-correlation functions are computed:
%          o  Cross-correlation between inputs and residuals
%          o  Cross-correlation between squared inputs and squared residuals
%          o  Cross-correlation between squared inputs and residuals
%          o  Cross-correlation between beta and residuals where
%             beta is a signal composed of the products of inputs and residuals
%          o  Cross-correlation between alpha and squared residuals where
%             alpha is a signal composed of the products of outputs and res.
%          o  Cross-correlation between alpha and squared inputs
%
%
%  References:
%  Billings, Jamaluddin, Chen: "Properties of Neural Networks with Applications
%            to Modelling Non-linear Dynamical Systems," Int. Journ. of Control
%            vol. 55, no. 1, pp. 193-224.
%  Billings, Zhu: "Nonlinear Model Validation Using Correlation Tests,"
%            Int. Journ. of Control, vol. 60, no. 6, pp. 1107-1120.
%
%
%  Call: 
%  Network generated by nnarx (or nnrarx):
%           xcorrel('nnarx',NetDef,NN,W1,W2,Y,U)
%
%  Network generated by nnarmax1 (or nnrarmx1):
%           xcorrel('nnarmax1',NetDef,NN,W1,W2,C,Y,U)
%
%  Network generated by nnarmax2 (or nnrarmx2):
%           xcorrel('nnarmax2',NetDef,NN,W1,W2,Y,U)
%
%  Network generated by nnoe:
%           xcorrel('nnoe',NetDef,NN,W1,W2,Y,U)
%
 
%  Programmed by : Magnus Norgaard, IAU/IMM
%  LastEditDate  : Oct. 17, 1997


% >>>>>>>>>>>>>>>>>>>>>>>>>>>>      GET PARAMETERS     <<<<<<<<<<<<<<<<<<<<<<<<<<<< 
skip = 1;
if strcmp(method,'nnarx') | strcmp(method,'nnrarx'),
  mflag=1;
  Y=par1;
  U=par2;
elseif strcmp(method,'nnarmax1') | strcmp(method,'nnrarmx1'),
  mflag=2;
  C=par1;
  Y=par2; 
  U=par3; 
elseif strcmp(method,'nnarmax2') | strcmp(method,'nnrarmx2'),
  mflag=3;
  Y=par1; 
  U=par2;
elseif strcmp(method,'nnoe'),
  mflag=4;
  Y=par1;
  U=par2;
else
  disp('Unknown method!!!!!!!!');
  break
end


% >>>>>>>>>>>>>>>>>>>>>>>>>>>>     INITIALIZATIONS     <<<<<<<<<<<<<<<<<<<<<<<<<<<< 
Ndat     = length(Y);                   % # of data
na = NN(1);

% ---------- NNARX model ----------
if mflag==1 | mflag==4,
  if length(NN)==1                      % nnar model
    nb = 0;
    nk = 0;
    nu = 0;
  else                                  % nnarx or nnoe model
    [nu,N] = size(U); 
    nb = NN(2:1+nu); 
    nk = NN(2+nu:1+2*nu);
  end
  nc = 0;

% --------- NNARMAX1 model --------
elseif mflag==2 | mflag==3,
  if length(NN)==2                      % nnarma model
    nc     = NN(2);
    nb     = 0;
    nk     = 0;
    nu     = 0;
  else                                  % nnarmax model
    [nu,Ndat]= size(U); 
    nb     = NN(2:1+nu);
    nc     = NN(2+nu);
    nk     = NN(2+nu+1:2+2*nu);
  end
end


% --------- Common initializations --------
nmax     = max([na,nb+nk-1,nc]);        % 'Oldest' signal used as input to the model
N        = Ndat - nmax;                 % Size of training set
nab      = na+sum(nb);                  % na+nb
nabc     = nab+nc;                      % na+nb+nc
outputs     = 1;                        % Only MISO models considered
L_hidden = find(NetDef(1,:)=='L')';     % Location of linear hidden neurons
H_hidden = find(NetDef(1,:)=='H')';     % Location of tanh hidden neurons
L_output = find(NetDef(2,:)=='L')';     % Location of linear output neurons
H_output = find(NetDef(2,:)=='H')';     % Location of tanh output neurons
[hidden,inputs] = size(W1);
inputs          = inputs-1;
E        = zeros(outputs,N);
y1       = zeros(hidden,N);
Yhat     = zeros(outputs,N);


% >>>>>>>>>>>>>>>>>>>>  CONSTRUCT THE REGRESSION MATRIX PHI   <<<<<<<<<<<<<<<<<<<<<
PHI = zeros(nab,N);
jj  = nmax+1:Ndat;
for k = 1:na, PHI(k,:)    = Y(jj-k); end
index = na;
for kk = 1:nu,
  for k = 1:nb(kk), PHI(k+index,:) = U(kk,jj-k-nk(kk)+1); end
  index = index + nb(kk);
end


% >>>>>>>>>>>>>>>>>>>>>>>>>>   COMPUTE NETWORK OUTPUT   <<<<<<<<<<<<<<<<<<<<<<<<<<<
% ---------- NNARX model ----------
if mflag==1,
  Y  = Y(nmax+1:Ndat);
  h1 = W1*[PHI;ones(1,N)];  
  y1(H_hidden,:) = pmntanh(h1(H_hidden,:));
  y1(L_hidden,:) = h1(L_hidden,:);
    
  h2 = W2*[y1;ones(1,N)];
  Yhat(H_output,:) = pmntanh(h2(H_output,:));
  Yhat(L_output,:) = h2(L_output,:);

  E     = Y - Yhat;                       % Error between Y and deterministic part
  SSE      = E*E';                        % Sum of squared errors (SSE)
  PI       = SSE/(2*N);                   % Performance index

% --------- NNARMAX1 model --------
elseif mflag==2,
  Y  = Y(nmax+1:Ndat);
  h1 = W1*[PHI;ones(1,N)];  
  y1(H_hidden,:) = pmntanh(h1(H_hidden,:));
  y1(L_hidden,:) = h1(L_hidden,:);
    
  h2 = W2*[y1;ones(1,N)];
  Yhat(H_output,:) = pmntanh(h2(H_output,:));
  Yhat(L_output,:) = h2(L_output,:);

  Ebar     = Y - Yhat;                    % Error between Y and deterministic part
  E        = filter(1,C,Ebar);            % Prediction error
  Yhat     = Y - E;                       % One step ahead prediction

  SSE      = E*E';                        % Sum of squared errors (SSE)
  PI       = SSE/(2*N);                   % Performance index


% --------- NNARMAX2 model --------
elseif mflag==3,
  Y  = Y(nmax+1:Ndat);
  PHI_aug=[PHI;zeros(nc,N);ones(1,N)];
  y1 = [y1;ones(1,N)];
  N2=N+1-skip;
  for t=1:N,
    h1 = W1*PHI_aug(:,t);  
    y1(H_hidden,t) = pmntanh(h1(H_hidden));
    y1(L_hidden,t) = h1(L_hidden);    

    h2 = W2*y1(:,t);
    Yhat(H_output,t) = pmntanh(h2(H_output,:));
    Yhat(L_output,t) = h2(L_output,:);

    E(:,t) = Y(:,t) - Yhat(:,t);          % Prediction error
    for d=1:min(nc,N-t),
      PHI_aug(nab+d,t+d) = E(:,t);
    end
  end
  SSE      = E(skip:N)*E(skip:N)';        % Sum of squared errors (SSE)
  PI       = SSE/(2*N2);                  % Performance index


% ---------- NNOE model ----------
elseif mflag==4,
  Y  = Y(nmax+1:Ndat);
  PHI_aug=[PHI;ones(1,N)];
  y1 = [y1;ones(1,N)];
  N2=N+1-skip;
  for t=1:N,
    h1 = W1*PHI_aug(:,t);;  
    y1(H_hidden,t) = pmntanh(h1(H_hidden));
    y1(L_hidden,t) = h1(L_hidden);    

    h2 = W2*y1(:,t);
    Yhat(H_output,t) = pmntanh(h2(H_output,:));
    Yhat(L_output,t) = h2(L_output,:);

    for d=1:min(na,N-t),
      PHI_aug(d,t+d) = Yhat(:,t);
    end
  end
  E     = Y - Yhat;                     % Error between Y and deterministic part
  SSE      = E(skip:N)*E(skip:N)';      % Sum of squared errors (SSE)
  PI       = SSE/(2*N2);                % Performance index
end


% >>>>>>>>>>>>>>>>>>>>>  CALCULATE CORRELATION FUNCTIONS   <<<<<<<<<<<<<<<<<<<<<<
M     = min(25,N-1);                    % Max. lag investigated
conf  = 1.96/sqrt(N);                   % Confidence bound
E2    = E(1,:).*E(1,:);                 % Signal composed of squared residuals
Ecov  = cov(E);                         % Sample variance of residual
E2cov = cov(E2);                        % Sample variance of squared residuals
alp   = Y(1,:).*E(1,:);                 % Mixture of output and residual
alpcov=cov(alp);                        % Sample variance of mixture signal

for i=1:nu,
  U2    = U(i,1:N).*U(i,1:N);           % Signal composed of squared inputs
  bet   = U(i,1:N).*E(1,:);             % Mixture of output and residual
  U2cov = cov(U2);                      % Sample variance of squared inputs
  Ucov  = cov(U(i,:));                  % Sample variance of input
  betcov=cov(bet);                      % Sample variance of mixture signal


  figure;
  titletext = '                   cross-correlation functions (input #';
  titletext = [titletext sprintf('%g',i) ', output #1)'];

  
  % --------- Correlation function 1 (u and e) ----------
  subplot(3,2,1)
  UEcross=xcov(E,U(i,1:N),'unbiased')/sqrt(Ecov*Ucov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 1') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);
  title(titletext)

  % --------- Correlation function 2 (u^2 and e^2) ----------
  subplot(3,2,2)
  UEcross=xcov(E2,U2,'unbiased')/sqrt(E2cov*U2cov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 2') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);


  % --------- Correlation function 3 (u^2 and e) ----------
  subplot(3,2,3)
  UEcross=xcov(E,U2,'unbiased')/sqrt(Ecov*U2cov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 3') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);


  % --------- Correlation function 4 (beta and e) ----------
  subplot(3,2,4)
  UEcross=xcov(E,bet,'unbiased')/sqrt(Ecov*betcov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 4') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);
  
  
  % --------- Correlation function 5 (alpha and e^2) ----------
  subplot(3,2,5)
  UEcross=xcov(E2,alp,'unbiased')/sqrt(E2cov*alpcov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 5') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);
  
  
  % --------- Correlation function 6 (alpha and u^2) ----------
  subplot(3,2,6)
  UEcross=xcov(U2,alp,'unbiased')/sqrt(U2cov*alpcov)';
  plot([-M:M], UEcross(N-M:N+M),'b-'); hold on
  plot([-M M],[conf -conf;conf -conf],'r--');hold off
  xlabel('lag')
  ylabel('Function 6') 
  ymax=min(5*conf,max([abs(UEcross)]));
  axis([-M M -ymax ymax]);

  subplot(111)
  drawnow
end