function [b,varb,R2]= olsnw(y,x,lags);

%   Computes heteroskedasticity consistent OLS of y on x
%   and returns the coefficients,the standard errors,
%   the t-statistics and the p-values.  The number of lagged
%   terms in the covariance matrix for Newey-West is k.


[T,k] = size(x);
k=k+1;              % count constant 
x = [ones(T,1) x];  % add constant to x

b = x\y;


u = y - x*b;

uz = x.*(u*ones(1,k)) ;
w = uz'*uz/T;

if lags ~= 0;
           j = 1;
           lend = lags + 1;
           while j <= lags
              
                 rzuj = uz(1+j:T,:)'*uz(1:T-j,:);
                 rzuj = rzuj./T;
                 w = w + (1 - j./lend).*(rzuj+rzuj');
					  j = j + 1;                 
           end
 end

		xx = inv(x'*x/T);

      varb =  xx'*w*xx/T;
      
      
      std = sqrt(diag(varb));
      
      %tstat = b./std;
      
  
     
% Calculate R^2 
        temp = y-mean(y);
        R2 = 1- u'*u/(temp'*temp);
  
    
  %disp([b';std']);
  %disp(R2);