arch.test tells you the residuals display autoregressive conditional heteroskedasticity. To properly account for it, you can model the conditional variance of your time series with a multivariate GARCH model, e.g. BEKK-GARCH or DCC-GARCH; the latter is available in the
rmgarch package in R. There you may
- specify (
dccspec) the conditional mean model as VAR (knowing that VECM has an equivalent representation as VAR) alongside the conditional variance model;
- fit the model (
- produce forecasts from them in an automated way (
Check out the
rmgarch vignette (p. 3-6) and the help files.
The tricky part will be specifying the VAR model so that it obeys the restrictions due to the cointegrating relationship, i.e. so that it corresponds to the VEC model. To avoid this, you may fit your model in two stages: first the conditional mean just the way you did it, and then the conditional variance of the model residuals. You will get the point forecasts just as before and you will have to add the confidence bounds based on the conditional variance forecasts from the conditional variance model. I am not aware of any R commands that would do this directly, so I am afraid you will have to do a little bit of programming yourself.