Getting understand HAC estimators Can you help me please with understanding HAC estimator ? I've searched whole internet about it and I didn't find any page which explains clearly algorithm of HAC. I would also see some mathematical explanation as well. As far as I understand we use HAC to solve problem with autocorrelation and heteroscedasticity, and to do so we use vcovHAC() function from sandwich package.
x <- sin(1:100)
y <- 1 + x + rnorm(100)
fm <- lm(y ~ x)
vcovHAC(fm)

what exactly is that matrix in output ? And what should I do next ? i.e. I'm not sure how can I include that in my model so that I can resolve issue with autocorrelation and heteroscedasticity. Can you please help me with understanding that ?
 A: In a linear model, we have $\hat\beta = (X^TX)^{-1}X^TY$.
A basic property of variances and matrices is that
$$\mathrm{var}[A^TY] = A^T\mathrm{var}[Y]A$$
So
$$\mathrm{var}[\hat\beta] = (X^TX)^{-1}X^T \mathrm{var}[Y] X(X^TX)^{-1}$$
It's usual when considering HAC  estimators to  break this into three pieces, two of which are the same, hence the name "sandwich"
$$\mathrm{var}[\hat\beta] = I^{-1}  H I^{-1}$$
I'm not going to do quite that. Instead, I'm going to write
$$\hat\beta-\beta = (X^TX)^{-1}X^T (Y-X\beta)$$
which works because $(X^TX)^{-1}(X^TX)\beta=\beta$
and note that
$$\hat\beta-\beta = \sum_{i=1}^n h_i(\beta)$$
where
$$h_i(\beta)=(X^TX)^{-1} x_i(y_i-x_i\beta).$$
These are the influence functions. They have mean zero and each one almost depends on only one $y_i$.  I say 'almost' because they all depend on $(X^TX)^{-1}$, but that is an average of $n$ observations and so is effectively constant for large $n$. If we were doing asymptotics, we'd replace it by its limiting value.
By basically the definition of covariance
$$\mathrm{var}\left[\sum_i h_i(\beta)\right] = \sum_{i,j} \mathrm{cov}[h_i(\beta),h_j(\beta)]$$
We know that $E[h_i(\beta)h_j(\beta)=\mathrm{cov}[h_i(\beta),h_j(\beta)]$ (because they have zero mean), and we might hope to estimate it by $h_i(\hat\beta)h_j(\hat\beta)$.
For any individual $(i,j)$ that's a  terrible  measurement, but it is approximately unbiased (it would be exactly unbiased if we evaluated it at the true $\beta$, but if we knew the true $\beta$ we wouldn't be doing any of this). Since it's (approximately) unbiased, we have a reasonable hope that the law of large numbers will turn the sum of these things into a good estimate of the sum of the variances.
Sadly, it doesn't. For a start, since we know $\sum_i h_i(\hat\beta)=0$ by construction, $\sum_{i,j} h_i(\hat\beta)h_j(\hat\beta)=0$.
However, we can rescue the estimator with a bias:variance tradeoff on the covariance terms.  Suppose we assume that $i$ indexes time and that observations well-separately in  time are very nearly independent. That's reasonable: an ARIMA model has exponentially decaying correlations.  We could then estimate $\mathrm{cov}[h_i(\beta),h_j(\beta)]$ by 0 if $|i-j|$ is large  enough, and use $h_i(\hat\beta)h_j(\hat\beta)$ when $|i-j|$ is small.
This does work.
It also works for spatial data and for various sparse correlation models. The proofs get a bit detailed, especially if you want nearly optimal conditions, because there is uniform convergence to  be proved. The general form of the result, though, is  fairly straightforward
Write ${\cal N}$ (neighbours) for the set of  $(i,j)$ such that $\mathrm{cov}[h_i(\beta),h_j(\beta)]$ is not small.  If

*

*$|{\cal N}|$ is much  smaller than $n^2$, (eg $O_p(n^{2-\delta})$)

*The sum of the true $\mathrm{cov}[h_i(\beta),h_j(\beta)]$ over pairs not in ${\cal N}$ is small (goes to zero)

then the HAC estimator is pretty good (is consistent). You can improve things a bit by not  taking a binary yes/no decision but instead taking $w_{ij} h_i(\hat\beta)h_j(\hat\beta)$ for some $0<w_{ij}<1$. Most of the HAC estimators do this.
That was all  for the linear model, but (apart from a bit of smoothness and moment assumptions) the only property of the  linear model we used was that each $h_i$ depends (approximately) only one observation, and  that the $h_i$ add up to $\hat\beta-\beta$. If you weaken the latter part of that to "add up to $\hat\beta-\beta$" plus error of smaller order, that's a definition of an influence function, and you can find them for generalised linear models, the Cox model,  and many parametric regression models and just use the same approach to get HC and HAC variance estimators.
The sandwich package also has some sophisticated improvements that give you slightly better performance (and noticeably better in small samples), but this is the basic idea.
A: The ideas the HAC estimators implemented in the sandwich package are explained in an accompanying paper that is also listed in the references in ?vcovHAC:

Zeileis A (2004).
"Econometric Computing with HC and HAC Covariance Matrix Estimators."
Journal of Statistical Software, 11(10), 1-17.
doi:10.18637/jss.v011.i10

The paper also provides further links to the relevant literature and explains what you can do with the estimated variance-covariance matrix in R. Typically, you plug it into functions that allow you to test the coefficients of your model based on Wald-type tests, e.g., coeftest(), coefci(), and waldtest() from the lmtest package or linearHypothesis(), Anova(), or deltaMethod() from the car package.
