function [bhat,se_beta] = hac_simple(x,y,m)
% Compute OLS coefficients and Newey-West HAC SEs.

[nobs,nvar]=size(x);
bhat=inv(x'*x)*x'*y;
e=y-x*bhat;

v=x.*(e*ones(1,nvar));

z=zeros(nvar);
for i=-m:m
    k=1-(abs(i)/(m+1));
    z=z+k*v(1+abs(i):nobs,:)'*v(1:nobs-abs(i),:);
end

v_beta=diag(inv(x'*x)*z*inv(x'*x)');

se_beta = sqrt(v_beta);

end