Lección 3.4.2 · 75 min

Árboles binomiales, Black-Scholes-Merton, Black-76 y Monte Carlo

Un árbol de un paso ya contiene toda la teoría: replicar la opción con el subyacente y un bono. De ahí salen el árbol de Cox-Ross-Rubinstein, la fórmula de Black-Scholes-Merton como su límite, Black-76 para opciones sobre futuros y Monte Carlo para lo que no tiene fórmula. Cada método se programa con NumPy y se contrasta con QuantLib, py_vollib y los precios de modelo de los datos.

Fig. — Convergencia del árbol binomial a Black-Scholes-Merton: call de ORX de strike 5.150 a 30 días, 18 de junio de 2025

02/07

Módulo: 3.4 · Teoría y valoración de opciones
Duración: 75 min
Requisitos: 3.4.1 · Payoffs, paridad put-call y límites de arbitraje, 0.7.2 · Probabilidad y distribuciones, correlación frente a causalidad y regresión simple, 1.6.1 · Rentabilidades, volatilidad, correlaciones móviles y drawdown

Lecciones del nivel 3

Convenciones antes de valorar nada

La lección 3.4.1 acotó los precios de las opciones sin ningún modelo. Para pasar de cotas a un número hace falta suponer cómo se mueve el subyacente. Antes de las fórmulas, las convenciones que se usan en todo el módulo, porque la mitad de los errores de valoración en la práctica son de unidades, no de modelo:

Magnitud	Convención del curso	Error típico
Plazo $T$	días naturales hasta el vencimiento / 365 (vencen al cierre)	días hábiles / 252 en la fórmula de precios
Tipo $r$	continuo, de la curva interpolada en $T$	tipo simple o del plazo equivocado
Dividendos	índice: rendimiento continuo $q$ (o el forward); acción: importes discretos en sus fechas ex	$q$ continuo para una acción con un dividendo la semana que viene
Volatilidad $\sigma$	anual, en tanto por uno (0,217 = 21,7 %)	pasar 21,7 en lugar de 0,217
Precio	por unidad del subyacente (puntos o USD por acción)	olvidar el multiplicador (100 en ORX y KRN; 1.000 barriles en el crudo)

Un paso del árbol: réplica y probabilidad neutral al riesgo

Un subyacente vale $S_0 = 100$ y dentro de un año valdrá $uS_0 = 110$ o $dS_0 = 90$ . Una call de strike 100 paga 10 o 0. Con un tipo continuo del 5 %, se puede replicar la call con $\Delta$ unidades del subyacente y un bono $B$ : basta pedir que la cartera pague lo mismo que la call en los dos estados.

Fuente Mermaid

flowchart LR
  S0["S = 100"] -- sube --> SU["S = 110, C = 10"]
  S0 -- baja --> SD["S = 90, C = 0"]

\Delta = \frac{C_u - C_d}{uS_0 - dS_0}, \qquad B = e^{-rT}\,(C_u - \Delta\,uS_0), \qquad C_0 = \Delta S_0 + B.

IN [01]python

import math
import numpy as np
import pandas as pd

S0, u, d, r, T, K = 100.0, 1.10, 0.90, 0.05, 1.0, 100.0
Cu, Cd = max(u * S0 - K, 0), max(d * S0 - K, 0)
delta = (Cu - Cd) / (u * S0 - d * S0)
B = math.exp(-r * T) * (Cu - delta * u * S0)
print(f"Réplica: {delta:.2f} unidades y bono {B:.2f} "
      f"-> C = {delta * S0 + B:.4f}")

p = (math.exp(r * T) - d) / (u - d)
print(f"p neutral al riesgo = {p:.4f} -> C = "
      f"{math.exp(-r * T) * (p * Cu + (1 - p) * Cd):.4f}")

OUT [01]Un paso: réplica frente a esperanza neutral al riesgo

Réplica: 0.50 unidades y bono -42.81 -> C = 7.1947
p neutral al riesgo = 0.7564 -> C = 7.1947

La cartera de media unidad del subyacente financiada con un préstamo de 42,81 paga exactamente lo que la call y cuesta 7,19: ese es el precio, porque cualquier otro permitiría comprar uno y vender el otro sin riesgo. La misma cifra sale como esperanza descontada con una probabilidad $p = (e^{rT} - d)/(u - d)$ = 0,756. Ese $p$ no es la probabilidad de que el subyacente suba: es el peso que hace que el subyacente, descontado al tipo libre de riesgo, no gane ni pierda de media. La probabilidad real de subida no aparece en ningún sitio, y por eso dos inversores que discrepan sobre ella están de acuerdo en el precio de la opción.

El árbol de Cox-Ross-Rubinstein

Repetir el paso $n$ veces con $u = e^{\sigma\sqrt{\Delta t}}$ , $d = 1/u$ y $\Delta t = T/n$ da el árbol de Cox, Ross y Rubinstein (1979). Con un rendimiento por dividendo continuo $q$ , la probabilidad es $p = (e^{(r-q)\Delta t} - d)/(u - d)$ . Se valora hacia atrás: en el último paso, el payoff; en cada nodo anterior, la esperanza descontada de sus dos hijos y, si la opción es americana, el máximo entre ese valor y ejercer ya. Con NumPy basta un vector por paso:

IN [02]python

def crr(S, K, T, r, q, sig, n, call=True, americana=False):
    dt = T / n
    u = math.exp(sig * math.sqrt(dt))
    p = (math.exp((r - q) * dt) - 1 / u) / (u - 1 / u)
    desc = math.exp(-r * dt)
    st = S * u ** (n - 2 * np.arange(n + 1))  # nodos finales
    v = np.maximum(st - K, 0) if call else np.maximum(K - st, 0)
    for _ in range(n):
        v = desc * (p * v[:-1] + (1 - p) * v[1:])
        if americana:
            st = st[:-1] / u
            v = np.maximum(v, st - K if call else K - st)
    return v[0]

