Sonrisa, skew y estructura temporal; superficie de volatilidad con SVI.

Una volatilidad por strike y vencimiento

Si Black-Scholes describiera bien el mercado, todas las opciones de un subyacente darían la misma volatilidad implícita. No es así, y la forma en que no lo es tiene estructura. Fijado un vencimiento, la volatilidad implícita cambia con el strike: es la sonrisa o, en los índices de renta variable, más bien una mueca, el skew, porque las puts fuera del dinero cotizan mucho más caras que las calls equidistantes. Fijado un strike relativo, cambia con el plazo: es la estructura temporal. Las dos juntas forman la superficie de volatilidad.

La causa del skew es económica antes que matemática. Los índices caen más deprisa de lo que suben y su volatilidad sube cuando caen (el efecto apalancamiento), de modo que la distribución de los rendimientos tiene la cola izquierda gruesa; además, hay una demanda estructural de puts como seguro de las carteras. Las dos cosas encarecen las puts lejanas frente a lo que diría una lognormal. En los datos del curso, el índice ORX tiene precisamente esas dos propiedades (correlación −0,72 entre rendimiento y varianza y saltos a la baja), y su superficie verdadera tiene skew.

Tres convenciones para toda la lección:

• Moneyness $k = \ln(K/F)$, con $F$ el forward del vencimiento: $k < 0$ son strikes por debajo del forward.
• Varianza total $w(k, T) = \sigma^2_{\text{imp}}(k, T)\,T$, con $T$ en días naturales / 365. Es la magnitud natural para comparar plazos: la varianza de un plazo largo es la suma de las de los tramos.
• Delta de 25 en la convención de contado con dividendo continuo, $e^{-qT}N(d_1)$ en la call; la put de «25 delta» es la de delta −0,25.

flowchart TB
  A[Cadena semanal] --> B[Vol implícita OTM]
  B --> C[SVI por corte]
  C --> D{Sin arbitraje}
  D -- no --> C
  D -- sí --> E[SSVI global]
  E --> F[Interpolación]
  F --> G[Superficie usable]

El diagrama es el orden de la lección. Cada caja se contrasta con la superficie verdadera, que en estos datos se conoce.

La superficie de un miércoles

opciones3_superficie.parquet es la foto de los miércoles de la cadena de ORX con todos los vencimientos hasta dos años: los seis mensuales y los trimestrales. Son datos sintéticos; su superficie verdadera es una SSVI de Gatheral y Jacquier (2014) cuyos parámetros por corte están en opciones3_verdad.parquet. Primero, la volatilidad implícita de cada opción fuera del dinero con bid positivo, con el Newton vectorizado con salvaguarda de la lección 3.4.4. El forward y el tipo de cada corte se toman de la verdad para separar el problema del ajuste del de la paridad, que la lección anterior ya resolvió:

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

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


def black(F, K, T, D, s, tipo):
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    c = D * (F * ndtr(d1) - K * ndtr(d1 - sd))
    return np.where(tipo == "C", c, c - D * (F - K))


def vega(F, K, T, D, s):
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    return D * F * np.exp(-d1 * d1 / 2) * np.sqrt(T / (2 * math.pi))


