Volatilidad implícita: Newton-Raphson, Brent, py_vollib y QuantLib.

Del precio a la volatilidad

Un modelo de valoración como Black-Scholes-Merton o Black-76 (lección 3.4.2) toma cinco datos y devuelve un precio: forward (o contado y dividendos), strike, plazo, tipo de interés y volatilidad. Los cuatro primeros se observan; la volatilidad, no. El mercado resuelve el problema al revés: cotiza el precio y deja que cada uno deduzca la volatilidad que lo justifica. Esa es la volatilidad implícita: la σ que, metida en el modelo, reproduce el precio de mercado.

Para una opción europea sobre un índice, con el forward $F$ al vencimiento y el factor de descuento $D = e^{-rT}$, Black-76 da

y la volatilidad implícita es la raíz de $f(\sigma) = C_{\text{modelo}}(\sigma) - C_{\text{mercado}} = 0$. No hay fórmula cerrada para despejarla: hay que buscarla con un método numérico. Dos propiedades garantizan que la búsqueda tiene sentido:

• Monotonía. La derivada del precio respecto a σ es la vega (lección 3.4.3), $\mathcal{V} = D\,F\,\varphi(d_1)\sqrt{T} > 0$. El precio crece estrictamente con σ, así que hay como mucho una raíz.
• Límites. Cuando σ tiende a 0, el precio tiende al intrínseco descontado, $D\max(F - K, 0)$ en la call; cuando σ crece sin límite, a $D\,F$. Un precio fuera de ese intervalo viola los límites de arbitraje de la lección 3.4.1 y no tiene volatilidad implícita. Entre los dos, existe y es única.

Las convenciones importan tanto como el método. En esta lección: precios en puntos del índice por opción (el contrato vale 100 USD por punto, pero eso no cambia la volatilidad), plazo en días naturales / 365 porque las opciones de ORX vencen al cierre, tipo continuo en tanto por uno y volatilidad anual en tanto por uno (0,25 = 25 %). Con el forward no hace falta el dividendo: ya está dentro de $F$.

flowchart TB
  A[Cadena bid y ask] --> B[Forward y tipo]
  B --> C{Filtros}
  C -- pasa --> D[Precio medio]
  C -- no pasa --> X[Descartar]
  D --> E{Entre límites}
  E -- sí --> F[Raíz en sigma]
  E -- no --> X
  F --> G[Sonrisa]

El orden del diagrama es el de la lección: primero el forward, después los métodos numéricos y, al final, los filtros, que es donde se decide si el número que sale significa algo.

La cadena de un día y su forward

Los datos son la cadena diaria de opciones sobre ORX, un índice ficticio que imita al S&P 500, con opciones europeas que imitan a las SPX (opciones3_cadena.parquet). Son sintéticos y llevan la respuesta conocida: cada fila tiene su iv_modelo, la volatilidad implícita verdadera sin ruido, y los bid y ask se generaron alrededor del precio de esa volatilidad con una horquilla realista y un poco de ruido. Se puede comprobar, por tanto, cuánto se equivoca cada paso. El día elegido es el 4 de junio de 2025, con el vencimiento del 18 de julio:

import math
import numpy as np
import pandas as pd

cad = pd.read_parquet("opciones3_cadena.parquet")
ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"])
ver = pd.read_parquet("opciones3_verdad.parquet")

hoy, vto = pd.Timestamp("2025-06-04"), pd.Timestamp("2025-07-18")
dia = cad[(cad["fecha"] == hoy) & (cad["vencimiento"] == vto)]
S = ind.loc[ind["fecha"] == hoy, "cierre"].item()
dias = (vto - hoy).days
T = dias / 365            # días naturales: vence al cierre
print(f"ORX {S:.2f} · {dias} días · T = {T:.4f}")
print(dia["tipo"].value_counts().to_dict())
cols = ["tipo", "strike", "bid", "ask", "iv_modelo"]
dia[dia["strike"].isin([4000, 5150, 6000])][cols]

ORX 5137.51 · 44 días · T = 0.1205
{'C': 45, 'P': 45}
       tipo  strike     bid     ask  iv_modelo
516554    C  4000.0  1153.9  1164.2     0.4130
516570    C  5150.0   160.2   161.7     0.2231
516583    C  6000.0     0.5     0.6     0.1709
516599    P  4000.0     9.8    10.1     0.4130
516615    P  5150.0   155.2   156.7     0.2231
516628    P  6000.0   838.2   844.7     0.1709

