--- title: "3.4.7 Heston, SABR y volatilidad local" description: "Heston con su función característica y QuantLib, SABR con la fórmula de Hagan y volatilidad local de Dupire, calibrados contra una superficie conocida." url: "https://pythonparatrading.com/curso-python-trading/nivel-3/heston-sabr-y-volatilidad-local/" lang: "es" image: "https://pythonparatrading.com/og/curso-python-trading/nivel-3/heston-sabr-y-volatilidad-local.png" --- # 3.4.7 · Modelos avanzados: Heston, SABR y volatilidad local > La superficie de volatilidad describe los precios de hoy; un modelo dice cómo se mueve. Esta lección programa los tres modelos de referencia, Heston, SABR y la volatilidad local de Dupire, los contrasta con QuantLib y los calibra sobre una superficie sintética cuya verdad se conoce, para ver qué ajusta cada uno, qué no y por qué los parámetros calibrados no son los de la serie histórica. - Curso: [Nivel 3: Derivados: futuros, spreads, forex y opciones](https://pythonparatrading.com/curso-python-trading/nivel-3/) - Módulo: 3.4 · Teoría y valoración de opciones - Duración: 75 min - Requisitos: [3.4.2 · Árboles binomiales, Black-Scholes-Merton, Black-76 y Monte Carlo](https://pythonparatrading.com/curso-python-trading/nivel-3/binomial-black-scholes-black-76-y-monte-carlo/), [3.4.4 · Volatilidad implícita: Newton-Raphson, Brent, py_vollib y QuantLib](https://pythonparatrading.com/curso-python-trading/nivel-3/volatilidad-implicita-newton-brent-py-vollib-y-quantlib/), [3.4.5 · Sonrisa, skew y estructura temporal; superficie de volatilidad con SVI](https://pythonparatrading.com/curso-python-trading/nivel-3/sonrisa-skew-y-superficie-de-volatilidad-svi/) ## Al terminar esta lección - Explicar qué añade cada modelo a Black-Scholes: varianza estocástica (Heston), volatilidad estocástica con elasticidad (SABR) o volatilidad función del precio y del tiempo (Dupire). - Valorar opciones europeas en Heston con su función característica y la fórmula de Lewis, y comprobarlo con el motor analítico de QuantLib. - Calibrar Heston a una superficie completa y SABR vencimiento a vencimiento, y medir su error frente a la volatilidad de mercado. - Obtener la volatilidad local de una superficie con la fórmula de Dupire en varianza total y comprobarla con una simulación de Monte Carlo. ## 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: ```mermaid 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 ``` | Modelo | Qué es aleatorio | Parámetros | Ajusta la superficie | Uso típico | |---|---|---|---|---| | Black-Scholes | Solo el precio | σ | No (una σ) | Cotizar en volatilidad, referencia | | Volatilidad local (Dupire) | Solo el precio; σ función de (S, t) | Una superficie entera | Exactamente, por construcción | Exóticas de renta variable, consistencia de precios | | Heston | Precio y varianza, correlacionados | v0, κ, θ, ξ, ρ | Aproximadamente, con 5 números | Exóticas, dinámica de la volatilidad | | SABR | Precio y volatilidad, correlacionados | α, β, ρ, ν por vencimiento | Muy bien cada vencimiento | Tipos de interés (swaptions), interpolar sonrisas | 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. ```python 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") ``` ```output 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$: $$ \frac{dS_t}{S_t} = (r - q)\,dt + \sqrt{v_t}\,dW^S_t, \qquad dv_t = \kappa(\theta - v_t)\,dt + \xi\sqrt{v_t}\,dW^v_t, \qquad dW^S dW^v = \rho\,dt . $$ - $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: $$ \varphi(u) = \exp\big(C(u, T) + D(u, T)\,v_0\big), \quad d = \sqrt{(\rho\xi iu - \kappa)^2 + \xi^2(iu + u^2)}, \quad g = \frac{\kappa - \rho\xi iu - d}{\kappa - \rho\xi iu + d}, $$ $$ C = \frac{\kappa\theta}{\xi^2}\left[(\kappa - \rho\xi iu - d)\,T - 2\ln\frac{1 - g e^{-dT}}{1 - g}\right], \qquad D = \frac{\kappa - \rho\xi iu - d}{\xi^2}\,\frac{1 - e^{-dT}}{1 - g e^{-dT}} . $$ 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$: $$ C = D_T\left[F - \frac{\sqrt{FK}}{\pi}\int_0^\infty \frac{\operatorname{Re}\left[e^{-iuk}\,\varphi(u - i/2)\right]}{u^2 + 1/4}\,du\right]. $$ ```python 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}") ``` ```output 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: ```python 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()) ``` ```output 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. ```python 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}") ``` ```output 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: $$ dF_t = \alpha_t F_t^{\beta}\,dW_t, \qquad d\alpha_t = \nu\,\alpha_t\,dZ_t, \qquad dW\,dZ = \rho\,dt . $$ $\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]$: $$ \sigma_{\text{imp}} \approx \frac{\alpha\,\big[1 + \big(\tfrac{(1-\beta)^2\alpha^2}{24 (FK)^{1-\beta}} + \tfrac{\rho\beta\nu\alpha}{4 (FK)^{(1-\beta)/2}} + \tfrac{2 - 3\rho^2}{24}\nu^2\big)T\big]} {(FK)^{(1-\beta)/2}\big[1 + \tfrac{(1-\beta)^2}{24}\ln^2\tfrac{F}{K} + \tfrac{(1-\beta)^4}{1920}\ln^4\tfrac{F}{K}\big]} \cdot\frac{z}{x(z)} . $$ ```python 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}") ``` ```output 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. ```python 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) ``` ```output 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). *Fig. — Error cuadrático medio de Heston y de SABR (β = 1) frente a la volatilidad de mercado, por vencimiento, ORX, 4 de junio de 2025, en puntos de volatilidad.* Dos curvas frente al plazo en días. Heston (destacado) tiene forma de U: 2,51 puntos a 16 días, por debajo de 0,9 entre 44 y 471 días y 1,61 a 653. SABR baja de forma continua, de 0,95 a 16 días a 0,13 a 653. 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): $$ \sigma_L^2(k, T) = \frac{\partial w / \partial T} {1 - \frac{k}{w}\frac{\partial w}{\partial k} + \frac{1}{4}\left(-\frac{1}{4} - \frac{1}{w} + \frac{k^2}{w^2}\right)\left(\frac{\partial w}{\partial k}\right)^2 + \frac{1}{2}\frac{\partial^2 w}{\partial k^2}} . $$ 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: ```python 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}") ``` ```output 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$). *Fig. — Volatilidad local de Dupire frente a la implícita a 30 días, superficie SSVI verdadera de ORX del 4 de junio de 2025: la local tiene casi el doble de pendiente en el dinero.* Dos curvas frente a k, de −0,4 a 0,25. Se cruzan en el dinero, en el 22 %. Hacia las puts, la local (destacada) sube mucho más deprisa: 65 % frente a 47 % en k = −0,3. Hacia las calls, la local baja menos y vuelve a subir antes: 23 % frente a 19 % en k = 0,2. La figura dibuja las dos curvas a 30 días: la local, destacada, es la implícita «desenrollada». ```python 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) ``` ```output 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? **Solución.** ```python 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}") ``` ```output calibrado: [38.54 22.03 16.66], Feller -2.97 ξ / 2: [34.2 24.87 19.65], Feller 0.02 ``` Con la mitad de $\xi$, la sonrisa se aplana: la volatilidad del ala de las puts baja de 38,5 % a 34,2 %, la del dinero sube de 22,0 % a 24,9 % y la del ala de las calls, de 16,7 % a 19,7 %. La volatilidad de la varianza da curvatura y, con $\rho$ negativa, pendiente; además, en Heston más $\xi$ abarata el dinero (la volatilidad implícita ATM es una función cóncava de la varianza). Feller pasa de −2,97 a 0,02: se cumple justo. 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? **Solución.** ```python 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})") ``` ```output 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) ``` Con $\beta = 0{,}5$ el error sube un poco (de 0,945 a 1,088 a 16 días y de 0,231 a 0,375 a 289) y $\rho$ se ajusta para compensar (−0,79 y −0,69). $\beta$ y $\rho$ generan los dos pendiente: $\beta < 1$ hace que la volatilidad suba cuando el precio baja, y $\rho$ negativa también. Con un solo vencimiento, los dos efectos casi se sustituyen y el ajuste no los separa bien; por eso $\beta$ se fija de antemano (1 en renta variable, a menudo 0,5 o menos en tipos) y se calibran $\alpha$, $\rho$ y $\nu$. 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*). **Solución.** Una opción con inicio diferido fija su strike dentro de tres meses al precio de ese día, así que su valor depende de la sonrisa que habrá entonces, condicionada al camino: del skew futuro. Los precios europeos de hoy solo fijan la distribución del precio en cada vencimiento por separado, no cómo se encadenan. Volatilidad local y volatilidad estocástica calibradas a la misma superficie dan distribuciones marginales iguales y sonrisas futuras distintas: la local, más plana (la función $\sigma_L$ se aplana con el plazo); la estocástica, regenerada. Valorar el producto exige elegir la dinámica, y contrastarla con precios de productos que la revelen (el propio *forward start*, cliquets u opciones sobre el VIX) o usar un modelo local estocástico. ## Autoevaluación 1. Heston calibrado a la superficie del 4 de junio de 2025 da una volatilidad de la varianza ξ ≈ 2,0, cuando la varianza del índice de los datos se simuló con ξ = 0,40. ¿Qué indica? - A. Que la calibración ha fallado - B. Que los parámetros calibrados son los que reproducen los precios (medida de valoración) con la estructura de Heston, no estimaciones estadísticas de la serie: para generar el skew de corto plazo, Heston necesita una volatilidad de la varianza mucho mayor - C. Que el índice tiene más volatilidad de la que se simuló - D. Que hay que usar más vencimientos *Respuesta B: Que los parámetros calibrados son los que reproducen los precios (medida de valoración) con la estructura de Heston, no estimaciones estadísticas de la serie: para generar el skew de corto plazo, Heston necesita una volatilidad de la varianza mucho mayor.* La superficie de los datos no es Heston: es una SSVI con prima de riesgo. Heston la aproxima con los parámetros que mejor reproducen los precios, y el skew empinado de los plazos cortos solo lo consigue con ξ grande. Calibrar y estimar son preguntas distintas. 2. ¿Qué hace la condición de Feller, 2κθ ≥ ξ², en el modelo de Heston? - A. Garantiza que no hay arbitraje en la superficie - B. Garantiza que la varianza no llega a cero; si no se cumple, puede tocar cero (y vuelve a subir), lo que complica la simulación pero no impide valorar con la función característica - C. Hace que la correlación sea negativa - D. Asegura que la calibración converge *Respuesta B: Garantiza que la varianza no llega a cero; si no se cumple, puede tocar cero (y vuelve a subir), lo que complica la simulación pero no impide valorar con la función característica.* Con 2κθ < ξ² el proceso de la varianza alcanza el cero. La fórmula semianalítica sigue siendo válida; los esquemas de Monte Carlo necesitan tratar el cero (truncado o reflejado). La calibración de la lección la incumple, como es habitual en índices. 3. Cerca del dinero y a plazo corto, ¿cómo se compara el skew de la volatilidad local con el de la implícita? - A. Son iguales - B. La local es aproximadamente el doble de empinada: la implícita es, en primera aproximación, una media de la local a lo largo del camino entre el contado y el strike - C. La local es la mitad - D. La local no tiene skew *Respuesta B: La local es aproximadamente el doble de empinada: la implícita es, en primera aproximación, una media de la local a lo largo del camino entre el contado y el strike.* Es la regla de Derman: la volatilidad implícita de un strike promedia la local entre el precio actual y ese strike. En los datos, la razón de pendientes en el dinero sale entre 1,80 y 1,86 a 8, 30 y 180 días. 4. En SABR con β = 1, el ν calibrado baja con el plazo de 3,0 (16 días) a 0,58 (653 días). ¿Qué patrón esconde? - A. Ninguno: es ruido de calibración - B. Que ν·√T cambia mucho menos (de 0,64 a 0,78): la superficie de los datos escala la sonrisa con la desviación típica, y un SABR por vencimiento lo reproduce con un ν que cae aproximadamente como 1/√T - C. Que la volatilidad de la volatilidad sube con el plazo - D. Que hay que usar β = 0 *Respuesta B: Que ν·√T cambia mucho menos (de 0,64 a 0,78): la superficie de los datos escala la sonrisa con la desviación típica, y un SABR por vencimiento lo reproduce con un ν que cae aproximadamente como 1/√T.* Con β = 1 y α parecido en todos los plazos, la curvatura y la pendiente de la sonrisa en k dependen de ν√T. La SSVI de los datos tiene la misma forma en desviaciones típicas en todos los plazos (lección 3.4.5), así que ν√T varía poco; lo que se aparta en los plazos cortos es el peso de las alas muy lejanas en el ajuste. ## Ideas clave - Heston hace estocástica la varianza (reversión κ a θ, volatilidad de la varianza ξ, correlación ρ con el precio); su función característica tiene forma cerrada y la fórmula de Lewis valora una europea con una integral. Coincide con QuantLib hasta la duodécima cifra. - Calibrado a toda la superficie, Heston deja errores cuadráticos de 0,5 a 2,5 puntos de volatilidad por vencimiento: no puede tener a la vez el skew empinado de 16 días y el suave de casi dos años. Sus parámetros calibrados (ξ ≈ 2,0) no son los de la serie (ξ = 0,40). - SABR se calibra vencimiento a vencimiento con la fórmula de Hagan, que coincide con QuantLib; con β = 1 ajusta cada corte con errores de 0,13 a 1,0 puntos, menores cuanto más largo es el plazo. - La volatilidad local de Dupire sale de la superficie implícita: su varianza es la derivada de w respecto al plazo dividida por g(k). Reproduce todos los precios de hoy (una simulación lo comprueba) y su skew cerca del dinero es 1,8 veces el implícito. - Ajustar los precios de hoy no basta: local y estocástica pueden dar la misma superficie y dinámicas distintas (skew futuro, cobertura), y ahí se diferencian para opciones exóticas y para el delta. ## Recursos ### Research - [Procesos Estocásticos - Paseos Aleatorios](https://pythonparatrading.com/procesos-estocasticos/) · Research-0005 · 12 min ### Documentación - [QuantLib 1.44: modelo de Heston (hestonmodel.hpp)](https://github.com/lballabio/QuantLib/blob/v1.44/ql/models/equity/hestonmodel.hpp) · github.com - [QuantLib-Python: motores de valoración (AnalyticHestonEngine)](https://quantlib-python-docs.readthedocs.io/en/latest/pricing_engines.html) · quantlib-python-docs.readthedocs.io ### Libros - [Heston, S. L. (1993): A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options. Review of Financial Studies 6(2), 327-343](https://doi.org/10.1093/rfs/6.2.327) · doi.org - [Hagan, P. S., Kumar, D., Lesniewski, A. S. y Woodward, D. E. (2002): Managing Smile Risk (artículo de Wilmott Magazine, copia en PDF)](https://www.next-finance.net/IMG/pdf/pdf_SABR.pdf) · next-finance.net - [Dupire, B. (1994): Pricing with a Smile. Risk, enero de 1994](https://www.risk.net/derivatives/equity-derivatives/1500211/pricing-with-a-smile) · risk.net ### Datos - [Datos del nivel 3: superficie semanal de opciones sobre ORX (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/opciones3_superficie.parquet) · Descarga · opciones3_superficie.parquet Lección anterior: [3.4.6 · Volatilidad realizada frente a implícita; prima de riesgo de volatilidad; VIX y VSTOXX](https://pythonparatrading.com/curso-python-trading/nivel-3/volatilidad-realizada-implicita-y-vix/) Lección siguiente: [3.5.1 · Estructuras direccionales y verticales: call, put, covered call, collar y spreads](https://pythonparatrading.com/curso-python-trading/nivel-3/opciones-direccionales-y-verticales/)