Modèles de volatilité GARCH

La volatilité se regroupe

Les rendements quotidiens sont presque imprévisibles, mais pas leur taille : un grand mouvement aujourd'hui rend un autre grand mouvement demain plus probable, et les jours calmes viennent par séries. Ce regroupement de la volatilité est l'un des faits les mieux établis en finance, et il met à mal l'hypothèse d'une volatilité constante derrière la VaR simple. GARCH (hétéroscédasticité conditionnelle autorégressive généralisée) modélise directement la variance conditionnelle, de sorte que le risque et les prix d'options suivent l'humeur du marché.

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

La récurrence GARCH(1,1)

Le modèle a trois paramètres : var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]. omega est un plancher, alpha la réaction au choc de la veille et beta la persistance de la variance de la veille. Chaque jour dépendant du précédent, la série se construit avec une boucle, à partir d'une variance initiale. Pour les actions au quotidien, alpha vaut typiquement 0,05 à 0,1 et beta 0,85 à 0,95 : un choc s'estompe lentement.

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

Où elle se stabilise, et la prévision

Si alpha + beta < 1, la variance est à retour à la moyenne et se stabilise au niveau de long terme omega / (1 - alpha - beta) ; si elle vaut 1 ou plus, elle ne le fait jamais. Pour les prévisions, la variance de demain est connue à partir du rendement et de la variance d'aujourd'hui. Au-delà, le rendement au carré attendu égale la variance, donc chaque pas vaut omega + (alpha + beta) * précédent : la prévision glisse du niveau d'aujourd'hui vers celui de long terme, à une vitesse fixée par la persistance 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])])

Choisir les paramètres : la vraisemblance

Les paramètres ne sont pas connus ; on les estime sur les données par maximum de vraisemblance. Si les rendements sont normaux de variance var[t], la log-vraisemblance d'un jour est -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]). Leur somme donne un nombre qui note un choix de (omega, alpha, beta) : les meilleurs paramètres rendent les rendements observés les plus probables. Un optimiseur (par exemple scipy.optimize.minimize) essaie de nombreux triplets et garde le meilleur.

# 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

• La volatilité se regroupe : GARCH(1,1) construit la variance par une boucle, var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1].
• Avec alpha + beta < 1, elle revient vers omega / (1 - alpha - beta) ; les prévisions utilisent omega + (alpha + beta) * précédent.
• Les paramètres se choisissent par maximum de vraisemblance, sum(-0.5 * (log(2 * pi * var) + r**2 / var)) ; utilisez des chocs à queues épaisses pour les queues.

À suivre : le mouvement brownien, la marche aléatoire en temps continu qui sous-tend les modèles de valorisation d'options.

Lab : ce qui est noté

• Dérouler la récurrence avec garch_variance25 pts
• Trouver où elle se stabilise avec long_run_variance25 pts
• Prévoir les jours suivants avec garch_forecast25 pts
• Noter l'ajustement avec garch_loglik25 pts