Hay una call y una put por strike, de 2.400 a 7.200 puntos. La de 4.000 es una put un 22 % por debajo del índice; la de 5.150, casi en el dinero; la de 6.000, una call un 17 % por encima. La columna iv_modelo es la respuesta que hay que recuperar.

Falta el forward. En vez de calcularlo con el dividendo (que en la realidad se estima), se deduce de los propios precios con la paridad put-call de la lección 3.4.1, $C - P = D\,(F - K)$: una recta en $K$ con pendiente $-D$ y ordenada $D\,F$. Con los strikes a menos de un 5 % del índice:

medio = dia.assign(medio=(dia["bid"] + dia["ask"]) / 2)
tabla = medio.pivot(index="strike", columns="tipo", values="medio")
cerca = tabla[np.abs(tabla.index / S - 1) < 0.05]
pend, corte = np.polyfit(cerca.index, cerca["C"] - cerca["P"], 1)
D = -pend
F = corte / D
r = -math.log(D) / T
v = ver[(ver["fecha"] == hoy) & (ver["vencimiento"] == vto)].iloc[0]
print(f"{len(cerca)} strikes")
print(f"Paridad: F = {F:.2f}, r = {r:.3%}")
print(f"Verdad:  F = {v.forward:.2f}, r = {v.tipo_libre:.3%}")

10 strikes
Paridad: F = 5155.15, r = 4.214%
Verdad:  F = 5155.25, r = 4.203%

Con diez strikes, la paridad da un forward de 5.155,15 puntos, a una décima del verdadero (5.155,25), y un tipo del 4,214 % frente al 4,203 %. La diferencia es el ruido de las cotizaciones. El forward está por encima del contado (5.137,51) porque el tipo supera a la rentabilidad por dividendo. Es la misma idea que usa Cboe para el VIX (lección 3.4.6): el mercado de opciones dice cuál es su forward, y ese forward ya lleva dentro el dividendo esperado, que así no hace falta estimar.

Newton-Raphson: rápido y frágil

Newton-Raphson aprovecha que se conoce la derivada: en cada paso sustituye la curva por su tangente y salta al punto donde la tangente corta el precio objetivo,

Cerca de la raíz converge de forma cuadrática: el número de decimales correctos se duplica en cada iteración. Primero, el precio y la vega de Black-76 vectorizados con numpy:

from scipy.special import ndtr


def black(F, K, T, D, s, tipo):
    """Precio de Black-76 (por punto de índice)."""
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    call = D * (F * ndtr(d1) - K * ndtr(d1 - sd))
    return np.where(tipo == "C", call, call - D * (F - K))


def vega(F, K, T, D, s):
    """Derivada del precio respecto a σ (por unidad de σ)."""
    sd = s * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    return D * F * np.exp(-d1 * d1 / 2) / math.sqrt(2 * math.pi) \
        * np.sqrt(T)


fila = dia[(dia["strike"] == 5150) & (dia["tipo"] == "C")].iloc[0]
p = black(F, 5150.0, T, D, fila.iv_modelo, "C")
print(f"Call 5150 con iv_modelo: {p:.2f}"
      f" (bid {fila.bid}, ask {fila.ask})")
print(f"Vega: {vega(F, 5150.0, T, D, fila.iv_modelo) / 100:.2f}"
      " puntos de precio por punto de volatilidad")

Call 5150 con iv_modelo: 160.96 (bid 160.2, ask 161.7)
Vega: 7.09 puntos de precio por punto de volatilidad

Con la volatilidad verdadera, el precio de Black-76 queda dentro de la horquilla. La vega dice cuánto vale un punto de volatilidad: 7,09 puntos de precio, 709 USD por contrato. Leída al revés, es la conversión que importa en esta lección: media horquilla de 0,75 puntos equivale a unos ±0,11 puntos de volatilidad.

def newton(precio, K, tipo, s0, tol=1e-10, maxit=50):
    s, camino = s0, [s0]
    for _ in range(maxit):
        dif = float(black(F, K, T, D, s, tipo)) - precio
        if abs(dif) < tol:
            break
        s = s - dif / float(vega(F, K, T, D, s))
        camino.append(s)
        if not 0 < s < 10:          # fuera de lo razonable
            break
    return s, camino


def medio_de(K, tipo):
    f = dia[(dia["strike"] == K) & (dia["tipo"] == tipo)].iloc[0]
    return float(f.bid + f.ask) / 2


