Forecasting Bayesian GARCH(1,1) volatilities As a beginner in Bayesian statistics, I was wondering how one can make a GARCH(1,1) volatility point forecast using a Bayesian approach in the following model:
$$\sigma^2_{t+1}=\alpha_0+\alpha_1\varepsilon^2_t+\beta\sigma^2_t$$
The $\varepsilon_t$ corresponds to the zero mean financial returns. The $\sigma^2_t$ corresponds to the conditional variance of these demeaned returns (i.e. $Var(\varepsilon_t|I_t)$ where $I_t$ denotes the information set available at time t).
Suppose you have a window of 100 observations and you simulate 1000 values from the posterior distribution of the parameters. So, you have 1000 simulated values for $\alpha_0$, $\alpha_1$ and $\beta$. 
How can one make a forecast using a Bayesian approach? Is one, for example, allowed to do the following: Calculate a 1000 values for $\sigma^2_{t+1}$ using the simulated values for $\alpha_0$, $\alpha_1$ and $\beta$ and take the average over the 1000 values for $\sigma^2_{t+1}$ and use this average over the 1000 values for $\sigma^2_{t+1}$ as your point forecast of the volatility?
And if this is allowed, how is making a Bayesian point forecast in this specific case related to the posterior predictive density (which I thought was the way to go in Bayes). 
So, if one assumes a Normal distribution for the residuals ($\varepsilon_t$) with mean zero and variance equal to $\sigma^2_{t}$, what is then the posterior predictive density for a GARCH(1,1) model (given that we use a truncated normal prior for our parameters $\alpha_0$ and $\alpha_1$ with hyperparameters $\mu_{\alpha}$ which is a $2\times1$ vector of zeros and $\Sigma_{\alpha}$. For the parameter $\beta$ we also use a truncated normal prior. Also, we assume prior independence between parameters $\alpha$ and $\beta$) ?
Many Thanks
 A: A Bayesian forecast that takes into account the parameter uncertainty is constructed as follows when we use the posterior predictive distribution shown below: 
$p(\sigma^2_{t+1}|\sigma^2)= \int_{\theta}p(\sigma^2_{t+1}|\theta,\sigma^2)p(\theta|\sigma^2)d\theta$
where $\theta = [\alpha_0, \alpha_1, 
\beta]$ and $p(\theta|\sigma^2)$ is the posterior distribution of the parameters. The $p(\sigma^2_{t+1}|\theta,\sigma^2)$ is the distribution of the predicted volatilities. 
According to Rachev et al. (2008) we should plug-in the estimated parameters in the GARCH equation to get (f.e. 1000) simulated values of $\sigma^2_{t+1}$. The thousand $\sigma^2_{t+1}$ are an estimation of this distribution of predicted volatilities: $p(\sigma^2_{t+1}|\theta,\sigma^2)$.
To approximate the integral above we use the following:
$p(\sigma^2_{t+1}|\sigma^2)= \int_{\theta}p(\sigma^2_{t+1}|\theta,\sigma^2)p(\theta|\sigma^2)d\theta = E_{\theta|\sigma^2}(p(\sigma^2_{t+1}|\theta,\sigma^2))$
The mean is approximated by taking the sample average of all the predicted volatilities. 
Basically, we approximate this integral (which is an epectation) by using simulated sample averages.  
Reference:
Rachev, S. T., Hsu, J. S., Bagasheva, B. S., & Fabozzi, F. J. (2008). Bayesian methods in finance (Vol. 153). John Wiley & Sons
