I have fitted a Bayesian GARCH(1,1) model with Student $t$ innovations to some time series data, $X_1,...,X_n$ and now want to estimate Value-at-Risk (VaR) (i.e., 5% quantiles) at each times $t=1,,...,n$. To do so, I'm using bayesGARCH in R, which returns an mcmc object from which I can obtain samples from the posterior distribution of the coefficients $\theta$, which in turn determine $\textrm{Var}(X_t)$.
More formally, I want to estimate the quantiles of a random variable $X_t$ which is constructed as follows: $$ X_t = \sigma_t(\theta) Z_t $$ $$ \theta \sim F $$ $$ Z_t \sim \sqrt{\frac{\nu-2}{\nu}} ~ t_\nu, $$ where $F$ is unknown but I have $M$ samples from it. ($t_\nu$ is scaled to that $\textrm{Var}(Z_t) = 1$ and hence $\textrm{Var}(X_t) = \sigma_t^2$.) Given $\theta$, I can readily compute $\sigma_t$.
If $\sigma_t$ were fixed, then clearly $X_t \sim \sigma_t \sqrt{\frac{\nu-2}{\nu}} t_\nu$. The stochasticity of $\sigma_t$ complicates estimating the quantiles of $X_t$.
I was thinking of 2 approaches:
For $j=1,...,M$, simulate $X_j \sim \sigma_t(\theta_j) \sqrt{\frac{\nu-2}{\nu}} t_\nu$ (perhaps multiple times). Then, take $\widehat{\textrm{VaR}}$ as the empirical 5% quantile of $X_1,...X_M$.
For $j=1,...,M$, let $\widehat{\textrm{VaR}}_j$ be the 5% quantile of $\sigma_t(\theta_j) \sqrt{\frac{\nu-2}{\nu}} t_\nu$. Then, take $\widehat{\textrm{VaR}} = \frac{1}{M} \sum_j \widehat{\textrm{VaR}}_j$.
My question is: which of the two approaches are valid? Does there exist a more appropriate method?
To me, approach (1) seems reasonable, yet perhaps inefficient. Approach (2), might be invalid, but I couldn't say why. Indeed, these two approaches yield quite different results, especially for more extreme quantiles).
