Documentation of iv4

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

Function Synopsis

TH=iv4(z,nn,maxsize,T,p)

Help text

IV4    Computes approximately optimal IV-estimates for ARX-models.
   TH = IV4(Z,NN)

   TH: returned as the estimate of the ARX model
   A(q) y(t) = B(q) u(t-nk) + v(t)
   along with estimated covariances and structure information.
   For the exact format of TH see also THETA.

   Z : the output-input data Z=[y u], with y and u as column vectors.
   For multivariable systems Z=[y1 y2 .. yp u1 u2 .. um].

   NN: NN = [na nb nk], the orders and delays of the above model.
   For multi-output systems, NN has as many rows as there are outputs
   na is then an ny|ny matrix whose i-j entry gives the order of the
   polynomial (in the delay operator) relating the j:th output to the
   i:th output. Similarly nb and nk are ny|nu matrices. (ny:# of outputs,
   nu:# of inputs).

   Some parameters associated with the algorithm are accessed by
   TH = IV4(Z,NN,maxsize,T)
   See also AUXVAR for an explanation of these and their default values.
   See also ARX, ARMAX, BJ, N4SID, OE, PEM.

Cross-Reference Information

Listing of function iv4

function TH=iv4(z,nn,maxsize,T,p)

%   L. Ljung 10-1-86,4-15-90
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.4 $  $Date: 1997/12/02 03:43:24 $
if nargin < 2
   disp('Usage: TH = IV4(Z,ORDERS)')
   disp('       TH = IV4(Z,ORDERS,MAXSIZE,T)')
   return
end

% *** Set up default values ***
[Ncap,nz]=size(z);
maxsdef=idmsize(Ncap);
if nargin<5, p=1;end
if nargin<4, T=1;end
if nargin<3, maxsize=maxsdef;end

if p<0,p=1;end, if T<0, T=1;end,if maxsize<0,maxsize=maxsdef;end
if isempty(T),T=1;end,if isempty(maxsize),maxsize=maxsdef;end
if any(any(nn<0))
   error('All orders must be non-negative.')
end
if nz>Ncap, error('The data should be organized in column vectors')
return,end
% *** Do some consistency tests ***
[nr,nc]=size(nn);if nr>1, eval('TH=iv4mv(z,nn,maxsize,T,p);'),return,end
[Ncap,nz]=size(z); nu=nz-1;

if nz==1, error('This routine does not make sense for a time series!')
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);
%
% *** First stage: compute an LS model ***
%
th=arx(z,[na nb nk],maxsize,T,0);
if na>0, a=fstab([1 th(1:na).']);
        else a=1;end
b=zeros(nu,max(nb+nk));
NBcum=cumsum([na nb]);
for k=1:nu, b(k,nk(k)+1:nk(k)+nb(k))=th(NBcum(k)+1:NBcum(k+1)).';,end
%
% *** Second stage: Compute the IV-estimates using the LS
%     model a, b for generating the instruments ***
%
th=iv(z,[na nb nk],a,b,maxsize,T,0);
if na>0,a=[1 th(1:na).'];
        else a=1;end
for k=1:nu, b(k,nk(k)+1:nk(k)+nb(k))=th(NBcum(k)+1:NBcum(k+1)).';,end
%
% *** Third stage: Compute the residuals, v, associated with the
%     current model, and determine an AR-model, A, for these ***
%
v=filter(a,1,z(:,1));
for k=1:nu, v=v-filter(b(k,:),1,z(:,k+1));,end
art=arx(v,na+sum(nb),maxsize,T,0);
Acap=[1 art.'];
%
% *** Fourth stage: Use the optimal instruments ***
%
yf=filter(Acap,[1],z(:,1));
for k=1:nu, uf(:,k)=filter(Acap,[1],z(:,k+1));,end
if p~=0, p=2;,end
if na>1,a=fstab(a);else a=1;end
TH=iv([yf uf],[na nb nk],a,b,maxsize,T,p);
TH(2,7)=4;
% *** Reference: Equations (15.221) - (15.26) in Ljung(1987) ***