function[beta, Cov, se, tstat, pbeta, R2, s2, e, r1, dw, pdw] = OLS(y, X)
% [beta, CovOLS, seOLS, pOLS, s2OLS, e, r1, dw, pdw] = OLS(y, X) 
% y = X * beta + e, y ~ n x 1.
% X ~ n x k, beta ~ k X 1.  beta(1) = const if X[:,1] = ones(n,1)
n = size(X,1);
k = size(X,2);
nmk = n - k;
XXi = inv(X'*X);
beta = XXi * X' * y;
e = y - X*beta;
ee = e'*e;
s2 = ee / nmk;     % unbiased if rho = 0
Cov = s2*XXi;           % OLS Cov matrix (assumes iid errors)
se = sqrt(diag(Cov));   % OLS se's (assumes iid errors)
tstat = beta./se;
pbeta = 2*(1-tcdf(abs(tstat),nmk));
ye = y - mean(y);
R2 = 1-ee/(ye'*ye);
r1 = e(1:n-1)'*e(2:n)/ee;
[pdw, dw] = dwtest(e, X);
return
% end of function OLS