for K, tipo in [(5150.0, "C"), (4000.0, "P")]:
    s, camino = newton(medio_de(K, tipo), K, tipo, s0=0.20)
    print(K, tipo, [round(x, 4) for x in camino])

5150.0 C [0.2, 0.2231, 0.2231]
4000.0 P [0.2, 12.7831]

La call casi en el dinero converge en dos pasos. La put de 4.000 salta en el primero a σ = 12,78, un 1.278 % de volatilidad anual. La razón está en la tangente: con σ = 20 % y 44 días, el strike 4.000 queda a 3,65 desviaciones del forward, el precio del modelo es 0,0099 puntos y la vega, 0,79 por unidad de σ. El paso es (0,0099 − 9,95) / 0,79 ≈ −12,6, y Newton suma 12,6 a σ. En las alas de la sonrisa, un arranque fijo es una mala idea.

def arranque_mk(K):
    """Punto de inflexión del precio en σ: sqrt(2·|ln(F/K)|/T)."""
    return math.sqrt(2 * abs(math.log(F / K)) / T)


for K, tipo in [(4000.0, "P"), (3400.0, "P"), (6000.0, "C")]:
    s0 = arranque_mk(K)
    s, camino = newton(medio_de(K, tipo), K, tipo, s0)
    print(f"{K:.0f} {tipo}: arranque {s0:.3f}, "
          f"{len(camino) - 1} iteraciones, σ = {s:.6f}")

4000 P: arranque 2.052, 7 iteraciones, σ = 0.412856
3400 P: arranque 2.628, 8 iteraciones, σ = 0.502269
6000 C: arranque 1.587, 9 iteraciones, σ = 0.170921

Manaster y Koehler (1982) propusieron arrancar en $\sigma^* = \sqrt{2\,|\ln(F/K)|/T}$, el punto donde la vega es máxima (allí $d_1 d_2 = 0$) y la curva del precio en σ cambia de convexa a cóncava. Si la raíz está por debajo, la función es convexa entre las dos y la tangente queda siempre bajo la curva: cada paso se acerca a la raíz sin pasarse. Si está por encima, es cóncava y ocurre lo mismo desde el otro lado. Con ese arranque, Newton llega a $10^{-10}$ en el precio en 7 a 9 iteraciones en las tres opciones, la put de 4.000 incluida.

La figura dibuja los dos arranques sobre la put de 4.000, con el precio en escala logarítmica. Desde 2,05, los pasos bajan por la curva sin pasarse de la raíz. La tangente en 0,20 es tan plana que, en todo el tramo dibujado, no llega a 1,6 puntos de precio.

Brent: lento y seguro

Si la derivada es traicionera, se puede prescindir de ella. Bisección parte de un intervalo $[\sigma_a, \sigma_b]$ en el que la función cambia de signo y lo divide por la mitad en cada paso: nunca falla, pero gana solo un bit por iteración. El método de Brent combina bisección con interpolación (secante e interpolación cuadrática inversa): cuando la interpolación promete, la usa; cuando no, biseca. Converge casi tan deprisa como la secante y conserva la garantía de la bisección. En SciPy es brentq:

from scipy.optimize import brentq


def iv_brent(precio, K, tipo, a=1e-4, b=5.0):
    f = lambda s: float(black(F, K, T, D, s, tipo)) - precio
    return brentq(f, a, b, xtol=1e-12, full_output=True)


for K, tipo in [(5150.0, "C"), (4000.0, "P"), (6000.0, "C")]:
    s, res = iv_brent(medio_de(K, tipo), K, tipo)
    print(f"{K:.0f} {tipo}: σ = {s:.6f}, {res.iterations} iteraciones,"
          f" {res.function_calls} evaluaciones")

5150 C: σ = 0.223091, 6 iteraciones, 7 evaluaciones
4000 P: σ = 0.412856, 15 iteraciones, 16 evaluaciones
6000 C: σ = 0.170921, 19 iteraciones, 20 evaluaciones

Brent llega a la misma σ con 6, 15 y 19 iteraciones, una evaluación del precio en cada una (más la de arranque). Newton desde Manaster-Koehler necesitaba 7 a 9, con precio y vega en cada paso. Para una opción la diferencia no importa; para millones de filas, sí, y entonces se usa un método vectorizado (al final de la lección) o el de Jäckel. El intervalo $[10^{-4}, 5]$ es la única hipótesis: si el precio está fuera de los límites, brentq no encuentra cambio de signo y lanza ValueError, que es justo lo que se quiere.