print(f"{crr(100, 100, 1, 0.05, 0, 0.2, 500):.4f}",
      f"{crr(100, 100, 1, 0.05, 0, 0.2, 500, call=False):.4f}",
      f"{crr(100, 100, 1, 0.05, 0, 0.2, 500, False, True):.4f}")

OUT [02]Árbol CRR europeo y americano

10.4466 5.5695 6.0888

Con 500 pasos, la call europea del ejemplo de libro ( $S = K = 100$ , un año, 5 %, 20 %) vale 10,4466, la put europea 5,5695 y la put americana 6,0888: medio punto más, el valor de poder cobrar el strike antes.

Black-Scholes-Merton como límite del árbol

Cuando $n \to \infty$ , el árbol converge a la fórmula de Black-Scholes-Merton (BSM), la de un subyacente cuyo logaritmo sigue un paseo aleatorio gaussiano con volatilidad constante:

C = S e^{-qT} N(d_1) - K e^{-rT} N(d_2), \qquad P = K e^{-rT} N(-d_2) - S e^{-qT} N(-d_1),

d_1 = \frac{\ln(S/K) + (r - q)T}{\sigma\sqrt{T}} + \frac{\sigma\sqrt{T}}{2}, \qquad d_2 = d_1 - \sigma\sqrt{T}.

Se aplica a una opción de los datos: la call de ORX de strike 5.150 (el más cercano al índice) del vencimiento del 18 de julio de 2025, el 18 de junio. El índice, el rendimiento por dividendo y la curva salen de opciones3_indice.csv; la volatilidad, de la columna iv_modelo, la volatilidad implícita del modelo con la que se generó la cadena (la lección 3.4.4 la recupera de los precios):

IN [03]python

from scipy.stats import norm

def bsm(S, K, T, r, q, sig, call=True):
    sd = sig * np.sqrt(T)
    d1 = (np.log(S / K) + (r - q) * T) / sd + sd / 2
    d2 = d1 - sd
    if call:
        return (S * np.exp(-q * T) * norm.cdf(d1)
                - K * np.exp(-r * T) * norm.cdf(d2))
    return (K * np.exp(-r * T) * norm.cdf(-d2)
            - S * np.exp(-q * T) * norm.cdf(-d1))

cad = pd.read_parquet("opciones3_cadena.parquet")
ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"],
                  index_col="fecha")
plazos = np.array([1 / 12, 0.25, 0.5, 1.0, 2.0])
cols = ["tipo_1m", "tipo_3m", "tipo_6m", "tipo_1a", "tipo_2a"]
hoy, venc = pd.Timestamp("2025-06-18"), pd.Timestamp("2025-07-18")
T = (venc - hoy).days / 365
r = float(np.interp(T, plazos, ind.loc[hoy, cols].to_numpy(float)
                    / 100))
q = float(ind.loc[hoy, "rentabilidad_dividendo"]) / 100
S = float(ind.loc[hoy, "cierre"])
corte = cad[(cad["fecha"] == hoy) & (cad["vencimiento"] == venc)]
K = 5150.0
op = corte[(corte["strike"] == K) & (corte["tipo"] == "C")].iloc[0]
sig = op["iv_modelo"]
c_bsm = bsm(S, K, T, r, q, sig)
print(f"S = {S}, T = {T:.4f} años ({(venc - hoy).days} días), "
      f"r = {r:.4%}, q = {q:.4%}, sigma = {sig}")
print(f"BSM: {c_bsm:.4f} | bid {op['bid']} / ask {op['ask']}"
      f" | por contrato: {100 * c_bsm:,.0f} USD")

OUT [03]BSM sobre una call de ORX

S = 5148.06, T = 0.0822 años (30 días), r = 4.1910%, q = 1.3438%, sigma = 0.217
BSM: 132.5896 | bid 131.8 / ask 133.1 | por contrato: 13,259 USD

La call vale 132,59 puntos, dentro de su horquilla de 131,80 a 133,10, y 13.259 USD por contrato. La figura de la cabecera muestra cómo llega a ese número el árbol CRR cuando crece el número de pasos, junto al árbol de Leisen-Reimer de QuantLib, que coloca los nodos para que el strike quede en el centro y converge mucho antes. La celda compara los dos árboles propios y los de QuantLib con BSM:

IN [04]python

import QuantLib as ql

hoy_ql = ql.Date(hoy.day, hoy.month, hoy.year)
ql.Settings.instance().evaluationDate = hoy_ql
dc = ql.Actual365Fixed()

def curva_plana(x):
    return ql.YieldTermStructureHandle(
        ql.FlatForward(hoy_ql, x, dc, ql.Continuous))

proceso = ql.BlackScholesMertonProcess(
    ql.QuoteHandle(ql.SimpleQuote(S)), curva_plana(q),
    curva_plana(r), ql.BlackVolTermStructureHandle(
        ql.BlackConstantVol(hoy_ql, ql.NullCalendar(), sig, dc)))
ejercicio = ql.EuropeanExercise(
    ql.Date(venc.day, venc.month, venc.year))
opcion = ql.VanillaOption(
    ql.PlainVanillaPayoff(ql.Option.Call, K), ejercicio)

filas = []
for n in [25, 50, 51, 100, 101, 500]:
    opcion.setPricingEngine(ql.BinomialVanillaEngine(proceso,
                                                     "crr", n))
    ql_crr = opcion.NPV()
    opcion.setPricingEngine(ql.BinomialVanillaEngine(proceso,
                                                     "lr", n))
    filas.append({"n": n,
                  "CRR propio": crr(S, K, T, r, q, sig, n) - c_bsm,
                  "CRR QuantLib": ql_crr - c_bsm,
                  "LR QuantLib": opcion.NPV() - c_bsm})
pd.DataFrame(filas).set_index("n").round(4)

OUT [04]Convergencia del árbol: propio y QuantLib

     CRR propio  CRR QuantLib  LR QuantLib
n
25       1.2747        1.2745      -0.0098
50      -0.5299       -0.5300      -2.6206
51       0.6222        0.6221      -0.0024
100     -0.2439       -0.2439      -1.3231
101      0.3127        0.3126      -0.0006
500     -0.0316       -0.0316      -0.2668

