Modelos avanzados: Heston, SABR y volatilidad local.

Por qué hacen falta otros modelos

Black-Scholes supone una volatilidad constante. La lección 3.4.5 mostró que el mercado no la cotiza así: hay sonrisa y estructura temporal. Una superficie paramétrica como SVI describe esos precios, pero no dice nada de cómo se moverán mañana, ni sirve para valorar una opción cuyo pago depende del camino (una barrera, una opción con inicio diferido) o para calcular coberturas coherentes. Para eso hace falta un modelo del subyacente que reproduzca la superficie de hoy. Hay dos familias:

flowchart TB
  B[Black-Scholes] --> L[Volatilidad local]
  B --> E[Vol. estocástica]
  L --> D[Dupire 1994]
  E --> H[Heston 1993]
  E --> S[SABR 2002]
  D --> SLV[Local estocástica]
  H --> SLV

Esta lección es de «concepto y código»: cada modelo se programa, se comprueba con QuantLib y se calibra sobre la superficie de ORX del 4 de junio de 2025 de los datos sintéticos del curso, cuya verdad es una SSVI conocida. Eso permite separar lo que es error del modelo de lo que es ruido de las cotizaciones.

La superficie de partida

Como en la lección 3.4.5: volatilidad implícita de cada opción fuera del dinero con bid positivo, con el forward y el tipo de cada corte de la verdad.

import math
import numpy as np
import pandas as pd
from scipy.special import ndtr
from scipy.optimize import least_squares

sup = pd.read_parquet("opciones3_superficie.parquet")
ver = pd.read_parquet("opciones3_verdad.parquet")


def bs_norm(k, T, s):
    """Call de Black-76 sin descontar y dividida por F."""
    sd = s * np.sqrt(T)
    d1 = -k / sd + sd / 2
    return ndtr(d1) - np.exp(k) * ndtr(d1 - sd)


def vega_norm(k, T, s):
    sd = s * np.sqrt(T)
    d1 = -k / sd + sd / 2
    return np.exp(-d1 * d1 / 2) * np.sqrt(T / (2 * math.pi))


def iv_norm(c, k, T):
    """Newton con salvaguarda (lección 3.4.4), call normalizada."""
    lo, hi = np.full(len(c), 1e-4), np.full(len(c), 5.0)
    s = np.full(len(c), 0.3)
    for _ in range(60):
        d = bs_norm(k, T, s) - c
        lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
        paso = s - d / np.maximum(vega_norm(k, T, s), 1e-12)
        mal = (paso <= lo) | (paso >= hi) | ~np.isfinite(paso)
        s = np.where(mal, (lo + hi) / 2, paso)
    return s


hoy = pd.Timestamp("2025-06-04")
x = sup[sup["fecha"] == hoy].merge(ver, on=["fecha", "vencimiento"])
x["T"] = x["dias"] / 365
x["k"] = np.log(x["strike"] / x["forward"])
x = x[np.where(x["tipo"] == "C", x["k"] >= 0, x["k"] < 0)
      & (x["bid"] > 0)].copy()
D = np.exp(-x["tipo_libre"] * x["T"])
medio = (x["bid"] + x["ask"]) / 2 / (D * x["forward"])
# put -> call por la paridad (en unidades de F, sin descontar)
x["c"] = np.where(x["tipo"] == "C", medio, medio + 1 - np.exp(x["k"]))
x["iv"] = iv_norm(x["c"].to_numpy(), x["k"].to_numpy(),
                  x["T"].to_numpy())
print(len(x), "opciones en", x["dias"].nunique(), "vencimientos")

341 opciones en 12 vencimientos

Son las mismas 341 opciones de los doce vencimientos de la lección 3.4.5. Para trabajar con una sola fórmula, cada precio se pasa a una call sin descontar en unidades del forward (las puts, por la paridad de la lección 3.4.1): así Heston, SABR y Black-76 se comparan sobre la misma magnitud.

Heston: la varianza también es aleatoria

Heston (1993) añade a Black-Scholes una segunda ecuación para la varianza instantánea $v_t$:

• $v_0$: varianza de hoy. $\theta$: varianza de largo plazo a la que vuelve. $\kappa$: velocidad de la vuelta (vida media $\ln 2/\kappa$ años).
• $\xi$: volatilidad de la varianza. Da curvatura a la sonrisa (alas más caras).
• $\rho$: correlación entre precio y varianza. Negativa, inclina la sonrisa: es el efecto apalancamiento, el skew.
• Condición de Feller, $2\kappa\theta \ge \xi^2$: si se cumple, la varianza nunca toca cero.