Librerías: vollib, py_vollib_vectorized y QuantLib

vollib (antes py_vollib)

py_vollib fue durante años la librería de referencia en Python para Black, Black-Scholes y Black-Scholes-Merton con el algoritmo Let's Be Rational de Jäckel (2015), que obtiene la volatilidad con precisión de máquina en dos iteraciones gracias a una aproximación racional del punto de partida. Su versión instalada, la 1.0.12, ya no trae código: es un paquete de transición que depende de vollib (el nombre nuevo) y avisa al importarlo:

import warnings

with warnings.catch_warnings(record=True) as avisos:
    warnings.simplefilter("always")
    import py_vollib.black.implied_volatility as pv_iv
print(avisos[0].category.__name__, "·", avisos[0].message)

from vollib.black.implied_volatility import implied_volatility
print(pv_iv.implied_volatility is implied_volatility)

K, tipo = 4000.0, "P"
iv = implied_volatility(medio_de(K, tipo), F, K, r, T, tipo.lower())
print(f"vollib, put 4000: σ = {iv:.6f}")

DeprecationWarning · py_vollib is deprecated and will be removed in a future release; please import from vollib instead.
True
vollib, put 4000: σ = 0.412856

El aviso sale una vez, al importar. py_vollib.black.implied_volatility es el mismo objeto que en vollib: el paquete viejo solo redirige. El código nuevo debe importar de vollib, como en el resto de la lección. La función recibe el precio descontado, F, K, r, t y 'c' o 'p', y da la misma σ que Newton y Brent: 0,412856.

py_vollib_vectorized

py_vollib_vectorized envolvía las funciones de py_vollib con numba para calcular cadenas enteras de una vez. Su última versión, la 0.1.1, es de febrero de 2021, y no compila con las versiones actuales de numba (0.68) y de py_lets_be_rational (1.1.2, de mayo de 2026): la primera llamada falla al tipar las funciones. Se comprueba en un proceso aparte, para no dejar numba cargado en este:

import subprocess
import sys

prueba = """
import numpy as np
import py_vollib_vectorized as pvv
pvv.vectorized_implied_volatility_black(
    np.array([3.9]), np.array([100.0]), np.array([100.0]),
    np.array([0.03]), np.array([0.25]), np.array(["c"]))
"""
res = subprocess.run([sys.executable, "-c", prueba],
                     capture_output=True, text=True)
error = [l for l in res.stderr.splitlines() if "Error" in l][-1]
print("Código de salida:", res.returncode)
print(error)

Código de salida: 1
numba.core.errors.TypingError: Failed in nopython mode pipeline (step: nopython frontend)

El proceso termina con código 1 y un TypingError de numba. Leyendo el código instalado se ve por qué: la función compilada de py_vollib_vectorized llama a normalised_black_call de py_lets_be_rational, que en la versión 1.1.2 es una función de Python normal (allí la compilación con numba pasó a ser opcional, con la variable de entorno PY_LETS_BE_RATIONAL_ENABLE_NUMBA), y numba no puede tipar una función sin compilar dentro de una compilada. Activar esa variable tampoco basta: la compilación falla más adelante.

import os
import numba

os.environ["NUMBA_DISABLE_JIT"] = "1"   # numba como Python normal
numba.config.reload_config()            # ya estaba importado
import py_vollib_vectorized as pvv

otm = dia[np.where(dia["tipo"] == "P", dia["strike"] < F,
                   dia["strike"] >= F) & (dia["bid"] > 0)].copy()
otm["medio"] = (otm["bid"] + otm["ask"]) / 2
n = len(otm)
iv_vec = pvv.vectorized_implied_volatility_black(
    otm["medio"].to_numpy(), np.full(n, F), otm["strike"].to_numpy(),
    np.full(n, r), np.full(n, T), otm["tipo"].str.lower().to_numpy(),
    return_as="numpy")
print(n, "opciones fuera del dinero con bid > 0")
print(np.round(iv_vec[:4], 6))

40 opciones fuera del dinero con bid > 0
[0.214444 0.205999 0.19825  0.190714]

Con NUMBA_DISABLE_JIT=1, numba deja las funciones decoradas como Python normal y la librería funciona, aunque sin la velocidad que era su razón de ser. Dos detalles del código. La variable se lee al importar numba, y aquí numba ya estaba cargado (lo importa vollib a través de py_lets_be_rational): por eso reload_config(). Y al importarse, py_vollib_vectorized sustituye las funciones de vollib por sus versiones vectorizadas; implied_volatility se importó antes y sigue siendo la original. Una dependencia sin mantenimiento puede romperse con cualquier actualización del entorno; para cadenas grandes, un Newton vectorizado propio (más abajo) o QuantLib son apuestas más seguras.

