GARCH-Volatilitätsmodelle

Volatilität tritt in Clustern auf

Tägliche Renditen sind nahezu unvorhersagbar, ihre Größe aber nicht: Eine große Bewegung heute macht eine weitere große morgen wahrscheinlicher, und ruhige Tage kommen in Serien. Dieses Volatilitäts-Clustering gehört zu den bestgesicherten Tatsachen der Finanzwelt und bricht die Annahme konstanter Volatilität hinter dem einfachen VaR. GARCH (verallgemeinerte autoregressive bedingte Heteroskedastizität) modelliert die bedingte Varianz direkt, sodass Risiko und Optionspreise der Stimmung des Marktes folgen können.

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

Die GARCH(1,1)-Rekursion

Das Modell hat drei Parameter: var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1]. omega ist ein Sockel, alpha die Reaktion auf den Schock von gestern und beta die Persistenz der gestrigen Varianz. Da jeder Tag vom vorigen abhängt, baut man die Reihe mit einer Schleife auf, ausgehend von einer Anfangsvarianz. Typisch für tägliche Aktienwerte sind alpha um 0,05 bis 0,1 und beta nahe 0,85 bis 0,95: Ein Schock klingt langsam ab.

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

Wo sie sich einpendelt, und die Prognose

Ist alpha + beta < 1, ist die Varianz mean-reverting und pendelt sich beim langfristigen Niveau omega / (1 - alpha - beta) ein; ist es 1 oder mehr, nie. Für Prognosen ist die Varianz von morgen aus Rendite und Varianz von heute bekannt. Danach entspricht das erwartete quadrierte Rendite der Varianz, jeder Schritt ist also omega + (alpha + beta) * vorheriger: Die Prognose gleitet vom heutigen Niveau zum langfristigen, mit einer durch die Persistenz alpha + beta bestimmten Geschwindigkeit.

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])])

Die Parameter wählen: die Likelihood

Die Parameter sind nicht bekannt; sie werden aus den Daten per Maximum Likelihood geschätzt. Sind die Renditen normal mit Varianz var[t], ist die Log-Likelihood eines Tages -0.5 * (log(2 * pi * var[t]) + r[t]**2 / var[t]). Ihre Summe ist eine Zahl, die eine Wahl von (omega, alpha, beta) bewertet: Die besten Parameter machen die beobachteten Renditen am wahrscheinlichsten. Ein Optimierer (z. B. scipy.optimize.minimize) probiert viele Tripel und behält das beste.

# 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

• Volatilität tritt in Clustern auf: GARCH(1,1) baut die Varianz per Schleife, var[t] = omega + alpha * r[t-1]**2 + beta * var[t-1].
• Bei alpha + beta < 1 kehrt sie zu omega / (1 - alpha - beta) zurück; Prognosen nutzen omega + (alpha + beta) * vorheriger.
• Die Parameter wählt man per Maximum Likelihood, sum(-0.5 * (log(2 * pi * var) + r**2 / var)); für die Ränder Schocks mit Fat Tails verwenden.

Als Nächstes: die Brownsche Bewegung, der Random Walk in stetiger Zeit, der Optionspreismodellen zugrunde liegt.

Lab: Was bewertet wird

• Mit garch_variance die Rekursion ausführen25 Pkt.
• Mit long_run_variance finden, wo sie sich einpendelt25 Pkt.
• Mit garch_forecast die nächsten Tage vorhersagen25 Pkt.
• Mit garch_loglik die Anpassung bewerten25 Pkt.