vcovHC, vcovHAC, NeweyWest – which function to use? I am trying to update my lm() based model to get correct standard errors and tests. I am really confused which VC matrix to use. The sandwich package offers vcovHC , vcovHAC and NeweyWest. While the former only accounts for heteroskedasticity the latter two account for both serial correlation and heteroskedasticity. Yet, the documentation does not tell much about the difference between the latter two (at least I don't get it). Looking to the function itself I realized that NeweyWest actually calls vcovHAC.
Empirically the results of coeftest(mymodel, vcov. = vcovHAC) and coeftest(mymodel, vcov. = NeweyWest)are mad different. While vcovHACis somewhat close to the naive lm results, using NeweyWest all coefficients turn insignificant (tests even close to 1). 
 A: The "sandwich" in question is two pieces of bread defined by the expected information enclosing a meat defined by the observed information. See my comments here and here. For a linear regression, the estimating equation is:
$$
U(\beta) = \mathbf{X}^T\left(Y - \mathbf{X}^T\beta\right)
$$
The expected information (bread) is:
$$
A = \frac{\partial U(\beta)}{\partial \beta} = -(\mathbf{X}^T\mathbf{X})
$$
The observed information (meat) is:
$$
B = E(U(\beta)U(\beta)^T) = \mathbf{X}^T(Y-\mathbf{X}^T\beta)(Y-\mathbf{X}^T\beta)^T\mathbf{X}
$$
Note the inner term is a diagonal of constant residuals when the homoscedasticity, independent data assumption is met, then the sandwich covariance estimator which is given by $A^{-1}BA^{-1}$ is the usual linear regression covariance matrix $\sigma^2 \left(\mathbf{X}^T\mathbf{X}\right)^{-1}$ where $\sigma^2$ is the variance of the residuals. However, that's rather strict. You get a considerably broader class of estimators by relaxing the assumptions involved around the $n \times n$ residual matrix: $$R = (Y-\mathbf{X}^T\beta)(Y-\mathbf{X}^T\beta)$$. 
The "HC0" vcovHC estimator is consistent even when the data are not independent. So I will not say that we "assume" the residuals are independent, but I will say that we use "a working independent covariance structure". Then the matrix $R$ is replaced by a diagonal of the residuals 
$$R_{ii} = (Y_i - \beta \mathbf{X}_{I.})^2, \quad 0\text{ elsewhere}$$
This estimator works really well except under small samples (<40 is often purported). The HC1-3 are various finite sample corrections. HC3 is generally the best performing.
However if there are autoregressive effects, the off-diagonal entries of $T$ are non-zero, so a scaled covariance matrix is produced based on commonly used autoregressive structures. This is the rationale for the "vcovHAC". Here, very flexible and general methods are produced to estimate the autoregressive effect: the details may be beyond the scope of your question. The "meatHAC" function is the general workhorse: the default method is Andrews'. Newey-West is a special case of the general autoregressive error estimator. These methods solve one of two problems: 1. at what rate does correlation decay between "adjacent" observations and 2. what is a reasonable distance between two observations? These If you have balanced panel data, this covariance estimator is overkill. You should use gee from the gee package instead specifying the covariance structure to AR-1 or similar.
As for which to use, it depends on the nature of the data analysis and the scientific question. I would not advise fitting all the types and picking the one that looks best, as it is a multiple testing issue. As I alluded to earlier, the vcovHC estimator is consistent even in the presence of an autoregressive effect, so you can use and justify a "working independence correlation model" in a variety of circumstances.