Su ventaja es que la función característica de $x_T = \ln(S_T/F)$ tiene forma cerrada. En la forma de Albrecher y otros (2007), que evita saltos de la rama del logaritmo complejo:

Con ella, la fórmula de Lewis (2001) da la call como una integral real, con $k = \ln(K/F)$ y el descuento $D_T$:

from scipy.integrate import quad
import QuantLib as ql


def heston_cf(u, T, v0, kappa, theta, xi, rho):
    """E[exp(i·u·ln(S_T/F))], forma de Albrecher y otros (2007)."""
    iu = 1j * u
    b = kappa - rho * xi * iu
    d = np.sqrt(b ** 2 + xi ** 2 * (iu + u * u))
    g = (b - d) / (b + d)
    e = np.exp(-d * T)
    C = kappa * theta / xi ** 2 * (
        (b - d) * T - 2 * np.log((1 - g * e) / (1 - g)))
    Dv = (b - d) / xi ** 2 * (1 - e) / (1 - g * e)
    return np.exp(C + Dv * v0)


def heston_call(F, K, T, DT, p):
    k = math.log(K / F)
    f = lambda u: (np.exp(-1j * u * k)
                   * heston_cf(u - 0.5j, T, *p)).real / (u * u + 0.25)
    integral = quad(f, 0, 200, limit=400)[0]
    return DT * (F - math.sqrt(F * K) / math.pi * integral)


p = (0.04, 2.0, 0.04, 0.5, -0.7)        # v0, κ, θ, ξ, ρ
F, r, dias = 100.0, 0.03, 183
T, DT = dias / 365, math.exp(-0.03 * dias / 365)

d0 = ql.Date(4, 6, 2025)
ql.Settings.instance().evaluationDate = d0
dc = ql.Actual365Fixed()
plana = lambda t: ql.YieldTermStructureHandle(
    ql.FlatForward(d0, t, dc, ql.Continuous))
proc = ql.HestonProcess(plana(r), plana(r),        # q = r: F = S
                        ql.QuoteHandle(ql.SimpleQuote(F)), *p)
motor = ql.AnalyticHestonEngine(ql.HestonModel(proc))
for K in (80.0, 100.0, 120.0):
    op = ql.VanillaOption(ql.PlainVanillaPayoff(ql.Option.Call, K),
                          ql.EuropeanExercise(d0 + dias))
    op.setPricingEngine(motor)
    propio = heston_call(F, K, T, DT, p)
    print(f"K = {K:.0f}: propio {propio:.6f}, QuantLib "
          f"{op.NPV():.6f}, diferencia {propio - op.NPV():.1e}")

K = 80: propio 20.458189, QuantLib 20.458189, diferencia -9.4e-13
K = 100: propio 5.109865, QuantLib 5.109865, diferencia -2.0e-12
K = 120: propio 0.158180, QuantLib 0.158180, diferencia -2.5e-12

Las dos implementaciones coinciden hasta la duodécima cifra: la diferencia es el error de la integración numérica. Una advertencia práctica: el plazo de QuantLib sale de las fechas (183 días naturales con Actual365Fixed), así que el propio tiene que usar exactamente el mismo $T$; con $T = 0{,}5$ la diferencia sería de milésimas.

Calibrar Heston a la superficie

Para calibrar hay que valorar cientos de opciones en cada iteración. quad es exacta pero lenta; con una cuadratura de Gauss-Legendre fija (160 nodos entre 0 y 150) la integral se vectoriza para todas las opciones a la vez. El objetivo es el error en volatilidad, aproximado por el error de precio dividido por la vega, que es mucho más barato que invertir la volatilidad en cada paso:

nodos, pesos = np.polynomial.legendre.leggauss(160)
U, W = (nodos + 1) * 75, pesos * 75          # de [-1, 1] a [0, 150]


def heston_norm(k, T, p):
    """Call de Heston sin descontar, en unidades de F."""
    f = (np.exp(-1j * U * k[:, None])
         * heston_cf(U - 0.5j, T[:, None], *p)).real / (U * U + 0.25)
    return 1 - np.exp(k / 2) / math.pi * (f * W).sum(axis=1)


