Documentation of inival

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

TH=inival(z,nn,maxsize,Tsamp);

Help text

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.

Cross-Reference Information

Listing of function inival

function TH=inival(z,nn,maxsize,Tsamp);

%   L. Ljung 10-1-86
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.4 $  $Date: 1997/12/02 03:42:43 $

% *** 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([naf 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 initialize 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 denominator
%    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
    if nb(k)>0
     t=iv([z(:,1) z(:,k+1)],[nf(k) nb(k) nk(k)],a,b(k,:), ...
           maxsize,Tsamp,0);
     if any(isinf(t)),TH=[];return,end
     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
     TH(3,na+s1:na+s1+nb(k)-1)=t(nf(k)+1:nf(k)+nb(k)).';
     s1=s1+nb(k);
    end
  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;