function [b, s2, se, vcv] = myols(y,x);
b = inv(x'*x)*x'*y;
e = y-x*b;
[N,K] = size(x);
s2 = e'*e/(N-K);
vcv = s2*inv(x'*x);
se = sqrt(diag(vcv));