function [beta,se,resid,r2,rb2,sereg,logl,dw,aic,sc,fstat]=ols(X,y);
%OLS   Ordinary least-square estimation
%      beta=ols(X,y) calculates the ordinary least squares 
%      estimates, beta, for each of n equations of the form:
% 
%           y=X*beta+u,     E[u]=0, E[u*u']=sigma^2 * I
%     
%      where y is a Txn matrix of T observations on the n dependent 
%      variables, X is a Txk matrix of T observations on k independent
%      variables, and u is a Txn matrix of residuals.  If 3 arguments:
%      [beta,se,resid]=ols(X,y) are given then the standard errors of
%      the estimates and residuals (u) are also returned.  If more than
%      three arguments are requested then the program returns:
%      [beta,se,resid,r2,rb2,sereg,logl,dw,aic,sc,fstat]=ols(X,y) 
%      where r2 is the R-squared statistic, rb2 is the adjusted 
%      R-squared statistic, sereg is the standard error of the regression,
%      logl is the log of the likelihood, dw is the Durbin-Watson 
%      statistic testing the serial correlation in the residuals,
%      aic is the Aikake information criterion, sc is the Schwarz
%      criterion, and fstat is the F-statistic testing the hypothesis
%      that all slope coefficients (other than the constant) are zero.
%     
%

%      Ellen R. McGrattan, 10-13-87
%      Revised, 12-2-05, ERM
%
beta = (X'*X)\(X'*y);
if nargout>1;
  resid = y-X*beta;
  [T,k] = size(X);
  df    = T-k;
  se    = sqrt(diag(inv(X'*X)/df))*sqrt(diag(resid'*resid))';
end;
if nargout>3;
  l     = eye(T)-ones(T)/T;
  r2    = 1-diag(resid'*resid)./diag(y'*l*y);
  rb2   = 1-(1-r2)*(T-1)/df;
  sereg = sqrt(diag(resid'*resid)/df);
  logl  = -T/2*(1+log(2*pi)+log(diag(resid'*resid)/T));
  dw    = (sum((resid(2:T,:)-resid(1:T-1,:)).^2)./sum(resid.^2))';
  aic   = -2*logl/T+2*k/T;
  sc    = -2*logl/T+k*log(T)/T;
  fstat = r2/(1-r2)*df/(k-1);
end;