Tres lecciones de la tabla:

El error de CRR oscila: con 50 pasos el árbol infravalora 0,53 puntos y con 51 sobrevalora 0,62. El signo depende de si el strike cae cerca de un nodo final o entre dos. Promediar dos árboles consecutivos o usar muchos pasos corrige casi todo; con 500 pasos el error es de 0,03.
El CRR de QuantLib no es idéntico al propio, aunque la diferencia es de diezmilésimas y se reduce con $n$ : QuantLib no construye el árbol con las mismas fórmulas exactas. Los dos convergen al mismo límite.
El árbol de Leisen-Reimer de QuantLib converge con muy pocos pasos si n es impar (0,0024 con 51); con n par, en QuantLib 1.44, da errores de hasta 2,6 puntos con 50 pasos y de 0,27 con 500. Antes de fiarse de un motor de una biblioteca hay que comprobarlo con algo conocido, que es lo que hace esta celda.

Valor antes del vencimiento

BSM da también el valor de la opción en cualquier momento anterior al vencimiento, y con él el resultado de una posición antes de que expire. La figura dibuja la call de la sección anterior en función del índice a 30, 10 y 1 días del vencimiento, con la misma volatilidad, junto a su payoff: lejos del vencimiento la curva es suave y está por encima del payoff (valor temporal); al acercarse, se pega a él, y la transición se concentra en el strike. Las pendientes y curvaturas de estas curvas son las griegas de la lección 3.4.3.

Fig. — Valor de una call de ORX de strike 5.150 a 30, 10 y 1 días del vencimiento, con Black-Scholes-Merton (volatilidad del 21,7 %), frente a su payoff

py_vollib y Black-76

py_vollib 1.0.12 es una implementación de las fórmulas de Black, Black-Scholes y Black-Scholes-Merton basada en el algoritmo «Let's Be Rational» de Peter Jäckel. En esta versión el código vive en el paquete vollib (1.0.11), que se instala con ella: import py_vollib sigue funcionando, pero emite un DeprecationWarning que pide importar vollib. Aquí se importa directamente vollib, que es la misma biblioteca, para no ensuciar la salida.

Las opciones sobre futuros se valoran con Black-76 (Black, 1976): la misma fórmula escrita sobre el precio del futuro $F$ , que no tiene coste de mantenimiento (comprarlo no cuesta nada hoy):

C = e^{-rT}\left[F N(d_1) - K N(d_2)\right], \qquad d_1 = \frac{\ln(F/K)}{\sigma\sqrt{T}} + \frac{\sigma\sqrt{T}}{2}.

Si $F = S e^{(r-q)T}$ es el forward del mismo vencimiento, Black-76 y BSM son la misma fórmula. Para un índice eso significa que se puede valorar sobre el forward implícito de la cadena (lección 3.4.1) sin conocer el dividendo:

IN [05]python

from vollib.black import black
from vollib.black_scholes_merton import black_scholes_merton

def black76(F, K, T, r, sig, call=True):
    sd = sig * np.sqrt(T)
    d1 = np.log(F / K) / sd + sd / 2
    d2 = d1 - sd
    if call:
        return np.exp(-r * T) * (F * norm.cdf(d1) - K * norm.cdf(d2))
    return np.exp(-r * T) * (K * norm.cdf(-d2) - F * norm.cdf(-d1))

F = S * math.exp((r - q) * T)
opcion.setPricingEngine(ql.AnalyticEuropeanEngine(proceso))
calc = ql.BlackCalculator(ql.PlainVanillaPayoff(ql.Option.Call, K),
                          F, sig * math.sqrt(T), math.exp(-r * T))
pd.Series({
    "BSM propio": c_bsm,
    "Black-76 propio": black76(F, K, T, r, sig),
    "vollib BSM": black_scholes_merton("c", S, K, T, r, sig, q),
    "vollib Black": black("c", F, K, T, r, sig),
    "QuantLib BSM": opcion.NPV(),
    "QuantLib Black": calc.value(),
}).map("{:.10f}".format)

OUT [05]BSM, Black-76, py_vollib y QuantLib

BSM propio         132.5896128016
Black-76 propio    132.5896128016
vollib BSM         132.5896128016
vollib Black       132.5896128016
QuantLib BSM       132.5896128016
QuantLib Black     132.5896128016
dtype: str

Las seis cifras coinciden en diez decimales. Lo que sí cambia el precio es el dato de entrada. La cadena entera de ese día, valorada con Black-76 sobre el forward de opciones3_verdad.parquet y la volatilidad de cada strike, cae dentro de su horquilla:

IN [06]python

verdad = pd.read_parquet("opciones3_verdad.parquet")
dia = (cad[cad["fecha"] == hoy]
       .merge(verdad, on=["fecha", "vencimiento"]))
Tv = dia["dias"] / 365
es_call = (dia["tipo"] == "C").to_numpy()
modelo = np.where(
    es_call,
    black76(dia["forward"], dia["strike"], Tv, dia["tipo_libre"],
            dia["iv_modelo"]),
    black76(dia["forward"], dia["strike"], Tv, dia["tipo_libre"],
            dia["iv_modelo"], call=False))
dentro = (modelo >= dia["bid"]) & (modelo <= dia["ask"])
medio = (dia["bid"] + dia["ask"]) / 2
print(f"{len(dia)} opciones; dentro de la horquilla: "
      f"{dentro.mean():.1%}")
print(f"|modelo - medio|: mediana {np.median(abs(modelo - medio)):.3f}"
      f", máximo {np.max(abs(modelo - medio)):.3f} puntos")

OUT [06]La cadena de un día con Black-76

266 opciones; dentro de la horquilla: 100.0%
|modelo - medio|: mediana 0.046, máximo 2.586 puntos

Una opción sobre un futuro de crudo

