function th=pem(z,nn,index,maxiter,tol,lim,maxsize,Tsamp)
%PEM	Computes the prediction error estimate of a general linear model
%	TH = pem(Z,THSTRUC)
%
%	TH : returned as the estimated parameters of the model
%	along with estimated covariances and structure information. 
%	For the exact format of TH see HELP THETA.
%
%	Z : the output input data [y u], with y and u as column vectors
%	For multi-variable systems, Z=[y1 y2 ... yp u1 u2 ... un]
%
% 	THSTRUC: A matrix that defines the model structure. Typically created
%	by POLY2TH, MS2TH or MF2TH or by some estimation routine. The minimiza-
%	tion is then initialized at the parameters contained in THSTRUC. 
%	A general, black-box multi-input structure 
%	A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
%	is obtained for	THSTRUC = [na nb nc nd nf nk], indicating the orders
%	in the above model. 
%	By TH=pem(Z,THSTRUC,INDEX) only parameters corresponding to the indices
%	in the row vector INDEX are estimated.
%	Some parameters associated with the algorithm are accessed by
%	TH = pem(Z,THSTRUC,index,maxiter,tol,lim,maxsize,T)
%	See HELP AUXVAR for an explanation of these and their default values.

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

%  *** Set up default values ***

[Ncap,nz]=size(z);
maxsdef=idmsize(Ncap);

if nargin<8, Tsamp=[];end
if nargin<7, maxsize=[];end
if nargin<6, lim=[];end
if nargin<5, tol=[];end
if nargin<4, maxiter=[];end
if nargin<3, index=[];end
[nr,cc]=size(nn);
if isempty(Tsamp),if nr>1,Tsamp=nn(1,2);else Tsamp=1;end,end
if isempty(maxsize),maxsize=maxsdef;end,if maxsize<0,maxsize=maxsdef;end
if isempty(lim),lim=1.6;end,if lim<0,lim=1.6;end
if isempty(tol),tol=0.01;end,if tol<0,tol=0.01;end
if isempty(maxiter),maxiter=10;end,if maxiter<0,maxiter=10;end
if nr>1,if isthss(nn)
	eval('th=pemss(z,nn,index,maxiter,tol,lim,maxsize,Tsamp);')
return,end,end 

if Tsamp<0, Tsamp=1;end

% *** Do some consistency tests ***

nu=nz-1;
if nz>Ncap, error('The data should be organized in column vectors!')
return,end
nb=nn(2:nz); if nz==1 
error('Use ARMAX instead for a time series!')
end
if nr==1 & cc~=3+3*nu,
    disp('Incorrect number of orders specified:')
    disp('You should have nn=[na nb nc nd nf nk],')
    disp('where nb, nf and nk are row vectors of length equal to')
    disp(' the number of inputs'),error('see above')
return,end


%*** if nn has more than 1 row, then it is a theta matrix
%    and we jump directly to mincrit ***
%
if nr>1,th=nn;end


if nr==1
   th=inival(z,nn,maxsize,Tsamp);
end
 
% *** Display the initial estimate ***

clc, disp([' INITIAL ESTIMATE'])
disp(['Current loss:' num2str(th(1,1))])
disp(['th-vector'])
nu=th(1,3); n=sum(th(1,4:6+2*nu));
disp(th(3,1:n)')

if maxiter==0,return,end

% *** Minimize the prediction error criterion ***

th=mincrit(z,th,index,maxiter,tol,lim,maxsize,Tsamp);