QuantLib

QuantLib ofrece dos caminos. El directo es blackFormulaImpliedStdDev, que invierte la fórmula de Black con el forward y el descuento y devuelve la desviación típica total $\sigma\sqrt{T}$. El completo monta el instrumento (EuropeanOption) sobre un proceso de Black-Scholes-Merton con curvas de tipos y de dividendos, y llama a impliedVolatility; es el que se usa cuando la opción forma parte de una cartera valorada con QuantLib. Sus valores por defecto son una tolerancia de $10^{-4}$ y σ entre $10^{-7}$ y 411Nota 1 · Firma de VanillaOption::impliedVolatility en QuantLib 1.44, ql/instruments/vanillaoption.hpp: https://github.com/lballabio/QuantLib/blob/v1.44/ql/instruments/vanillaoption.hpp; aquí se le pasan $10^{-10}$, 200 evaluaciones y σ entre $10^{-4}$ y 5:

import QuantLib as ql

K, precio = 4000.0, medio_de(4000.0, "P")
sd = ql.blackFormulaImpliedStdDev(ql.Option.Put, K, F, precio, D)
print(f"blackFormulaImpliedStdDev: σ = {sd / math.sqrt(T):.6f}")
# desplazamiento 0, arranque 0,3, tolerancia 1e-12, 200 iteraciones
sd = ql.blackFormulaImpliedStdDev(ql.Option.Put, K, F, precio, D,
                                  0.0, 0.3, 1e-12, 200)
print(f"... con tolerancia 1e-12: σ = {sd / math.sqrt(T):.6f}")

ql.Settings.instance().evaluationDate = ql.Date(4, 6, 2025)
dc = ql.Actual365Fixed()
hoy_ql = ql.Settings.instance().evaluationDate
q = r - math.log(F / S) / T     # dividendo implícito en el forward
curva = lambda x: ql.YieldTermStructureHandle(
    ql.FlatForward(hoy_ql, x, dc, ql.Continuous))
vol = ql.BlackVolTermStructureHandle(
    ql.BlackConstantVol(hoy_ql, ql.TARGET(), 0.20, dc))
proc = ql.BlackScholesMertonProcess(
    ql.QuoteHandle(ql.SimpleQuote(S)), curva(q), curva(r), vol)
put = ql.EuropeanOption(ql.PlainVanillaPayoff(ql.Option.Put, K),
                        ql.EuropeanExercise(ql.Date(18, 7, 2025)))
put.setPricingEngine(ql.AnalyticEuropeanEngine(proc))
iv_ql = put.impliedVolatility(precio, proc, 1e-10, 200, 1e-4, 5.0)
print(f"EuropeanOption.impliedVolatility: σ = {iv_ql:.6f}")
print(f"Rentabilidad por dividendo implícita: {q:.3%}")

blackFormulaImpliedStdDev: σ = 0.412854
... con tolerancia 1e-12: σ = 0.412856
EuropeanOption.impliedVolatility: σ = 0.412856
Rentabilidad por dividendo implícita: 1.370%

Tres detalles. blackFormulaImpliedStdDev devuelve σ√T, no σ. Su tolerancia por defecto, $10^{-6}$ (se ve en la cabecera blackformula.hpp de QuantLib), deja mal la sexta cifra decimal: 0,412854; con $10^{-12}$ coincide con los demás. Y el instrumento completo necesita un proceso con contado, dividendo y tipo; la volatilidad constante de 0,20 con la que se monta no influye, porque impliedVolatility la sustituye mientras busca. El dividendo implícito en el forward de la paridad, 1,370 %, queda cerca del 1,3435 % del fichero del índice: la diferencia es el ruido con el que se estimaron forward y tipo.

Los cinco métodos, opción a opción

Sobre las opciones fuera del dinero del día, los cinco métodos deberían dar lo mismo. La comparación sirve de prueba de que el código propio está bien:

