Documentation of ivx

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

TH=ivx(z,nn,x,maxsize,Tsamp,p)

Help text

IVX    Computes instrumental variable estimates for ARX-models.

   TH = ivx(Z,NN,X)

   TH: returned as the IV-estimate of the ARX-model
   A(q) y(t) = B(q) u(t-nk) + v(t)
   along with relevant structure information. See HELP THETA for
   the exact structure of TH.

   Z : the output-input data Z=[y u], with y and u as column vectors.
   For multi-input systems u=[u1 u2 ... un]

   NN: NN=[na nb nk] gives the orders and delays associated with the
   above model.
   X : is the vector of instrumental variables. This should be of the
   same size as the y-vector (i.e. the first column of Z).

   See IV4 for good, automatic choices of instruments.
   A multioutput variant ig given by IVXMV
   TH=ivx(Z,NN,X,maxsize,T)
   allows access to some parameters associated with the algorithm.
   See HELP AUXVAR for an explanation of these.

Cross-Reference Information

Listing of function ivx

function TH=ivx(z,nn,x,maxsize,Tsamp,p)

%   L. Ljung 10-1-86, 4-12-87
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.3 $  $Date: 1997/12/02 03:42:56 $

[Ncap,nz]=size(z); nu=nz-1;[Nc,nx]=size(x);[nnr,nnc]=size(nn);
maxsdef=idmsize(Ncap);
%
% *** Some initial tests on the input arguments ***
%

if nargin<6, p=1;end
if nargin<5, Tsamp=1;end
if nargin<4, maxsize=maxsdef;end
if maxsize<0,maxsize=maxsdef;end
if Tsamp<0,Tsamp=1;end
if p<0,p=1;end
if isempty(Tsamp),Tsamp=1;end,if isempty(maxsize),maxsize=maxsdef;end
if nnr>1,eval('TH=ivxmv(z,nn,x,maxsize,Tsamp,p);'),return,end
if Ncap~=Nc | nx~=1,error('The x-vector should be a column vector with the same number of elements as z!'),return,end
if length(nn)~=1+2*nu,disp('Incorrect number of orders specified!'),
        disp('nn should be nn=[na nb nk]'),
        disp('where nb and nk are row vectors of length equal to the number of inputs'),error('see above'),return,end
na=nn(1);nb=nn(2:1+nu);nk=nn(2+nu:1+2*nu);n=na+sum(nb);
%
% construct regression matrix
%
nmax=max([na+1 nb+nk])-1;
M=floor(maxsize/n);
Rxx=zeros(na);Ruu=zeros(sum(nb));Rxu=zeros(na,sum(nb));Rxy=zeros(na);
Ruy=zeros(sum(nb),na); F=zeros(n,1);
for k=nmax:M:Ncap-1
     jj=(k+1:min(Ncap,k+M));
     phix=zeros(length(jj),na); phiy=phix; phiu=zeros(length(jj),sum(nb));
     for kl=1:na, phiy(:,kl)=-z(jj-kl,1); phix(:,kl)=-x(jj-kl); end
     ss=0;
     for ku=1:nu
            for kl=1:nb(ku), phiu(:,ss+kl)=z(jj-kl-nk(ku)+1,ku+1);end
             ss=ss+nb(ku);
     end
     Rxy=Rxy+phix'*phiy;
     if nu>0,Ruy=Ruy+phiu'*phiy;
        Rxu=Rxu+phix'*phiu;
        Ruu=Ruu+phiu'*phiu;
     end
     Rxx=Rxx+phix'*phix;
     if na>0, F(1:na)=F(1:na)+phix'*z(jj,1);end
      F(na+1:n)=F(na+1:n)+phiu'*z(jj,1);
end
clear phiu, clear phix, clear phiy,
%
% compute estimate
%
if nu==0,TH=Rxy\F;end
if nu>0,TH=[Rxy Rxu;Ruy Ruu]\F;end
if p==0, return,end
%
% proceed to build up THETA-matrix
%
t=TH;, clear TH;
I=[0 Tsamp nu na nb 0 0 zeros(1,nu) nk];
n=na+sum(nb);
TH=zeros(n+3,max(length(I),n));
TH(1,1:length(I))=I;
ti=fix(clock); ti(1)=ti(1)/100;
TH(2,2:6)=ti(1:5);
TH(2,7)=11;
TH(3,1:n)=t.';
if p==2
e=pe(z,TH); TH(1,1)=e'*e/(length(e)-max(na,sum(nb)));
TH(2,1)=TH(1,1)*(1+n/Ncap)/(1-n/Ncap);
 TH(4:3+n,1:n)=TH(1,1)*inv([Rxx Rxu;Rxu.' Ruu]);end