k, T, iv = (x[c].to_numpy() for c in ("k", "T", "iv"))
objetivo, vg = bs_norm(k, T, iv), vega_norm(k, T, iv)
res = least_squares(
    lambda p: (heston_norm(k, T, p) - objetivo) / vg,
    [0.05, 3.0, 0.05, 0.8, -0.7],
    bounds=([1e-4, 0.1, 1e-4, 0.05, -0.99], [1, 30, 1, 5, 0.99]))
v0, kappa, theta, xi, rho_h = ph = res.x
print(f"v0 = {v0:.4f} (vol {math.sqrt(v0):.1%}), κ = {kappa:.2f}, "
      f"θ = {theta:.4f} (vol {math.sqrt(theta):.1%})")
print(f"ξ = {xi:.3f}, ρ = {rho_h:.3f}, "
      f"Feller 2κθ − ξ² = {2 * kappa * theta - xi ** 2:.2f}")
x["iv_heston"] = iv_norm(heston_norm(k, T, ph), k, T)
err = 100 * (x["iv_heston"] - x["iv"])
print(err.groupby(x["dias"]).agg(
    lambda z: np.sqrt(np.mean(z ** 2))).round(2).to_dict())

v0 = 0.0691 (vol 26.3%), κ = 7.55, θ = 0.0669 (vol 25.9%)
ξ = 1.995, ρ = -0.750, Feller 2κθ − ξ² = -2.97
{16: 2.51, 44: 0.52, 72: 0.76, 107: 0.87, 135: 0.86, 170: 0.78, 198: 0.71, 289: 0.49, 379: 0.6, 471: 0.88, 562: 1.28, 653: 1.61}

Heston calibrado arranca con una volatilidad del 26,3 % ($v_0$) que vuelve con $\kappa = 7{,}55$ (vida media de unas cinco semanas) a un nivel de largo plazo del 25,9 %, con $\xi = 1{,}995$ y $\rho = -0{,}750$. La condición de Feller no se cumple ($-2{,}97$), lo habitual al calibrar índices: para generar el skew empinado de los plazos cortos, el modelo necesita una volatilidad de la varianza muy alta. El error por vencimiento dibuja una U: 2,51 puntos a 16 días, entre 0,49 y 0,88 de 44 a 471 días y 1,61 a 653. Cinco parámetros no bastan para el skew de 16 días y el de casi dos años a la vez.

S = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"]) \
    .set_index("fecha").loc[hoy, "cierre"]
c44 = x[x["dias"] == 44].iloc[0]
r44 = c44["tipo_libre"]
q44 = r44 - math.log(c44["forward"] / S) / c44["T"]
proc = ql.HestonProcess(plana(r44), plana(q44),
                        ql.QuoteHandle(ql.SimpleQuote(S)),
                        0.05, 3.0, 0.05, 0.8, -0.7)
modelo = ql.HestonModel(proc)
motor = ql.AnalyticHestonEngine(modelo)
ayudas = []
for f in x.itertuples():
    K = S * math.exp((r44 - q44) * f.T) * math.exp(f.k)
    a = ql.HestonModelHelper(
        ql.Period(int(f.dias), ql.Days), ql.NullCalendar(), S, K,
        ql.QuoteHandle(ql.SimpleQuote(f.iv)), plana(r44), plana(q44),
        ql.BlackCalibrationHelper.ImpliedVolError)
    a.setPricingEngine(motor)
    ayudas.append(a)
modelo.calibrate(ayudas, ql.LevenbergMarquardt(),
                 ql.EndCriteria(500, 50, 1e-8, 1e-8, 1e-8))
th_q, ka_q, xi_q, rho_q, v0_q = modelo.params()
print(f"QuantLib: v0 = {v0_q:.4f}, κ = {ka_q:.2f}, θ = {th_q:.4f}, "
      f"ξ = {xi_q:.3f}, ρ = {rho_q:.3f}")
iv_q = iv_norm(heston_norm(k, T, (v0_q, ka_q, th_q, xi_q, rho_q)),
               k, T)
for nombre, ivm in [("propia", x["iv_heston"]), ("QuantLib", iv_q)]:
    e = 100 * np.sqrt(np.mean((ivm - x["iv"]) ** 2))
    print(f"Error cuadrático medio, calibración {nombre}: {e:.3f}")