filas = []
for f in otm.itertuples():
    K, tipo, pm = f.strike, f.tipo, f.medio
    filas.append({
        "newton": newton(pm, K, tipo, arranque_mk(K))[0],
        "brent": iv_brent(pm, K, tipo)[0],
        "vollib": implied_volatility(pm, F, K, r, T, tipo.lower()),
        "quantlib": ql.blackFormulaImpliedStdDev(
            ql.Option.Call if tipo == "C" else ql.Option.Put,
            K, F, pm, D, 0.0, 0.3, 1e-12, 200) / math.sqrt(T)})
ivs = pd.DataFrame(filas, index=otm.index)
ivs["vectorized"] = iv_vec
dif = ivs.sub(ivs["brent"], axis=0).abs().max()
print("Diferencia máxima con Brent:")
print(dif.map(lambda x: f"{x:.1e}").to_string())

Diferencia máxima con Brent:
newton        4.6e-12
brent         0.0e+00
vollib        3.2e-13
quantlib      2.4e-12
vectorized    3.2e-13

Newton, vollib y la versión vectorizada coinciden con Brent por debajo de $10^{-11}$; QuantLib con tolerancia estricta, en $2{,}4 \cdot 10^{-12}$. Son diferencias de redondeo. El método numérico no es lo que separa una volatilidad buena de una mala: lo que la separa es el precio que se le da.

Lo que cuesta la horquilla

Con el método resuelto, la pregunta que importa es otra: ¿cuánto se parece la volatilidad que sale del precio medio a la verdadera? El medio no es el precio del modelo: está dentro de una horquilla, con ruido, y redondeado al tick (0,05 puntos por debajo de 3 y 0,10 por encima). Un error de precio $\Delta C$ se convierte en un error de volatilidad de aproximadamente $\Delta C / \mathcal{V}$, y la vega cae muy deprisa al alejarse del dinero. Por eso conviene medir también la media horquilla en volatilidad, $(a - b) / (2\mathcal{V})$: la banda dentro de la cual cualquier volatilidad es compatible con las cotizaciones.

otm["iv"] = ivs["brent"]
otm["k"] = np.log(otm["strike"] / F)
vg = vega(F, otm["strike"].to_numpy(), T, D,
          otm["iv_modelo"].to_numpy())
otm["error_pts"] = 100 * (otm["iv"] - otm["iv_modelo"])
otm["semi_pts"] = 100 * (otm["ask"] - otm["bid"]) / 2 / vg
sel = otm["strike"].isin([3400, 4000, 4600, 5150, 5600, 6000, 6200])
print(otm.loc[sel, ["tipo", "strike", "bid", "ask", "iv_modelo",
                    "iv", "error_pts", "semi_pts"]]
      .round(4).to_string(index=False))
print(f"Error medio absoluto: {otm['error_pts'].abs().mean():.3f}"
      " puntos")

tipo  strike   bid    ask  iv_modelo     iv  error_pts  semi_pts
   C  5600.0  11.1  11.50     0.1690 0.1692     0.0247    0.0731
   C  6000.0   0.5   0.60     0.1709 0.1709     0.0021    0.1719
   C  6200.0   0.1   0.15     0.1772 0.1763    -0.0867    0.2887
   P  3400.0   2.0   2.10     0.5018 0.5023     0.0469    0.1509
   P  4000.0   9.8  10.10     0.4130 0.4129    -0.0144    0.1150
   P  4600.0  41.9  42.70     0.3188 0.3191     0.0269    0.1014
   P  5150.0 155.2 156.70     0.2231 0.2233     0.0170    0.1057
Error medio absoluto: 0.032 puntos

Cerca del dinero el error es de centésimas de punto: la put de 5.150 recupera un 22,33 % frente al 22,31 % verdadero, con una media horquilla de 0,11 puntos. En la call de 6.200, con bid 0,10 y ask 0,15, la media horquilla sube a 0,29 puntos y el error a −0,09. En las 40 opciones, el error medio absoluto es de 0,032 puntos, y en todas es menor que la media horquilla: el medio está dentro de la banda de volatilidades compatibles con el bid y el ask (el ejercicio 2 lo comprueba).

La figura de la cabecera dibuja la sonrisa completa del día. Los medios caen sobre la línea de la verdad salvo las cinco calls lejanas con bid nulo, de 6.400 a 7.200, que la celda anterior había dejado fuera. Su medio es la mitad del ask (0,025 puntos) y su volatilidad se separa hasta 5,8 puntos de la verdadera.

Dónde falla: bids nulos, dentro del dinero y cotizaciones congeladas

