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_i]=0, E[u_i*u_i']
%         =sigma_i^2 * I.   AUTOREG returns the estimates, A=[A0,..Ap],
%         their standard errors, the R^2 statistic for each regression,
%         and the Rbar^2 statistic for each regression.
%         If p is a n-dimensional vector, p(i) is the lag length assumed
%         for the ith variable in X and A=[A0,..Am], m=max(p).
% 
%         [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);
m=max(p);
if n>length(p); 
  p=ones(n,1)*p(1);
end;
if  T<=m; error('Too few observations in X for AUTOREG'); end;
if nocnst; 
  XX=[ ];
else;
  XX=ones(T-m,1);
end;
for i=1:m;
  XX=[XX,X(m-i+1:T-i,p>=i)];
end;
y=X(m+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;
  m=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'*m*l*m*y));
  rb2=(T*r2-k)/df;
end;
A=A';
se=se';
XXa=XX*A';
