In a GARCH(1,1) model

$$ x_t = \sigma_tz_t$$ $$\sigma_{t+1}^2=a_0 + a_1x_t^2 + b_1\sigma_t^2$$

the kurtosis (when it exists) can be shown to be equal to

$$ \kappa_x = \kappa_z \frac{1-(a_1+b_1)^2}{1 - (a_1+b_1)^2 - a_1^2 (\kappa_z - 1) }$$

where $\kappa_z$ is the kurtosis of $z_t$. For standard normal innovations, i.e. when $z_t \sim N(0,1)$, $\kappa_z = 3$, and with $$ a_0 = 0.01, a_1 = 0.09, b_1=0.9$$

this gives

$$ \kappa_x = 3 \frac{1-(0.99)^2}{1 - (0.99)^2 - 0.09^2 (3 - 1) } \approx 16.14$$

However, when I run simulation of this GARCH process, I find that the sample kurtosis is somewhere around 7-8. For example, the plot below shows sample kurtosis calculated on 1,000 simulations of the above GARCH process with 10,000 time steps each. I don't understand where this mismatch between the theoretical value above and estimates is coming from.

Sample kurtosis calculated on 1000 GARCH(1,1) simulations of 10,000 time-steps each.

Add 1

I have run simulations with different parameter settings to move the process away from IGARCH, but it doesn't seem to have improved much the results. For example, with $$ a_0 = 0.2, a_1 = 0.383, b_1 = 0.417$$, the theoretical kurtosis is about $16.2$ but in the simulations the sample estimates produce a mean kurtosis of about $9.7$ and median of about $6.7$ (plot below).Sample kurtosis for 1,000 GARCH simulations with 10,000 time steps each.

Add 2

I also run simulations using simple t-distributed random number with 4.45 degrees of freedom (which also gives a kurtosis of around 16.3) and calculated sample kurtosis on these. The plot below shows the results, which are also "bunched up" around 8 with the median of 8.78, but the mean is driven up to 20 by the few of the very large outliers (actually, it is mainly the most extreme one which pushes it up to 20, without it the average kurtosis is 12.3). So it is qualitatively similar to GARCH results and supports the arguments by Matthew Gunn.

Sample kurtosis from 1,000 samples of t distributed RVs of size 10,000 each

  • 1
    $\begingroup$ I believe this is caused by your parameters: your model specification is almost an IGARCH model ($\alpha + \beta = 1$). The IGARCH model has an infinite variance. $\endgroup$
    – Malick
    Commented Oct 1, 2018 at 17:28
  • $\begingroup$ @Malick Thank you for your comment. The sample variance seems to remain fairly close to the theoretical (which is 1 in this case). I thought it might be convergence issue, but even with $10^8$ time steps it is still far away from the theoretical value. $\endgroup$
    – Confounded
    Commented Oct 1, 2018 at 17:34

1 Answer 1


You've found parameterizations where fantastically long samples are required for sample 4th moments to converge on population 4th moments.

Quick evidence of imprecise estimation

Let $k_i$ denote your estimated kurtosis in simulation $i$. Looking across your $i = (1,\ldots, 1000)$ simulations, your $k_i$ estimates are all over the place. What's your standard error for $\frac{1}{n} \sum_i k_i$? It's huge.

You're saying your sample kurtosis is around 7-8, but if your standard error is huge, how can you say your result is inconsistent with an actual kurtosis of 16? You can't.

Quick Theory: $\{x^2_t\}$ is close to non-stationary

A GARCH(1,1) implies an ARMA(1,1) in the squared process. If your GARCH model is: $$ x_t = \sigma_tz_t$$ \begin{align*} \sigma_{t}^2&=\omega + a_1x_{t-1}^2 + b_1\sigma_{t-1}^2 \end{align*}

It implies an ARMA(1,1) representation for $\{x^2_t\}$. Observe that you can write $x_t^2 = \operatorname{E}_{t-1}[\sigma^2_tz_t^2] + u_t = \sigma_t^2 + u_t$. Using the lag operator $L$ we can write $(1 - b_1L) \sigma_t^2 = \omega + a_1 L x_{t}^2$. Combining those equations you get:

$$ x_t^2 = \omega + (a_1 + b_1) x_{t-1}^2 + u_t - b_1 u_{t-1} $$

You can see here that if $a_1 + b_1 = 1$, the model is non-stationary. For you, $a_1+b_1 = .99$ and $\{x_t^2\}$ is incredibly persistent.

Another requirement for the existence of a 4th moment of a GARCH(1,1) is $b_1^2 + 2a_1b_1 + 3a_1^2 < 1$. (See Theorem (2.3) and Example 2.4 from Petra Posedel.) Part of what's happening is your 2nd example is getting close to that constraint. Speaking loosely, you've found parameterizations where there's a small probability of very large volatility shocks.

Simple example with slow convergence of sample mean on population mean

Let $X$ denote a random variability with a $p=.000001$ (one in a million) probability of 10,000,000 and a $q = 1 - p$ probability of 0.

  • Trivially, $\operatorname{E}[X] = 10$.

  • In an IID sample of 1000 observations, there's a 99.9 percent probability the sample mean is 0 and approximately .0999 percent probability of a sample mean of 10,000.

The distribution of the sample mean is not symmetric around the expected value: far more observations will be below than above. Also observe that you need extreme sample sizes for the sample mean to converge on the population mean.

More broadly, small measure events that you're unlikely to sample via Monte-Carlo methods can have significant impact on overall population moments if the random variables takes on huge values in those events.

  • $\begingroup$ Thank you for your answer. The standard deviation might be large indeed, but the estimates are not even centered on the expected value; with would suggest a strong bias. Also, I made further simulations (see the addendum 1 to my post) where I moved $a_1+b_1 = 0.8$, but the results are still biased and with large standard deviation (about 18). $\endgroup$
    – Confounded
    Commented Oct 1, 2018 at 18:21
  • 1
    $\begingroup$ @Confounded Why should they be centered on the expected value? Is the sample distribution required to be symmetric? No. If your estimating the expected value of the lottery, most every small subsample of tickets will fall below the actual expectation. $\endgroup$ Commented Oct 1, 2018 at 18:26
  • 1
    $\begingroup$ Also, you state that "This makes the variance of $x_t^2$ quite large." But $Var[x_t^2] = E[x_t^4] - (E[x_t^2])^2 = E[x_t^4] - 1$ since for the parameters used $E[x_t^2] = 1$, so this would give large kurtosis, but it doesn't. $\endgroup$
    – Confounded
    Commented Oct 1, 2018 at 18:29
  • 1
    $\begingroup$ On the parameters used: I believe that the values of the parameters are not too far from what is estimated for various financial time series. If it is indeed the case that convergence for these parametrisation is so slow, then what can we then possibly infer from the estimates based on the real observations where available time series are typically rather short? (It's really more of a rhetorical question). $\endgroup$
    – Confounded
    Commented Oct 1, 2018 at 18:49
  • 1
    $\begingroup$ Thank you for your answer. It seems to me then that kurtosis is not an appropriate way to assess "flat-tailness" of the GARCH process since in the sample sizes it is used for there is effectively 0 probability of observing these extreme vols that drive up the kurtosis. What other measures/statistics one can calculate in order to measure "flat-tailness" of GARCH and set the parameters appropriately? Thanks $\endgroup$
    – Confounded
    Commented Oct 2, 2018 at 8:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.