home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IVXMV.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.2 KB  |  135 lines
  1. function eta=ivxmv(z,nn,x,maxsize,Tsamp,p);
  2. %IVXMV    Estimates (multivariate) ARX-models using instrumental variables
  3. %
  4. %    TH = ivxmv(Z,[NA NB NK],X)
  5. %
  6. %    Z: The output-input data Z = [Y U] with the outputs Y
  7. %    and the inputs U arranged in column vectors
  8. %    [NA NB NK]: defines the model structure as follows:
  9. %    NA: an ny|ny matrix whose i-j entry is the order of the
  10. %        polynomial (in the delay operator) that relates the
  11. %        j:th output to the i:th output
  12. %    NB: an ny|nu matrix whose i-j entry is the order of the 
  13. %        polynomial that related the j:th input to the i:th output
  14. %    NK: an ny|nu matrix whose i-j entry is the delay from
  15. %        the j:th input to the i:th output
  16. %    (ny: the number of outputs, nu: the number of inputs)
  17. %    X:  The instrumental variables. Same format as Y.
  18. %    
  19. %    TH: the resulting model, given in the THETA-format (see HELP THETA)
  20. %
  21. %    With TH=ivxmv(Z,[NA NB NK],X,MAXSIZE,T) some additional
  22. %    options can be achieved. See HELP AUXVAR for details
  23.  
  24. %    
  25. %    L. Ljung 10-3-90
  26. %    Copyright (c) 1987-90 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. [nnr,nnc]=size(nn);
  30. ny=nnr;nu=(nnc-ny)/2;
  31. [Ncap,nz]=size(z);[Nc,nx]=size(x); 
  32. if nz>Ncap,error(' Data should be organized in column vectors!'),return,end
  33. if Nc~=Ncap | nx~=ny, error('The matrix x shall have the same number of rows as z, and the number of columns must be equal to the number of outputs'),end
  34. 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]'),
  35. error('see above'),return,end
  36. na=nn(:,1:ny);
  37. if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end
  38. if nu>0,nk=nn(:,ny+nu+1:nnc);else nk=zeros(ny,1);end
  39.  
  40.  
  41. nma=max(max(na)');nbkm=max(max(nb+nk)')-1;nkm=min(min(nk)');
  42. nd=sum(sum(na)')+sum(sum(nb)');
  43. n=nma*ny+(nbkm-nkm+1)*nu;
  44.  
  45. % *** Set up default values **
  46. maxsdef=idmsize(Ncap,nd);
  47. if nargin<6, p=1;end
  48. if nargin<4, maxsize=maxsdef;end
  49.  
  50. if maxsize<0, maxsize=maxsdef; end
  51. if nargin<5, Tsamp=1;end
  52. if isempty(Tsamp),Tsamp=1;end,if isempty(maxsize),maxsize=maxsdef;end
  53. % *** construct regression matrix ***
  54.  
  55.  nmax=max([nma nbkm]');
  56. M=floor(maxsize/n);
  57. R=zeros(n);F=zeros(n,ny);
  58. for k=nmax:M:Ncap
  59.       if min(Ncap,k+M)<k+1,ntz=0;else ntz=1;end
  60.     
  61.       if ntz,jj=(k+1:min(Ncap,k+M));phi=zeros(length(jj),n);phix=phi;end
  62.       
  63.       for kl=1:nma, 
  64.       if ntz,
  65.     phi(:,(kl-1)*ny+1:kl*ny)=z(jj-kl,1:ny);
  66.     phix(:,(kl-1)*ny+1:kl*ny)=x(jj-kl,1:ny);end
  67.       end
  68.       ss=nma*ny;
  69.   
  70.       for kl=nkm:nbkm,
  71.       if ntz,phi(:,ss+(kl-nkm)*nu+1:ss+(kl-nkm+1)*nu)=z(jj-kl,ny+1:nz);
  72.             phix(:,ss+(kl-nkm)*nu+1:ss+(kl-nkm+1)*nu)=z(jj-kl,ny+1:nz);    end
  73.       end    
  74.       if ntz,R=R+phix'*phi; F=F+phix'*z(jj,1:ny);end
  75.       
  76. end
  77. %
  78. % *** compute loss functions ***
  79. %
  80. par=[];parb=[];p11=[];p12=[];p22=[];lamm=[];
  81. for outp=1:ny
  82.     rowind=[];
  83.     for kk=1:ny
  84.     rowind=[rowind kk:ny:ny*na(outp,kk)];
  85.     end
  86.     for kk=1:nu
  87.     baslev=nma*ny+(nk(outp,kk)-nkm)*nu;
  88.     rowind=[rowind baslev+kk:nu:baslev+nu*nb(outp,kk)];
  89.     end
  90.    rowind=sort(rowind);
  91.  
  92.  
  93.     th=R(rowind,rowind)\F(rowind,outp); 
  94.  
  95. k00=sum((nk(outp,:)==0).*(nb(outp,:)>0));
  96. naux=sum(na(outp,:));lth(outp)=length(th);
  97. ind1=1:naux;ind2=naux+k00+1:length(th);ind3=naux+1:naux+k00;
  98. par=[par th(ind1)' th(ind2)'];parb=[parb th(ind3)'];
  99.  
  100. end
  101. eta=mketaarx(nn,[par parb],eye(ny),Tsamp);
  102. e=pe(z,eta);lamm=e'*e/(Ncap-nmax);eta=mketaarx(nn,[par parb],lamm,Tsamp);
  103. eta(2,7)=38;
  104. if p,
  105. for outp=1:ny
  106.     rowind=[];
  107.     for kk=1:ny
  108.     rowind=[rowind kk:ny:ny*na(outp,kk)];
  109.     end
  110.     for kk=1:nu
  111.     baslev=nma*ny+(nk(outp,kk)-nkm)*nu;
  112.     rowind=[rowind baslev+kk:nu:baslev+nu*nb(outp,kk)];
  113.     end
  114.    rowind=sort(rowind);
  115.  
  116.    pth=inv(R(rowind,rowind));
  117.    LAM=lamm(outp,outp);
  118.    pth=LAM*pth;
  119. naux=sum(na(outp,:));
  120. ind1=1:naux;ind2=naux+k00+1:lth(outp);ind3=naux+1:naux+k00;
  121.  
  122. [sp11,ssp11]=size(p11);lni=length(ind1)+length(ind2);
  123. p11=[[p11,zeros(sp11,lni)];[zeros(lni,sp11),[[pth(ind1,ind1),pth(ind1,ind2)];...
  124. [pth(ind2,ind1),pth(ind2,ind2)]]]];
  125. [sp12,ssp12]=size(p12);
  126. p12=[[p12,zeros(sp11,length(ind3))];[zeros(lni,ssp12),[pth(ind1,ind3);pth(ind2,ind3)]]];
  127. [sp22,ssp22]=size(p22);
  128. p22=[[p22,zeros(sp22,length(ind3))];[zeros(length(ind3),ssp22),pth(ind3,ind3)]];
  129.  
  130. end
  131. pvar=[[p11,p12];[p12',p22]];
  132. d=length([par parb]);
  133. eta(4:3+d,1:d)=pvar;
  134. end
  135.