# VAR(p) Model in R with HAC estimator

I'm running a VAR model in R and found with several tests (arch.test, serial.test) that my model still contains heteroscedasticity and autocorrelation. I took the AIC lag order selection criteria, since I got least heteroscedasticity and autocorrelation compared to HQ, SC, FPE (I'm running a great number of VAR Models) I want to correct the heteroscedasticity and autocorrelation, if possible. Is the HAC estimator possible to use for a VAR Model. And if the answer is yes, how can I implement it in R? If the answer is no, is there any other possibility to get better estimations (other than higher lag order)?

Thank you in advance! Appreciate every help!

I think it is possible to use HAC with VAR.

Regarding implementation in R:

• There are multiple R packages offering VAR estimation, perhaps one of them has the relevant functionality?
• Otherwise, see packages for multivariate models (seemingly unrelated regression SUR and the like), perhaps they have the functionality. You would have to specify the model manually.
• Finally, you could do it entirely by hand. VAR can be estimated equation by equation by OLS, so you could just apply HAC equation by equation. This might be the easiest option.

Also, note that arch.test is about autoregressive conditional heteroskedasticity which is not the kind of heteroskedasticity that HAC is designed to deal with. In ARCH, the variance of the error $$\varepsilon$$ is not a function of the right-hand-side variables (lags of the dependent variable ($$X$$), so I think a HAC variance estimator would not account for it. You could model autoregressive conditional heteroskedasticity explicitly using a multivariate GARCH model (of which there are many varieties) e.g. using the rmgarch package in R (it allows specifying and estimating the VAR and the multivariate GARCH in one step). You would get more efficient estimates of the VAR coefficients and if you are not directly interested in the time-varying volatility, you could just ignore the fitted multivariate GARCH model afterwards.

• Thank you, Richard, for the answer! What is the difference between the heteroscedasticity that the HAC estimator tries to rule out and the heteroscedasticity that the arch.test is testing for? Also I tried to model the VAR by hand, because I did it with the vars package (VAR()) before. If I use dynlm() for the OLS I can use vcovHAC, and I get a correlation matrix. I am just not sure how to use this matrix then...
– Anna
Commented Mar 1, 2021 at 7:55
• @Anna, standard concept of heteroskedasticity (which HAC deals with): variance of the error $\varepsilon$ is a function of the right-hand-side variables so that $\mathbb{E}(\varepsilon^\top X^\top X\varepsilon)\neq\sigma^2 X^\top X$. Autoregressive conditional heteroskedasticity: variance of the error is a function of past errors where $\mathbb{E}(\varepsilon^\top X^\top X\varepsilon)=\sigma^2 X^\top X$ is not violated. Commented Mar 1, 2021 at 8:09
• Assume a bivariate VAR(1) and assume you are interested in testing $H_0\colon x\xrightarrow{Granger}y$. You need to test whether the coefficient in front of $x_{t-1}$ in the equation for $y_t$, call it $\varphi_{2,1}$ is $=0$. You have to find the diagonal element of the HAC covariance matrix corresponding to this coefficient, say $\hat\sigma_{j,j}$. You can then use it for constructing the test statistic $t=\frac{\hat\varphi_{2,1}}{\sqrt{\hat\sigma_{j,j}}}$. $t$ will be asymptotically $N(0,1)$ under $H_0$. Commented Mar 8, 2021 at 20:31
• Assuming a bivariate VAR(p), you need to test a joint hypothesis $H_0\colon \varphi_{2,1}=...=0$ where the dots contain coefficients in front of all the lags of $x_t$ (except for the first lag which is already represented by $\varphi_{2,1}$) in the equation of $y_t$. You can do an $F$-test using the HAC covariance matrix just like you would using the normal covariance matrix. You can look up the construction of the $F$ statistic in a textbook. Also, there is absolutely nothing to be sorry about :) Commented Mar 8, 2021 at 20:35
• Assuming a VAR with dimension above 2 ($x,y,z,\dots$) makes it harder, and your equation-by-equation HAC covariance matrices might not suffice -- unless you care about $x\xrightarrow{Granger}y$ where $x$ and $y$ are univariate time series just like above. But outside this simple case, you might need a system-wide HAC covariance matrix (one that represents all equations at once). Commented Mar 8, 2021 at 20:38