[R-SIG-Finance] How are errors terms calculated in GARCH model by rugarch package?

[Xie Yijun]

Hi,

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 at fit$z
###################

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}^2$

$z_t = X_t / \sigma_t$

The code is:

###########################
     omega1 = myfit at fit$coef[1]
     alpha1 = myfit at fit$coef[2]
     beta1 = myfit at fit$coef[3]

     z2 = rep(0,(N-1000))
     sighat1 = rep(0,(N-1000))
     sighat1[2] = 1
     sighat1[1] = 1
     for(i in 2:(N - 1000)){
       sighat1[i] = sqrt(omega1 + alpha1 * X1[i-1] + 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!
Patrick