GARCH volatility models

Volatility clusters

Daily returns are nearly unpredictable, but their size is not: a big move today makes another big move tomorrow more likely, and quiet days come in runs. This volatility clustering is one of the best-established facts in finance, and it breaks the assumption of a constant volatility behind simple VaR. GARCH (generalised autoregressive conditional heteroskedasticity) models the conditional variance directly, so risk and option prices can follow the market's mood.

import numpy as np

rng = np.random.default_rng(4)
n = 3000
omega, alpha, beta = 2e-6, 0.10, 0.85
var = np.empty(n)
returns = np.empty(n)
var[0] = omega / (1 - alpha - beta)
returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
    var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
    returns[t] = np.sqrt(var[t]) * rng.normal()

# Clustering: the size of a return today predicts the size tomorrow, though the sign is a coin flip
print(round(float(np.corrcoef(returns[:-1], returns[1:])[0, 1]), 3))                 # about 0
print(round(float(np.corrcoef(np.abs(returns[:-1]), np.abs(returns[1:]))[0, 1]), 3))  # clearly positive

The GARCH(1,1) recursion

The model has three parameters: var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]. omega is a floor, alpha the reaction to yesterday's shock, and beta the persistence of yesterday's variance. Because each day depends on the one before, the series is built with a loop, starting from an initial variance. Typical daily equity values are alpha around 0.05 to 0.1 and beta near 0.85 to 0.95: a shock fades slowly.

import numpy as np

def garch_variance(returns, omega, alpha, beta, initial_var):
    var = np.empty(len(returns))
    var[0] = initial_var
    for t in range(1, len(returns)):
        var[t] = omega + alpha * returns[t - 1] ** 2 + beta * var[t - 1]
    return var

returns = np.array([0.01, -0.02, 0.015, 0.0, -0.03, 0.01, 0.0, 0.0])
var = garch_variance(returns, 1e-6, 0.08, 0.9, 1e-4)
print(np.round(np.sqrt(var) * 100, 3))             # daily volatility in %: it jumps after the -3% day, then fades

Where it settles, and the forecast

If alpha + beta < 1 the variance is mean-reverting and settles to the long-run level omega / (1 - alpha - beta); if it is 1 or more it never does. For forecasts, tomorrow's variance is known from today's return and variance. Beyond that, the expected squared return equals the variance, so each step is omega + (alpha + beta) * previous: the forecast glides from today's level towards the long-run one, at a speed set by the persistence alpha + beta.

def long_run_variance(omega, alpha, beta):
    if alpha + beta >= 1:
        raise ValueError("alpha + beta must be below 1")
    return omega / (1 - alpha - beta)

def garch_forecast(omega, alpha, beta, last_var, last_return, steps):
    out = [omega + alpha * last_return ** 2 + beta * last_var]
    for _ in range(steps - 1):
        out.append(omega + (alpha + beta) * out[-1])
    return out

omega, alpha, beta = 1e-6, 0.08, 0.9
print(long_run_variance(omega, alpha, beta) ** 0.5)                        # about 0.00707, 0.71% a day
forecast = garch_forecast(omega, alpha, beta, 4e-4, -0.03, 60)             # after a -3% day in a volatile spell
print([round(v ** 0.5 * 100, 3) for v in (forecast[0], forecast[9], forecast[29], forecast[59])])

Choosing the parameters: the likelihood

The parameters are not known; they are estimated from the data by maximum likelihood. If returns are normal with variance var[t], the log-likelihood of a day is -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]). Summing it gives one number that scores a choice of (omega, alpha, beta): the best parameters make the observed returns most probable. An optimiser (for example scipy.optimize.minimize) tries many triples and keeps the highest.

# Score three parameter choices on data simulated with alpha = 0.10, beta = 0.85
import numpy as np

rng = np.random.default_rng(4)
n = 3000
true = (2e-6, 0.10, 0.85)
var = np.empty(n); returns = np.empty(n)
var[0] = true[0] / (1 - true[1] - true[2]); returns[0] = np.sqrt(var[0]) * rng.normal()
for t in range(1, n):
    var[t] = true[0] + true[1] * returns[t - 1] ** 2 + true[2] * var[t - 1]
    returns[t] = np.sqrt(var[t]) * rng.normal()

def loglik(params):
    omega, alpha, beta = params
    v = np.empty(n)
    v[0] = returns.var()
    for t in range(1, n):
        v[t] = omega + alpha * returns[t - 1] ** 2 + beta * v[t - 1]
    return float(np.sum(-0.5 * (np.log(2 * np.pi * v) + returns ** 2 / v)))

for params in ((2e-6, 0.10, 0.85), (2e-6, 0.30, 0.50), (1e-5, 0.0, 0.0)):
    print(params, round(loglik(params), 1))                       # the true parameters score highest

• Volatility clusters: GARCH(1,1) builds the variance by a loop, var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1].
• With alpha + beta < 1 it reverts to omega / (1 - alpha - beta); forecasts use omega + (alpha + beta) * previous.
• Parameters are chosen by maximum likelihood, sum(-0.5 * (log(2 * pi * var) + r**2 / var)); use fat-tailed shocks for tails.

Next up: Brownian motion, the random walk in continuous time that underlies option pricing models.

Lab · what's graded

• Run the recursion with garch_variance25 pts
• Find where it settles with long_run_variance25 pts
• Forecast the next days with garch_forecast25 pts
• Score the fit with garch_loglik25 pts