How are errors terms calculated in GARCH model by rugarch package?

I am fitting a GARCH(1,1) model to the data and want to look at the innovation distribution.

#generate the data
set.seed(1)
N = 5000
omega = 0.5
alpha = 0.08
beta = 0.91
X1 = rep(0,N)
X2 = rep(0,N)
sig1 = rep(0,N)
sig2 = rep(0,N)

for(i in 2:N){
sig1[i] = sqrt(omega + alpha * X1[i-1] + beta * sig1[i-1]^2)
X1[i] = sig1[i] * rnorm(1)
sig2[i] = sqrt(omega + alpha * X2[i-1] + beta * sig2[i-1]^2)
X2[i] = sig2[i] * rt(1, df = 8)
}

X1 = X1[-c(1:1000)]

I first generate the data and fit it to GARCH(1,1) model with t innovation.

spec = ugarchspec(mean.model=list(armaOrder=c(0,0),include.mean=F),
distribution.model="std") # GARCH(1,1) model
myfit = ugarchfit(spec, X1)

Suppose the fitted model is called myfit then I can get the error terms by z1 = myfit@fit$z which is calculated as$(z_1,...z_n) = (\frac{X_1}{\hat{\sigma}_1},...,\frac{X_n}{\hat{\sigma}_n})$However, if I only extract the parameters$\omega,\alpha$, and$\beta$estimated by the rugarch package and calculate the error terms manually as$\sigma_t^2 = \omega + \alpha X_{t-1}^2 + \beta \sigma_{t-1}^2z_t = X_t / \sigma_t$The code is: omega1 = myfit@fit$coef[1]
alpha1 = myfit@fit$coef[2] beta1 = myfit@fit$coef[3]

z2 = rep(0,(N-1000))
sighat1 = rep(0,(N-1000))
sighat1[1] = 1
z2[1] = X1[1]/sighat1[1]
for(i in 2:(N - 1000)){
sighat1[i] = sqrt(omega1 + alpha1 * X1[i-1]^2 + beta1 * sighat1[i-1]^2)
z2[i] = X1[i]/sighat1[i]
}

I got very different results between these two approaches by comparing the Q-Q plot of z1 and z2.

qqnorm(z1)
qqline(z1)
qqnorm(z2)
qqline(z2)

z1 seems to be normally distributed following the data generating process, while z2 has a heavy tail following the model specification. I was wondering why they are so different? I arbitrarily chose 1 for for first two terms of $\hat{\sigma}$, so does the difference come from the initial values?

And more generally, how does rugarch package fit the GARCH model and choose initial values? My understanding is that we need to find parameters using either QMLE or MLE and then find error terms iteratively using my second approach. But I am not sure how is the initial value chosen.

Thanks!

• I have no idea how to answer this, but I do know where you are likely to find an answer. The R-SIG-Finance mailing list is populated by a good percentage of the authors of important R finance packages. Your question will be seen, and undoubtedly answered by an expert on this package.
– meh
Jun 8, 2016 at 12:12
• You could also examine - cran.r-project.org/web/packages/rugarch/vignettes/…
– meh
Jun 8, 2016 at 12:14
• @aginensky Thank you so much! I would try to seek help from the mailing list. Jun 8, 2016 at 22:05

I also posted this question on the R-SIG-Finance mailing list (thanks to @aginensky's suggestion), and the author of rugarch replied it. The initial value of $\hat{\sigma}^2_1$ by default is mean of squared residuals, and in this case it's sighat1[1] = sqrt(mean(X1^2)).