Un día no basta para ver los fallos raros. Con la cadena entera, 544.548 filas en doce años, los métodos escalares son lentos, pero Newton con salvaguarda de bisección se vectoriza en pocas líneas: cada opción mantiene su intervalo $[\sigma_a, \sigma_b]$ y, si el paso de Newton se sale de él, biseca. El forward y el tipo de cada corte se toman de opciones3_verdad.parquet para aislar el error de las cotizaciones del de la paridad:

x = cad.merge(ver[["fecha", "vencimiento", "dias", "tipo_libre",
                   "forward"]], on=["fecha", "vencimiento"])
Tx = x["dias"].to_numpy() / 365
Dx = np.exp(-x["tipo_libre"].to_numpy() * Tx)
Fx, Kx = x["forward"].to_numpy(), x["strike"].to_numpy()
tx = x["tipo"].to_numpy()
px = ((x["bid"] + x["ask"]) / 2).to_numpy()
intr = Dx * np.maximum(np.where(tx == "C", Fx - Kx, Kx - Fx), 0)
techo = Dx * np.where(tx == "C", Fx, Kx)
valida = (px > intr) & (px < techo)

lo, hi = np.full(len(x), 1e-4), np.full(len(x), 5.0)
s = np.full(len(x), 0.3)
for _ in range(60):
    d = black(Fx, Kx, Tx, Dx, s, tx) - px
    lo, hi = np.where(d < 0, s, lo), np.where(d > 0, s, hi)
    paso = s - d / np.maximum(vega(Fx, Kx, Tx, Dx, s), 1e-12)
    fuera = (paso <= lo) | (paso >= hi) | ~np.isfinite(paso)
    s = np.where(fuera, (lo + hi) / 2, paso)
x["iv"] = np.where(valida, s, np.nan)
x["error_pts"] = 100 * (x["iv"] - x["iv_modelo"])
x["otm"] = np.where(tx == "C", Kx >= Fx, Kx < Fx)
print(f"{len(x)} filas; sin solución: {(~valida).sum()}"
      f" ({(~valida).mean():.1%})")

544548 filas; sin solución: 37662 (6.9%)

Las 37.662 filas sin solución (6,9 % de la cadena) son opciones dentro del dinero cuyo medio queda por debajo del intrínseco descontado. No es un error del método: con esos precios no existe volatilidad implícita.

inc = pd.read_csv("opciones3_incidencias.csv",
                  parse_dates=["fecha", "vencimiento"])
x = x.merge(inc, on=["fecha", "vencimiento", "tipo", "strike"],
            how="left")
x["clase"] = x["incidencia"]
x.loc[x["clase"].isna() & (x["bid"] <= 0), "clase"] = "bid_nulo"
x.loc[x["clase"].isna() & ~x["otm"], "clase"] = "dentro_del_dinero"
x["clase"] = x["clase"].fillna("fuera_del_dinero")
res = x.groupby("clase").agg(
    filas=("iv", "size"),
    sin_solucion=("iv", lambda z: z.isna().mean()),
    mediana_abs=("error_pts", lambda z: z.abs().median()),
    p95_abs=("error_pts", lambda z: z.abs().quantile(0.95)))
res.sort_values("mediana_abs").round(3)

                       filas  sin_solucion  mediana_abs  p95_abs
clase
fuera_del_dinero      210949         0.000        0.013    0.170
dentro_del_dinero     272261         0.138        0.074   16.265
cotizacion_cruzada        12         0.000        0.127   15.770
paridad_violada            8         0.000        0.494    2.834
cotizacion_congelada    1500         0.000        2.491   27.635
bid_nulo               59818         0.000        4.817   59.862

Fuera del dinero y sin incidencias, el error mediano es de 0,013 puntos y el percentil 95, de 0,17: la volatilidad se recupera casi exacta. Dentro del dinero, el 13,8 % de las filas no tiene solución y el percentil 95 llega a 16 puntos, porque el valor temporal es pequeño frente a la horquilla. Los bids nulos son 59.818 filas, el 11 % de la cadena, con una mediana de 4,8 puntos y un percentil 95 de 60. Las congeladas, con 2,5 puntos de mediana: el índice se movió más de un 1 % y la cotización no. Las cruzadas y las de paridad violada son pocas; la paridad violada se descubre comparando la put con la call del mismo strike (lección 3.4.1).

La figura ordena las clases por su error mediano en escala logarítmica: entre la fila limpia y la de bid nulo hay un factor de 370.

