home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IV4.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.0 KB  |  88 lines
  1. function TH=iv4(z,nn,maxsize,T,p)
  2. %IV4    Computes approximately optimal IV-estimates for ARX-models.
  3. %
  4. %    TH = iv4(Z,NN)
  5. %
  6. %    TH: returned as the estimate of the ARX model
  7. %    A(q) y(t) = B(q) u(t-nk) + v(t)
  8. %    along with estimted covariances and structure information.
  9. %    For the exact format of TH see HELP THETA.
  10. %
  11. %    Z : the output-input data Z=[y u], with y and u as column vectors.
  12. %    For multivariable systems Z=[y1 y2 .. yp u1 u2 .. um].
  13. %
  14. %    NN: NN = [na nb nk], the orders and delays of the above model.
  15. %    For multi-output systems, NN has as many rows as there are outputs
  16. %    na is then an ny|ny matrix whose i-j entry gives the order of the
  17. %    polynomial (in the delay operator) relating the j:th output to the
  18. %    i:th output. Similarily nb and nk are ny|nu matrices. (ny:# of outputs,
  19. %    nu:# of inputs).
  20. %
  21. %    Some parameters associated with the algorithm are accessed by
  22. %    TH = iv4(Z,NN,maxsize,T)
  23. %    See HELP AUXVAR for an explanation of these and their default values.
  24.  
  25. %    L. Ljung 10-1-86,4-15-90
  26. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. % *** Set up default values ***
  30. [Ncap,nz]=size(z);
  31. maxsdef=idmsize(Ncap);
  32. if nargin<5, p=1;end
  33. if nargin<4, T=1;end
  34. if nargin<3, maxsize=maxsdef;end
  35.  
  36. if p<0,p=1;end, if T<0, T=1;end,if maxsize<0,maxsize=maxsdef;end
  37. if isempty(T),T=1;end,if isempty(maxsize),maxsize=maxsdef;end
  38.  
  39. % *** Do some consistency tests ***
  40. [nr,nc]=size(nn);if nr>1, eval('TH=iv4mv(z,nn,maxsize,T,p);'),return,end
  41. [Ncap,nz]=size(z); nu=nz-1;
  42.  
  43. if nz>Ncap, error('The data should be organized in column vectors')
  44. return,end
  45. if nz==1, error('This routine does not make sense for a time series!')
  46. return,end
  47. if length(nn)~=1+2*nu, disp('Incorrect number of orders specified:')
  48. disp(' nn should be  nn=[na nb nk]')
  49. disp(' where nb anb nk are row vectors of length equal to the number of inputs'),error('see above')
  50. return,end
  51.  
  52. na=nn(1);nb=nn(2:1+nu);nk=nn(2+nu:1+2*nu); n=na+sum(nb);
  53. %
  54. % *** First stage: compute an LS model ***
  55. %
  56. th=arx(z,[na nb nk],maxsize,T,0);
  57. if na>0, a=fstab([1 th(1:na)']);
  58.     else a=1;end
  59. b=zeros(nu,max(nb+nk));
  60. NBcum=cumsum([na nb]);
  61. for k=1:nu, b(k,nk(k)+1:nk(k)+nb(k))=th(NBcum(k)+1:NBcum(k+1))';,end
  62. %
  63. % *** Second stage: Compute the IV-estimates using the LS
  64. %     model a, b for generating the instruments ***
  65. %
  66. th=iv(z,[na nb nk],a,b,maxsize,T,0); 
  67. if na>0,a=[1 th(1:na)'];
  68.     else a=1;end
  69. for k=1:nu, b(k,nk(k)+1:nk(k)+nb(k))=th(NBcum(k)+1:NBcum(k+1))';,end
  70. %
  71. % *** Third stage: Compute the residuals, v, associated with the
  72. %     current model, and determine an AR-model, A, for these ***
  73. %
  74. v=filter(a,1,z(:,1));
  75. for k=1:nu, v=v-filter(b(k,:),1,z(:,k+1));,end
  76. art=arx(v,na+sum(nb),maxsize,T,0);
  77. Acap=[1 art'];
  78. %
  79. % *** Fourth stage: Use the optimal instruments ***
  80. %
  81. yf=filter(Acap,[1],z(:,1));
  82. for k=1:nu, uf(:,k)=filter(Acap,[1],z(:,k+1));,end
  83. if p~=0, p=2;,end
  84. if na>1,a=fstab(a);else a=1;end
  85. TH=iv([yf uf],[na nb nk],a,b,maxsize,T,p);
  86. TH(2,7)=4;
  87. % *** Reference: Equations (15.221) - (15.26) in Ljung(1987) ***
  88.