Los datos no traen opciones sobre futuros, pero sí la curva completa del crudo ficticio PEC (imita al WTI de NYMEX: 1.000 barriles, en USD por barril). Se valora una opción europea hipotética sobre el contrato de febrero de 2026 (PECG26) el 15 de octubre de 2025, con 60 días de plazo y la volatilidad realizada de las 63 sesiones anteriores del propio contrato. Si la opción fuera americana (lo dice la ficha de cada contrato), Black-76 ya no basta: QuantLib valora la americana con un árbol sobre el futuro (BlackProcess):

IN [07]python

curvas = pd.read_parquet("futuros3_curvas.parquet",
                         filters=[("contrato", "==", "PECG26")])
fut = curvas.set_index("fecha")["liquidacion"]
dia_f = pd.Timestamp("2025-10-15")
F_c = float(fut.loc[dia_f])
rv = np.log(fut).diff().loc[:dia_f].tail(63).std() * math.sqrt(252)
tipos = pd.read_csv("tipos3_curva.csv", parse_dates=["fecha"],
                    index_col="fecha")
r_c = math.log(1 + tipos.loc[dia_f, "usd_1m"] / 100)
T_c = 60 / 365
print(f"F = {F_c}, vol. realizada = {rv:.3f}, r = {r_c:.4f}")

hoy_c = ql.Date(15, 10, 2025)
ql.Settings.instance().evaluationDate = hoy_c
vto_c = hoy_c + 60
proc_f = ql.BlackProcess(
    ql.QuoteHandle(ql.SimpleQuote(F_c)),
    ql.YieldTermStructureHandle(
        ql.FlatForward(hoy_c, r_c, dc, ql.Continuous)),
    ql.BlackVolTermStructureHandle(
        ql.BlackConstantVol(hoy_c, ql.NullCalendar(), rv, dc)))
filas = []
for tipo, Kc in [("C", 60.0), ("C", 66.0), ("P", 66.0), ("P", 75.0)]:
    pago = ql.PlainVanillaPayoff(
        ql.Option.Call if tipo == "C" else ql.Option.Put, Kc)
    eu = black76(F_c, Kc, T_c, r_c, rv, call=tipo == "C")
    vl = black(tipo.lower(), F_c, Kc, T_c, r_c, rv)
    am = ql.VanillaOption(pago, ql.AmericanExercise(hoy_c, vto_c))
    am.setPricingEngine(ql.BinomialVanillaEngine(proc_f, "lr", 801))
    filas.append({"tipo": tipo, "K": Kc, "Black-76": eu,
                  "vollib": vl, "americana QL": am.NPV(),
                  "USD/contrato (eu.)": 1000 * eu})
pd.DataFrame(filas).round(4)

OUT [07]Black-76 sobre el futuro PECG26, europea y americana

F = 66.11, vol. realizada = 0.206, r = 0.0398
  tipo     K  Black-76  vollib  americana QL  USD/contrato (eu.)
0    C  60.0    6.3875  6.3875        6.4024           6387.5328
1    C  66.0    2.2415  2.2415        2.2438           2241.4764
2    P  66.0    2.1322  2.1322        2.1343           2132.1940
3    P  75.0    8.9992  8.9992        9.0261           8999.1781

Una put de strike 66 sobre un futuro a 66,11 vale 2,13 USD por barril, 2.132 USD por contrato de 1.000 barriles. El ejercicio anticipado de una opción sobre futuros vale poco salvo muy dentro del dinero: lo único que se gana ejerciendo antes son los intereses de la prima intrínseca, y en la put de 75, la más dentro del dinero, son 0,027 USD por barril (27 USD por contrato).

Opciones americanas con dividendos discretos

Un índice paga dividendos casi a diario y un rendimiento continuo $q$ lo describe bien. Una acción paga cuatro al año, en fechas conocidas, y eso cambia dos cosas: el precio cae de golpe en la fecha ex, y una call americana puede valer más ejerciéndola la víspera. Los datos de KRN (acción ficticia, opciones americanas, dividendo trimestral) traen la cadena con su precio de modelo, calculado con un árbol CRR de 150 pasos con el modelo «escrowed»: el árbol se construye sobre $S$ menos el valor presente de los dividendos hasta el vencimiento, y en cada nodo se suma el valor presente de los que faltan para comparar con el ejercicio. El objetivo es recuperar esos precios con un árbol propio y con QuantLib.

Fuente Mermaid

flowchart TB
  O[Opción] --> E{¿Europea?}
  E -- sí --> P{¿Payoff simple?}
  P -- sí --> BS[BSM o Black-76]
  P -- no --> MC[Monte Carlo]
  E -- no --> D{¿Dividendos?}
  D -- discretos --> AR[Árbol escrowed]
  D -- continuos --> AR2[Árbol o FD]

IN [08]python

krn = pd.read_parquet("accion3_opciones.parquet")
acc = pd.read_csv("accion3_diario.csv", parse_dates=["fecha"],
                  index_col="fecha")
ev = pd.read_csv("accion3_eventos.csv", parse_dates=["fecha"])
divs = ev[ev["evento"] == "dividendo"].set_index("fecha")[
    "dividendo"]

f_k, v_k = pd.Timestamp("2025-11-03"), pd.Timestamp("2025-11-21")
S_k = float(acc.loc[f_k, "cierre"])
T_k = (v_k - f_k).days / 365
r_k = float(np.interp(T_k, plazos,
                      ind.loc[f_k, cols].to_numpy(float) / 100))
d_k = divs[(divs.index > f_k) & (divs.index <= v_k)]
t_div = ((d_k.index - f_k).days / 365).to_numpy()
print(f"S = {S_k}, T = {(v_k - f_k).days} días, r = {r_k:.4%};",
      "dividendo", d_k.to_dict())

