--- title: "3.4.2 Binomial, Black-Scholes, Black-76 y Monte Carlo" description: "Cuatro formas de valorar una opción, programadas desde cero y contrastadas con QuantLib y py_vollib: árbol CRR, BSM, Black-76 y Monte Carlo." url: "https://pythonparatrading.com/curso-python-trading/nivel-3/binomial-black-scholes-black-76-y-monte-carlo/" lang: "es" image: "https://pythonparatrading.com/og/curso-python-trading/nivel-3/binomial-black-scholes-black-76-y-monte-carlo.png" --- # 3.4.2 · Á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. - Curso: [Nivel 3: Derivados: futuros, spreads, forex y opciones](https://pythonparatrading.com/curso-python-trading/nivel-3/) - Módulo: 3.4 · Teoría y valoración de opciones - Duración: 75 min - Requisitos: [3.4.1 · Payoffs, paridad put-call y límites de arbitraje](https://pythonparatrading.com/curso-python-trading/nivel-3/payoffs-paridad-put-call-y-limites-de-arbitraje/), [0.7.2 · Probabilidad y distribuciones, correlación frente a causalidad y regresión simple](https://pythonparatrading.com/curso-python-trading/nivel-0/probabilidad-distribuciones-y-regresion/), [1.6.1 · Rentabilidades, volatilidad, correlaciones móviles y drawdown](https://pythonparatrading.com/curso-python-trading/nivel-1/volatilidad-correlaciones-y-drawdown/) ## Al terminar esta lección - Valorar una opción por réplica en un árbol de un paso y explicar por qué la probabilidad neutral al riesgo no es una previsión. - Programar un árbol CRR europeo y americano, medir su convergencia a Black-Scholes-Merton y compararlo con los árboles de QuantLib. - Aplicar Black-Scholes-Merton y Black-76 con las convenciones correctas (plazo, tipo continuo, dividendos, multiplicador) y contrastarlos con py_vollib y QuantLib. - Valorar opciones americanas con dividendos discretos con un árbol «escrowed» y recuperar los precios de modelo de la acción KRN. - Estimar precios por Monte Carlo con su error típico y reducirlo con variables antitéticas y de control, también en una opción asiática. ## 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. ```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. $$ ```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}") ``` ```output 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: ```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}") ``` ```output 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): ```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") ``` ```output 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: ```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) ``` ```output 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.* Cuatro curvas en función del índice, de 4.700 a 5.600: el payoff al vencimiento (cero hasta 5.150 y recta de pendiente 1 después) y el valor de la call a 30 días (destacada), 10 y 1 día. A 30 días la curva es suave y queda por encima del payoff (en el strike, unos 132 puntos de valor temporal); a 1 día casi se confunde con el payoff salvo cerca del strike. ## 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: ```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) ``` ```output 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: ```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") ``` ```output 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`): ```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) ``` ```output 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. ```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] ``` ```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") ``` ```output 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`): ```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) ``` ```output 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.* Diferencia entre el precio americano y el europeo de cada strike, de 145 a 225 USD. En las calls (destacadas) crece al bajar el strike hasta 0,88 USD en la de 145, cerca del dividendo de 1,04 (línea de trazos): compensa ejercer la víspera del ex. En las puts es pequeña, de 0,19 USD como mucho en las más dentro del dinero. - 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. ```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) ``` ```output 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.* En escala logarítmica, el error típico de los tres estimadores cae en rectas paralelas de pendiente −1/2: simple, de 6,44 puntos con 1.000 trayectorias a 0,20 con un millón; antitética, de 3,29 a 0,10; y con el índice como variable de control (destacada), de 3,07 a 0,096. ```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})") ``` ```output 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: ```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) ``` ```output 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 1. 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.** ```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}") ``` ```output 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í. 2. 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.** ```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}") ``` ```output 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). 3. 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. ## Autoevaluación 1. 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? - A. Que el mercado espera que el subyacente suba con un 75,6 % de probabilidad - 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 - C. Que la opción vale el 75,6 % de su valor intrínseco - D. Que hay arbitraje, porque debería ser 0,5 *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. 2. 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? - A. Nada: el árbol está mal programado - 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 - C. Usar siempre 50 pasos, porque el error es menor - D. Cambiar el tipo de interés *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. 3. ¿Cuándo coinciden Black-Scholes-Merton y Black-76? - A. Nunca: una es para acciones y otra para futuros - 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 - C. Solo con tipos de interés nulos - D. Solo para opciones americanas *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. 4. 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? - A. Ninguno: el dividendo ya está en el precio de la acción - 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 - C. Solo un error de redondeo - D. Se sobrevalora la call *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. 5. 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? - A. 1.000.000 - B. 10.000.000 - C. 200.000 - D. Es imposible *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)](https://www.quantlib.org/docs.shtml) · quantlib.org - [py_vollib (vollib): documentación](https://vollib.org/) · vollib.org ### Libros - [Cox, J. C., Ross, S. A. y Rubinstein, M. (1979): Option pricing: A simplified approach. Journal of Financial Economics 7(3), 229-263](https://doi.org/10.1016/0304-405X(79)90015-1) · doi.org - [Black, F. y Scholes, M. (1973): The Pricing of Options and Corporate Liabilities. Journal of Political Economy 81(3), 637-654](https://doi.org/10.1086/260062) · doi.org - [Black, F. (1976): The pricing of commodity contracts. Journal of Financial Economics 3(1-2), 167-179](https://doi.org/10.1016/0304-405X(76)90024-6) · doi.org - [Boyle, P. P. (1977): Options: A Monte Carlo approach. Journal of Financial Economics 4(3), 323-338](https://doi.org/10.1016/0304-405X(77)90005-8) · doi.org ### Datos - [Datos del nivel 3: opciones americanas de la acción KRN con precio de modelo (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/accion3_opciones.parquet) · Descarga · accion3_opciones.parquet Lección anterior: [3.4.1 · Payoffs, paridad put-call y límites de arbitraje](https://pythonparatrading.com/curso-python-trading/nivel-3/payoffs-paridad-put-call-y-limites-de-arbitraje/) Lección siguiente: [3.4.3 · Griegas: delta, gamma, vega, theta, rho, vanna, volga y charm](https://pythonparatrading.com/curso-python-trading/nivel-3/griegas-de-primer-y-segundo-orden/)