Bayesian GARCH: forecasting volatility and returns

Question

I am interested in forecasting conditional volatility, $h_t$, and returns, $y_t$, in a Bayesian GARCH framework. I am using the bayesGARCH package by Ardia in R (https://cran.r-project.org/web/packages/bayesGARCH/index.html) in order to estimate the Bayesian GARCH. The model setup implemented by the bayesGARCH package can be seen below.

$r_t$ = $\varepsilon_t(\frac{v-2}{v}\omega_th_t)^{1/2}$ $\quad$ with $\quad$ $t=1,...,T$

$\varepsilon_t \overset{iid}{\sim}N(0,1)$

$\omega_t \overset{iid}{\sim}IG(\frac{v}{2},\frac{v}{2})$

$h_t = \alpha_0 + \alpha_1r^{2}_{t-1}+\beta h_{t-1}$

This setup represents a GARCH(1,1) model with student-t innovations. Currently I am able to forecast $h_t$. I have however only managed to do so using one estimation window. Ideally, I would like to create a rolling window forecast which essentially allows me to update the posterior distribution of the parameters say every 50 observations instead of using the same posterior parameter estimates for the entire test sample. Does anyone have any ideas as to how to go about this?

I set out my existing R code below. Any advice/critiques/comments would be greatly appreciated.

1. Fit bayesGARCH on training data

   MCMC.1 = bayesGARCH(train_data, mu.alpha = c(0,0), Sigma.alpha = 10000 * diag(1,2),mu.beta = 0, Sigma.beta = 10000, lambda = 0.01, delta = 4, control = list(n.chain=2,l.chain=nr_iterations,refresh=1000))
   

2. Extract parameter posterior distribution

   smpl <- formSmpl(MCMC.1, l.bi = 10000, batch.size = 2)        
   alpha0 = smpl[, 1]
   alpha1 = smpl[, 2]
   beta   = smpl[, 3]
   nu     = smpl[, 4]
   

3. Create an empty conditional volatility series, with 100 rows (one entry for every percentile) and where the amount of columns equals the length of your test(forecast) sample - i.e. you'll have a sigma distribution for each day in your forecast sample

   ind  = seq(1,nr_iterations - burn-in,by=100)
   nind = length(ind)
   sigma2   = matrix(0,nind, test_length - 1)
   

4. For the conditional volatility forecast merely substitute the estimated parameters into the equation and use the training data as your return series for period t and then the test data for every period thereafter.

   l    = 0 
   for (i in ind){ 
    sigma2_temp = rep(0, test_length - 1)
    sigma2_temp[1] = alpha0[i]+alpha1[i]*train_data[train_length - 1]^2 
     for (t in 2:(test_length - 1))
      sigma2_temp[t] = alpha0[i]+alpha1[i]*test_data[t-1]^2+beta[i]*sigma2_temp[t-1]
      l = l + 1 
    sigma2[l,] = sigma2_temp
    }
   

I am also interested in forecasting $y_t$. While I have seen the math and intuition for forecasting $y_t$ (https://stats.stackexchange.com/a/152427/98939) I am yet to find any code. Any ideas?

Answer 1

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.

Answer 2

Look into pymc3 deterministic and gaussianrandomwalk features.

From: https://pymc-devs.github.io/pymc3/getting_started/

from pymc3 import Exponential, StudentT, exp, Deterministic
from pymc3.distributions.timeseries import GaussianRandomWalk

with Model() as sp500_model:

nu = Exponential('nu', 1./10, testval=5.)

sigma = Exponential('sigma', 1./.02, testval=.1)

s = GaussianRandomWalk('s', sigma**-2, shape=len(returns))

volatility_process = Deterministic('volatility_process', exp(-2*s))

r = StudentT('r', nu, lam=1/volatility_process, observed=returns['S&P500'])