Global Index (short | long) | Local contents | Local Index (short | long)
TH=inival(z,nn,maxsize,Tsamp);
INIVAL Computes initial values for the prediction error method. TH = inival(Z,NN,MAXSIZE,T) The routine is intended primarily as a subroutine to PEM. See PEM for an explanation of the arguments.
This function calls | This function is called by |
---|---|
function TH=inival(z,nn,maxsize,Tsamp); % L. Ljung 10-1-86 % Copyright (c) 1986-98 by The MathWorks, Inc. % $Revision: 2.4 $ $Date: 1997/12/02 03:42:43 $ % *** Set up the relevant orders *** [Ncap,nz]=size(z); nu=nz-1; na=nn(1);nb=nn(2:1+nu); nc=nn(2+nu); nd=nn(3+nu); nf=nn(4+nu:3+2*nu); nk=nn(4+2*nu:3+3*nu); Ndcap=sum(nn(1:3+nu)); n=Ndcap+sum(nf); % *** Begin construction the TH-matrix *** TH=zeros(3+n,max([n 6+3*nu 7])); a=[1];, b=[0]; TH(1,2:6+3*nu)=[Tsamp nu nn]; % *** step 1: compute the LS-estimate of A and B naf=na+max(nf); t=arx(z,[naf nb nk],maxsize,Tsamp,0); % ** Stabilize the a-estimate ** a=fstab([1 t(1:naf).']); b=zeros(nu,max(nb+nk)); NBcap=cumsum([naf nb]); for k=1:nu if nb(k)>0 b(k,nk(k)+1:nk(k)+nb(k))=t(NBcap(k)+1:NBcap(k+1)).'; end end % *** step 2: compute IV-estimates if naf>0 % if na>0 we compute the IV-estimate of A and B in a % straightforward fashion and initialize F (if any) % as F=1. % if na>0 t=iv(z,[na nb nk],a,b,maxsize,Tsamp,0); TH(3,1:na+sum(nb))=t.'; end % if na=0 we initialize each of the F as the denominator % dynamics of each of the corresponding single input systems % estimated by the IV-method. The estimate is first made stable if na==0 s1=1; s=1; for k=1:nu if nb(k)>0 t=iv([z(:,1) z(:,k+1)],[nf(k) nb(k) nk(k)],a,b(k,:), ... maxsize,Tsamp,0); if any(isinf(t)),TH=[];return,end if nf(k)>0, f=fstab([1 t(1:nf(k)).']); TH(3,Ndcap+s:Ndcap+s+nf(k)-1)=f(2:length(f)); s=s+nf(k); end TH(3,na+s1:na+s1+nb(k)-1)=t(nf(k)+1:nf(k)+nb(k)).'; s1=s1+nb(k); end end end end % *** step 3: calculate the residuals associated with the % current model if nc+nd>0 v=pe(z,TH); % *** step 4: build an ARMA(nd,nc) model of these residuals. if nc==0 t=arx(v,nd,maxsize,Tsamp,0); TH(3,na+sum(nb)+1:na+sum(nb)+nd)=t.'; end if nc>0 art=arx(v,n,maxsize,Tsamp,0); ep=filter([1 art.'],1,v); t=arx([v ep],[nd nc 1],maxsize,Tsamp,0); % ** check stability of C-polynomial ** c=fstab([1 t(nd+1:nd+nc).']); TH(3,na+sum(nb)+1:na+sum(nb)+nc)=c(2:length(c)); if nd>0 TH(3,na+sum(nb)+nc+1:na+sum(nb)+nc+nd)=t(1:nd).'; end end end % *** Build up the TH-matrix *** e=pe(z,TH); TH(1,1)=e'*e/length(e); % The estimate of lambda TH(2,1)=(1+n/Ncap)/(1-n/Ncap)*TH(1,1); % Akaike's FPE ti=fix(clock); ti(1)=ti(1)/100; TH(2,2:6)=ti(1:5); TH(2,7)=9;