function TH=inival(z,nn,maxsize,Tsamp);
%INIVAL   Computes initial values for the prediction error method.
%	
%	TH = inival(Z,NN,MAXSIZE,T)
%
%	The routine is intended primarily as a subroutine to PEM.
%	See PEM for an explanation of the arguments.

%	L. Ljung 10-1-86
%	Copyright (c) 1986 by the MathWorks, Inc.
%	All Rights Reserved.

% *** Set up the relevant orders ***

[Ncap,nz]=size(z); nu=nz-1;
na=nn(1);nb=nn(2:1+nu); nc=nn(2+nu); nd=nn(3+nu);
nf=nn(4+nu:3+2*nu); nk=nn(4+2*nu:3+3*nu);
Ndcap=sum(nn(1:3+nu)); n=Ndcap+sum(nf);

% *** Begin construction the TH-matrix ***

TH=zeros(3+n,max([n 6+3*nu 7])); a=[1];, b=[0];
TH(1,2:6+3*nu)=[Tsamp nu nn];

% *** step 1: compute the LS-estimate of A and B
 
 naf=na+max(nf);
 t=arx(z,[naf nb nk],maxsize,Tsamp,0);
% ** Stabilize the a-estimate **
   a=fstab([1 t(1:naf)']); 
   b=zeros(nu,max(nb+nk));
  NBcap=cumsum([na nb]);
  for k=1:nu
        if nb(k)>0
           b(k,nk(k)+1:nk(k)+nb(k))=t(NBcap(k)+1:NBcap(k+1))';
        end
  end
 
% *** step 2: compute IV-estimates

if naf>0
%     if na>0 we compute the IV-estimate of A and B in a
%     straightforward fashion and initilize F (if any)
%     as F=1.
%
  if na>0
    t=iv(z,[na nb nk],a,b,maxsize,Tsamp,0);
    TH(3,1:na+sum(nb))=t';
  end

%    if na=0 we initialize each of the F as the donominator
%    dynamics of each of the corresponding single input systems
%    estimated by the IV-method. The estimate is first made stable
 
 if na==0
  s1=1; s=1;
  for k=1:nu
     t=iv([z(:,1) z(:,k+1)],[nf(k) nb(k) nk(k)],a,b(k,:), ...
           maxsize,Tsamp,0);
     if nf(k)>0, f=fstab([1 t(1:nf(k))']);
        TH(3,Ndcap+s:Ndcap+s+nf(k)-1)=f(2:length(f));
        s=s+nf(k);
     end
     if nb(k)>0, TH(3,na+s1:na+s1+nb(k)-1)=t(nf(k)+1:nf(k)+nb(k))';,end
     s1=s1+nb(k);
  end
 end
end

%  *** step 3: calculate the residuals associated with the
%              current model

if nc+nd>0
   v=pe(z,TH);

%   *** step 4: build an ARMA(nd,nc) model of these residuals.

   if nc==0
      t=arx(v,nd,maxsize,Tsamp,0);
      TH(3,na+sum(nb)+1:na+sum(nb)+nd)=t';
   end

   if nc>0
      art=arx(v,n,maxsize,Tsamp,0);
      ep=filter([1 art'],1,v);
      t=arx([v ep],[nd nc 1],maxsize,Tsamp,0);
%     ** check stability of C-polynomial **
      c=fstab([1 t(nd+1:nd+nc)']);
      TH(3,na+sum(nb)+1:na+sum(nb)+nc)=c(2:length(c));
      if nd>0
         TH(3,na+sum(nb)+nc+1:na+sum(nb)+nc+nd)=t(1:nd)';
      end
   end
end

% *** Build up the TH-matrix ***

e=pe(z,TH);
TH(1,1)=e'*e/length(e); % The estimate of lambda
TH(2,1)=(1+n/Ncap)/(1-n/Ncap)*TH(1,1); % Akaike's FPE
ti=fix(clock); ti(1)=ti(1)/100;
TH(2,2:6)=ti(1:5); TH(2,7)=9;