def crr_escrowed(S, K, T, r, sig, t_div, divs, call, n=150,
                 americana=True):
    dt = T / n
    u = math.exp(sig * math.sqrt(dt))
    p = (math.exp(r * dt) - 1 / u) / (u - 1 / u)
    desc = math.exp(-r * dt)
    t = np.arange(n + 1) * dt
    vp = np.zeros(n + 1)       # VP en cada paso de los que faltan
    for td, D in zip(t_div, divs):
        vp += np.where(td > t + 1e-12, D * np.exp(-r * (td - t)), 0)
    st = (S - vp[0]) * u ** (n - 2 * np.arange(n + 1))
    v = np.maximum(st - K, 0) if call else np.maximum(K - st, 0)
    for i in range(n - 1, -1, -1):
        v = desc * (p * v[:-1] + (1 - p) * v[1:])
        st = st[:-1] / u
        if americana:
            s_total = st + vp[i]
            v = np.maximum(v, s_total - K if call else K - s_total)
    return v[0]

c_k = krn[(krn["fecha"] == f_k) & (krn["vencimiento"] == v_k)].copy()
c_k["arbol"] = [crr_escrowed(S_k, k, T_k, r_k, s, t_div, d_k.values,
                             tp == "C")
                for k, s, tp in zip(c_k["strike"], c_k["iv_modelo"],
                                    c_k["tipo"])]
err = (c_k["arbol"] - c_k["precio_modelo"]).abs()
print(f"{len(c_k)} opciones; error máximo frente al dato: "
      f"{err.max():.4f} USD")

OUT [08]KRN: el árbol con dividendos discretos

S = 184.27, T = 18 días, r = 3.8990%; dividendo {Timestamp('2025-11-13 00:00:00'): 1.04}
50 opciones; error máximo frente al dato: 0.0011 USD

El árbol propio reproduce los 50 precios de modelo de ese corte con un error máximo de 0,0011 USD, el que dejan los redondeos del cierre (al céntimo), de la curva y de la volatilidad (a cuatro decimales). QuantLib resuelve el mismo problema con diferencias finitas sobre la ecuación de Black-Scholes, con los dividendos como importes en sus fechas y el mismo modelo «escrowed» (FdBlackScholesVanillaEngine.Escrowed):

IN [09]python

hoy_k = ql.Date(f_k.day, f_k.month, f_k.year)
ql.Settings.instance().evaluationDate = hoy_k
fechas_div = [hoy_k + int(round(x * 365)) for x in t_div]
dividendos = ql.DividendVector(fechas_div, list(d_k.values))
spot_k = ql.QuoteHandle(ql.SimpleQuote(S_k))
r_ts = ql.YieldTermStructureHandle(
    ql.FlatForward(hoy_k, r_k, dc, ql.Continuous))
q_ts = ql.YieldTermStructureHandle(
    ql.FlatForward(hoy_k, 0.0, dc, ql.Continuous))
vto_k = ql.Date(v_k.day, v_k.month, v_k.year)

def ql_americana(K, sig, call):
    vol = ql.BlackVolTermStructureHandle(
        ql.BlackConstantVol(hoy_k, ql.NullCalendar(), sig, dc))
    proc = ql.BlackScholesMertonProcess(spot_k, q_ts, r_ts, vol)
    motor = ql.FdBlackScholesVanillaEngine(
        proc, dividendos, 400, 400, 0, ql.FdmSchemeDesc.Douglas(),
        False, -ql.nullDouble(),
        ql.FdBlackScholesVanillaEngine.Escrowed)
    o = ql.VanillaOption(
        ql.PlainVanillaPayoff(ql.Option.Call if call
                              else ql.Option.Put, K),
        ql.AmericanExercise(hoy_k, vto_k))
    o.setPricingEngine(motor)
    return o.NPV()

c_k["quantlib"] = [ql_americana(k, s, tp == "C")
                   for k, s, tp in zip(c_k["strike"],
                                       c_k["iv_modelo"], c_k["tipo"])]
c_k["prima_ej"] = c_k["precio_modelo"] - c_k["europea_modelo"]
print("QuantLib FD frente al dato: error máximo "
      f"{(c_k['quantlib'] - c_k['precio_modelo']).abs().max():.4f}")
tabla = c_k[c_k["strike"].isin([165.0, 175.0, 185.0, 195.0, 205.0])]
tabla.pivot(index="strike", columns="tipo",
            values=["precio_modelo", "europea_modelo", "prima_ej"]
            ).round(4)

OUT [09]KRN con diferencias finitas de QuantLib

QuantLib FD frente al dato: error máximo 0.0076
       precio_modelo          europea_modelo          prima_ej
tipo               C        P              C        P        C       P
strike
165.0        19.6928   0.4962        19.0424   0.4950   0.6504  0.0012
175.0        10.5096   1.5681        10.1290   1.5625   0.3806  0.0056
185.0         3.4471   4.7965         3.3529   4.7671   0.0942  0.0294
195.0         0.5779  12.0636         0.5740  11.9690   0.0039  0.0946
205.0         0.0783  21.6055         0.0783  21.4541   0.0000  0.1514

QuantLib, con una rejilla de 400 × 400, queda a menos de 0,008 USD de los precios de los datos; las diferencias son de discretización (de su rejilla y del árbol de 150 pasos), muy por debajo de un tick de 0,01. La tabla y la figura muestran el valor del ejercicio anticipado:

Fig. — Valor del ejercicio anticipado de las opciones americanas de KRN por strike, 3 de noviembre de 2025, vencimiento del 21 con el dividendo ex el día 13

En las calls dentro del dinero el ejercicio anticipado vale buena parte del dividendo: con el ex a 10 días, la call de 165 vale 0,65 USD más que la europea (el 63 % del dividendo de 1,04) y la de 175, 0,38. Quien la tenga y no la ejerza la víspera del ex regala ese importe; quien la haya vendido debe contar con la asignación.
En las puts el premio es pequeño a 18 días: los intereses de cobrar antes el strike son de céntimos y además el dividendo, que hará caer la acción, aconseja esperar.

Valorar estas calls con BSM europeo y un rendimiento continuo equivalente infravalora justo las opciones que más se negocian alrededor del ex. Es el error de convenciones de la tabla inicial con nombre propio.

Monte Carlo

