function [out,R2,Epsilon]=Regression(Y,X,n_lag)
if nargin==2,
n_lag=5;
end
A=[ones(length(Y),1) X];
b=inv(A'*A)*(A'*Y);
Epsilon=Y-A*b;
n_obs=length(Y);
H=A.*kron(ones(1,size(A,2)),Epsilon);
W=H'*H/n_obs;
for j=1:n_lag,
S_j=H(1:end-j,:)'*H(j+1:end,:)/(n_obs-j);
W=W+(1-j/(1+n_lag))*(S_j+S_j');
end
SE_b=sqrt(n_obs*diag(inv(A'*A)*W*inv(A'*A)));
t_Stat=b./SE_b;
out=[b'; SE_b'; t_Stat'];
R2=1-var(Epsilon)/var(Y);
%R2=1-(size(Y,1)-1)/(size(Y,1)-1-size(X,2))*(1-R2);