home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IVAR.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.2 KB  |  77 lines
  1. function TH=ivar(y,na,nc,maxsize,T)
  2. %IVAR    Computes IV-estimates for the AR-part of a scalar time series
  3. %
  4. %    TH = ivar(Y,NA)
  5. %
  6. %    TH: returned as the IV-estimate of the AR-part of a time series
  7. %    model
  8. %    A(q) y(t) =  v(t)
  9. %    along with relevant structure information. (See HELP THETA for
  10. %    the exact structure of TH.
  11. %
  12. %    Y : the time series as a column vector.
  13. %
  14. %    NA: the order of the AR-part (the A(q)-polynomial)
  15. %    
  16. %    In the above model v(t) is an arbitrary process, assumed to be a
  17. %    (possibly time-varying) moving average process of order NC. The
  18. %    default value of NC is NA. Other values of NC are obtained by
  19. %    TH = ivar(Y,NA,NC)
  20. %     
  21. %    With
  22. %    TH = ivar(Y,NA,NC,maxsize,T)
  23. %    access to some parameters associated with the algorithm is obtained.
  24. %    See HELP AUXVAR for an explanation of these.
  25.  
  26. %    L. Ljung 4-21-87
  27. %    Copyright (c) 1987-90 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. [Ncap,nz]=size(y); 
  31.  
  32. % 
  33. % *** Some initial tests on the input arguments ***
  34. %
  35. maxsdef=idmsize(Ncap);
  36.  
  37. if nargin<5, T=1;end
  38. if nargin<4, maxsize=maxsdef;end
  39. if nargin<3, nc=na;end
  40. if T<0,T=1;end, if maxsize<0,maxsize=maxsdef;,end,if nc<0,nc=na;end
  41. if isempty(T),T=1;end,if isempty(maxsize),maxsize=maxsdef;end,
  42. if isempty(nc),nc=na;end
  43. if nz>Ncap,y=y';end
  44. if nz>1,error('The data should be a single time series!'),return,end
  45. if length(na)>1,error(' NA should be a single number!'),return,end
  46.  
  47. % *** Estimate a C-polynomial by first building a high order AR-model and
  48. %     then using the corresponding residuals as inputs in an ARX-model ***
  49.  
  50. ar=arx(y,na+nc,maxsize,T,0);
  51. e=filter([1 ar'],1,y);
  52. ac=arx([y e],[na nc 1],maxsize,T,0);
  53. c=fstab([1 ac(na+1:na+nc)']);
  54.  
  55. % *** construct instruments ***
  56.  
  57.      x=[zeros(nc,1);y(1:Ncap-nc)];
  58.      cc=conv(c,c);
  59.      x=filter(1,cc,x);
  60.      TH=ivx(y,na,x,maxsize,T);
  61.      V=e'*e/Ncap;
  62.      yf=filter(1,c,y);
  63.      ef=filter(1,c,e);
  64.      M=floor(maxsize/(na+nc));
  65.      R=zeros(na+nc);
  66.      for k1=max(na,nc):M:Ncap-1
  67.          jj=(k1+1:min(Ncap,k1+M));
  68.     psi=zeros(length(jj),na+nc);
  69.     for k=1:na, psi(:,k)=yf(jj-k);end
  70.     for k=1:nc, psi(:,na+k)=ef(jj-k);end
  71.     R=R+psi'*psi;
  72.      end
  73.      P=V*inv(R);
  74.      TH(1,1)=V; TH(2,1)=V*(1+na/Ncap)/(1-na/Ncap);
  75.      TH(4:3+na,1:na)=P(1:na,1:na);
  76.     TH(2,7)=10;
  77.