Generating GARCH time series with Python

The GARCH model is commonly employed in modeling financial time series. I show how to implement this model in python and generate data that have the properties of the stock market.


The GARCH process was subject for the Nobel Prize in Economics in 2003, rewarded to Robert F. Engle and Clive Granger, "for methods of analyzing economic time series with time-varying volatility (ARCH)". In fact, GARCH is a generalized version of ARCH.


The GARCH(p, q) model is defined by \begin{align*} y_k &= \sigma_k \varepsilon_k \\ \sigma_k^2 &= \omega + \sum_{i=1}^p \alpha_i y_{k-i}^2 + \sum_{j=1}^q \beta_j \sigma_{k-j}^2 \end{align*} where $\omega > 0$, $\alpha_i \ge 0$, $\beta_j \ge 0$ and $\varepsilon_i$ are IID. See Lecture notes: Financial time series, ARCH and GARCH models for more details.

Remark: Sometimes the definition is written as \[ y_k = x_k \theta + \sigma_k \varepsilon_k \] where $(x_k)$ is some input time series, used as a baseline for the generated data, where $\theta$ is a parameter deciding "how much" to re-use the information from $(x_k)$ time series in the resulting $(y_k)$ time series.

Implementing the GARCH model

There are many theoretical sources on the web discussing the GARCH process. Code examples are less common, which is why I share some code here.

Python code

import numpy as np
import matplotlib.pyplot as plt

def garch(ω, α, β, n_out=1000):
    p = len(α)
    q = len(β)

    # Since the first max(p, q) number of points are not generated from the garch
    # process, the first points are garbage (possibly extending beyond max(p, q)),
    # so we drop n_pre > max(p, q) number of points.
    n_pre = 1000
    n = n_pre + n_out

    # Sample noise
    ɛ = np.random.normal(0, 1, n)

    y = np.zeros(n)
    σ = np.zeros(n)

    # Pre-populate first max(p, q) values, because they are needed in the iteration.
    for k in range(np.max([p, q])):
        σ[k] = np.random.normal(0, 1)
        y[k] = σ[k] * ɛ[k]

    # Run the garch process, notation from
    for k in range(np.max([p, q]), n):
        α_term = sum([α[i] * y[k-i]**2 for i in range(p)])
        β_term = sum([β[i] * σ[k-i]**2 for i in range(q)])
        σ[k] = np.sqrt(ω + α_term + β_term)
        y[k] = σ[k] * ɛ[k]

    return y[n_pre:]

# Make a cumulative series from a "delta series".
def delta_to_cum(ys):
    ys_cum = []
    y_cum = 0
    for y in ys:
        y_cum += y
    return ys_cum

# Define a garch(1,1) process
ω = 0.1
α = [0.3]
β = [0.2]
y = garch(ω, α, β)
x = range(len(y))

# The resulting series looks like stock returns and
# the cumulative series looks like stock prices.
plt.plot(x, y)
plt.plot(x, delta_to_cum(y))

Fitting data to a GARCH model

The python package arch can be used to fit a GARCH model to existing data.