function [w,A,C,sbc,msc,fpe,th]=ar(v,pmin,pmax,selector,no_const)
%AR	Least squares estimation for multivariate AR process.
%
%	[w,A,C,SBC,MSC,FPE,th]=AR(v,pmin,pmax) produces parameter 
%	estimates for an AR process of order p, 
%
%        v(k,:)' = w' + A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C),
%
%	where p lies between pmin and pmax and is chosen as the 
%	minimizer of a modified Schwarz order selection criterion.
%	The input matrix v must contain the time series data in form 
%	of row vectors.  AR returns least squares estimates for
%	the intercept vector w, for the coefficient matrices
%	A1,...,Ap (as A=[A1 ... Ap]), and for the noise covariance 
%	matrix C.
%
%	As order selection criteria, AR computes approximations to
%	Schwarz' Bayesian Criterion, to a modification of Schwarz' 
%	criterion, and to the logarithm of Akaike's final prediction 
%	error estimate. The order selection criteria for models of 
%	order pmin:pmax are returned as the vectors SBC, MSC, and FPE, 
%	respectively. 
%
%	On output, the matrix th contains information needed in the 
%	computation of confidence intervals for the estimated
%	parameters. ARMODE and ARCONF require th as input argument.
%       
%	If the optional argument SELECTOR is included in the function
% 	call as in AR(v,pmin,pmax,SELECTOR), SELECTOR is used as the
%	order selection criterion to determine the optimum process
%	order. The three letter string SELECTOR must have one of 
%	the three values 'sbc', 'msc', or 'fpe'. (By default, the
%	modified Schwarz criterion MSC is used.) Choosing the bounds
%       for the order to coincide, i.e., pmin=pmax, amounts to 
%	circumventing the order selection procedure and fixing the 
%	order of the estimated model.     
%
%	If the function call contains the optional argument 'zero' as
%	fourth or fifth argument, then a model of the form
%
%           v(k,:)' = A1*v(k-1,:)' +...+ Ap*v(k-p,:)' + noise(C). 
%
%	is fitted to the data. That is, the intercept vector w is
%	taken to be zero, which amounts to assuming that the AR(p)
%	process has zero mean.

%	Modified 6-Sep-97
%	Author: Tapio Schneider
%		AOS Program, Princeton University 
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU
%
%               and  
%
%               Arnold Neumaier
%               Institut fuer Mathematik, Universitaet Wien
%               Strudlhofgasse 4, A-1090 Wien, Austria
%               neum@cma.univie.ac.at

% n: number of time steps; m: dimension of state vectors
[n,m] = size(v);     

if (pmin ~= round(pmin) | pmax ~= round(pmax))
   error('Order must be integer.');
end
if (pmax < pmin)
   error('PMAX must be greater than or equal to PMIN.')
end

% set defaults and check for optional arguments
if (nargin == 3)        % no optional arguments => set default values
   mcor     = 1; 	% fit intercept vector
   selector = 'msc';	% use MSC as order selection criterion
elseif (nargin == 4)    % one optional argument present
   selector = lower(selector);
   if strcmp(selector, 'zero')
      mcor     = 0;     % no intercept vector fitted
      selector = 'msc';	% default order selection 
   else
      mcor     = 1;
   end
elseif (nargin == 5)
   selector = lower(selector);
   if strcmp(lower(no_const), 'zero')
      mcor   = 0;       % no intercept vector fitted
   else
      error(['Bad argument. Usage: ', ...
              '[w,A,C,SBC,MSC,FPE,th]=AR(v,pmin,pmax,SELECTOR,''zero'')'])
   end
end

ne    = n-pmax;		% number of block equations of size m
npmax = m*pmax+mcor;	% maximum number of parameter vectors of length m

if (ne <= npmax)
    error('Time series too short.')
end

% compute QR factorization for model of order pmax
[R, scale]=arqr(v,pmax,mcor);

% compute approximate order selection criteria for models 
% of order pmin:pmax
[sbc,msc,fpe]= arord(R,m,mcor,ne,pmin,pmax);

% get index iopt of order that minimizes the order selection 
% criterion specified by selector
[val, iopt] = min(eval(selector)); 

% select order of process
p     = pmin : pmax;
popt  = p(iopt);                % optimum order 
np    = m*popt + mcor;          % number of parameter vectors of length m

% decompose R for the optimal model order popt according to 
%
%   | R11  R12 |
% R=|          |
%   | 0    R22 |
%
R11 = R(1:np, 1:np);
R12 = R(1:np, npmax+1:npmax+m);    
R22 = R(np+1:npmax+m, npmax+1:npmax+m);

% get augmented parameter matrix Aaug=[w A] if mcor=1 and Aaug=A if mcor=0
if (np > 0)   
   if (mcor == 1)
      % improve condition of R11 by re-scaling first column
      con = max(scale(2:npmax+m))/scale(1); 
      R11(:,1) = R11(:,1)*con; 
   end;
   Aaug = (R11\R12)';

%  return coefficient matrix A and intercept vector w separately
   if (mcor == 1)
%     intercept vector w is first column of Aaug, rest of Aaug is 
%     coefficient matrix A
      w = Aaug(:,1)*con; 	% undo condition-improving scaling
      A = Aaug(:,2:np);
   else
%     return zeros as intercept vector 
      w = zeros(m,1);
      A = Aaug;
   end
else
   % no parameters have been estimated 
   % => return only covariance matrix estimate and order selection 
   % criteria for ``zeroth order model''  
   w = zeros(m,1);
   A = [];
end

% return covariance matrix
dof = ne-np;            % number of block degrees of freedom
C   = R22'*R22./dof; 	% bias-corrected estimate 

% for later computation of confidence intervals return in th: 
% (i)  the inverse of U=R11'*R11, which appears in the estimator's 
%      asymptotic covariance matrix
% (ii) the number of degrees of freedom in the noise covariance matrix 
%      estimate
invR11 = inv(R11);
Uinv   = invR11*invR11';
th     = [dof zeros(1,size(Uinv,2)-1); Uinv];



