function TH=iv4(z,nn,maxsize,T,p)
%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 estimted covariances and structure information.
%	For the exact format of TH see HELP 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. Similarily 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 HELP AUXVAR for an explanation of these and their default values.

%	L. Ljung 10-1-86,4-15-90
%	Copyright (c) 1986-90 by the MathWorks, Inc.
%	All Rights Reserved.

% *** 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

% *** 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>Ncap, error('The data should be organized in column vectors')
return,end
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 anb 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) ***
