function [sbc,msc,fpe,logdp,np]=arord(R,m,mcor,ne,pmin,pmax)
%ARORD	Order selection criteria for AR model.
%
%	[SBC,MSC,FPE]=ARORD(R,m,mcor,ne,pmin,pmax) returns approximate
%	values for the order selection criteria SBC, MSC, and FPE for
%	models of order pmin:pmax. The input matrix R is the upper
%	triangular factor in the QR factorization of the AR model, m
%	is the dimension of the state space, the flag mcor indicates
%	whether or not an intercept vector is being fitted, and ne is
%	the number of block equations of size m that are used in the
%	estimation. 
%
%	As a lower-level routine, ARORD is called by AR. 

%	For testing purposes, ARORD also returns the vectors logdp and
%	np, containing the logarithms of the determinants of the (scaled)  
%	covariance matrix estimates and the number of parameters vectors 
%	at each order pmin:pmax.

%	Modified 6-Sep-97
%	Author: Tapio Schneider
%		AOS Program, Princeton University 
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU

imax = pmax-pmin+1;		% maximum index of output vectors

% initialize output vectors
sbc     = zeros(1, imax); 	% Schwarz' Bayesian Criterion
msc     = zeros(1, imax);	% Modified Schwarz Criterion
fpe     = zeros(1, imax); 	% log of Akaike's Final Prediction Error
logdp   = zeros(1, imax);	% determinant of (scaled) covariance matrix
np      = zeros(1, imax); 	% number of parameter vectors of length m
np(imax)= m*pmax+mcor;

% Get lower right triangle R22 of R: 
%
%   | R11  R12 |
% R=|          |
%   | 0    R22 |
%
R22 = R(np(imax)+1 : np(imax)+m, np(imax)+1 : np(imax)+m);

% From R22, get (scaled) inverse of noise covariance matrix for model
% of order pmax 
invR22 = inv(R22);
Mp     = invR22*invR22';

% For order selection: determinant of covariance matrix (not bias corrected) 
% scaled by ne:
%            logdp = log det(noise covariance matrix * ne)
logdp(imax) = 2.*log(abs(prod(diag(R22))));

% Compute approximate order selection criteria for models of 
% order pmin:pmax
i = imax;
for p = pmax:-1:pmin
   np(i) = m*p + mcor; 	% number of parameter vectors of length m
   if p < pmax
      % Update determinant of noise covariance matrix
      % Rp: Part of R to be newly added to covariance matrix
      Rp = R(np(i)+1:np(i)+m, np(imax)+1:np(imax)+m);

      % Get Mp, the updated (scaled) inverse of the noise covariance
      % matrix, using the Sherman-Morrison-Woodbury formula
      L    = chol(eye(m) + Rp*Mp*Rp')';
      N    = L \ Rp*Mp;
      Mp   = Mp - N'*N;

      % Get updated logarithm of determinant
      logdp(i) = logdp(i+1) + 2.* log(abs(prod(diag(L))));
   end

   % Schwarz' Bayesian criterion
   sbc(i) = logdp(i)/m - log(ne) * (ne-np(i))/ne;

   % Modification of Schwarz' criterion (MSC):
   msc(i) = logdp(i)/m - (log(ne) - 2.5) * (1 - 2.5*np(i)/(ne-np(i)));

   % logarithm of Akaike's final prediction error estimate
   fpe(i) = logdp(i)/m - log(ne*(ne-np(i))/(ne+np(i)));
   i = i-1;
end

