A note on simulation in the rugarch package.

In the rugarch package there are two main ways to simulate a GARCH process:

  1. The ugarchsim method which takes an already estimated object of class uGARCHfit.
  2. The ugarchpath method which takes a specification of class uGARCHspec with fixed parameters.

Before proceeding to a demonstration to show how to obtain equivalence between the 2, it is important to say a few words about the initialization of the simulation which has led to some confusion. Consider the typical GARCH(1,1) process:
\[ \begin{gathered}
{y_{t + 1}} = {\mu _{t + 1}} + {\varepsilon _{t + 1}},{\text{ }}{\varepsilon _{t + 1}} = {z_{t + 1}}\sqrt {{h_{t + 1}}}\\
{h_{t + 1}} = \omega + {\alpha _1}\varepsilon _t^2 + {\beta _1}{h_t}\\
\end{gathered} \]
Now assume that we wish to simulate \( t+1 \) given estimated information upto time \( t \). As can be seen from the equation, we require a value for \( \varepsilon_{t}^2 \). In the rugarch package, the ugarchsim method has an option (‘startMethod’) of using either the estimated object’s last known values at time \( t \) to initialize the simulation (“sample”) or their unconditional long run values (“unconditional”). In the former case, \( \varepsilon_{t}^2 \) is based on the last value from the estimated object, whilst in the latter case \( z_{t} \) is set to zero which means that \( \varepsilon_{t}^2 \) is also zero. In either case, it is clear that there is no distributional ‘uncertainty’ in the value of \( h_{t+1} \). At \( t+2 \) however, we need a value for \( \varepsilon_{t+1}^2 \) which we do not have from the estimated object which stops at time \( t \) so that we need to use simulated \( z_{t+1} \) from the standardized distribution used in conjunction with the value of \( h_{t+1} \) calculated in the previous step, and hence the distributional uncertainty now enters the simulation.

The following code shows how to obtain equivalence between the 2 simulation methods.

data(sp500ret)
spec1 = ugarchspec()
fit = ugarchfit(spec1, sp500ret)
sim1 = ugarchsim(fit, n.sim = 1, m.sim = 5, startMethod = 'sample', rseed = 10)
spec2 = spec1
setfixed(spec2) < - as.list(coef(fit))
sim2 = ugarchpath(spec2, presigma = as.numeric(tail(sigma(fit), 1)), preresiduals = as.numeric(tail(residuals(fit),
    1)), prereturns = tail(sp500ret[, 1], 1), n.sim = 1, m.sim = 5, rseed = 10)
sigma(sim1)
##       [,1]   [,2]   [,3]   [,4]   [,5]
## T+1 0.0248 0.0248 0.0248 0.0248 0.0248
# simulated sigma
sigma(sim1) == sigma(sim2)
##     [,1] [,2] [,3] [,4] [,5]
## T+1 TRUE TRUE TRUE TRUE TRUE
# simulated returns
fitted(sim1) == fitted(sim2)
##     [,1] [,2] [,3] [,4] [,5]
## T+1 TRUE TRUE TRUE TRUE TRUE

This short note described the simulation methods used in the rugarch package and its rational based on distributional uncertainty which does not affect \( h_{t+1} \). There is also a method called ugarchboot which can take into account parameter uncertainty (via simulation of the parameter distribution) which then leads to uncertainty in the \( h_{t+1} \) value.