function [b,se] = ols_het(y,X);
% ols_het.m
% OLS with heteroskedasticity-consistent standard errors
  
[n,k] = size(X);
XXinv = inv(X'*X);
b = XXinv*X'*y;
e = y - X*b;
Xtilde = (e*ones(1,k)).*X;
Delta = (1/n)*Xtilde'*Xtilde;
Var = n*XXinv*Delta*XXinv;
se = sqrt(diag(Var));