% proc for OLS regression with Newey-West standard errors @
function [b_vee, s_vee, sig_vee] = nw(y_vee,x_vee,lag_vee)
% 	@ inputs:
% 		y_vee = (capt_vee x 1) vector of observations on dependent variable
% 		x_vee = (capt_vee x k_vee) matrix of observations on explanatory variables
% 		lag_vee = positive integer indicating lag length 
% 			     (j_vee = 0 means White standard errors)
% 	   outputs:
% 		b_vee = (k_vee x 1) vector of OLS regression coefficients
% 		s_vee = (k_vee x 1) vector of Newey-West standard errors
% 		sig_vee = (k_vee x k_vee) matrix of Newey-West coefficient variances @

% do OLS regression and calculate OLS resid matrix @
	b_vee = lsq(y_vee,x_vee);
	eps_vee = y_vee - x_vee*b_vee;
    
    %eps_vee = eps_vee .* x_vee;  does not work in Matlab
    eps_vee = matmul(eps_vee, x_vee);
	
% calculate the S matrix @
	sig_vee = autocov(eps_vee,0);
	j_vee = 1;
	while j_vee <= lag_vee
		sig_j = autocov(eps_vee,j_vee);
		sig_j = (sig_j + sig_j')*(1 - (j_vee/(lag_vee + 1)));
		sig_vee = sig_vee + sig_j;
        j_vee = j_vee + 1;
    end

% get the variance matrix @
	sig_j = inv(x_vee'*x_vee);
	sig_vee = length(y_vee)*sig_j * sig_vee * sig_j;
	s_vee = sqrt(diag(sig_vee));
end