Cuando no hay fórmula ni árbol cómodo (opciones que dependen de la trayectoria, varios subyacentes, modelos con volatilidad estocástica), se simula. Bajo la medida neutral al riesgo, el índice al vencimiento es $S_T = F\,e^{-\sigma^2 T/2 + \sigma\sqrt{T} Z}$ con $Z$ normal estándar, y el precio es la media de los payoffs descontados (Boyle, 1977). La simulación es exacta en $T$ : no hace falta discretizar el camino. Lo que no es exacto es la media, que tiene un error típico $s/\sqrt{N}$ . Dos técnicas lo reducen sin más trayectorias:

Variables antitéticas: cada $Z$ se usa también con signo contrario; los dos payoffs están correlados negativamente y su media varía menos.
Variable de control: se simula a la vez algo cuyo valor exacto se conoce (aquí, el propio $e^{-rT}S_T$ , que vale $e^{-rT}F$ ) y se corrige la estimación con su error: $\hat C_{cv} = \hat C - \beta(\hat X - E[X])$ , con $\beta$ la pendiente de la regresión del payoff sobre el control.

IN [10]python

rng = np.random.default_rng(2029)
DF = math.exp(-r * T)

def mc_call(N, z):
    ST = F * np.exp(-0.5 * sig**2 * T + sig * math.sqrt(T) * z)
    return DF * np.maximum(ST - K, 0), DF * ST

filas = []
for N in [1_000, 10_000, 100_000, 1_000_000]:
    z = rng.standard_normal(N)
    pago, ctrl = mc_call(N, z)
    pago_a, _ = mc_call(N, -z)
    anti = (pago + pago_a) / 2                # N/2 pares bastarían
    beta = np.cov(pago, ctrl)[0, 1] / ctrl.var(ddof=1)
    cv = pago - beta * (ctrl - DF * F)
    for nombre, x in [("simple", pago), ("antitética", anti),
                      ("control", cv)]:
        filas.append({"N": N, "estimador": nombre,
                      "error": x.mean() - c_bsm,
                      "error_tipico": x.std(ddof=1) / math.sqrt(N)})
mc = pd.DataFrame(filas)
mc.pivot(index="N", columns="estimador",
         values="error_tipico").round(4)

OUT [10]Monte Carlo de la call de ORX: tres estimadores

estimador  antitética  control  simple
N
1000           3.2873   3.0652  6.4429
10000          1.0361   0.9643  1.9480
100000         0.3268   0.3046  0.6244
1000000        0.1030   0.0961  0.1972

El error típico cae como $1/\sqrt{N}$ en los tres estimadores: multiplicar las trayectorias por 100 lo divide por 10. La variable de control lo divide además por algo más de dos con cualquier $N$ (de 0,62 a 0,30 con 100.000), porque el payoff de una call en el dinero está muy correlado con el índice; las antitéticas, por 1,9 (y cuestan el doble de evaluaciones del payoff). La figura lo muestra en escala logarítmica. Los errores realizados frente a BSM son del orden de esos errores típicos:

Fig. — Error típico de Monte Carlo de la call de ORX según el número de trayectorias, con tres estimadores

IN [11]python

print(mc.pivot(index="N", columns="estimador",
               values="error").round(3))
ql.Settings.instance().evaluationDate = hoy_ql
opcion.setPricingEngine(ql.MCEuropeanEngine(
    proceso, "pseudorandom", timeSteps=1, requiredSamples=200_000,
    seed=42, antitheticVariate=True))
print(f"QuantLib MC: {opcion.NPV():.3f} ± "
      f"{opcion.errorEstimate():.3f} (BSM {c_bsm:.3f})")

OUT [11]Errores realizados y QuantLib

estimador  antitética  control  simple
N
1000            2.334    2.494   0.004
10000           0.535    0.478  -2.391
100000          0.233    0.200   0.572
1000000         0.064    0.063   0.171
QuantLib MC: 132.368 ± 0.230 (BSM 132.590)

QuantLib, con 200.000 trayectorias antitéticas y su propio generador, da 132,37 ± 0,23: a menos de un error típico de BSM.

Una opción asiática: donde Monte Carlo es la herramienta

Una call asiática paga $\max(\overline{S} - K, 0)$ , con $\overline{S}$ la media aritmética de los cierres de un periodo. No tiene fórmula cerrada, pero su hermana con media geométrica sí (el logaritmo de una media geométrica de lognormales es normal). Es la variable de control ideal: casi idéntica a la aritmética y con valor exacto, que QuantLib calcula analíticamente. Se valora una call asiática sobre ORX con la media de las 20 sesiones de la NYSE que quedan hasta el vencimiento de julio:

IN [12]python

import exchange_calendars as xcals

xnys = xcals.get_calendar("XNYS")
ses = xnys.sessions_in_range(hoy + pd.Timedelta(days=1), venc)
t_fix = ((ses - hoy).days / 365).to_numpy()
print(len(t_fix), "fijaciones, de", ses[0].date(), "a",
      ses[-1].date())

N = 100_000
dt = np.diff(np.r_[0, t_fix])
z = rng.standard_normal((N, len(t_fix)))
logS = np.log(S) + np.cumsum((r - q - sig**2 / 2) * dt
                             + sig * np.sqrt(dt) * z, axis=1)
arit = DF * np.maximum(np.exp(logS).mean(axis=1) - K, 0)
geom = DF * np.maximum(np.exp(logS.mean(axis=1)) - K, 0)

fix_ql = [ql.Date(x.day, x.month, x.year) for x in ses]
asi = ql.DiscreteAveragingAsianOption(
    ql.Average.Geometric, fix_ql,
    ql.PlainVanillaPayoff(ql.Option.Call, K), ejercicio)
asi.setPricingEngine(
    ql.AnalyticDiscreteGeometricAveragePriceAsianEngine(proceso))
geo_exacta = asi.NPV()

beta = np.cov(arit, geom)[0, 1] / geom.var(ddof=1)
cv = arit - beta * (geom - geo_exacta)
asi_a = ql.DiscreteAveragingAsianOption(
    ql.Average.Arithmetic, fix_ql,
    ql.PlainVanillaPayoff(ql.Option.Call, K), ejercicio)
