How to calculate the kurtosis implied by a GARCH model? I have estimatate a GARCH model using the "rugarch" package in R. The data are the returns on the S&P500 named xst and stored as timeSeries object.
I would like to compare the kurtosis implied by the model specification vs the kurtosis observed in the real data used to specify the model. I searched in the documentation of the "rugarch" package but could not find how to do it.
Questions:


*

*If I extract the fitted values of the GARCH model and use the kurtosis function (from package "moments") to calculate the kurtosis implied by the model, is this procedure the appropriate one to find the implied kurtosis?  

*If I were to compare it with the observed kurtosis, what should I expect? 

*How would I tell if it was the correct solution or not?


Here are the R functions:
model <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder =                
c(1,1)), mean.model = list(armaOrder = c(1,1), include.mean = T), 
                      distribution.model = "norm")

modelfit <- ugarchfit(model, xst)

s <-modelfit@fit$fitted.values

kurtosis(s)

 A: Ruey S. Tsay derives the excess kurtosis of a GARCH(1,1) model in his "Analysis of Financial Time Series" (2010) textbook, section 3.16, p. 165. (He also notes that the same idea applies to other GARCH models.)
The model under consideration is
$$
\begin{align}
 a_t &= \sigma_t \epsilon_t, \\ 
\sigma_t^2 &= \alpha_0 + \alpha_1 a_{t-1}^2 + \beta_1 \sigma_{t-1}^2, 
\end{align}
$$
where $\alpha_0 > 0$, $\alpha_1 \geqslant 0$, $\beta_1 \geqslant 0$, $\alpha_1 + \beta_1 < 1$, $\epsilon_t$ is $i.i.d.$ with $\mathbb{E}(\epsilon_t) = 0$, $\text{Var}(\epsilon_t)=1$, $\mathbb{E}(\epsilon_t^4)=\text{K}_{\epsilon}+3$ where $\text{K}_{\epsilon}$ is the excess kurtosis of the standardized errors.
If the excess kurtosis of $a_t$, $\text{K}_{a}$, exists and if $1 - \alpha_1^2 (\text{K}_{\epsilon} + 2) - (\alpha_1 + \beta_1)^2 > 0$, then
$$ \text{K}_{a} = \frac{ (\text{K}_{\epsilon} + 3)(1 - (\alpha_1 + \beta_1)^2 ) }{ 1 - 2 \alpha_1^2 - (\alpha_1 + \beta_1)^2 - \text{K}_{\epsilon} \alpha_1^2 } - 3. $$
[Non-excess] kurtosis of $a_t$ will be $\text{K}_{a} + 3$.
In case the standardized errors are Gaussian, $\text{K}_{\epsilon} = 3$ and the excess kurtosis of $a_t$ becomes
$$ \text{K}_{a}^{Gaussian} = \frac{ 6 \alpha_1^2 }{ 1 - 2 \alpha_1^2 - (\alpha_1 + \beta_1)^2 }. $$
The textbook section contains a few more related results.
