function [VLs,S,VS]=estorderiv(y,u,nmax,L);
//function [V,S,VS]=estorderiv(y,u,nmax,L);
//function for model order estimation. The ordre is estimated using Instrumental Variable Method
//with delayed inputs as the instrumental variables. The criterion has the following 
//forme:
//                            1						      
//                      V(n)=---Min{|Y(0) - Rn(n).theta|^2}  n=1,2,3,...nmax
//                            N
//where
//                      Y(0)= Zm.y(0)
//                      Rn(n)=Zm.R(n)
//inputs:
//y ... 	vector of acquired system outputs
//u ... 	vector of applied inputs
//nmax ...maximum order to be considered
//L ...	size of instrumental matrix. Recommended value for L is 70, otherwise L should be >2nmax.	
//otputs:
//V...		normalized vector of the criterion (V=[V(0) V(1) V(2), ...]^T, where V(0)=1)
//S...		vector of penalty coefficients. Model complexity n is penalised.
//VS...	normalized vector of penalised criterion VS=V+S
//programmed by: Hynek Prochazka 
//
//written by:  H. Prochazka, I.D. Landau
//7th june 2002

sz=size(y);
if(sz(2)~=1), y=y';end;//column
sz=size(u);
if(sz(2)~=1), u=u';end;

VLs=[];
S=0;
VS=[];
L=2*nmax;
for nn=1:nmax,
	N=length(y)+1-L;   
   //construction of R(nn) and Z
   Z=[];
   uz=u;
   for k=1:2*nmax,
      Z=[Z uz(2*nmax-k+1:2*nmax-k+N-1)];
   end;
   Rnny=[];
   Rnnu=[];
   for k=1:nn,
      yp=y(L-k+1:L-k+N-1);
      up=u(L-k+1:L-k+N-1);
      Rnny=[Rnny yp];
      Rnnu=[Rnnu up];
   end;
   y0=y(L+1:length(y));
   //QR-triangularization
   [Q,R]=qr(Z);
	Z1=R(1:L,:);
	Zm=inv(Z1')*Z';
   Y0=Zm*y0;
   R1nn=Zm*[Rnny Rnnu];
	//theta computation
   theta_min=inv(R1nn'*R1nn)*R1nn'*Y0;//inv(R1nn)*Y0;
   criter=sum((Y0-R1nn*theta_min).^2)/N;
   VLs=[VLs;criter];
   penal=0.4*nn*log10(N)/N;
   S=[S;penal];
end;
V0=sum(Y0.^2)/N;
VLs=[V0;VLs];
VLs=VLs/VLs(1);
VS=VLs+S;
VS=VS/VS(1);