home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IVX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.1 KB  |  96 lines
  1. function TH=ivx(z,nn,x,maxsize,Tsamp,p)
  2. %IVX    Computes instrumental variable estimates for ARX-models.
  3. %
  4. %    TH = ivx(Z,NN,X)
  5. %
  6. %    TH: returned as the IV-estimate of the ARX-model
  7. %    A(q) y(t) = B(q) u(t-nk) + v(t)
  8. %    along with relevant structure information. See HELP THETA for
  9. %    the exact structure of TH.
  10. %
  11. %    Z : the output-input data Z=[y u], with y and u as column vectors.
  12. %    For multi-input systems u=[u1 u2 ... un]
  13. %
  14. %    NN: NN=[na nb nk] gives the orders and delays associated with the
  15. %    above model.
  16. %    X : is the vector of instrumental variables. This should be of the
  17. %    same size as the y-vector (i.e. the first column of Z).
  18. %
  19. %     See IV4 for good, automatic choices of instruments.
  20. %    A multioutput variant ig given by IVXMV
  21. %    TH=ivx(Z,NN,X,maxsize,T)
  22. %    allows access to some parameters associated with the algorithm.
  23. %    See HELP AUXVAR for an explanation of these.
  24.  
  25. %    L. Ljung 10-1-86, 4-12-87
  26. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. [Ncap,nz]=size(z); nu=nz-1;[Nc,nx]=size(x);[nnr,nnc]=size(nn);
  30. maxsdef=idmsize(Ncap);
  31. % 
  32. % *** Some initial tests on the input arguments ***
  33. %
  34.  
  35. if nargin<6, p=1;end
  36. if nargin<5, Tsamp=1;end
  37. if nargin<4, maxsize=maxsdef;end
  38. if maxsize<0,maxsize=maxsdef;end
  39. if Tsamp<0,Tsamp=1;end
  40. if p<0,p=1;end
  41. if isempty(Tsamp),Tsamp=1;end,if isempty(maxsize),maxsize=maxsdef;end
  42. if nnr>1,eval('TH=ivxmv(z,nn,x,maxsize,Tsamp,p);'),return,end
  43. if Ncap~=Nc | nx~=1,error('The x-vector should be a column vector with the same number of elements as z!'),return,end
  44. if length(nn)~=1+2*nu,disp('Incorrect number of orders specified!'),
  45.     disp('nn should be nn=[na nb nk]'),
  46.     disp('where nb and nk are row vectors of length equal to the number of inputs'),error('see above'),return,end
  47. na=nn(1);nb=nn(2:1+nu);nk=nn(2+nu:1+2*nu);n=na+sum(nb);
  48. %
  49. % construct regression matrix
  50. %
  51. nmax=max([na+1 nb+nk])-1;
  52. M=floor(maxsize/n);
  53. Rxx=zeros(na);Ruu=zeros(sum(nb));Rxu=zeros(na,sum(nb));Rxy=zeros(na);
  54. Ruy=zeros(sum(nb),na); F=zeros(n,1);
  55. for k=nmax:M:Ncap-1
  56.      jj=(k+1:min(Ncap,k+M));
  57.      phix=zeros(length(jj),na); phiy=phix; phiu=zeros(length(jj),sum(nb));
  58.      for kl=1:na, phiy(:,kl)=-z(jj-kl,1); phix(:,kl)=-x(jj-kl); end
  59.      ss=0;
  60.      for ku=1:nu
  61.             for kl=1:nb(ku), phiu(:,ss+kl)=z(jj-kl-nk(ku)+1,ku+1);end
  62.              ss=ss+nb(ku);
  63.      end
  64.      Rxy=Rxy+phix'*phiy;
  65.      if nu>0,Ruy=Ruy+phiu'*phiy;
  66.          Rxu=Rxu+phix'*phiu;
  67.     Ruu=Ruu+phiu'*phiu;
  68.      end    
  69.      Rxx=Rxx+phix'*phix;
  70.      if na>0, F(1:na)=F(1:na)+phix'*z(jj,1);end
  71.       F(na+1:n)=F(na+1:n)+phiu'*z(jj,1);
  72. end
  73. clear phiu, clear phix, clear phiy, 
  74. %
  75. % compute estimate
  76. %
  77. if nu==0,TH=Rxy\F;end
  78. if nu>0,TH=[Rxy Rxu;Ruy Ruu]\F;end
  79. if p==0, return,end
  80. %
  81. % proceed to build up THETA-matrix
  82. %
  83. t=TH;, clear TH;
  84. I=[0 Tsamp nu na nb 0 0 zeros(1,nu) nk];
  85. n=na+sum(nb);
  86. TH=zeros(n+3,max(length(I),n));
  87. TH(1,1:length(I))=I;
  88. ti=fix(clock); ti(1)=ti(1)/100;
  89. TH(2,2:6)=ti(1:5);
  90. TH(2,7)=11;
  91. TH(3,1:n)=t';
  92. if p==2
  93. e=pe(z,TH); TH(1,1)=e'*e/(length(e)-max(na,sum(nb)));
  94. TH(2,1)=TH(1,1)*(1+n/Ncap)/(1-n/Ncap);
  95.  TH(4:3+n,1:n)=TH(1,1)*inv([Rxx Rxu;Rxu' Ruu]);end
  96.