def iv_vec(p, F, K, T, D, tipo):
    """Newton con salvaguarda de bisección (lección 3.4.4)."""
    lo, hi = np.full(len(p), 1e-4), np.full(len(p), 5.0)
    s = np.full(len(p), 0.3)
    for _ in range(60):
        d = black(F, K, T, D, s, tipo) - p
        lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
        paso = s - d / np.maximum(vega(F, K, T, D, 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"])
otm = np.where(x["tipo"] == "C", x["k"] >= 0, x["k"] < 0)
x = x[otm & (x["bid"] > 0)].copy()
x["iv"] = iv_vec(((x["bid"] + x["ask"]) / 2).to_numpy(),
                 x["forward"].to_numpy(), x["strike"].to_numpy(),
                 x["T"].to_numpy(),
                 np.exp(-x["tipo_libre"] * x["T"]).to_numpy(),
                 x["tipo"].to_numpy())
x["w"] = x["iv"] ** 2 * x["T"]
x.groupby("dias").agg(opciones=("iv", "size"), k_min=("k", "min"),
                      k_max=("k", "max")).round(3)

      opciones  k_min  k_max
dias
16          32 -0.414  0.103
44          40 -0.765  0.185
72          39 -0.767  0.245
107         37 -0.689  0.329
135         39 -0.692  0.354
170         38 -0.694  0.352
198         18 -0.622  0.322
289         19 -0.629  0.315
379         20 -0.636  0.363
471         20 -0.642  0.356
562         20 -0.649  0.350
653         19 -0.655  0.343

Doce vencimientos, de 16 a 653 días, con entre 18 y 40 opciones útiles por corte. Los cortes cortos llegan a $k = -0{,}77$, puts un 53 % por debajo del forward: tras la caída de abril de 2025 siguen listados los strikes bajos (un strike no desaparece hasta su vencimiento). Los de más de seis meses tienen la mitad de strikes porque su separación es doble.

Sonrisa, skew y estructura temporal

Antes de ajustar nada conviene saber leer una sonrisa con las tres cifras que usa el mercado:

• Volatilidad ATM: la del forward, $k = 0$. Mide el nivel.
• Risk reversal de 25 delta: $RR_{25} = \sigma_{25C} - \sigma_{25P}$. Mide la inclinación; en un índice, negativo.
• Butterfly de 25 delta: $BF_{25} = (\sigma_{25C} + \sigma_{25P})/2 - \sigma_{ATM}$. Mide la curvatura.

En el ATM, la varianza total verdadera es el parámetro theta de la verdad: $\sigma_{ATM} = \sqrt{\theta/T}$. La estructura temporal es esa volatilidad para cada vencimiento, y se ve mejor comparando días distintos:

def theta_en(g, T):
    """θ interpolado linealmente en T, con θ(0) = 0."""
    return np.interp(T, np.r_[0, g["dias"] / 365], np.r_[0, g["theta"]])


def atm(fecha, T):
    """Vol ATM verdadera a plazo T."""
    g = ver[ver["fecha"] == pd.Timestamp(fecha)]
    return np.sqrt(theta_en(g, T) / T)


vx = pd.read_csv("vol3_indices.csv", parse_dates=["fecha"])
plazos = np.array([30, 91, 182, 365, 548]) / 365
filas = {}
for f in ["2017-06-14", "2025-06-04", "2020-03-18"]:
    vorx = vx.loc[vx["fecha"] == f, "vorx"].item()
    filas[f"{f} (VORX {vorx:.1f})"] = 100 * atm(f, plazos)
pd.DataFrame(filas, index=["1m", "3m", "6m", "1a", "18m"]).round(1)

     2017-06-14 (VORX 11.1)  2025-06-04 (VORX 25.9)  2020-03-18 (VORX 72.8)
1m                     10.3                    22.1                    63.0
3m                     12.4                    22.4                    48.3
6m                     14.0                    22.4                    38.1
1a                     15.5                    22.0                    30.3
18m                    16.4                    21.6                    27.0

En calma, en junio de 2017 con el VORX (el índice tipo VIX de los datos, lección 3.4.6) en 11,1, la estructura sube del 10,3 % a un mes al 16,4 % a 18 meses: el mercado cobra que la calma no dure. En marzo de 2020 se invierte, del 63,0 % a un mes al 27,0 % a 18 meses: se espera que la crisis pase. En junio de 2025, dos meses después de la caída de abril, es casi plana, en torno al 22 %. Detrás está la reversión a la media de la varianza: los plazos cortos reflejan la varianza de hoy y los largos, su nivel de largo plazo. En el modelo de los datos, la varianza total ATM combina la varianza instantánea, un nivel lento y una media de largo plazo, multiplicadas por una prima de riesgo de varianza.

La figura dibuja las tres curvas de 10 a 630 días; la de marzo de 2020, destacada, cae a menos de la mitad entre un mes y un año.

mi = ver[ver["fecha"].dt.weekday == 2]
tabla = mi.groupby("fecha").apply(
    lambda g: pd.Series({"un_mes": np.sqrt(theta_en(g, 30 / 365)
                                           * 365 / 30),
                         "un_ano": np.sqrt(theta_en(g, 1.0))}))
invertida = tabla["un_mes"] > tabla["un_ano"]
print(f"{len(tabla)} miércoles; invertida: {invertida.mean():.1%}")
print(tabla[invertida].index.year.value_counts().sort_index()
      .to_dict())

619 miércoles; invertida: 23.4%
{2014: 2, 2015: 12, 2016: 17, 2017: 6, 2018: 19, 2019: 6, 2020: 20, 2021: 8, 2022: 28, 2023: 2, 2024: 7, 2025: 18}

El 23,4 % de los 619 miércoles tiene la volatilidad ATM a un mes por encima de la de un año (con θ interpolado linealmente en el plazo). Se concentran en los años de estrés del guion: 28 en 2022, 20 en 2020, 19 en 2018 y 18 en 2025. Una estrategia que vende volatilidad corta y compra larga (un calendario, lección 3.5.3) gana con el contango y pierde de golpe cuando la curva se invierte.

SVI: cinco parámetros por vencimiento

La parametrización SVI (stochastic volatility inspired), ideada en Merrill Lynch en 1999 y presentada en público por Gatheral en 2004 (así lo cuentan Gatheral y Jacquier, 2014), describe la varianza total de un corte con cinco parámetros:

Sus alas son rectas en varianza total, como exige la fórmula del momento de Lee: lejos del dinero, $w(k)$ no puede crecer más deprisa que $2|k|$. El ajuste se hace por mínimos cuadrados sobre la volatilidad (no sobre la varianza, que pesaría demasiado las alas), con límites para que los parámetros no salgan de su dominio:

from scipy.optimize import least_squares


def svi(k, a, b, rho, m, sig):
    return a + b * (rho * (k - m) + np.sqrt((k - m) ** 2 + sig ** 2))


LIM = ([-1, 1e-6, -0.999, -1, 1e-4], [1, 5, 0.999, 1, 2])


def ajusta_svi(k, iv, T):
    def resid(p):
        return np.sqrt(np.maximum(svi(k, *p), 1e-10) / T) - iv
    p0 = [0.5 * (iv ** 2 * T).min(), 0.1, -0.5, 0.0, 0.1]
    return least_squares(resid, p0, bounds=LIM).x


cols = ["svi_a", "svi_b", "svi_rho", "svi_m", "svi_sigma"]
ajustes, filas = {}, []
for dias, g in x.groupby("dias"):
    k, T = g["k"].to_numpy(), g["T"].iloc[0]
    p = ajusta_svi(k, g["iv"].to_numpy(), T)
    ajustes[dias] = p
    verdad = g[cols].iloc[0].to_numpy()
    iv_v = np.sqrt(svi(k, *verdad) / T)
    filas.append({"dias": dias, "rho": p[2], "rho_v": verdad[2],
                  "b": p[1], "b_v": verdad[1],
                  "ruido_max": 100 * np.abs(g["iv"] - iv_v).max(),
                  "error_max": 100 * np.abs(
                      np.sqrt(svi(k, *p) / T) - iv_v).max()})
pd.DataFrame(filas).set_index("dias").round(3)

        rho  rho_v      b    b_v  ruido_max  error_max
dias
16   -0.673 -0.662  0.021  0.022      0.338      0.091
44   -0.665 -0.662  0.037  0.037      0.412      0.087
72   -0.678 -0.662  0.046  0.047      0.132      0.098
107  -0.656 -0.662  0.058  0.057      0.103      0.044
135  -0.660 -0.662  0.064  0.064      0.054      0.030
170  -0.661 -0.662  0.072  0.072      0.044      0.021
198  -0.657 -0.662  0.078  0.077      0.027      0.013
289  -0.661 -0.662  0.092  0.092      0.026      0.008
379  -0.657 -0.662  0.105  0.104      0.023      0.019
471  -0.688 -0.662  0.112  0.114      0.030      0.014
562  -0.655 -0.662  0.123  0.123      0.018      0.012
653  -0.634 -0.662  0.133  0.131      0.041      0.007

Las columnas _v son la verdad. El ρ ajustado va de −0,634 a −0,688, alrededor del −0,662 verdadero, que es el mismo en todos los cortes porque la superficie es una SSVI; $b$ sale casi exacto. ruido_max es la mayor distancia entre la volatilidad de una opción y la verdad: hasta 0,41 puntos en el corte de 44 días, el ruido de las cotizaciones que se vio en la lección 3.4.4. error_max es la del SVI ajustado: 0,087 en ese mismo corte y menos de 0,1 en todos. Un modelo de cinco parámetros no puede seguir el ruido de cuarenta cotizaciones: actúa como un filtro.

La figura de la cabecera dibuja cinco de los doce cortes ajustados, con la volatilidad de cada opción del corte de 16 días como punto. La de 16 días, destacada, es la más empinada en $k$: baja del 62 % al 16 % en menos de medio punto de moneyness, mientras que la de 653 días tarda casi un punto entero en bajar del 35 % al 16 %. Medida en desviaciones típicas, la diferencia casi desaparece: $k = -0{,}41$ a 16 días son 9,0 desviaciones ($\sqrt{\theta} = 0{,}046$), y $k = -0{,}66$ a 653 días, solo 2,3 ($\sqrt{\theta} = 0{,}286$). En SSVI, la pendiente ATM de la volatilidad es $\partial\sigma/\partial k = \rho\,\varphi(\theta)\sqrt{\theta}/(2\sqrt{T}) \approx \rho\,\eta/(2\sqrt{T})$: crece como $1/\sqrt{T}$. A 16 días son −1,5 puntos de volatilidad por cada 1 % de $k$.

Con un SVI por corte ya se pueden calcular las medidas de mercado en cualquier delta. El strike de 25 delta se busca con Brent sobre $k$, porque la volatilidad de cada strike depende del propio strike:

from scipy.optimize import brentq

S = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"]) \
    .set_index("fecha").loc[hoy, "cierre"]


def strike_delta(p, F, T, r, objetivo):
    """k con delta de contado e^(-qT)·N(d1) (call) o −e^(-qT)·N(−d1)."""
    eq = math.exp(-r * T) * F / S          # e^(-qT)
    def f(k):
        sd = math.sqrt(svi(k, *p))
        d1 = -k / sd + sd / 2
        delta = eq * (ndtr(d1) if objetivo > 0 else ndtr(d1) - 1)
        return delta - objetivo
    return brentq(f, -1.5, 1.0)


g = x[x["dias"] == 44].iloc[0]
T44 = g["T"]
for nombre, p in [("ajuste", ajustes[44]),
                  ("verdad", g[cols].to_numpy(float))]:
    vol = lambda k: 100 * math.sqrt(svi(k, *p) / T44)
    k25p = strike_delta(p, g["forward"], T44, g["tipo_libre"], -0.25)
    k25c = strike_delta(p, g["forward"], T44, g["tipo_libre"], 0.25)
    print(f"{nombre}: ATM {vol(0):.2f}, 25P {vol(k25p):.2f} "
          f"(K = {g['forward'] * math.exp(k25p):.0f}), "
          f"25C {vol(k25c):.2f}, RR {vol(k25c) - vol(k25p):.2f}, "
          f"BF {(vol(k25c) + vol(k25p)) / 2 - vol(0):.2f}")

ajuste: ATM 22.21, 25P 27.52 (K = 4856), 25C 18.52, RR -9.01, BF 0.81
verdad: ATM 22.22, 25P 27.52 (K = 4856), 25C 18.51, RR -9.01, BF 0.80

La volatilidad ATM sale en el 22,21 %, frente al 22,22 % verdadero. La put de 25 delta, en el strike 4.856 (un 5,8 % por debajo del forward), está en el 27,52 %, y la call de 25 delta, en el 18,52 %: un risk reversal de −9,01 puntos y un butterfly de 0,81, frente a −9,01 y 0,80 de la verdad. El skew de un índice se lee en el risk reversal: la put de 25 delta cuesta 9 puntos de volatilidad más que la call de 25 delta.

Arbitraje estático: mariposa y calendario

Una superficie puede ajustar muy bien los precios observados y aun así ser imposible. Dos condiciones la hacen coherente con precios sin arbitraje:

• Mariposa (por corte). La segunda derivada del precio de la call respecto al strike es la densidad riesgo neutral descontada (Breeden y Litzenberger): una mariposa de strikes $K - h, K, K + h$ no puede valer menos de cero. Escrita con $w(k)$, la condición es que la función de Gatheral sea no negativa:

• Calendario (entre cortes). Con el forward como referencia, la varianza total tiene que crecer con el plazo a $k$ fijo: $\partial w / \partial T \ge 0$. Si no, una call larga valdría menos que una corta equivalente.

def g_svi(k, a, b, rho, m, sig):
    r = np.sqrt((k - m) ** 2 + sig ** 2)
    w = svi(k, a, b, rho, m, sig)
    w1 = b * (rho + (k - m) / r)
    w2 = b * sig ** 2 / r ** 3
    return (1 - k * w1 / (2 * w)) ** 2 \
        - w1 ** 2 / 4 * (1 / w + 0.25) + w2 / 2


def controles(ajustes, kk=np.linspace(-1.5, 0.6, 421)):
    dias = sorted(ajustes)
    w = np.array([svi(kk, *ajustes[d]) for d in dias])
    a, b, rho, m, sig = np.array([ajustes[d] for d in dias]).T
    return pd.DataFrame({
        "g_min": [g_svi(kk, *ajustes[d]).min() for d in dias],
        "w_min": a + b * sig * np.sqrt(1 - rho ** 2),
        "lee": b * (1 + np.abs(rho)),
        "dw_min_sig": np.r_[np.diff(w, axis=0).min(axis=1), np.nan]},
        index=pd.Index(dias, name="dias"))


controles(ajustes).round(4)

       g_min   w_min     lee  dw_min_sig
dias
16    0.2604  0.0012  0.0360      0.0019
44    0.2668  0.0033  0.0609      0.0019
72    0.2722  0.0055  0.0780      0.0025
107   0.2771  0.0083  0.0952      0.0019
135   0.2807  0.0105  0.1068      0.0023
170   0.2842  0.0132  0.1196      0.0018
198   0.2870  0.0152  0.1285      0.0058
289   0.2950  0.0219  0.1528      0.0055
379   0.3005  0.0282  0.1732      0.0048
471   0.3084  0.0343  0.1887      0.0058
562   0.3130  0.0403  0.2038      0.0056
653   0.3177  0.0463  0.2171         NaN

Los doce cortes ajustados con todos los strikes pasan los controles en todo el tramo de $k = -1{,}5$ a $0{,}6$, que va más allá de los strikes listados: $g$ mínima entre 0,26 y 0,32, varianza mínima positiva, condición de Lee muy por debajo de 2 y cada corte por debajo del siguiente (dw_min_sig, la menor diferencia de varianza total con el corte siguiente, siempre positiva).

¿Qué pasa si se ajusta con pocos datos? Un caso habitual es quedarse solo con los strikes cercanos al dinero, los más líquidos, porque las alas «tienen mucho ruido»:

cerca = {}
for dias, g in x[np.abs(x["k"]) < 0.06].groupby("dias"):
    cerca[dias] = ajusta_svi(g["k"].to_numpy(), g["iv"].to_numpy(),
                             g["T"].iloc[0])
c = controles(cerca)
print(c.round(4).to_string())
print("Cortes con g < 0:", list(c.index[c["g_min"] < 0]))
print("Con varianza mínima negativa:", list(c.index[c["w_min"] < 0]))
print("Cruzados con el siguiente:", list(c.index[c["dw_min_sig"] < 0]))

       g_min   w_min     lee  dw_min_sig
dias
16    0.2606  0.0012  0.0358      0.0025
44   -0.2289  0.0036  0.1778     -0.0798
72    0.2738  0.0044  0.0771      0.0031
107   0.2832  0.0075  0.0905      0.0021
135   0.2868  0.0098  0.1017      0.0029
170   0.2940  0.0127  0.1109      0.0035
198   0.3046  0.0220  0.1136      0.0100
289   0.1305  0.0320  0.2377     -0.1163
379   0.3270  0.0362  0.1485     -0.0152
471   0.2501  0.0106  0.2583     -0.0051
562   0.2320 -0.0132  0.2964     -0.0243
653   0.2779 -0.0058  0.2581         NaN
Cortes con g < 0: [44]
Con varianza mínima negativa: [562, 653]
Cruzados con el siguiente: [44, 289, 379, 471, 562]

Con solo 5 a 13 strikes por corte, el de 44 días tiene $g < 0$ (−0,23), los dos más largos tienen varianza mínima negativa (una volatilidad imaginaria en algún strike) y cinco se cruzan con el siguiente en algún punto del tramo. En el centro ajustan casi igual de bien: el problema está donde no hay datos, que es justo donde la superficie se usa para valorar opciones lejanas o para calcular un índice de volatilidad.

La figura traduce $g$ a densidad riesgo neutral de $k$ en el corte de 44 días. Con todos los strikes es una campana asimétrica, positiva en todo el tramo y con la cola izquierda larga (el skew). Con los strikes cercanos es casi igual en el centro, pero se hace negativa entre $k = 0{,}17$ y $0{,}35$, con un mínimo de −0,16: una probabilidad negativa de que el índice termine entre un 19 % y un 42 % por encima del forward. Una mariposa centrada ahí tendría precio negativo.

SSVI: la superficie entera con dos parámetros

SVI ajusta cada corte por separado, y nada obliga a que los cortes sean coherentes entre sí. La SSVI (surface SVI) de Gatheral y Jacquier (2014) parametriza la superficie entera con la varianza total ATM de cada plazo, $\theta_T$, y una función de forma común:

Con $\theta_T$ creciente en $T$ no hay arbitraje de calendario, y con $\eta\,(1 + |\rho|) \le 2$ (y $\gamma = 1/2$) tampoco de mariposa: es una condición suficiente que se comprueba con una línea. Cada corte de una SSVI es un SVI crudo, por eso los parámetros de la verdad vienen en forma SVI. El ajuste usa como $\theta_T$ la varianza ATM de los SVI de cada corte y busca $\rho$, $\eta$ y $\gamma$ comunes:

def w_ssvi(k, th, rho, eta, gam=0.5):
    phi = eta / (th ** gam * (1 + th) ** (1 - gam))
    return th / 2 * (1 + rho * phi * k
                     + np.sqrt((phi * k + rho) ** 2 + 1 - rho ** 2))


theta = {d: svi(0.0, *p) for d, p in ajustes.items()}
grupos = [(theta[d], g["k"].to_numpy(), g["iv"].to_numpy(),
           g["T"].iloc[0]) for d, g in x.groupby("dias")]


def resid(p):
    return np.concatenate([np.sqrt(w_ssvi(k, th, *p) / T) - iv
                           for th, k, iv, T in grupos])


fit = least_squares(resid, [-0.5, 0.5, 0.4],
                    bounds=([-0.999, 0.01, 0.01], [0.999, 5, 0.99]))
rho, eta, gam = fit.x
print(f"ρ = {rho:.4f} (verdad {x['svi_rho'].iloc[0]:.4f}), "
      f"η = {eta:.3f} (0,95), γ = {gam:.3f} (0,5)")
print(f"η·(1 + |ρ|) = {eta * (1 + abs(rho)):.3f} <= 2")
ecm = 100 * np.sqrt(np.mean(fit.fun ** 2))
print(f"Error cuadrático medio: {ecm:.3f} puntos")

ρ = -0.6617 (verdad -0.6616), η = 0.953 (0,95), γ = 0.499 (0,5)
η·(1 + |ρ|) = 1.584 <= 2
Error cuadrático medio: 0.037 puntos

SSVI recupera ρ = −0,6617 frente a −0,6616, η = 0,953 frente a 0,95 y γ = 0,499 frente a 0,5 (los valores del modelo de los datos, en su README), con un error cuadrático medio de 0,037 puntos. Tres parámetros de forma más los doce $\theta$ describen la superficie entera. $\eta(1 + |\rho|) = 1{,}584 < 2$: sin arbitraje de mariposa por construcción. Y los $\theta$ ajustados crecen con el plazo en los doce cortes (es la varianza ATM de cada SVI, ya comprobada al hacer los controles), así que tampoco hay arbitraje de calendario en el dinero.

Interpolar entre vencimientos

Los vencimientos listados son unos pocos; las estrategias y los índices piden plazos constantes (30 días para el VIX, lección 3.4.6) o el vencimiento de una opción a medida. Hay que interpolar, y la forma de hacerlo importa. La prueba más limpia es quitar un corte que se conoce, el de 44 días, reconstruirlo a partir de sus vecinos de 16 y 72 y compararlo con su verdad. Cuatro métodos, con el peso $\alpha = (T - T_1)/(T_2 - T_1)$:

1. Volatilidad implícita lineal en $T$, a $k$ fijo.
2. Varianza total lineal en $T$, a $k$ fijo.
3. Varianza total lineal en $T$, a desviaciones típicas fijas $z = k/\sqrt{\theta_T}$: cada vecino se evalúa en su propio $k = z\sqrt{\theta_{T_i}}$.
4. SSVI con $\theta_T$ interpolado linealmente y los $\rho$, $\eta$, $\gamma$ ajustados.

g = x[x["dias"] == 44]
k, T = g["k"].to_numpy(), g["T"].iloc[0]
verdad = np.sqrt(svi(k, *g[cols].iloc[0].to_numpy(float)) / T)
T1, T2 = 16 / 365, 72 / 365
p1, p2 = ajustes[16], ajustes[72]
al = (T - T1) / (T2 - T1)
th = (1 - al) * theta[16] + al * theta[72]
z = k / np.sqrt(th)

metodos = {
    "vol a k fijo": (1 - al) * np.sqrt(svi(k, *p1) / T1)
                    + al * np.sqrt(svi(k, *p2) / T2),
    "varianza a k fijo": np.sqrt(((1 - al) * svi(k, *p1)
                                  + al * svi(k, *p2)) / T),
    "varianza a k/raiz(θ) fijo": np.sqrt((
        (1 - al) * svi(z * np.sqrt(theta[16]), *p1)
        + al * svi(z * np.sqrt(theta[72]), *p2)) / T),
    "SSVI": np.sqrt(w_ssvi(k, th, rho, eta, gam) / T)}
err = pd.DataFrame({m: 100 * np.abs(v - verdad)
                    for m, v in metodos.items()})
err.agg(["mean", "max"]).T.round(3)

                            mean    max
vol a k fijo               1.499  4.914
varianza a k fijo          0.546  1.786
varianza a k/raiz(θ) fijo  0.062  0.075
SSVI                       0.076  0.116

Interpolar la volatilidad a $k$ fijo se equivoca hasta en 4,9 puntos, y la varianza total a $k$ fijo, hasta en 1,8. La razón es la de la figura de la cabecera: la forma de la sonrisa escala con la desviación típica, y entre 16 y 72 días la escala cambia mucho. A $k$ fijo se mezclan puntos de las dos sonrisas que no se corresponden: el mismo $k$ es el ala lejana a 16 días y casi el centro a 72. A $z = k/\sqrt{\theta}$ fijo, el error máximo baja a 0,075 puntos, y con SSVI a 0,12; el error que queda en SSVI viene de interpolar $\theta$ linealmente en $T$ (0,005990 frente a 0,005950 verdadero), porque la varianza ATM verdadera no es lineal en el plazo.

La figura ordena los cuatro métodos por su error máximo, en escala logarítmica: entre el peor y el mejor hay un factor de 65.

Lo que cambia en un sistema real

• El forward manda. Toda la superficie se construye en $k = \ln(K/F)$: con un forward mal estimado (dividendos, tipos) la sonrisa se desplaza y aparecen falsas violaciones de paridad. La paridad sobre los medios (lección 3.4.4) es la forma estándar de obtenerlo.
• Pesos del ajuste. Ajustar en volatilidad con pesos iguales da tanta importancia a una put lejana con horquilla de 2 puntos de volatilidad como a una ATM con 0,1. Se suele ponderar por la inversa de la horquilla en volatilidad o por la vega.
• Calibrar con restricciones. En producción no basta con comprobar el arbitraje a posteriori: se impone en el ajuste (límites como los de Lee y Gatheral-Jacquier, o un ajuste SSVI y después un SVI por corte que no se aleje de él).
• La superficie se mueve. Entre dos días, la sonrisa se desplaza con el índice de forma distinta según el régimen (sticky strike o sticky delta); una estrategia de skew (lección 3.6.2) depende de cuál se cumple.

Ejercicios

1. Comprueba la condición de Lee, $b(1 + |\rho|) \le 2$, y la de varianza mínima no negativa en los doce cortes ajustados con todos los strikes. ¿Qué corte se acerca más al límite de Lee y por qué?

c = controles(ajustes)
print(c[["lee", "w_min"]].round(4).T.to_string())
print("Máximo de Lee:", round(c["lee"].max(), 4), "en",
      c["lee"].idxmax(), "días")

dias      16      44      72      107     135     170     198     289     379     471     562     653
lee    0.0360  0.0609  0.0780  0.0952  0.1068  0.1196  0.1285  0.1528  0.1732  0.1887  0.2038  0.2171
w_min  0.0012  0.0033  0.0055  0.0083  0.0105  0.0132  0.0152  0.0219  0.0282  0.0343  0.0403  0.0463
Máximo de Lee: 0.2171 en 653 días

2. Repite el ajuste SVI del corte de 44 días ponderando cada residuo por la inversa de la media horquilla en volatilidad (calcúlala como en la lección 3.4.4) y compara el error máximo frente a la verdad con el del ajuste sin pesos.

g = x[x["dias"] == 44]
k, T, iv = g["k"].to_numpy(), g["T"].iloc[0], g["iv"].to_numpy()
D = np.exp(-g["tipo_libre"].iloc[0] * T)
vg = vega(g["forward"].to_numpy(), g["strike"].to_numpy(), T, D, iv)
semi = ((g["ask"] - g["bid"]) / 2).to_numpy() / vg
peso = 1 / np.maximum(semi, 1e-4)


def resid_p(p):
    return peso * (np.sqrt(np.maximum(svi(k, *p), 1e-10) / T) - iv)


pp = least_squares(resid_p, [0.001, 0.1, -0.5, 0.0, 0.1],
                   bounds=LIM).x
verdad = np.sqrt(svi(k, *g[cols].iloc[0].to_numpy(float)) / T)
for nombre, p in [("sin pesos", ajustes[44]), ("con pesos", pp)]:
    e = 100 * np.abs(np.sqrt(svi(k, *p) / T) - verdad)
    print(f"{nombre}: error máximo {e.max():.3f}, medio {e.mean():.3f}")

sin pesos: error máximo 0.087, medio 0.020
con pesos: error máximo 0.026, medio 0.010

3. Un compañero interpola la volatilidad a 30 días de la put de strike 4.500 tomando la volatilidad de esa misma put en los vencimientos de 16 y 44 días y haciendo la media ponderada por los días. Explica los dos errores del método y cómo lo harías.