Le modèle de Heston

Quand la volatilité bouge aussi

Black-Scholes suppose une volatilité constante, mais la vraie volatilité monte dans les krachs et s'éteint dans le calme, et le marché des options le montre : la volatilité implicite varie selon le strike (le sourire ou skew). Le modèle de Heston garde le processus géométrique du prix mais rend aussi la *variance* v aléatoire : elle est rappelée vers un niveau de long terme theta à la vitesse kappa, secouée par un bruit de taille xi proportionnel à sqrt(v), et les chocs sur le prix et la variance ont une corrélation rho. Un rho négatif (le prix baisse quand la volatilité monte, l'effet de levier) produit le skew descendant observé sur les options d'actions.

# The Heston equations, in words:  dS = mu S dt + sqrt(v) S dW      dv = kappa (theta - v) dt + xi sqrt(v) dZ      corr(dW, dZ) = rho
kappa, theta, xi = 2.0, 0.04, 0.3             # reversion speed, long-run variance (a 20% volatility), volatility of variance
print(round(theta ** 0.5, 2))                 # 0.2    the long-run volatility
print(2 * kappa * theta, xi ** 2)             # 0.16 0.09   the Feller condition 2 kappa theta >= xi^2 holds: the variance stays positive

Le simuler : troncature et corrélation

Pour simuler, on avance par petits intervalles dt. La mise à jour de la variance est v + kappa * (theta - v) * dt + xi * sqrt(v * dt) * dz, mais un simple pas d'Euler peut pousser v sous zéro, où la racine carrée échoue. Le schéma de troncature complète utilise max(v, 0) partout où v figure au membre de droite (une valeur négative agit comme zéro) et conserve le résultat brut pour le pas suivant. La condition de Feller 2 * kappa * theta >= xi^2 est le régime de paramètres où le vrai processus n'atteint jamais zéro ; si elle échoue, la troncature compte encore plus. Les deux chocs doivent être corrélés : tirez des normales indépendantes z1, z2, posez dz = z1 et dw = rho * z1 + sqrt(1 - rho^2) * z2.

# One variance step with full truncation, and the correlated shocks for the price
import math

def variance_step(v, kappa, theta, xi, dz, dt):
    vp = max(v, 0.0)
    return v + kappa * (theta - vp) * dt + xi * math.sqrt(vp * dt) * dz

def correlated_normals(z1, z2, rho):
    return z1, rho * z1 + math.sqrt(1 - rho ** 2) * z2

print(round(variance_step(0.04, 2.0, 0.04, 0.3, 0.5, 1 / 252), 6))      # 0.04189   at the long-run level only the noise moves it
dz, dw = correlated_normals(1.0, 1.0, -0.7)
print(round(dz, 4), round(dw, 4))             # 1.0 0.0141   a negative rho: the price shock pulls against the variance shock
print(variance_step(-0.01, 2.0, 0.04, 0.3, 1.0, 1 / 252))               # -0.00968...   still negative, but it moved up by kappa * theta * dt: the negative value acted as 0

Une trajectoire complète et une vérification

Un pas met à jour ensemble les deux variables d'état : le prix utilise la variance au *début* du pas, S * exp((mu - v/2) * dt + sqrt(v * dt) * dw), puis la variance bouge. Répétez pour le nombre de pas de l'horizon et vous avez une trajectoire ; de nombreuses trajectoires donnent une distribution de prix terminaux, et moyenner un gain actualisé sur elles valorise une option. Une bonne vérification : avec la dérive mu, le prix terminal moyen doit être proche de S0 * exp(mu * T), quels que soient les paramètres de volatilité. Sinon, il y a un bug dans le pas de prix. La calibration de kappa, theta, xi, rho sur des prix d'options viendra dans les sujets suivants.

# A short Monte Carlo: simulate many paths, and check the average terminal price against S0 * exp(mu * T)
import math
import random

def simulate(s0, v0, mu, kappa, theta, xi, rho, steps, years, rng):
    dt = years / steps
    s, v = s0, v0
    for _ in range(steps):
        z1, z2 = rng.gauss(0, 1), rng.gauss(0, 1)
        dz, dw = z1, rho * z1 + math.sqrt(1 - rho ** 2) * z2
        vp = max(v, 0.0)
        s *= math.exp((mu - vp / 2) * dt + math.sqrt(vp * dt) * dw)
        v = v + kappa * (theta - vp) * dt + xi * math.sqrt(vp * dt) * dz
    return s

rng = random.Random(11)
finals = [simulate(100.0, 0.04, 0.05, 2.0, 0.04, 0.3, -0.7, 50, 1.0, rng) for _ in range(4000)]
mean = sum(finals) / len(finals)
print(abs(mean / (100.0 * math.exp(0.05)) - 1) < 0.02)       # True   the average terminal price follows S0 * exp(mu * T)
print(min(finals) > 0)                                       # True   prices stay positive

# A check that costs nothing: with xi = 0 and v0 = theta the variance never moves (Black-Scholes with volatility sqrt(theta))
import math

kappa, theta, xi = 2.0, 0.04, 0.0
v = theta
for z in (1.5, -2.0, 0.3, 0.9):
    vp = max(v, 0.0)
    v = v + kappa * (theta - vp) * (1 / 252) + xi * math.sqrt(vp / 252) * z
print(v == theta)                                    # True   no noise and already at the long-run level
print(round(math.sqrt(v), 2))                        # 0.2    a constant 20% volatility

• Heston rend la variance aléatoire : elle revient vers theta à la vitesse kappa, a un bruit xi * sqrt(v) et est corrélée au prix par rho (négatif pour les actions).
• Simulez avec troncature complète (max(v, 0) au membre de droite), corrélez les chocs par dw = rho z1 + sqrt(1 - rho^2) z2 et vérifiez la condition de Feller 2 kappa theta >= xi^2.
• Le pas de prix utilise la variance au début du pas ; vérifiez que le prix terminal moyen suit S0 * exp(mu * T) et que xi = 0 redonne Black-Scholes.

À suivre : la volatilité locale, une surface de volatilité qui dépend du prix et du temps.

Lab : ce qui est noté

• Garder la variance positive avec feller_condition25 pts
• Faire avancer la variance avec variance_step25 pts
• Corréler les chocs avec correlated_normals25 pts
• Faire avancer prix et variance avec heston_step25 pts