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.
The behavior that you see is due to the presample variance option in EViews.
rugarch uses the variance of all data points and EViews uses backcasting using a parameter of 0.7.
More precisely, EViews uses this formula for initialization of the variance:
$\sigma_0^2 = \lambda \hat\sigma^2 + (1-\lambda) \sum_{t=0}^T \lambda^{T-t-1} \cdot \varepsilon_{T-t}^2$
where $\lambda$ is the backcast parameter (default in EViews: 0.7, default in rugarch, fGarch, and gretl: 1.0) and $\hat\sigma^2$ is the unconditional variance of all residuals $\varepsilon_1, \ldots, \varepsilon_T$.
This is the explanation in the EViews manual regarding this choice of the variance initialization (whatever outperform means for them):
Our experience has been that GARCH models initialized using backcast
exponential smoothing often outperform models initialized
using the unconditional variance.
Using the same approach as rugarch, we get:
Dependent Variable: RETURN
Method: ML - ARCH (Marquardt) - Normal distribution
Date: 10/30/17 Time: 20:26
Sample: 1 438
Included observations: 438
Convergence achieved after 11 iterations
Presample variance: unconditional
GARCH = C(2) + C(3)*RESID(-1)^2 + C(4)*GARCH(-1)
Variable Coefficient Std. Error z-Statistic Prob.
C 0.001093 0.000547 1.996164 0.0459
Variance Equation
C 1.33E-06 9.92E-07 1.345676 0.1784
RESID(-1)^2 0.026633 0.007118 3.741463 0.0002
GARCH(-1) 0.963334 0.011717 82.21705 0.0000
which is quite close to the estimates of rugarch.
However, we still see some differences in the standard errors: they are lower for rugarch resulting in differences regarding hypothesis testing of the GARCH parameters.
For comparison, these are the estimates using the fGarch library for R:
Estimate Std. Error t value Pr(>|t|)
mu 1.088e-03 5.198e-04 2.094 0.03627 *
omega 1.333e-06 1.060e-06 1.257 0.20864
alpha1 2.663e-02 9.728e-03 2.738 0.00618 **
beta1 9.633e-01 1.283e-02 75.078 < 2e-16 ***
and gretl:
coefficient std. error z p-value
-------------------------------------------------------
const 0.00108849 0.000519879 2.094 0.0363 **
alpha(0) 1.33330e-06 1.12351e-06 1.187 0.2353
alpha(1) 0.0266348 0.0101192 2.632 0.0085 ***
beta(1) 0.963349 0.0138129 69.74 0.0000 ***
which seem to be more similar to EViews than to rugarch.
What astonishes me most is the drop in the standard errors of the rugarch estimates: e.g. from 57 to 28, or from 1.7 to 0.9 - they appear as if they were halved.
Maybe this finding is worth posting on the R-SIG-Finance mailing list as we don't observe such a drop in the estimated standard errors for the other packages / programs.
fGarch without mean:
Estimate Std. Error t value Pr(>|t|)
omega 1.488e-06 1.219e-06 1.221 0.22227
alpha1 2.517e-02 9.602e-03 2.621 0.00875 **
beta1 9.635e-01 1.424e-02 67.681 < 2e-16 ***
gretl without mean:
coefficient std. error z p-value
-------------------------------------------------------
alpha(0) 1.48801e-06 1.31591e-06 1.131 0.2581
alpha(1) 0.0251710 0.0101023 2.492 0.0127 **
beta(1) 0.963509 0.0156173 61.69 0.0000 ***
EViews without mean:
Variable Coefficient Std. Error z-Statistic Prob.
Variance Equation
C 1.49E-06 1.13E-06 1.314825 0.1886
RESID(-1)^2 0.025179 0.006580 3.826449 0.0001
GARCH(-1) 0.963505 0.012718 75.75811 0.0000
Best Answer
Here is an example of implementation using the
rugarchpackage and with to some fake data. The functionugarchfitallows for the inclusion of external regressors in the mean equation (note the use ofexternal.regressorsinfit.specin the code below).To fix notations, the model is \begin{align*} y_t &= \lambda_0 + \lambda_1 x_{t,1} + \lambda_2 x_{t,2} + \epsilon_t, \\ \epsilon_t &= \sigma_t Z_t , \\ \sigma_t^2 &= \omega + \alpha \epsilon_{t-1}^2 + \beta \sigma_{t-1}^2 , \end{align*} where $x_{t,1}$ and $x_{t,2}$ denote the covariate at time $t$, and with the "usual" assumptions/requirements on parameters and the innovation process $Z_t$.
The parameter values used in the example are as follows.
The image below shows the series of covariate $x_{t,1}$ and $x_{t,2}$ as well as the series $y_t$. The
Rcode used to generate them is provided below.A fit is done with
ugarchfitas follows.Parameter estimates are
And corresponding true values are
The following figure shows a parameter estimates with 95% confidence intervals, and the true values. The
Rcode used to generate it is provided is below.