function [A,se,resid,r2,rb2,y,XXa]=autoreg(X,p,nocnst);
%AUTOREG  Ordinary least-square estimation for an autoreggressive process.
%         [A,se,resid,r2,rb2]=autoreg(X,p,nocnst) calculates the 
%         ordinary least squares estimates, A, for the equation(s):
% 
%              X(t)=A0+A1*X(t-1)+A2*X(t-2)+... Ap*X(t-p)+u(t),
%             n x 1                   
%
%         where the input X is a T x n matrix and E[u]=0, E[u*u']=Sigma.
%         AUTOREG returns the estimates, A=[A0,..Ap], their standard 
%         errors, equation residuals, the R^2 statistic for each regression,
%         and the Rbar^2 statistic for each regression.
% 
%         [A,se,resid,r2,rb2]=autoreg(X,p) can be used for the default
%         case in which a constant term is desired.
% 
%         See also OLS.

%         Ellen R. McGrattan, 7-1-89
%
if nargin<3; nocnst=0; end;
[T,n]=size(X);
if  T<=p; error('Too few observations in X for AUTOREG'); end;
if nocnst; 
  XX=[ ];
else;
  XX=ones(T-p,1);
end;
for i=1:p;
  XX=[XX,X(p-i+1:T-i,:)];
end;
y=X(p+1:T,:);
A=(XX'*XX)\(XX'*y);
if nargout>1;
  resid=y-XX*A;
  [T,k]=size(XX);
  df=T-k;
  se=sqrt(diag(inv(XX'*XX)/df))*sqrt(diag(resid'*resid))';
end;
if nargout>3;
  p=XX*((XX'*XX)\XX');
  l=eye(T)-ones(T)/T;
  r2=((diag(y'*XX*A)-(sum(y)').^2/T).^2)./(diag(y'*l*y).*diag(y'*p*l*p*y));
  rb2=(T*r2-k)/df;
end;
A=A';
se=se';
XXa=XX*A';