QuantLib: v0 = 0.0710, κ = 9.16, θ = 0.0661, ξ = 2.261, ρ = -0.737
Error cuadrático medio, calibración propia: 1.126
Error cuadrático medio, calibración QuantLib: 1.075

QuantLib llega a parámetros distintos ($\kappa = 9{,}16$, $\xi = 2{,}261$) con un error ligeramente menor (1,075 frente a 1,126 puntos). Las dos cosas tienen explicación. QuantLib minimiza el error exacto en volatilidad, y la calibración propia lo aproxima con precio dividido por vega. Y la función objetivo es casi plana a lo largo de una cresta en la que $\kappa$ y $\xi$ suben juntos: una reversión más rápida compensa una varianza más volátil. Dos juegos de parámetros muy distintos dan casi los mismos precios, y por eso los parámetros calibrados de Heston saltan de un día a otro aunque la superficie apenas cambie.

SABR: un modelo por vencimiento

SABR (Hagan, Kumar, Lesniewski y Woodward, 2002) modela el forward de un vencimiento y su volatilidad:

$\beta$ fija cómo depende la volatilidad del nivel (1 lognormal, 0 normal), $\alpha$ es el nivel, $\rho$ la inclinación y $\nu$ la volatilidad de la volatilidad (la curvatura). No tiene reversión a la media ni enlaza unos vencimientos con otros: se calibra corte a corte. Su atractivo es que Hagan y otros dieron una aproximación cerrada de la volatilidad implícita; con $z = (\nu/\alpha)(FK)^{(1-\beta)/2}\ln(F/K)$ y $x(z) = \ln\big[(\sqrt{1 - 2\rho z + z^2} + z - \rho)/(1 - \rho)\big]$:

def hagan(K, F, T, a, b, rho, nu):
    K = np.asarray(K, float)
    lfk = np.log(F / K)
    fk = (F * K) ** ((1 - b) / 2)
    z = nu / a * fk * lfk
    xz = np.log((np.sqrt(1 - 2 * rho * z + z * z) + z - rho)
                / (1 - rho))
    zx = np.where(np.abs(z) < 1e-8, 1.0, z / np.where(xz == 0, 1, xz))
    den = fk * (1 + (1 - b) ** 2 / 24 * lfk ** 2
                + (1 - b) ** 4 / 1920 * lfk ** 4)
    corr = 1 + ((1 - b) ** 2 / 24 * a * a / fk ** 2
                + rho * b * nu * a / (4 * fk)
                + (2 - 3 * rho * rho) / 24 * nu * nu) * T
    return a / den * zx * corr


for K, b in [(90.0, 1.0), (105.0, 1.0), (105.0, 0.5)]:
    a = 0.2 if b == 1 else 2.0
    propio = float(hagan(K, 100.0, 0.5, a, b, -0.5, 0.6))
    # argumentos de QuantLib: strike, forward, T, α, β, ν, ρ
    qlv = ql.sabrVolatility(K, 100.0, 0.5, a, b, 0.6, -0.5)
    print(f"K = {K:.0f}, β = {b}: propio {propio:.8f}, "
          f"QuantLib {qlv:.8f}")

K = 90, β = 1.0: propio 0.21787385, QuantLib 0.21787385
K = 105, β = 1.0: propio 0.19353697, QuantLib 0.19353697
K = 105, β = 0.5: propio 0.19187243, QuantLib 0.19187243

La implementación propia coincide con ql.sabrVolatility en las ocho cifras mostradas, también con $\beta = 0{,}5$. Cuidado con el orden de los argumentos en QuantLib: strike, forward, plazo, $\alpha$, $\beta$, $\nu$ y $\rho$, con $\nu$ antes que $\rho$, al revés que en la notación habitual.

filas = []
for dias, g in x.groupby("dias"):
    F, T = g["forward"].iloc[0], g["T"].iloc[0]
    K, ivm = g["strike"].to_numpy(), g["iv"].to_numpy()
    s = least_squares(
        lambda q: hagan(K, F, T, q[0], 1.0, q[1], q[2]) - ivm,
        [0.2, -0.5, 1.0], bounds=([1e-3, -0.999, 1e-3], [3, 0.999, 20]))
    a, rho_s, nu = s.x
    x.loc[g.index, "iv_sabr"] = hagan(K, F, T, a, 1.0, rho_s, nu)
    filas.append({"dias": dias, "alfa": a, "rho": rho_s, "nu": nu,
                  "nu_raizT": nu * math.sqrt(T),
                  "rmse_sabr": 100 * np.sqrt(np.mean(s.fun ** 2)),
                  "rmse_heston": 100 * np.sqrt(np.mean(
                      (g["iv"] - x.loc[g.index, "iv_heston"]) ** 2))})
