Over rolling windows you can update your posteriors by re-estimating the model every 50 observations (or however often you want to do it). This can easily be accomplished in a for loop.
Posterior samples, forecasts, and other quantities can get re-generate with each re-estimation of the model. At the end of the day it is just a lot of loops. You may also want to consider running the models in parallel (like with the foreach and doParallel libraries) to speed up computation time if that is an issue for you.
Unfortunately, I do not know of any R packages that provide functions for automatically computing rolling window forecasts in a Bayesian setting. For the most part I think you have to program it yourself. It is a little bit tedious figuring out how index and store everything but the good news is that it is not that hard to do.
Though I do not think you need it, I give more intuition below (which does not differ much from the intuition I gave in the previous post you referenced).
I Assume your test set consists of $r_1,r_2,r_3,...,r_t$ and that you want to generate forecasts for an out of sample window $r_{t+1},r_{t+2},...,r_{t+K}$.
In a Bayesian setting the forecasts are treated as random samples with their own joint posterior distribution. So instead of producing point forecasts, you will draw a sample from the forecast distribution;
$$P(r_{t+1},r_{t+2},...,r_{t+K}|r_1,r_2,...,r_t,\boldsymbol{\theta}_t)$$ where $\boldsymbol{\theta}_t$ is just my shorthand notation for a vector containing the model parameters.
You can sample form the above forecast distribution with a Gibbs sampler. If your posterior is made up of $G$ samples, then for each $g=1,2,3,...,G$ you can draw
$r^{(g)}_{t+1} \sim P(r_{t+1}|r_1,...,r_t,\boldsymbol{\theta}^{(g)}_t)$
then
$r^{(g)}_{t+2}\sim P(r_{t+2}|r_1,...,r_t,r^{(g)}_{t+1},\boldsymbol{\theta}^{(g)}_t)$
and keep going until you sample
$r^{(g)}_{t+K} \sim P(r_{t+K}|r_1,...,r_t,r^{(g)}_{t+1},r^{(g)}_{t+K-1},\boldsymbol{\theta}^{(g)}_t)$
for any $k=1,2,...,K$ the conditional posterior $P(r_{t+k}|r_1,...,r_t,r^{(g)}_{t+1},r^{(g)}_{t+k-1},\boldsymbol{\theta}^{(g)}_t)$ can be sampled from in the following manner
$$
\omega_{t+k}^{(g)} \sim IG\bigg(\frac{v^{(g)}}{2},\frac{v^{(g)}}{2}\bigg)\;\;\;\;\;\varepsilon_{t+k}^{(g)} \sim N(0,1)
$$
$$
r_{t+k}^{(g)} = \varepsilon_{t+k}^{(g)}\times \bigg(\frac{v^{(g)}-2}{v^{(g)}} \omega_{t+k}^{(g)}h_{t+k}^{(g)}\bigg)^{1/2}
$$
$h^{(g)}_{t+k+1}$, which is needed to sample $r^{(g)}_{t+k+1}$, can be calculated deterministically once $r^{(g)}_{t+k}$ is given as shown below
$$
h_{t+k+1}^{(g)} = \alpha_0^{(g)} + \alpha_1^{(g)}r_{t+k}^{(g)} + \beta^{(g)}h_{t+k}^{(g)}
$$
Also, I see that you simply fixed $h_0=0$. I do not know whether or not that is a good assumption. However, I do not know off the top of my head a better way of going about it.
Do you by any chance "predict" raw residuals rather than standardized residuals? (Then you would get exactly the same result of the ARCH-LM test when "pre-testing" and when "post-testing" in the case of no conditional mean model $r_t=\epsilon_t$.) The raw residuals will contain ARCH effects and that is why you want to apply a GARCH model and obtain standardized residuals that do not contain ARCH effects.
On the other hand, if you find ARCH effects in the standardized residuals, the GARCH model you are using may be inappropriate. Try different variants of the GARCH model (EGARCH, APARCH and whatever else) and different lag orders.
Also note that the original ARCH-LM test is inappropriate for testing for remaining ARCH effects in the standardized residuals of a GARCH model; Li-Mak test should be used instead. (I do not know whether the Li-Mak test is available in Stata. Also, use of the original ARCH-LM test seems to be relatively widespread in the applied literature, even though it is inappropriate.)
References:
- Li, W. K. and Mak, T. K. (1994) On the squared residual autocorrelations in non-linear time series with conditional heteroscedasticity. Journal of Time Series Analysis 15, 627–36.
Best Answer
It sounds like your question is about the properties of the QMLE estimates of the $GARCH$ coefficients when the $GARCH$ process is heavy tailed. In general, you will still get consistent and asymptotically normal estimates via the QMLE, but they won't be efficient or unbiased.
If $y_t$ are the innovations in the financial returns, and you assume $y_t$ is a standard $GARCH(p,q)$, so $y_t = h_t \eta_t$, where $\eta_t$ are the innovations to the $GARCH$ process (normalized to $E[\eta_t] = 0$, $E[\eta_t^2] = 1$) and $h_t$ is the volatility process, then even if $\eta_t$ has finite moments (for example a $\eta_t \sim N(0,1)$), then $y_t$ can still have an infinite fourth moment (and hence a fat tail) since $y_t^2$ is regularly varying with index $\alpha \in (1, \infty)$. The index $\alpha$ depends on the $GARCH$ coefficients and the distribution of $\eta$. For example, in a $GARCH(1,1)$, $\alpha$ is the solution to $E[a_1\eta^2 + b_1] = 1$.
But the QMLE will still be consistent and asymptotically normal if $\eta$ is not heavy tailed. If $\eta$ is heavy tailed itself, then the QMLE will be converge slower than $\sqrt{n}$, and instead will converge to a stable distribution.