home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ARX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.9 KB  |  96 lines
  1. function th=arx(z,nn,maxsize,Tsamp,p);
  2. %ARX    Computes LS-estimates of ARX-models
  3. %
  4. %    TH = arx(Z,NN)
  5. %
  6. %    TH: returned as the estimated parameters of the ARX model
  7. %    A(q) y(t) = B(q) u(t-nk) + e(t)
  8. %    along with estimated covariances and structure information.
  9. %    For the exact structure 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]. For time series
  13. %    Z=y only.
  14. %
  15. %    NN: NN = [na nb nk], the orders and delays of the above model.
  16. %    For multi-output systems, NN has as many rows as there are outputs
  17. %    na is then an ny|ny matrix whose i-j entry gives the order of the
  18. %    polynomial (in the delay operator) relating the j:th output to the
  19. %    i:th output. Similarily nb and nk are ny|nu matrices. (ny:# of outputs,
  20. %    nu:# of inputs). For a time series, NN=na only.
  21. %    Some parameters associated with the algorithm are accessed by
  22. %    TH = arx(Z,NN,maxsize,T)
  23. %    See HELP AUXVAR for an explanation of these and their default values.
  24.  
  25. %    L. Ljung 10-1-86,12-8-91
  26. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. % *** Set up default values **
  30. [Ncap,nz]=size(z); [nr,nc]=size(nn);
  31. maxsdef=idmsize(Ncap); % Mod 891007
  32. if nargin<5, p=1;end
  33. if nargin<4, Tsamp=[];end
  34. if nargin<3, maxsize=[];end
  35. if isempty(Tsamp),Tsamp=1;end,if Tsamp<0,Tsamp=1;end
  36. if isempty(maxsize),maxsize=maxsdef;end,if maxsize<0,maxsize=maxsdef;end
  37. if p<0,p=1;end
  38.  
  39. if nz>Ncap,error('Data should be organized in column vectors!'),return,end
  40. if nr>1, eval('th=mvarx(z,nn,maxsize,Tsamp);'),return,end
  41. nu=nz-1;
  42. if length(nn)~=1+2*(nz-1)
  43. disp(' Incorrect number of orders specified!')
  44. disp('For an AR-model nn=na'),disp('For an ARX-model, nn=[na nb nk]')
  45. error('see above'),return,end
  46. if nz>1, na=nn(1);nb=nn(2:1+nu);nk=nn(2+nu:1+2*nu);
  47. else na=nn(1); nb=0;,nk=0;end
  48. n=na+sum(nb);
  49. if nn==0,return,end
  50.  
  51. % *** construct regression matrix ***
  52.  
  53. nmax=max([na+1 nb+nk])-1;
  54. M=floor(maxsize/n);
  55. R=zeros(n);F=zeros(n,1);
  56. for k=nmax:M:Ncap-1
  57.       jj=(k+1:min(Ncap,k+M));
  58.       phi=zeros(length(jj),n);
  59.       for kl=1:na, phi(:,kl)=-z(jj-kl,1);end
  60.       ss=na;
  61.       for ku=1:nu
  62.            for kl=1:nb(ku), phi(:,ss+kl)=z(jj-kl-nk(ku)+1,ku+1);end
  63.            ss=ss+nb(ku);
  64.       end
  65.       if Ncap>M | p~=0, R=R+phi'*phi; F=F+phi'*z(jj,1);end
  66. end
  67. %
  68. % *** compute estimate ***
  69. %
  70. if Ncap>M, th=R\F; else th=phi\z(jj,1);end
  71. if p==0, return,end
  72. %
  73. % proceed to compute loss function and covariance matrix
  74. %
  75. t=th;, clear th;
  76.  
  77. %
  78. % build up the theta-matrix
  79. %
  80. if nu>0
  81.      I=[1 Tsamp nu na nb 0 0 zeros(1,nu) nk];
  82. else I=[1 Tsamp 0 na 0 0];end
  83.  
  84. n=na+sum(nb);
  85. th=zeros(n+3,max([length(I) n 7]));
  86. th(1,1:length(I))=I;
  87. ti=fix(clock); ti(1)=ti(1)/100;
  88. th(2,2:6)=ti(1:5);
  89. th(2,7)=1;
  90. th(3,1:n)=t.';
  91. e=pe(z,th);V=e'*e/(Ncap-nmax);
  92. th(4:n+3,1:n)=V*inv(R); % mod 911209
  93. th(1,1)=V;
  94. th(2,1)=V*(1+n/Ncap)/(1-n/Ncap); %Akaike's FPE
  95.  
  96.