function [b,s2,Vb,ess,tss] = ols(X,y)
% [b,s2,Vb,ess,tss] = ols(X,y) computes basic ordinary least squares statistics
%
%    X  is a matrix of explanatory variables, observations by row, variables by column
%    y  is a vector of dependent variable observations
%
%    b  is the fitted coefficient vector for the variables in X
%    s2 is the fitted variance
%    Vb is the estimated variance-covariance matrix of  b
%    ess is the error sum of squares
%    tss is the total sum of squares
%
%    collinear variables are assigned coefficients and variances equal to zero

% Paul A. Ruud, 2/25/00
% Department of Economics
% University of California, Berkeley
%
% This file supports An Introduction to Classical Econometric Theory, Oxford Univ. Press.
% Page references point to parts of this book.  It is meant for instruction, not as a
% substitute for specialized software.

[n,k] = size(X);

% find a basis for Col(X); see Section 2.5, Exact Multicollineary, and
% Section 2.6.1, Graham-Schmidt Orthonormalization

basis = logical([1;zeros(k-1,1)]);  % always include the first RHS (assuming it is nonzero)
for i=2:k
   Z = X(:,basis);
   Q = eye(n) - Z*inv(Z'*Z)*(Z');   % orthogonal projector onto orthogonal complement of Col(Z)
   chk = norm(Q*X(:,i))/norm(X(:,i)) >= 1e-10;  % this is a crude test criterion
   basis(i,1) = logical(chk);       % exclude any column that is
                                    %    (numerically) linearly dependent on Z
end

% fit coefficients to basis vectors for Col(X); see Section 2.4, Orthogonal Projection

Z     = X(:,basis);                 % final basis for Col(X)
c     = inv(Z'*Z)*Z'*y;             % OLS fitted coefficients (this could also be  c = Z\y)

b     = zeros(k,1);                 % expand coefficient vector to include zeros for
b(basis,1) = c;                     %    collinear RHS variables                    
                                    
mu    = X*b;                        % fitted values
ess   = (y-mu)'*(y-mu);             % error sum of squares
avgy  = mean(y);                    % sample mean of y
tss   = (y-avgy)'*(y-avgy);         % total sum of squares

% compute estimator of the conditional variance of each element of y conditional on X;
% see Section 8.2, Second-Moment Properties, Section 8.3, Variance and Covariance Matrices,
% and Section 8.4, Estimation of the Variance Parameter

s2    = ess/(n-k);                  % (sum of squared residuals)/(degrees of freedom)
Vb    = zeros(k,k);
Vb(basis,basis) = s2*inv(Z'*Z);     % estimated variance of fitted coefficients padded with
                                    %    zeros for coefficients of collinear RHS variables

return
