home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / THSS2TH.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  100 lines
  1. function th=thss2th(eta,iy)
  2. %THSS2TH Converts an internal TH(SS)-format to the standard multi-input-
  3. %    singe-output THETA-format. Auxiliary routine to th2ff and th2zp
  4. %
  5. %    TH = thss2th(TH_SS,IY)
  6. %
  7. %    TH_SS: The model defined in the TH(SS) format (See help th(ss))
  8. %    IY:  The output number in TH_SS to be considered as the
  9. %         output in TH (The THETA-format handles only single-
  10. %         output models). Default IY=1
  11. %    TH:  The corresponding model in the THETA-format including
  12. %         the transformed covariance matrix (see help theta)
  13. %    
  14. %    The tranformation of the covariance matrix is carried out using
  15. %    the Gauss' approximation formula with numerical differentiation.
  16. %    The stepsize in the differentiation is decided by the
  17. %    m-file nuderst.
  18.  
  19. %    L. Ljung 10-2-90
  20. %    Copyright (c) 1990 by the MathWorks, Inc.
  21. %    All Rights Reserved.
  22.  
  23. if nargin<2, iy=1;end
  24. [nr,nc]=size(eta);
  25. T=eta(1,2);nu=eta(1,3);ny=eta(1,4);
  26.  
  27. sspm=getmfth(eta);
  28. nd=eta(1,5);
  29. [etapar,p,lambda]=th2par(eta);if length(etapar)==0,etapar=0;end
  30. dt=nuderst(etapar);
  31. arg=getargth(eta);
  32. if any(eta(2,8)==[2 3 4]) & T<0,flag=1;else flag=0;end
  33. if eta(2,8)==3 & T<0,flaga=1;else flaga=0;end
  34. if any(eta(2,8)==[2 3]),Tmod=-1;else Tmod=T;end
  35. [a,b,c,d,k,x0]=feval(sspm,etapar,Tmod,arg);
  36. if flag, [a,bk,c,d]=contth(a,[b k],c,d,abs(T),flaga);
  37. b=bk(:,1:nu);k=bk(:,nu+1:nu+ny);end
  38. [ny,nx]=size(c);
  39. c=c(iy,:);
  40. if nu>0,d=d(iy,:);else d=[];end
  41. Apol=poly(a);
  42.  
  43. [nx,nu]=size(b);
  44. for kk=1:nu
  45.   Bpol(kk,:)=poly(a-b(:,kk)*c)+(d(kk)-1)*Apol;
  46.   bs=(abs(Bpol(kk,:))>eps);
  47.   Bpol(kk,:)=bs.*Bpol(kk,:);
  48. end
  49. Cpol=poly(a-k(:,iy)*c);
  50. if T>0,ia=max(find(abs(Cpol)>eps));Cpol=Cpol(1:ia);end
  51. if T>0,ia=max(find(abs(Apol)>eps));Apol=Apol(1:ia);end
  52.  
  53. Ncap=getncap(eta);lambb=lambda(iy,iy);
  54. th=poly2th(Apol,Bpol,Cpol,1,ones(nu,1),lambb,T);
  55. if ~isempty(Ncap),
  56.    ndtemp=length(th2par(th));th(2,1)=lambb*(1+ndtemp/Ncap)/(1-ndtemp/Ncap);
  57. else th(2,1)=lambb;
  58. end
  59. if nu>0
  60.     na=th(1,4);nb=th(1,5:4+nu);nc=th(1,5+nu);
  61.     ndd=th(1,6+nu);nf=th(1,7+nu:6+2*nu); nk=th(1,7+2*nu:6+3*nu);
  62. else na=th(1,4);nb=0;nc=th(1,5);ndd=th(1,6);nf=0;nk=0;
  63. end
  64.  
  65. if ~isempty(p), if norm(p)>0
  66.   for kl=1:nd
  67.       th1=etapar;th1(kl)=th1(kl)+dt(kl)/2;
  68.       th2=etapar;th2(kl)=th2(kl)-dt(kl)/2;
  69.       [A1,B1,C1,D1,K1,X1]=feval(sspm,th1,Tmod,arg);
  70.       [A2,B2,C2,D2,K2,X2]=feval(sspm,th2,Tmod,arg);
  71.       if flag [A1,BK1,C1,D1]=contth(A1,[B1 K1],C1,D1,abs(T),flaga);
  72.      [A2,BK2,C2,D2]=contth(A2,[B2 K2],C2,D2,abs(T),flaga);
  73.      if nu>0,B1=BK1(:,1:nu);B2=BK2(:,1:nu);end
  74.      K1=BK1(:,nu+1:nu+ny);K2=BK2(:,nu+1:nu+ny);
  75.       end
  76.       Apol1=poly(A1);Apol2=poly(A2);C1=C1(iy,:);C2=C2(iy,:);
  77.       if nu>0,D1=D1(iy,:);D2=D2(iy,:);end
  78.       for kk=1:nu
  79.     B1pol(kk,:)=poly(A1-B1(:,kk)*C1)+(D1(kk)-1)*Apol1;
  80.     B2pol(kk,:)=poly(A2-B2(:,kk)*C2)+(D2(kk)-1)*Apol2;
  81.       end
  82.       C1pol=poly(A1-K1(:,iy)*C1);C2pol=poly(A2-K2(:,iy)*C2);
  83.       if na>0, diff(kl,1:na)=(Apol1(2:na+1)-Apol2(2:na+1))/dt(kl);end
  84.       if nc>0, diff(kl,na+sum(nb)+1:na+nc+sum(nb))=... 
  85.       (C1pol(2:nc+1)-C2pol(2:nc+1))/dt(kl);end
  86.  
  87.       sb=na;sf=na+sum(nb)+nc+ndd;
  88.       for k=1:nu
  89.         if nb(k)>0,diff(kl,sb+1:sb+nb(k))=...
  90.            (B1pol(k,nk(k)+1:nk(k)+nb(k))-B2pol(k,nk(k)+1:nk(k)+nb(k)))/dt(kl);
  91.         end
  92.         sb=sb+nb(k);
  93.    
  94.       end
  95.    end
  96.  
  97.    PTH=diff'*p*diff;
  98.    [d,nd]=size(diff');
  99.    if nd>0,th(4:d+3,1:d)=PTH;end
  100. end,end