function [zeta,J] = whiteness_res_val(Yiu1,N,ordmax)

[Lt,nc] = size(Yiu1); % Lunghezza effettiva della sequenza.

L = floor(Lt/2); % Set di dati per l'identificazione

N = L - ordmax-8;

zeta = zeros(ordmax,1);
J = zeros(ordmax,1);
chi = 20.1 * ones(ordmax,1); % Chi^2 per 8 gradi di liberta'.

for i = 1:ordmax,

 [Theta,Ji,sigma_e,cov_e,err_e,yprev] = minquad3(Yiu1(1:L,:),N,i,0);
 
 J(i,1) = Ji;
 zeta(i,1) = bianres2( Yiu1(L+1:2*L,:) , N , 8 , i , Theta' , 0 );
 
 clear Theta;

end;

figure, plot([1:ordmax],zeta,'-',[1:ordmax],chi,':',[1:ordmax],zeta,'o')
title('Residual Whiteness Test'), xlabel('Identified ARX Model Order'), ylabel('\chi^{2}_{\alpha}(M) and \zeta_{N,M}')

figure, plot([1:ordmax],J,'-',[1:ordmax],J,'o')
title('Loss Function'), xlabel('Identified ARX Model Order'), ylabel('J(\theta_N)')

return

%%%
%%% Local functions
%%%

function Hn=myhank(Yvet,Nseq,Ord,Off)

   for i=1:Ord,

      Hn(:,i) = Yvet( Off+i: Nseq-Ord+Off+i-1 );

  end;

return  
  
%%%%%%%

function [Theta,J,sigma_e,cov,err,Yprev] = minquad3( Yvet , Nseq , Ord , graph )

  [ nr , nc ] = size( Yvet );    % ci sono ingressi ed uscite;

  if( nc >= 2 ) , %  Prime colonne ingressi e ultima col. uscita.

  H=[ myhank(Yvet(:,nc),Nseq+Ord,Ord,0), myhank(Yvet(:,1),Nseq+Ord , Ord ,0 )];

  else % Nessun ingresso, solo una uscita. 
  
  H= myhank( Yvet(:,1) , Nseq+Ord , Ord , 0 ); % Yvet e' un vettore.

  end;

  if( nc > 2 ),

    for input = 2:nc-1,
     
    H = [ H , myhank( Yvet(:,input) , Nseq+Ord , Ord , 0 )];

    end;

  end;
   
     Theta = pinv( H ) * Yvet( Ord+1:Nseq+Ord , nc) ;
     
   %  Theta = ( inv( H'* H ) )*( H')*Yvet( Ord+1:Nseq+Ord , nc) ;

if (graph == 1),


figure
plot([1:Nseq],H*Theta,'-',[1:Nseq],Yvet( Ord+1:Nseq+Ord ,nc ),'--');
title('Simulazione linea cont. e uscita misurata linea tratt.')
ylabel('Y (-) simulata e (--) misurata'),xlabel('Campioni'), pause;

%print -dps simPe.ps

end;

  
  Yprev = H*Theta;
  errps = Yvet( Ord+1:Nseq+Ord ,nc )-Yprev;

if (graph == 1),

figure
plot([1:Nseq],errps,'-');
title('Errore di previsione')
ylabel('Yv-Yp'),xlabel('Campioni'), pause;

%print -dps errPe.ps

end;


   J = ( errps'* errps )/Nseq;

   sigma_e = ( J * Nseq )/( Nseq - nc*Ord );

   cov = sigma_e * pinv( H' * H );

   err = errps;


return

%%%%

function zeta = bianres2( Y , Ncov , Mgdl , Ord , Theta , graph )

  [ L , nc ] = size( Y ); %L numero di campioni;

  % Ncov lunghezza della sequenza considerata;

  Yo = Y(Ord+1:L,nc);
  
        if( nc >= 2 )

    H = [ myhank( Y(:,nc) , L , Ord , 0 ) , myhank( Y(:,1),L,Ord,0 ) ];

        else H = myhank(  Y(:,1) , L , Ord , 0 );

        end;


        if( nc > 2)

          for input = 2:nc-1,

              H = [ H , myhank(  Y(:,input) , L , Ord , 0 ) ];

          end;

        end;

  Yp = H*Theta'; %vettore delle uscite previste;

if( graph == 1 )
  figure, subplot(211)
  plot([0:Ncov-1],Yo(1:Ncov),'-',[0:Ncov-1],Yp(1:Ncov),':')
% plot([0:L-Ord-1],Yo,'-',[0:L-Ord-1],Yp(1:L-Ord),':')
  title('Prediction (:) and observations (-)')
  xlabel('Samples'), ylabel('y(t)')
end;

  epserr=[]; %vettore dei residui;

  res=[]; %vettore delle covarianze campionarie;

  epserr=Yo-Yp(1:L-Ord);

if( graph == 1 )
  figure, subplot(212)
  plot([0:Ncov-1],epserr(1:Ncov)),title('Identified model residuals')
  xlabel('Samples'), ylabel('\epsilon(t)'), pause;


end;

           for i=1:Mgdl+1,

             res(i) = ( (epserr(1:Ncov))' * epserr(i:Ncov+i-1) )/Ncov;

           end;

if( graph == 1 )
     figure, plot([1:Mgdl+1],res,'-',[1:Mgdl+1],res,'o');
     title('Residual whiteness test'), xlabel('Correlation (\tau)'), pause;
end;

  zeta = Ncov *( res(2:Mgdl+1)*( res(2:Mgdl+1)' ) )/(res(1)*res(1));

return