home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IV4MV.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.9 KB  |  96 lines
  1. function TH=iv4mv(z,nn,maxsize,T,p)
  2. %IV4MV    Computes approximately optimal IV-estimates for multivariate ARX-models
  3. %    Auxiliary routine to iv4
  4. %    TH = iv4mv(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 each output, input as column
  12. %    vectors: Z=[y1 y2 ... yp u1 u2 ... un]
  13. %
  14. %    NN: NN=[na nb nk], giving the orders and delays of the above model.
  15. %
  16. %    Some parameters associated with the algorithm are accessed by
  17. %    TH = iv4mv(Z,NN,maxsize,T)
  18. %    See HELP AUXVAR for an explanation of these and their default values.
  19.  
  20. %    L. Ljung 10-3-90
  21. %    Copyright (c) 1986-90 by the MathWorks, Inc.
  22. %    All Rights Reserved.
  23.  
  24. % *** Set up default values ***
  25. [nnr,nnc]=size(nn);
  26. ny=nnr;nu=(nnc-ny)/2;
  27. if nu<=0,error('This routine is not intended for time series. Use ARX or AR instead!'),end
  28. [Ncap,nz]=size(z); 
  29. if nz>Ncap,error(' Data should be organized in column vectors!'),return,end
  30. if nz~=(nu+ny),disp(' Incorrect number of orders specified!'),disp('For an AR-model nn=na'),disp('For an ARX-model, nn=[na nb nk]'),
  31. error('see above'),return,end
  32. na=nn(:,1:ny);
  33. if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end
  34. if nu>0,nk=nn(:,ny+nu+1:nnc);else nk=zeros(ny,1);end
  35.  
  36. maxsdef=idmsize(Ncap);
  37. if nargin<5, p=1;end
  38. if nargin<4, T=1;end
  39. if nargin<3, maxsize=maxsdef;end
  40.  
  41. if p<0,p=1;end, if T<0, T=1;end,if maxsize<0,maxsize=maxsdef;end
  42.  
  43. % *** Do some consistency tests ***
  44.  
  45.  
  46. if nz>Ncap, error('The data should be organized in column vectors')
  47. return,end
  48. %
  49. % *** First stage: compute an LS model ***
  50. %
  51. th=mvarx(z,nn,maxsize,T);
  52. a=th2ss(th);
  53. maxa=max(abs(eig(a)));
  54. if maxa>1, %Stabilize! 
  55.     [a,b]=th2arx(th);[nar,nac]=size(a);
  56.     for ka=1:(nac/ny)-1
  57.     a(:,ka*ny+1:(ka+1)*ny)=((0.9/maxa)^ka)*a(:,ka*ny+1:(ka+1)*ny);
  58.     end
  59.   th=arx2th(a,b,ny,nu);
  60. end  
  61. %
  62. % *** Second stage: Compute the IV-estimates using the LS
  63. %     model a, b for generating the instruments ***
  64. %
  65. yh=idsimss(z(:,ny+1:nz),th);
  66. for outp=1:ny
  67. if norm(yh(:,outp),inf)<eps,
  68. error(['Output number ' int2str(outp) ' is decoupled from all inputs. IV4 cannot be used!']),end,end
  69. th=ivxmv(z,nn,yh,maxsize,T,0); 
  70. a=th2ss(th);th1=th;
  71. maxa=max(abs(eig(a)));
  72. if maxa>1, %Stabilize! 
  73.     [a,b]=th2arx(th);[nar,nac]=size(a);
  74.     for ka=1:(nac/ny)-1
  75.     a(:,ka*ny+1:(ka+1)*ny)=((0.9/maxa)^ka)*a(:,ka*ny+1:(ka+1)*ny);
  76.     end
  77.   th1=arx2th(a,b,ny,nu);
  78. end  
  79.  
  80. %
  81. % *** Third stage: Compute the residuals, v, associated with the
  82. %     current model, and determine an AR-model, A, for these ***
  83. %
  84. v=pe(z,th);if ny>1, v=sum(v')';end
  85. art=arx(v,(max(sum(na))+max(sum(nb))),maxsize,T,0);
  86.  
  87. %
  88. % *** Fourth stage: Use the optimal instruments ***
  89. %
  90. for kz=1:nz
  91. zf(:,kz)=filter([1 art'],1,z(:,kz));
  92. end
  93. xf=idsim(zf(:,ny+1:nz),th1);
  94. TH=ivxmv(zf,nn,xf,maxsize,T,1);
  95. TH(2,7)=34;
  96.