sabr = pd.DataFrame(filas).set_index("dias")
sabr.round(3)

       alfa    rho     nu  nu_raizT  rmse_sabr  rmse_heston
dias
16    0.228 -0.797  3.045     0.638      0.945        2.506
44    0.230 -0.786  1.834     0.637      1.001        0.522
72    0.231 -0.778  1.501     0.667      0.746        0.759
107   0.229 -0.760  1.309     0.709      0.499        0.875
135   0.229 -0.758  1.189     0.723      0.425        0.855
170   0.228 -0.756  1.077     0.735      0.365        0.780
198   0.228 -0.755  1.014     0.747      0.302        0.707
289   0.227 -0.753  0.856     0.762      0.231        0.489
379   0.225 -0.750  0.758     0.773      0.183        0.596
471   0.225 -0.752  0.682     0.775      0.166        0.882
562   0.223 -0.752  0.628     0.779      0.144        1.281
653   0.222 -0.754  0.583     0.780      0.128        1.609

$\alpha$ queda entre 0,222 y 0,231 en todos los plazos: con $\beta = 1$ es aproximadamente la volatilidad ATM. $\rho$ va de −0,80 a −0,75 y $\nu$ cae de 3,0 a 0,58, mientras $\nu\sqrt{T}$ se mueve mucho menos, de 0,64 a 0,78: la superficie escala con la desviación típica (lección 3.4.5). El error de SABR baja con el plazo, de 0,94 y 1,00 puntos en los dos primeros cortes a 0,13 en el último. En los cortes cortos, los strikes llegan a 10 desviaciones típicas por debajo del forward, y la forma de Hagan con $\beta = 1$ no sigue a la vez el ala y el centro: deja errores que cambian de signo a lo largo de la sonrisa, de ±1,4 puntos en el centro y hasta 2,9 en el ala más lejana. Heston, con un solo juego de parámetros para todo, gana a SABR solo a 44 días (0,52 frente a 1,00).

La figura compara los dos errores por plazo: la U de Heston, destacada, frente a la caída continua de SABR. Un modelo con parámetros por vencimiento ajusta mejor cada corte; uno con parámetros globales es coherente entre plazos y sirve para valorar productos que dependen de varios. La figura de la cabecera enseña el detalle: a 16 días Heston se queda 6,9 puntos corto en el ala de las puts y 2 puntos alto en el dinero; a 289 días, los dos modelos casi tocan los puntos.

Volatilidad local: la superficie de Dupire

El camino contrario es no añadir ninguna fuente de aleatoriedad y dejar que la volatilidad sea una función determinista del precio y del tiempo, $dS_t/S_t = (r - q)\,dt + \sigma_L(S_t, t)\,dW_t$. Dupire (1994) demostró que hay una única función $\sigma_L$ que reproduce todos los precios europeos de hoy, y que se lee de la propia superficie. Escrita en varianza total $w(k, T)$ y moneyness respecto al forward de cada plazo (Gatheral, 2006):

El denominador es la $g(k)$ de la lección 3.4.5 y el numerador, la pendiente de calendario: sin arbitraje, los dos son positivos y la volatilidad local existe. Con arbitraje, sale negativa o infinita, y por eso se calcula sobre una superficie limpia. Aquí, la SSVI verdadera del día (su ρ, η = 0,95 y γ = 0,5, que la lección 3.4.5 recuperó con el ajuste) con $\theta_T$ interpolado linealmente entre vencimientos:

v = ver[ver["fecha"] == hoy]
RHO, ETA, GAM = float(v["svi_rho"].iloc[0]), 0.95, 0.5
T_n = np.r_[0, v["dias"] / 365]
TH_n = np.r_[0, v["theta"]]


def w_sup(k, T):
    th = np.interp(T, T_n, TH_n)
    phi = ETA / (th ** GAM * (1 + th) ** (1 - GAM))
    return th / 2 * (1 + RHO * phi * k
                     + np.sqrt((phi * k + RHO) ** 2 + 1 - RHO ** 2))