asi_a.setPricingEngine(ql.MCDiscreteArithmeticAPEngine(
    proceso, "pseudorandom", requiredSamples=100_000, seed=42,
    controlVariate=True))
pd.DataFrame({
    "precio": [geom.mean(), geo_exacta, arit.mean(), cv.mean(),
               asi_a.NPV()],
    "error_tipico": [geom.std() / N**0.5, 0, arit.std() / N**0.5,
                     cv.std() / N**0.5, asi_a.errorEstimate()]},
    index=["geométrica MC", "geométrica exacta (QL)",
           "aritmética MC simple", "aritmética MC control",
           "aritmética QL MC control"]).round(4)

OUT [12]Asiática aritmética con control geométrico

20 fijaciones, de 2025-06-20 a 2025-07-18
                           precio  error_tipico
geométrica MC             81.8746        0.3860
geométrica exacta (QL)    81.6093        0.0000
aritmética MC simple      82.7353        0.3893
aritmética MC control     82.4677        0.0032
aritmética QL MC control  82.4743        0.0046

La geométrica simulada coincide con la exacta dentro de su error típico, lo que valida la simulación. Con ella como control, el error típico de la aritmética baja de 0,39 a 0,003 puntos: unas 120 veces menos, el equivalente a multiplicar las trayectorias por casi 15.000. QuantLib, con el mismo control, da el mismo precio. La asiática vale bastante menos que la call europea (82,47 frente a 132,59) porque la media de un mes varía menos que el cierre del último día.

Errores típicos

Unidades. Días hábiles / 252 en el plazo, volatilidad en porcentaje, un tipo simple como si fuera continuo o el multiplicador olvidado. Cada uno, por separado, mueve el precio más que la diferencia entre dos modelos.
Contado y dividendo supuesto en un índice. El forward implícito de la cadena ya descuenta el dividendo que el mercado espera; con Black-76 sobre ese forward se evita adivinarlo.
Fórmulas europeas para americanas. Con dividendos discretos o puts muy dentro del dinero, BSM infravalora; hay que usar un árbol o diferencias finitas.
Árboles con pocos pasos. CRR oscila con la paridad de $n$ ; con una biblioteca, se comprueba el motor contra un caso conocido (como el LR con $n$ par en la celda de convergencia).
Monte Carlo sin error típico. Una estimación sin su error no se puede comparar con nada. Se fija la semilla para que sea reproducible y se informa siempre el error.
Volatilidad. Todo lo anterior supone una volatilidad dada. Cuál usar (implícita por strike, realizada, de una superficie) es el tema de las lecciones 3.4.4 a 3.4.7; aquí se usó la del modelo de los datos para aislar el método.

Ejercicios

Valora con crr la put americana de ORX del ejemplo (strike 5.150, vencimiento de julio, el 18 de junio de 2025) como si fuera americana, con 100, 500 y 1.000 pasos. ¿Cuánto vale el ejercicio anticipado frente a la europea?

Solución

IN [13]python

put_eu = bsm(S, K, T, r, q, sig, call=False)
for n in [100, 500, 1000]:
    am = crr(S, K, T, r, q, sig, n, call=False, americana=True)
    print(n, f"americana {am:.4f} | europea BSM {put_eu:.4f} | "
          f"prima {am - put_eu:.4f}")

OUT [13]Put americana sobre ORX

100 americana 123.1397 | europea BSM 122.5030 | prima 0.6367
500 americana 123.2939 | europea BSM 122.5030 | prima 0.7908
1000 americana 123.3085 | europea BSM 122.5030 | prima 0.8055

El ejercicio anticipado de una put en el dinero a 30 días vale 0,81 puntos con 1.000 pasos: los intereses de cobrar el strike unos días antes, menos el dividendo del índice que se deja de ganar. Es poco porque la put en el dinero todavía tiene mucho valor temporal. Las opciones sobre índices como SPX son europeas, así que esta prima no existe en ellas; en las opciones americanas sobre ETF del mismo índice sí.

Estima por Monte Carlo una call digital de ORX (paga 100 puntos si el índice acaba por encima de 5.150) con 100.000 trayectorias y compárala con la fórmula $100\,e^{-rT}N(d_2)$ .

Solución

IN [14]python

z = rng.standard_normal(100_000)
ST = F * np.exp(-0.5 * sig**2 * T + sig * math.sqrt(T) * z)
dig = DF * 100 * (ST > K)
sd = sig * math.sqrt(T)
d2 = (math.log(S / K) + (r - q) * T) / sd - sd / 2
print(f"MC {dig.mean():.3f} ± {dig.std() / 1e5**0.5:.3f} | "
      f"fórmula {100 * DF * norm.cdf(d2):.3f}")

OUT [14]Digital por Monte Carlo

MC 49.688 ± 0.158 | fórmula 49.846

La digital vale lo que la probabilidad neutral al riesgo de acabar dentro, descontada y por 100. Monte Carlo la recupera dentro de su error típico. Su error relativo es mayor que el de la call porque el payoff es discontinuo; una call digital cotizada se replica en la práctica con un call spread estrecho (lección 3.4.1).

Un compañero valora las calls de KRN del 3 de noviembre de 2025 con BSM y un rendimiento continuo $q = 4 \times 1{,}04 / S$ . Explica qué dos errores comete y en qué strikes se notan más.

Solución

El primero es repartir en todo el año un dividendo que cae en un día concreto, dentro de 10 días: el forward a 18 días lleva el dividendo entero (1,04 USD), no 18/365 de cuatro dividendos, así que el rendimiento continuo infravalora la caída que descuenta el mercado y sobrevalora las calls europeas. El segundo es tratar como europeas unas opciones americanas: ignora el ejercicio anticipado la víspera del ex, que en las calls dentro del dinero vale casi todo el dividendo (0,65 USD en la de 165). Los dos errores tienen signo contrario y no se compensan por igual en todos los strikes: el primero afecta a toda la cadena; el segundo, sobre todo a las calls dentro del dinero.

Respuesta plegada en cada pregunta

Autoevaluación

