function [b,se,s2] = ols(y,X)
[n K] = size(X);
Q = inv(X'*X);
b = Q*X'*y;
e = y - X*b;
ess = e'*e;
s2 = ess/(n-K);
Vb = s2*Q;
se = sqrt(diag(Vb));