El modelo de Heston

Cuando la volatilidad también se mueve

Black-Scholes supone volatilidad constante, pero la volatilidad real sube en los desplomes y se apaga en la calma, y el mercado de opciones lo muestra: la volatilidad implícita varía según el strike (la sonrisa o skew). El modelo de Heston mantiene el proceso geométrico del precio pero hace aleatoria también la *varianza* v: la atrae hacia un nivel de largo plazo theta a velocidad kappa, la sacude un ruido de tamaño xi que escala con sqrt(v), y los choques del precio y de la varianza tienen correlación rho. Un rho negativo (el precio cae cuando sube la volatilidad, el efecto apalancamiento) produce el skew descendente de las opciones sobre acciones.

# 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

Simularlo: truncamiento y correlación

Para simular, se avanza en intervalos pequeños dt. La actualización de la varianza es v + kappa * (theta - v) * dt + xi * sqrt(v * dt) * dz, pero un paso de Euler simple puede llevar v por debajo de cero, donde falla la raíz cuadrada. El esquema de truncamiento completo usa max(v, 0) donde v aparece en el lado derecho (un valor negativo actúa como cero) y conserva el resultado bruto para el paso siguiente. La condición de Feller 2 * kappa * theta >= xi^2 es el régimen de parámetros donde el proceso verdadero nunca toca cero; si falla, el truncamiento importa aún más. Los dos choques deben estar correlacionados: saca normales independientes z1, z2, pon dz = z1 y 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

Una trayectoria completa y una comprobación

Un paso actualiza juntas las dos variables de estado: el precio usa la varianza al *inicio* del paso, S * exp((mu - v/2) * dt + sqrt(v * dt) * dw), y luego se mueve la varianza. Repite durante el número de pasos del horizonte y tienes una trayectoria; muchas trayectorias dan una distribución de precios finales, y promediar un pago descontado sobre ellas valora una opción. Una buena comprobación: con la deriva mu, el precio final medio debe estar cerca de S0 * exp(mu * T), sean cuales sean los parámetros de volatilidad. Si no, hay un error en el paso del precio. Calibrar kappa, theta, xi, rho a precios de opciones llegará en los temas siguientes.

# 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 hace aleatoria la varianza: revierte a theta a velocidad kappa, tiene ruido xi * sqrt(v) y se correlaciona con el precio por rho (negativo en acciones).
• Simula con truncamiento completo (max(v, 0) en el lado derecho), correlaciona los choques con dw = rho z1 + sqrt(1 - rho^2) z2 y comprueba la condición de Feller 2 kappa theta >= xi^2.
• El paso del precio usa la varianza al inicio del paso; comprueba que el precio final medio sigue S0 * exp(mu * T) y que xi = 0 reproduce Black-Scholes.

A continuación: volatilidad local, una superficie de volatilidad que depende del precio y del tiempo.

Lab: qué se califica

• Mantener la varianza positiva con feller_condition25 pts
• Avanzar la varianza con variance_step25 pts
• Correlacionar los choques con correlated_normals25 pts
• Mover precio y varianza con heston_step25 pts