En un árbol de un paso con u = 1,10, d = 0,90 y un tipo del 5 %, la probabilidad neutral al riesgo de subida es 0,756. ¿Qué significa?
- AQue el mercado espera que el subyacente suba con un 75,6 % de probabilidad
- BQue es el peso que hace que el subyacente, descontado, sea una martingala; el precio de la opción sale de la réplica y no depende de la probabilidad real de subida
- CQue la opción vale el 75,6 % de su valor intrínseco
- DQue hay arbitraje, porque debería ser 0,5
Ver respuesta
Respuesta B. Que es el peso que hace que el subyacente, descontado, sea una martingala; el precio de la opción sale de la réplica y no depende de la probabilidad real de subida
El precio se fija por réplica: la cartera de 0,5 unidades del subyacente y un préstamo de 42,81 paga lo mismo que la call en los dos estados y cuesta 7,19. La probabilidad neutral al riesgo es solo la forma de escribir ese precio como esperanza descontada.
El árbol CRR de una call en el dinero da errores de +0,62 con 51 pasos y de −0,53 con 50. ¿Qué se hace?
- ANada: el árbol está mal programado
- BUsar muchos más pasos, promediar dos árboles consecutivos o un árbol que converja de forma monótona, como el de Leisen-Reimer con un número impar de pasos
- CUsar siempre 50 pasos, porque el error es menor
- DCambiar el tipo de interés
Ver respuesta
Respuesta B. Usar muchos más pasos, promediar dos árboles consecutivos o un árbol que converja de forma monótona, como el de Leisen-Reimer con un número impar de pasos
La oscilación par-impar viene de dónde cae el strike respecto a los nodos finales. El error de CRR se reduce como 1/n pero oscila; el árbol de Leisen-Reimer de QuantLib con 51 pasos ya está a 0,002 puntos de Black-Scholes-Merton.
¿Cuándo coinciden Black-Scholes-Merton y Black-76?
- ANunca: una es para acciones y otra para futuros
- BSiempre que el forward que se pasa a Black-76 sea F = S·e^((r − q)T) para el mismo vencimiento: entonces son la misma fórmula
- CSolo con tipos de interés nulos
- DSolo para opciones americanas
Ver respuesta
Respuesta B. Siempre que el forward que se pasa a Black-76 sea F = S·e^((r − q)T) para el mismo vencimiento: entonces son la misma fórmula
Black-Scholes-Merton con rendimiento continuo q es Black-76 aplicado al forward de ese vencimiento. En la opción de ORX de la lección, las seis implementaciones coinciden en diez decimales. Por eso con índices es más robusto valorar sobre el forward implícito de la cadena que sobre el contado con un dividendo supuesto.
Una acción paga un dividendo de 1,04 USD dentro de 10 días. ¿Qué error se comete al valorar sus calls americanas con Black-Scholes-Merton europeo?
- ANinguno: el dividendo ya está en el precio de la acción
- BSe ignora el valor del ejercicio anticipado justo antes de la fecha ex, que en las calls dentro del dinero de KRN llega a 0,59 USD por acción
- CSolo un error de redondeo
- DSe sobrevalora la call
Ver respuesta
Respuesta B. Se ignora el valor del ejercicio anticipado justo antes de la fecha ex, que en las calls dentro del dinero de KRN llega a 0,59 USD por acción
El poseedor de una call americana puede ejercer la víspera del ex y cobrar el dividendo. La europea no puede. En la cadena de KRN del 3 de noviembre de 2025 la americana de strike 167,5 vale 17,31 y la europea 16,71: 0,59 USD por acción, 59 USD por contrato.
Una estimación de Monte Carlo con 100.000 trayectorias tiene un error típico de 0,62 puntos. ¿Cuántas hacen falta para bajarlo a 0,062 sin cambiar de método?
- A1.000.000
- B10.000.000
- C200.000
- DEs imposible
Ver respuesta
Respuesta B. 10.000.000
El error típico cae con la raíz del número de trayectorias: dividirlo por 10 exige 100 veces más. Por eso se usan reducciones de varianza: con el subyacente como variable de control, el error con 100.000 trayectorias baja de 0,62 a 0,30: lo mismo que cuadruplicar las trayectorias.

Ideas clave

El precio de una opción sale de la réplica (subyacente más bono); la probabilidad neutral al riesgo es una forma de escribirlo, no una previsión.
El árbol CRR converge a Black-Scholes-Merton con un error que cae como 1/n y oscila con la paridad de n; el de Leisen-Reimer con n impar converge mucho antes.
Black-Scholes-Merton con rendimiento continuo es Black-76 sobre el forward del mismo vencimiento: con índices y futuros conviene valorar sobre el forward.
Convenciones explícitas: plazo en días naturales / 365, tipo continuo de la curva en ese plazo, volatilidad anual en tanto por uno, precio por punto y multiplicador aparte.
Las opciones americanas con dividendos discretos piden un árbol con dividendos («escrowed»): recupera los precios de KRN con un error de milésimas y QuantLib lo confirma.
El error de Monte Carlo cae con la raíz del número de trayectorias; las variables antitéticas y de control lo aceleran, y es el método natural para opciones dependientes de la trayectoria, como las asiáticas.

Recursos

Documentación

QuantLib: documentación de la biblioteca (motores de valoración de opciones)quantlib.org
py_vollib (vollib): documentaciónvollib.org

Libros

Datos

Datos del nivel 3: opciones americanas de la acción KRN con precio de modelo (Parquet)Descarga · accion3_opciones.parquet

Volver al nivel 3

Convenciones antes de valorar nada

Un paso del árbol: réplica y probabilidad neutral al riesgo

El árbol de Cox-Ross-Rubinstein

Black-Scholes-Merton como límite del árbol

Valor antes del vencimiento

py_vollib y Black-76

Una opción sobre un futuro de crudo

Opciones americanas con dividendos discretos

Monte Carlo

Una opción asiática: donde Monte Carlo es la herramienta

Errores típicos

Ejercicios

Autoevaluación.

Ideas clave.

Recursos.

Documentación

Libros

Datos

Autoevaluación

Ideas clave

Recursos