function [bhat,vcov]=newey(y,x,l);
%OLS with HAC standard errors
if l==0; [bhat,vcov,rsq]=olswhite(y,x);
else;
bhat=inv(x'*x)*x'*y;
res=y-(x*bhat);
[n,p]=size(x);
for i=1:p; z(:,i)=x(:,i).*res; end;
g=(z'*z)/n;
for j=1:l;
h=(z(1:n-j,:)'*z(j+1:n,:))/n;  h=h*(l+1-j)/(l+1);
g=g+h+(h');
end;
vcov=n*n/(n-p)*inv(x'*x)*g*inv(x'*x);
end;
