Global Index (short | long) | Local contents | Local Index (short | long)
eta=mvarx(z,nn,maxsize,Tsamp);
MVARX Estimates (multivariate) ARX-models Auxiliary routine to ARX TH = mvarx(Z,[NA NB NK]) Z: The output-input data Z = [Y U] with the outputs Y and the inputs U arranged in column vectors [NA NB NK]: defines the model structure as follows: NA: an ny|ny matrix whose i-j entry is the order of the polynomial (in the delay operator) that relates the j:th output to the i:th output NB: an ny|nu matrix whose i-j entry is the order of the polynomial that related the j:th input to the i:th output NK: an ny|nu matrix whose i-j entry is the delay from the j:th input to the i:th output (ny: the number of outputs, nu: the number of inputs) Each row of [NA NB NK] is thus consistent with the way single output ARX structures is defined (See Help ARX) TH: the resulting model, given in the THETA-format (see HELP THETA) With TH=mvarx(Z,[NA NB NK],MAXSIZE,T) some additional options can be achieved. See HELP AUXVAR for details
This function calls | This function is called by |
---|---|
function eta=mvarx(z,nn,maxsize,Tsamp); % % L. Ljung 10-2-90 % Copyright (c) 1986-98 by The MathWorks, Inc. % $Revision: 3.3 $ $Date: 1997/12/02 03:42:38 $ [nnr,nnc]=size(nn); ny=nnr;nu=(nnc-ny)/2; [Ncap,nz]=size(z); if nz>Ncap,error(' Data should be organized in column vectors!'),return,end 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]'), error('see above'),return end na=nn(:,1:ny); if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end if nu>0,nk=nn(:,ny+nu+1:nnc);else nk=zeros(ny,1);end nma=max(max(na)');nbkm=max(max(nb+nk)')-1;nkm=min(min(nk)'); nd=sum(sum(na)')+sum(sum(nb)'); n=nma*ny+(nbkm-nkm+1)*nu; % *** Set up default values ** maxsdef=idmsize(Ncap,nd); if nargin<4, Tsamp=1;end if nargin<3, maxsize=maxsdef;end if isempty(Tsamp),Tsamp=1;end,if Tsamp<0,Tsamp=1;end % Changed T to Tsamp, JNL 9/19/91 if isempty(maxsize),maxsize=maxsdef;end,if maxsize<0, maxsize=maxsdef; end % *** construct regression matrix *** nmax=max([nma nbkm]'); M=floor(maxsize/n); R=zeros(n);F=zeros(n,ny); for k=nmax:M:Ncap if min(Ncap,k+M)<k+1,ntz=0;else ntz=1;end if ntz,jj=(k+1:min(Ncap,k+M));phi=zeros(length(jj),n);end for kl=1:nma, if ntz, phi(:,(kl-1)*ny+1:kl*ny)=z(jj-kl,1:ny); end end ss=nma*ny; for kl=nkm:nbkm, if ntz phi(:,ss+(kl-nkm)*nu+1:ss+(kl-nkm+1)*nu)=z(jj-kl,ny+1:nz); end end if ntz,R=R+phi'*phi; F=F+phi'*z(jj,1:ny);end end % % *** compute loss functions *** % par=[];parb=[];p11=[];p12=[];p22=[];lamm=[]; for outp=1:ny rowind=[]; for kk=1:ny rowind=[rowind kk:ny:ny*na(outp,kk)]; end for kk=1:nu baslev=nma*ny+(nk(outp,kk)-nkm)*nu; rowind=[rowind baslev+kk:nu:baslev+nu*nb(outp,kk)]; end rowind=sort(rowind); if Ncap>M, th=R(rowind,rowind)\F(rowind,outp); else th=phi(:,rowind)\z(jj,outp); end pth=inv(R(rowind,rowind)); LAM=(z(nmax+1:Ncap,outp)'*z(nmax+1:Ncap,outp)-F(rowind,outp)'*pth*... F(rowind,outp))/(Ncap-nmax); pth=LAM*pth; lamm=[[lamm,zeros(length(lamm),1)];[zeros(1,length(lamm)),LAM]]; k00=sum((nk(outp,:)==0).*(nb(outp,:)>0)); naux=sum(na(outp,:)); ind1=1:naux;ind2=naux+k00+1:length(th);ind3=naux+1:naux+k00; par=[par th(ind1).' th(ind2).'];parb=[parb th(ind3).']; [sp11,ssp11]=size(p11);lni=length(ind1)+length(ind2); p11=[[p11,zeros(sp11,lni)];[zeros(lni,sp11),... [[pth(ind1,ind1),pth(ind1,ind2)];... [pth(ind2,ind1),pth(ind2,ind2)]]]]; [sp12,ssp12]=size(p12); [sp22,ssp22]=size(p22); p12=[[p12,zeros(sp11,length(ind3))];[zeros(lni,ssp12),... [pth(ind1,ind3);pth(ind2,ind3)]]]; p22=[[p22,zeros(sp22,length(ind3))];... [zeros(length(ind3),ssp22),pth(ind3,ind3)]]; end eta=mketaarx(nn,[par parb],lamm,Tsamp); eta(2,7)=31;e=pe(z,eta); lamm=e'*e/(Ncap-nmax); pvar=[[p11,p12];[p12',p22]]; d=length([par parb]); eta(4:3+d,1:d)=pvar;rarg=eta(1,6); eta(4+d+rarg:3+ny+d+rarg,1:ny)=lamm;eta(1,1)=det(lamm); eta(2,1)=eta(1,1)*(1+nd/Ncap)/(1-nd/Ncap);