home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / INIVAL.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.7 KB  |  103 lines
  1. function TH=inival(z,nn,maxsize,Tsamp);
  2. %INIVAL   Computes initial values for the prediction error method.
  3. %    
  4. %    TH = inival(Z,NN,MAXSIZE,T)
  5. %
  6. %    The routine is intended primarily as a subroutine to PEM.
  7. %    See PEM for an explanation of the arguments.
  8.  
  9. %    L. Ljung 10-1-86
  10. %    Copyright (c) 1986 by the MathWorks, Inc.
  11. %    All Rights Reserved.
  12.  
  13. % *** Set up the relevant orders ***
  14.  
  15. [Ncap,nz]=size(z); nu=nz-1;
  16. na=nn(1);nb=nn(2:1+nu); nc=nn(2+nu); nd=nn(3+nu);
  17. nf=nn(4+nu:3+2*nu); nk=nn(4+2*nu:3+3*nu);
  18. Ndcap=sum(nn(1:3+nu)); n=Ndcap+sum(nf);
  19.  
  20. % *** Begin construction the TH-matrix ***
  21.  
  22. TH=zeros(3+n,max([n 6+3*nu 7])); a=[1];, b=[0];
  23. TH(1,2:6+3*nu)=[Tsamp nu nn];
  24.  
  25. % *** step 1: compute the LS-estimate of A and B
  26.  
  27.  naf=na+max(nf);
  28.  t=arx(z,[naf nb nk],maxsize,Tsamp,0);
  29. % ** Stabilize the a-estimate **
  30.    a=fstab([1 t(1:naf)']); 
  31.    b=zeros(nu,max(nb+nk));
  32.   NBcap=cumsum([na nb]);
  33.   for k=1:nu
  34.         if nb(k)>0
  35.            b(k,nk(k)+1:nk(k)+nb(k))=t(NBcap(k)+1:NBcap(k+1))';
  36.         end
  37.   end
  38.  
  39. % *** step 2: compute IV-estimates
  40.  
  41. if naf>0
  42. %     if na>0 we compute the IV-estimate of A and B in a
  43. %     straightforward fashion and initilize F (if any)
  44. %     as F=1.
  45. %
  46.   if na>0
  47.     t=iv(z,[na nb nk],a,b,maxsize,Tsamp,0);
  48.     TH(3,1:na+sum(nb))=t';
  49.   end
  50.  
  51. %    if na=0 we initialize each of the F as the donominator
  52. %    dynamics of each of the corresponding single input systems
  53. %    estimated by the IV-method. The estimate is first made stable
  54.  
  55.  if na==0
  56.   s1=1; s=1;
  57.   for k=1:nu
  58.      t=iv([z(:,1) z(:,k+1)],[nf(k) nb(k) nk(k)],a,b(k,:), ...
  59.            maxsize,Tsamp,0);
  60.      if nf(k)>0, f=fstab([1 t(1:nf(k))']);
  61.         TH(3,Ndcap+s:Ndcap+s+nf(k)-1)=f(2:length(f));
  62.         s=s+nf(k);
  63.      end
  64.      if nb(k)>0, TH(3,na+s1:na+s1+nb(k)-1)=t(nf(k)+1:nf(k)+nb(k))';,end
  65.      s1=s1+nb(k);
  66.   end
  67.  end
  68. end
  69.  
  70. %  *** step 3: calculate the residuals associated with the
  71. %              current model
  72.  
  73. if nc+nd>0
  74.    v=pe(z,TH);
  75.  
  76. %   *** step 4: build an ARMA(nd,nc) model of these residuals.
  77.  
  78.    if nc==0
  79.       t=arx(v,nd,maxsize,Tsamp,0);
  80.       TH(3,na+sum(nb)+1:na+sum(nb)+nd)=t';
  81.    end
  82.  
  83.    if nc>0
  84.       art=arx(v,n,maxsize,Tsamp,0);
  85.       ep=filter([1 art'],1,v);
  86.       t=arx([v ep],[nd nc 1],maxsize,Tsamp,0);
  87. %     ** check stability of C-polynomial **
  88.       c=fstab([1 t(nd+1:nd+nc)']);
  89.       TH(3,na+sum(nb)+1:na+sum(nb)+nc)=c(2:length(c));
  90.       if nd>0
  91.          TH(3,na+sum(nb)+nc+1:na+sum(nb)+nc+nd)=t(1:nd)';
  92.       end
  93.    end
  94. end
  95.  
  96. % *** Build up the TH-matrix ***
  97.  
  98. e=pe(z,TH);
  99. TH(1,1)=e'*e/length(e); % The estimate of lambda
  100. TH(2,1)=(1+n/Ncap)/(1-n/Ncap)*TH(1,1); % Akaike's FPE
  101. ti=fix(clock); ti(1)=ti(1)/100;
  102. TH(2,2:6)=ti(1:5); TH(2,7)=9;
  103.