I am using this package in R to do Bayesian estimation of GARCH models. I want to forecast $y_t$ (i.e. the mean equation), but it seems that the package has no built-in function for this. The model setup in the package is as follows: enter image description here

Now I was thinking how I should implement a forecasting procedure. I thought of doing this:

  1. Estimate model to obtain parameter posteriors
  2. Take mean of each posterior to get a point estimate
  3. Obtain $h_t$
  4. Simulate from normal and gamma to obtain $y_t$.

My questions/problems:

  1. Is this a proper way of doing it?
  2. This might take very long, because every forecast means that I have to simulate three times, first to get the posteriors, than for $\epsilon_t$ and then for $\bar{\omega}_t$
  3. Are there any better options (preferably already coded)?
  • $\begingroup$ I've been busy with a similar bayesGARCH analysis. I've posted the code I've come up with so far as a response to a similar question - have a look and see if it helps. Also any comments/improvements to my code would be appreciated. See here: stats.stackexchange.com/a/191756/98939 $\endgroup$ Jan 21, 2016 at 12:10

1 Answer 1


I am sure you already know $E[y_t]=0$ as it follows directly from the first equation you have listed. So when you say "I want to forecast $y_t$", I interpret it as "I would like to find the forecast distribution of $y_t$". Note that because $\varepsilon \stackrel{iid}{\sim} \mathcal{N}(0,1)$ we know that $$y_t \stackrel{iid}{\sim} \mathcal{N}(0,\frac{\nu-2}{\nu}\varpi_t h_t)$$

where $\frac{\nu-2}{\nu}\varpi_t h_t$ is the variance of $y_t$. For notation purpose define $$f(y_t\,|\,\nu,\varpi_t,h_t):=\mathcal{N}(0,\frac{\nu-2}{\nu}\varpi_t h_t)$$

Since the variance of $y_t$ is a random variable, the marginal distribution of $y_t$, $f(y_t)$, is not directly apparent. In this case, we can use a Gibbs sampler approach to draw from $f(y_t\,|\,\nu,\varpi_t,h_t)$ many times, each for different values of $h_t$ and $\varpi_t$ of whom were previously drawn from their own distributions. The resulting sample will in effect be drawn from $f(y_t)$ .

I am not familiar with this R function. Posting the output of the function would be helpful. I am assuming it returns a posterior sample of the $\alpha_0,\alpha_1,\beta$ and $\varpi$. Suppose the sample has $G$ draws. Then

  1. For each posterior draw $g=1,2,\cdots,G$, compute $h^{(g)}_t$
  2. simulate $y^{(g)}_t \sim \mathcal{N}(0,\frac{\nu-2}{\nu}\varpi^{(g)} h^{(g)}_t)$. Equivalently, you could simulate $G$ $\varepsilon_t$'s and multiply each one by a respective $(\frac{\nu-2}{\nu}\varpi^{(g)}_t h^{(g)}_t)^{1/2}$

The second step can be done multiple times each draw.

I would not take the sample mean of the posterior and use it to generate an estimate for $h_t$. This is because there is a non-linearity in the recursion between $\beta$ and the other variables. So $h_t$ evaluated at the mean of the posterior for the other variables may not be equal to $E[h_t]$ which could inturn lead to bias.

As far as something already coded goes...I do not know of anything off hand. If you are unsatisfied or would like help with coding, then editing your post to include the code you already have and the resulting output would help users answer the question more effectively.

  • $\begingroup$ Thank you very much! This is extremely helpful to me! I will try to program it myself first, but may I get back to you through commenting on this question if I don't succeed? Thanks again! $\endgroup$
    – rbm
    May 15, 2015 at 6:35
  • $\begingroup$ Actually one small question came up: "the second step can be done multiple times each draw". Does that mean that for instance I would do that 10000 and there I would take the mean to get $h_t^{(g)}$? $\endgroup$
    – rbm
    May 15, 2015 at 6:44
  • $\begingroup$ No. $h^{(g)}_t$ comes from the $g^{th}$ posterior draw. In fact you can probably forget that sentence altogether. What I meant was for each $g$ you could potentially simulate more than one $y^{(g)}_t$ though it is not necessary and people often do not do that. For example, if $G=10,000$ you can simulate 100,000 $y_t$ by simulating 10 $y^{(g)}_t$ for each $g$. All it does is increase the sample size for your forecast distribution. In general, I don't think it makes much difference when $G$ is sufficiently large. $\endgroup$ May 15, 2015 at 7:15
  • $\begingroup$ But won't my draws of $y_t$ become very much dependent on the outcome of $\epsilon_t$ if I don't simulate for instance 10000 times and take the mean of that? $\endgroup$
    – rbm
    May 15, 2015 at 7:18
  • $\begingroup$ Because what I mean is that if I have just one value for $\epsilon_t$, and this turns out to be far off from the mean of $\epsilon$, won't I get values for $y_t$ that are off even if I do this for all G? $\endgroup$
    – rbm
    May 15, 2015 at 11:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.