function TH=iv4mv(z,nn,maxsize,T,p)
%IV4MV	Computes approximately optimal IV-estimates for multivariate ARX-models
%	Auxiliary routine to iv4
%	TH = iv4mv(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 each output, input as column
%	vectors: Z=[y1 y2 ... yp u1 u2 ... un]
%
%	NN: NN=[na nb nk], giving the orders and delays of the above model.
%
%	Some parameters associated with the algorithm are accessed by
%	TH = iv4mv(Z,NN,maxsize,T)
%	See HELP AUXVAR for an explanation of these and their default values.

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

% *** Set up default values ***
[nnr,nnc]=size(nn);
ny=nnr;nu=(nnc-ny)/2;
if nu<=0,error('This routine is not intended for time series. Use ARX or AR instead!'),end
[Ncap,nz]=size(z); 
if nz>Ncap,error(' Data should be organized in column vectors!'),return,end
if nz~=(nu+ny),disp(' Incorrect number of orders specified!'),disp('For an AR-model nn=na'),disp('For an ARX-model, nn=[na nb nk]'),
error('see above'),return,end
na=nn(:,1:ny);
if nu>0,nb=nn(:,ny+1:ny+nu);else nb=zeros(ny,1);end
if nu>0,nk=nn(:,ny+nu+1:nnc);else nk=zeros(ny,1);end

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

% *** Do some consistency tests ***


if nz>Ncap, error('The data should be organized in column vectors')
return,end
%
% *** First stage: compute an LS model ***
%
th=mvarx(z,nn,maxsize,T);
a=th2ss(th);
maxa=max(abs(eig(a)));
if maxa>1, %Stabilize! 
	[a,b]=th2arx(th);[nar,nac]=size(a);
	for ka=1:(nac/ny)-1
	a(:,ka*ny+1:(ka+1)*ny)=((0.9/maxa)^ka)*a(:,ka*ny+1:(ka+1)*ny);
	end
  th=arx2th(a,b,ny,nu);
end  
%
% *** Second stage: Compute the IV-estimates using the LS
%     model a, b for generating the instruments ***
%
yh=idsimss(z(:,ny+1:nz),th);
for outp=1:ny
if norm(yh(:,outp),inf)<eps,
error(['Output number ' int2str(outp) ' is decoupled from all inputs. IV4 cannot be used!']),end,end
th=ivxmv(z,nn,yh,maxsize,T,0); 
a=th2ss(th);th1=th;
maxa=max(abs(eig(a)));
if maxa>1, %Stabilize! 
	[a,b]=th2arx(th);[nar,nac]=size(a);
	for ka=1:(nac/ny)-1
	a(:,ka*ny+1:(ka+1)*ny)=((0.9/maxa)^ka)*a(:,ka*ny+1:(ka+1)*ny);
	end
  th1=arx2th(a,b,ny,nu);
end  

%
% *** Third stage: Compute the residuals, v, associated with the
%     current model, and determine an AR-model, A, for these ***
%
v=pe(z,th);if ny>1, v=sum(v')';end
art=arx(v,(max(sum(na))+max(sum(nb))),maxsize,T,0);

%
% *** Fourth stage: Use the optimal instruments ***
%
for kz=1:nz
zf(:,kz)=filter([1 art'],1,z(:,kz));
end
xf=idsim(zf(:,ny+1:nz),th1);
TH=ivxmv(zf,nn,xf,maxsize,T,1);
TH(2,7)=34;