def var_local(k, T, h=1e-4):
    w = w_sup(k, T)
    w1 = (w_sup(k + h, T) - w_sup(k - h, T)) / (2 * h)
    w2 = (w_sup(k + h, T) - 2 * w + w_sup(k - h, T)) / h ** 2
    wt = (w_sup(k, T + h) - w_sup(k, T - h)) / (2 * h)
    g = (1 - k * w1 / (2 * w)) ** 2 - w1 ** 2 / 4 * (1 / w + 0.25) \
        + w2 / 2
    return wt / g


kk = np.array([-0.3, -0.2, -0.1, 0.0, 0.1, 0.2])
T30 = 30 / 365
print(pd.DataFrame({
    "k": kk, "implícita": 100 * np.sqrt(w_sup(kk, T30) / T30),
    "local": 100 * np.sqrt(var_local(kk, T30))}).round(2)
    .to_string(index=False))
for dias in (8, 30, 180):
    T, h = dias / 365, 0.005
    p_loc = (np.sqrt(var_local(h, T))
             - np.sqrt(var_local(-h, T))) / (2 * h)
    p_imp = (np.sqrt(w_sup(h, T) / T) - np.sqrt(w_sup(-h, T) / T)) \
        / (2 * h)
    print(f"{dias} días: pendiente ATM local {p_loc:.2f}, implícita "
          f"{p_imp:.2f}, razón {p_loc / p_imp:.2f}")

   k  implícita  local
-0.3      47.37  65.08
-0.2      40.52  54.41
-0.1      32.34  40.91
 0.0      22.13  22.09
 0.1      16.65  16.52
 0.2      18.96  22.79
8 días: pendiente ATM local -3.89, implícita -2.11, razón 1.84
30 días: pendiente ATM local -2.04, implícita -1.09, razón 1.86
180 días: pendiente ATM local -0.80, implícita -0.44, razón 1.80

En el dinero las dos coinciden (22,1 %), pero hacia las puts la local sube mucho más deprisa: 65,1 % frente a 47,4 % en $k = -0{,}3$. La pendiente en el dinero de la local es 1,84, 1,86 y 1,80 veces la implícita a 8, 30 y 180 días: es la regla de Derman, que dice que la implícita de un strike es aproximadamente la media de la local en el camino entre el precio actual y ese strike, y por eso tiene la mitad de pendiente. Hacia las calls, la local baja por debajo de la implícita y vuelve a subir antes (22,8 % frente a 19,0 % en $k = 0{,}2$).

La figura dibuja las dos curvas a 30 días: la local, destacada, es la implícita «desenrollada».

rng = np.random.default_rng(347)
N, T44, pasos = 100_000, 44 / 365, 176      # 4 pasos por día
dt = T44 / pasos
rejilla = np.linspace(-2.0, 1.5, 701)
X = np.zeros(N)                               # ln(S_t / F_t)
for i in range(pasos):
    s2 = np.interp(X, rejilla, var_local(rejilla, (i + 0.5) * dt))
    X += -0.5 * s2 * dt + np.sqrt(s2 * dt) * rng.standard_normal(N)
ST = np.exp(X) / np.exp(X).mean()          # martingala exacta
filas = []
for k_ in (-0.2, -0.1, -0.05, 0.0, 0.05):
    c = np.maximum(ST - math.exp(k_), 0).mean()
    filas.append({"k": k_,
                  "MC": 100 * iv_norm(np.array([c]), np.array([k_]),
                                      np.array([T44]))[0],
                  "superficie": 100 * math.sqrt(w_sup(k_, T44) / T44)})
mc = pd.DataFrame(filas)
mc["dif"] = mc["MC"] - mc["superficie"]
print(f"E[S_T/F] antes de normalizar: {np.exp(X).mean():.5f}")
mc.round(2)

E[S_T/F] antes de normalizar: 0.99983
      k     MC  superficie   dif
0 -0.20  37.78       37.93 -0.15
1 -0.10  30.81       30.81  0.00
2 -0.05  26.72       26.69  0.03
3  0.00  22.29       22.22  0.07
4  0.05  18.29       18.22  0.07

