function [beta,se,R2,mdl,pval]=ols_robust(Y,X)
mdl=fitlm(X,Y);
tmp= mdl.Coefficients.Variables;
beta=tmp(:,1);
R2=mdl.Rsquared.Adjusted;
T=rows(Y);
maxLag = floor(4*(T/100)^(2/9));
[EstCov,se] = hac(mdl,'bandwidth',maxLag+1,'display','full');
trat=beta./se;
pval= 2*(1-tcdf(abs(trat),T));