OLS regression with Newey-West error term I need to estimate a linear regression with the OLS method. Since I assume the error terms to be correlated, I would like to account for heteroskedasticity and autocorrelation in the error terms. I already wrote some code, but I wonder whether it is correct. 
So far, my code is:
myregression <- lm(x ~ y + z)
coeftest(myregression, vcov=NeweyWest(myregression))

Any help is highly appreciated!
 A: This question is about code but seeing as I've been looking at HAC estimates recently in R I will "answer".


*

*I have not checked the R implementation of Newey-West is exactly as in their original paper.  However, your code does indeed calculate R's NeweyWest HAC estimate using the default bandwidth selection/lag method. (You can view this parameter with the "verbose=T" option.)

*If you know the form of the correlations in your data then you can take less of a "sledgehammer" approach than Newey-West.  E.g. Prais-Winsten or Cochrane-Orcutt.

*Be aware that serial correlation is being examined here and so the order that your observations are sorted in does matter.
If you are interested you might consider a toy example where you generate correlated residuals on purpose to see how the Newey-West std error estimates/p-values differ.  In the example below the correlated residuals make it "look like" the response can be fitted against a straight line. The lag parameter is automatically (and correctly) chosen = 1 (seen with verbose=T option).  If you play with the process generating the residuals you can see how it changes.
set.seed(04012017)
n<-50
correlated_residuals<-arima.sim(list(ar = .9), n)
y<-correlated_residuals
x<-1:n
plot(x,correlated_residuals)
fit<-lm(y~x)
abline(fit)
summary(fit) # standard estimates
coeftest(fit,vcov=NeweyWest(fit,verbose=T))

