function [bhat,vcov,rsq]=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);
rsq=1-(sum(res.^2)/sum((y-mean(y)).^2));
[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-j);  h=h*(l+1-j)/(l+1);
g=g+h+(h');
end;
vcov=n*inv(x'*x)*g*inv(x'*x);
end;

function [bhat,vcov,rsq]=olswhite(y,x);
%OLS with White standard errors
bhat=inv(x'*x)*x'*y;
res=y-(x*bhat);
rsq=1-(sum(res.^2)/sum((y-mean(y)).^2));
[n,p]=size(x);
m1=zeros(p,p);
for i=1:n; m1=m1+((res(i)^2)*(x(i,:)'*x(i,:))); end;
vcov=inv(x'*x)*m1*inv(x'*x);