orden = ["vencimiento", "tipo", "strike", "fecha"]
y = x.sort_values(orden).copy()
grupo = y.groupby(["vencimiento", "tipo", "strike"])
igual = ((grupo["bid"].shift() == y["bid"])
         & (grupo["ask"].shift() == y["ask"]))
rend = ind.set_index("fecha")["cierre"].pct_change()
movido = y["fecha"].map(rend).abs() > 0.01
y["sospechosa"] = igual & movido & (y["bid"] > 0)
real = y["clase"] == "cotizacion_congelada"
marca = y["sospechosa"]
print(f"Congeladas: {real.sum()}, con bid > 0: "
      f"{(real & (y['bid'] > 0)).sum()}")
print(f"Marcadas: {marca.sum()}; aciertos: {(marca & real).sum()}")
err = y["error_pts"].abs()
print(f"Mediana del error, aciertos: {err[marca & real].median():.2f}"
      f" puntos; falsas alarmas: {err[marca & ~real].median():.2f}")

Congeladas: 1500, con bid > 0: 1090
Marcadas: 1923; aciertos: 1081
Mediana del error, aciertos: 2.07 puntos; falsas alarmas: 0.16

La regla (misma cotización que la sesión anterior, índice movido más de un 1 % y bid positivo) marca 1.923 filas y acierta 1.081 de las 1.090 congeladas con bid positivo; las otras 410 congeladas tienen bid nulo y ya las quita el filtro anterior. Las 842 falsas alarmas son opciones baratas cuyo precio, redondeado al tick, no tenía por qué cambiar: su error mediano es de 0,16 puntos, así que descartarlas cuesta poco. El de los aciertos es de 2,07.

flowchart TB
  A[Fila de la cadena] --> B{Bid mayor que 0}
  B -- no --> X[Descartar]
  B -- sí --> C{Fuera del dinero}
  C -- no --> P[Usar la otra pata]
  C -- sí --> D{Cotización igual}
  D -- sí, índice movido --> X
  D -- no --> E[Calcular sigma]

El diagrama resume los filtros en el orden en que conviene aplicarlos. «Usar la otra pata» quiere decir tomar la put del mismo strike en vez de la call dentro del dinero (o al revés): por la paridad tienen la misma volatilidad.

Ejercicios

1. Programa la bisección pura sobre $[10^{-4}, 5]$ con tolerancia $10^{-10}$ en el precio y cuenta las iteraciones que necesita para la put de 4.000 y la call de 6.000 del día. Compáralas con Newton desde Manaster-Koehler y con Brent.

def biseccion(precio, K, tipo, a=1e-4, b=5.0, tol=1e-10):
    n = 0
    while True:
        m = (a + b) / 2
        dif = float(black(F, K, T, D, m, tipo)) - precio
        n += 1
        if abs(dif) < tol or n == 200:
            return m, n
        a, b = (m, b) if dif < 0 else (a, m)


for K, tipo in [(4000.0, "P"), (6000.0, "C")]:
    s, n = biseccion(medio_de(K, tipo), K, tipo)
    print(f"{K:.0f} {tipo}: σ = {s:.6f} en {n} iteraciones")

4000 P: σ = 0.412856 en 41 iteraciones
6000 C: σ = 0.170921 en 39 iteraciones

2. Calcula la volatilidad implícita del bid y del ask de cada opción fuera del dinero del día (con bid positivo) y mide en qué fracción la iv_modelo queda dentro de esa banda. ¿Dónde se abre más la banda?

banda = []
for f in otm.itertuples():
    lado = [iv_brent(p, f.strike, f.tipo)[0] for p in (f.bid, f.ask)]
    banda.append(lado)
otm[["iv_bid", "iv_ask"]] = banda
dentro = otm["iv_modelo"].between(otm["iv_bid"], otm["iv_ask"])
otm["ancho_pts"] = 100 * (otm["iv_ask"] - otm["iv_bid"])
print(f"iv_modelo dentro de la banda: {dentro.mean():.1%}")
print(otm.groupby(pd.cut(otm["k"], [-0.5, -0.1, 0.1, 0.25]),
                  observed=True)["ancho_pts"].median().round(2))

iv_modelo dentro de la banda: 100.0%
k
(-0.5, -0.1]    0.21
(-0.1, 0.1]     0.18
(0.1, 0.25]     0.25
Name: ancho_pts, dtype: float64

3. Una call de strike 2.400 del día cotiza bid 2.730,40 y ask 2.756,10. Explica por qué no se debe calcular su volatilidad implícita directamente y cómo obtendrías la volatilidad de ese strike.