Simular el índice con la volatilidad local de cada instante y valorar las calls a 44 días reproduce la sonrisa de la superficie con diferencias de −0,15 a 0,07 puntos, el error de Monte Carlo con 100.000 caminos y de la discretización en 176 pasos. La media de $S_T/F$ sale en 0,99983 antes de normalizar: el término $-\sigma^2/2$ del drift mantiene la martingala salvo por el error de muestreo. Esa es la propiedad que define el modelo de Dupire: por construcción, devuelve todos los precios europeos de hoy.

Lo que la superficie de hoy no decide

Volatilidad local y Heston calibrados a la misma superficie dan los mismos precios europeos de hoy (local exactamente, Heston aproximadamente) y, sin embargo, son modelos muy distintos:

• Skew futuro. En volatilidad local la sonrisa futura se aplana: la función $\sigma_L(S, t)$ está fijada hoy y, a plazos largos, su dependencia del precio es suave. En los modelos de volatilidad estocástica, la sonrisa se regenera. Una opción con inicio diferido o un cliquet valen distinto en los dos.
• Delta. Con volatilidad local, cuando el índice baja, la volatilidad sube por el propio camino de la función; el delta incluye ese efecto. Con Heston, la varianza se mueve por su cuenta, correlacionada. La cobertura de una opción (lección 3.6.2) cambia con el modelo.
• Calibración frente a estimación. Ningún parámetro calibrado es una medida estadística. El índice de los datos se simuló con una varianza tipo Heston con κ = 8, ξ = 0,40 y ρ = −0,72 bajo la medida real (README de los datos); calibrado a los precios, Heston pide ξ ≈ 2,0 y ρ = −0,75, porque los precios llevan una prima de riesgo y la superficie no tiene forma de Heston.

En la práctica de los bancos se combinan: los modelos de volatilidad local estocástica (SLV) usan una dinámica de Heston multiplicada por una función local que se calibra para que la superficie de hoy salga exacta.

Ejercicios

1. Con los parámetros de Heston calibrados, divide ξ por dos y calcula la volatilidad implícita del corte de 44 días en $k = -0{,}2$, $0$ y $0{,}1$. ¿Qué cambia en la sonrisa? ¿Se cumple ahora la condición de Feller?

kq = np.array([-0.2, 0.0, 0.1])
Tq = np.full(3, 44 / 365)
for nombre, pp in [("calibrado", ph),
                   ("ξ / 2", (v0, kappa, theta, xi / 2, rho_h))]:
    vol = 100 * iv_norm(heston_norm(kq, Tq, pp), kq, Tq)
    feller = 2 * pp[1] * pp[2] - pp[3] ** 2
    print(f"{nombre}: {np.round(vol, 2)}, Feller {feller:.2f}")

calibrado: [38.54 22.03 16.66], Feller -2.97
ξ / 2: [34.2  24.87 19.65], Feller 0.02

2. Calibra SABR con β = 0,5 en vez de 1 en los cortes de 16 y 289 días y compara el error con el de β = 1. ¿Por qué en renta variable el valor de β importa poco para ajustar un solo vencimiento?

for dias in (16, 289):
    g = x[x["dias"] == dias]
    F, T = g["forward"].iloc[0], g["T"].iloc[0]
    K, ivm = g["strike"].to_numpy(), g["iv"].to_numpy()
    s = least_squares(
        lambda q: hagan(K, F, T, q[0], 0.5, q[1], q[2]) - ivm,
        [0.2 * F ** 0.5, -0.5, 1.0],
        bounds=([1e-3, -0.999, 1e-3], [500, 0.999, 20]))
    print(f"{dias} días, β = 0,5: ρ = {s.x[1]:.3f}, ν = {s.x[2]:.3f}, "
          f"error {100 * np.sqrt(np.mean(s.fun ** 2)):.3f} "
          f"(β = 1: {sabr.loc[dias, 'rmse_sabr']:.3f})")

16 días, β = 0,5: ρ = -0.792, ν = 2.818, error 1.088 (β = 1: 0.945)
289 días, β = 0,5: ρ = -0.693, ν = 0.703, error 0.375 (β = 1: 0.231)

3. Un modelo de volatilidad local calibrado al 4 de junio de 2025 reproduce exactamente todos los precios europeos de ese día. Explica por qué eso no basta para valorar una opción que empieza dentro de tres meses con el strike igual al precio de ese día (forward start).