--- title: "3.4.2 Binomial, Black-Scholes, Black-76 and Monte Carlo" description: "Four ways to price an option, coded from scratch and checked against QuantLib and py_vollib: CRR tree, BSM, Black-76 and Monte Carlo." url: "https://pythonparatrading.com/en/python-trading-course/level-3/binomial-black-scholes-black-76-monte-carlo/" lang: "en" image: "https://pythonparatrading.com/og/en/python-trading-course/level-3/binomial-black-scholes-black-76-monte-carlo.png" --- # 3.4.2 · Binomial trees, Black-Scholes-Merton, Black-76 and Monte Carlo > A one-step tree already contains the whole theory: replicating the option with the underlying and a bond. From there come the Cox-Ross-Rubinstein tree, the Black-Scholes-Merton formula as its limit, Black-76 for options on futures and Monte Carlo for what has no formula. Each method is coded with NumPy and checked against QuantLib, py_vollib and the model prices in the data. - Course: [Level 3: Derivatives: futures, spreads, forex and options](https://pythonparatrading.com/en/python-trading-course/level-3/) - Module: 3.4 · Option theory and pricing - Duration: 75 min - Prerequisites: [3.4.1 · Payoffs, put-call parity and arbitrage bounds](https://pythonparatrading.com/en/python-trading-course/level-3/option-payoffs-put-call-parity-arbitrage/), [0.7.2 · Probability and distributions, correlation versus causation and simple regression](https://pythonparatrading.com/en/python-trading-course/level-0/probability-distributions-and-regression/), [1.6.1 · Returns, volatility, rolling correlations and drawdown](https://pythonparatrading.com/en/python-trading-course/level-1/returns-volatility-correlation-and-drawdown/) ## By the end of this lesson - Price an option by replication on a one-step tree and explain why the risk-neutral probability is not a forecast. - Code a European and American CRR tree, measure its convergence to Black-Scholes-Merton and compare it with QuantLib's trees. - Apply Black-Scholes-Merton and Black-76 with the correct conventions (time to expiry, continuous rate, dividends, multiplier) and check them against py_vollib and QuantLib. - Price American options with discrete dividends with an “escrowed” tree and recover the model prices of the KRN stock. - Estimate prices by Monte Carlo with their standard error and reduce it with antithetic and control variates, also on an Asian option. ## Conventions before pricing anything Lesson 3.4.1 bounded option prices without any model. Going from bounds to a number requires assuming how the underlying moves. Before the formulas, the conventions used throughout the module, because half of the pricing errors in practice are about units, not about the model: | Quantity | Course convention | Typical error | |---|---|---| | Time $T$ | calendar days to expiry / 365 (they expire at the close) | business days / 252 in the pricing formula | | Rate $r$ | continuous, from the curve interpolated at $T$ | simple rate or the wrong maturity | | Dividends | index: continuous yield $q$ (or the forward); stock: discrete amounts on their ex-dates | continuous $q$ for a stock with a dividend next week | | Volatility $\sigma$ | annual, as a decimal (0.217 = 21.7%) | passing 21.7 instead of 0.217 | | Price | per unit of the underlying (points or USD per share) | forgetting the multiplier (100 for ORX and KRN; 1,000 barrels for crude) | ## One step of the tree: replication and risk-neutral probability An underlying is worth $S_0 = 100$ and in one year it will be worth $uS_0 = 110$ or $dS_0 = 90$. A call with strike 100 pays 10 or 0. With a continuous rate of 5%, the call can be **replicated** with $\Delta$ units of the underlying and a bond $B$: it is enough to require that the portfolio pays the same as the call in both states. ```mermaid flowchart LR S0["S = 100"] -- up --> SU["S = 110, C = 10"] S0 -- down --> 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"Replication: {delta:.2f} units and bond {B:.2f} " f"-> C = {delta * S0 + B:.4f}") p = (math.exp(r * T) - d) / (u - d) print(f"risk-neutral p = {p:.4f} -> C = " f"{math.exp(-r * T) * (p * Cu + (1 - p) * Cd):.4f}") ``` ```output Replication: 0.50 units and bond -42.81 -> C = 7.1947 risk-neutral p = 0.7564 -> C = 7.1947 ``` The portfolio of half a unit of the underlying financed with a loan of 42.81 pays exactly what the call does and costs 7.19: that is the price, because any other would allow buying one and selling the other without risk. The same figure comes out as a discounted expectation with a probability $p = (e^{rT} - d)/(u - d)$ = 0.756. That $p$ is not the probability that the underlying rises: it is the weight that makes the underlying, discounted at the risk-free rate, neither gain nor lose on average. The real probability of an up move appears nowhere, and that is why two investors who disagree about it agree on the price of the option. ## The Cox-Ross-Rubinstein tree Repeating the step $n$ times with $u = e^{\sigma\sqrt{\Delta t}}$, $d = 1/u$ and $\Delta t = T/n$ gives the **Cox, Ross and Rubinstein** (1979) tree. With a continuous dividend yield $q$, the probability is $p = (e^{(r-q)\Delta t} - d)/(u - d)$. It is priced backward: at the last step, the payoff; at each earlier node, the discounted expectation of its two children and, if the option is American, the maximum between that value and exercising now. With NumPy one vector per step is enough: ```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)) # final nodes 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 ``` With 500 steps, the European call of the textbook example ($S = K = 100$, one year, 5%, 20%) is worth 10.4466, the European put 5.5695 and the American put 6.0888: half a point more, the value of being able to collect the strike earlier. ## Black-Scholes-Merton as the limit of the tree When $n \to \infty$, the tree converges to the **Black-Scholes-Merton** (BSM) formula, that of an underlying whose logarithm follows a Gaussian random walk with constant volatility: $$ 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}. $$ It is applied to an option from the data: the ORX call with strike 5,150 (the closest to the index) of the July 18, 2025 expiry, on June 18. The index, the dividend yield and the curve come from `opciones3_indice.csv`; the volatility, from the `iv_modelo` column, the model implied volatility with which the chain was generated (lesson 3.4.4 recovers it from the prices): ```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} years ({(venc - hoy).days} days), " f"r = {r:.4%}, q = {q:.4%}, sigma = {sig}") print(f"BSM: {c_bsm:.4f} | bid {op['bid']} / ask {op['ask']}" f" | per contract: {100 * c_bsm:,.0f} USD") ``` ```output S = 5148.06, T = 0.0822 years (30 days), r = 4.1910%, q = 1.3438%, sigma = 0.217 BSM: 132.5896 | bid 131.8 / ask 133.1 | per contract: 13,259 USD ``` The call is worth 132.59 points, within its spread of 131.80 to 133.10, and 13,259 USD per contract. The header figure shows how the CRR tree reaches that number as the number of steps grows, next to QuantLib's Leisen-Reimer tree, which places the nodes so that the strike sits at the center and converges much sooner. The cell compares our own two trees and QuantLib's with 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, "own CRR": 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 own CRR 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 ``` Three lessons from the table: - The CRR error **oscillates**: with 50 steps the tree undervalues by 0.53 points and with 51 it overvalues by 0.62. The sign depends on whether the strike falls near a final node or between two. Averaging two consecutive trees or using many steps corrects almost all of it; with 500 steps the error is 0.03. - QuantLib's CRR is not identical to ours, although the difference is in the ten-thousandths and shrinks with $n$: QuantLib does not build the tree with exactly the same formulas. Both converge to the same limit. - QuantLib's Leisen-Reimer tree converges with very few steps **if n is odd** (0.0024 with 51); with even n, in QuantLib 1.44, it gives errors of up to 2.6 points with 50 steps and 0.27 with 500. Before trusting a library engine it has to be checked against something known, which is what this cell does. ## Value before expiry BSM also gives the value of the option at any moment before expiry, and with it the result of a position before it expires. The figure draws the call from the previous section as a function of the index at 30, 10 and 1 days to expiry, with the same volatility, next to its payoff: far from expiry the curve is smooth and lies above the payoff (time value); as expiry approaches, it hugs the payoff, and the transition concentrates at the strike. The slopes and curvatures of these curves are the Greeks of lesson 3.4.3. *Fig. — Value of an ORX 5,150 call 30, 10 and 1 days before expiry, with Black-Scholes-Merton (21.7% volatility), against its payoff.* Four curves as a function of the index, from 4,700 to 5,600: the payoff at expiry (zero up to 5,150, then a slope-1 line) and the value of the call 30 days (highlighted), 10 and 1 day before expiry. At 30 days the curve is smooth and above the payoff (about 132 points of time value at the strike); at 1 day it almost merges with the payoff except near the strike. ## py_vollib and Black-76 **py_vollib** 1.0.12 is an implementation of the Black, Black-Scholes and Black-Scholes-Merton formulas based on Peter Jäckel's “Let's Be Rational” algorithm. In this version the code lives in the `vollib` package (1.0.11), which is installed with it: `import py_vollib` still works, but it emits a `DeprecationWarning` asking to import `vollib`. Here `vollib`, which is the same library, is imported directly, so as not to clutter the output. Options on futures are priced with **Black-76** (Black, 1976): the same formula written on the futures price $F$, which has no carry cost (buying it costs nothing today): $$ 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}. $$ If $F = S e^{(r-q)T}$ is the forward of the same expiry, Black-76 and BSM are the same formula. For an index this means that it can be priced on the chain's implied forward (lesson 3.4.1) without knowing the dividend: ```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({ "own BSM": c_bsm, "own Black-76": 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 own BSM 132.5896128016 own Black-76 132.5896128016 vollib BSM 132.5896128016 vollib Black 132.5896128016 QuantLib BSM 132.5896128016 QuantLib Black 132.5896128016 dtype: str ``` The six figures agree to ten decimals. What does change the price is the input. The whole chain of that day, priced with Black-76 on the forward from `opciones3_verdad.parquet` and the volatility of each strike, falls within its spread: ```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)} options; within the spread: " f"{dentro.mean():.1%}") print(f"|model - mid|: median {np.median(abs(modelo - medio)):.3f}" f", max {np.max(abs(modelo - medio)):.3f} points") ``` ```output 266 options; within the spread: 100.0% |model - mid|: median 0.046, max 2.586 points ``` ### An option on a crude oil future The data has no options on futures, but it does have the full curve of the fictitious crude PEC (it mimics NYMEX WTI: 1,000 barrels, in USD per barrel). A hypothetical European option on the February 2026 contract (`PECG26`) is priced on October 15, 2025, with 60 days to expiry and the realized volatility of the contract's own previous 63 sessions. If the option were American (each contract's spec sheet says so), Black-76 is no longer enough: QuantLib prices the American one with a tree on the future (`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}, realized vol. = {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, "American QL": am.NPV(), "USD/contract (eu.)": 1000 * eu}) pd.DataFrame(filas).round(4) ``` ```output F = 66.11, realized vol. = 0.206, r = 0.0398 tipo K Black-76 vollib American QL USD/contract (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 ``` A put with strike 66 on a future at 66.11 is worth 2.13 USD per barrel, 2,132 USD per 1,000-barrel contract. Early exercise of an option on futures is worth little except deep in the money: the only thing gained by exercising early is the interest on the intrinsic premium, and on the 75 put, the deepest in the money, it is 0.027 USD per barrel (27 USD per contract). ## American options with discrete dividends An index pays dividends almost daily and a continuous yield $q$ describes it well. A stock pays four a year, on known dates, and that changes two things: the price drops at once on the ex-date, and an American call can be worth more by exercising it the day before. The **KRN** data (fictitious stock, American options, quarterly dividend) comes with the chain and its model price, computed with a 150-step CRR tree with the “escrowed” model: the tree is built on $S$ minus the present value of the dividends up to expiry, and at each node the present value of the remaining ones is added back to compare with exercise. The goal is to recover those prices with our own tree and with QuantLib. ```mermaid flowchart TB O[Option] --> E{European?} E -- yes --> P{Simple payoff?} P -- yes --> BS[BSM or Black-76] P -- no --> MC[Monte Carlo] E -- no --> D{Dividends?} D -- discrete --> AR[Escrowed tree] D -- continuous --> AR2[Tree or 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} days, r = {r_k:.4%};", "dividend", 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) # PV at each step of the remaining ones 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)} options; maximum error against the data: " f"{err.max():.4f} USD") ``` ```output S = 184.27, T = 18 days, r = 3.8990%; dividend {Timestamp('2025-11-13 00:00:00'): 1.04} 50 options; maximum error against the data: 0.0011 USD ``` Our own tree reproduces the 50 model prices of that slice with a maximum error of 0.0011 USD, the one left by the roundings of the close (to the cent), the curve and the volatility (to four decimals). QuantLib solves the same problem with finite differences on the Black-Scholes equation, with the dividends as amounts on their dates and the same “escrowed” model (`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 against the data: maximum error " 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 against the data: maximum error 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, with a 400 × 400 grid, stays within 0.008 USD of the prices in the data; the differences are from discretization (of its grid and of the 150-step tree), far below a 0.01 tick. The table and the figure show the value of early exercise: *Fig. — Early-exercise value of KRN American options by strike, 3 November 2025, 21 November expiry with the dividend going ex on the 13th.* Difference between the American and the European price of each strike, from 145 to 225 USD. For calls (highlighted) it grows as the strike falls, up to 0.88 USD at 145, close to the 1.04 dividend (dashed line): it pays to exercise the day before the ex-date. For puts it is small, at most 0.19 USD for the deepest in the money. - For in-the-money **calls**, early exercise is worth a good part of the dividend: with the ex-date 10 days away, the 165 call is worth 0.65 USD more than the European one (63% of the 1.04 dividend) and the 175 one, 0.38. Whoever holds it and does not exercise it the day before the ex-date gives that amount away; whoever has sold it must count on assignment. - For **puts** the premium is small at 18 days: the interest from collecting the strike earlier is a matter of cents and, moreover, the dividend, which will make the stock fall, argues for waiting. Pricing these calls with European BSM and an equivalent continuous yield undervalues precisely the options that trade the most around the ex-date. It is the conventions error from the opening table, with a name of its own. ## Monte Carlo When there is no formula or convenient tree (path-dependent options, several underlyings, stochastic volatility models), one simulates. Under the risk-neutral measure, the index at expiry is $S_T = F\,e^{-\sigma^2 T/2 + \sigma\sqrt{T} Z}$ with $Z$ standard normal, and the price is the mean of the discounted payoffs (Boyle, 1977). The simulation is exact at $T$: there is no need to discretize the path. What is not exact is the mean, which has a **standard error** $s/\sqrt{N}$. Two techniques reduce it without more paths: - **Antithetic variates**: each $Z$ is also used with the opposite sign; the two payoffs are negatively correlated and their mean varies less. - **Control variate**: something whose exact value is known is simulated at the same time (here, $e^{-rT}S_T$ itself, which is worth $e^{-rT}F$) and the estimate is corrected with its error: $\hat C_{cv} = \hat C - \beta(\hat X - E[X])$, with $\beta$ the slope of the regression of the payoff on the 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 pairs would suffice beta = np.cov(pago, ctrl)[0, 1] / ctrl.var(ddof=1) cv = pago - beta * (ctrl - DF * F) for nombre, x in [("simple", pago), ("antithetic", 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 antithetic 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 ``` The standard error falls as $1/\sqrt{N}$ for all three estimators: multiplying the paths by 100 divides it by 10. The control variate also divides it by somewhat more than two for any $N$ (from 0.62 to 0.30 with 100,000), because the payoff of an at-the-money call is highly correlated with the index; the antithetic variates, by 1.9 (and they cost twice the payoff evaluations). The figure shows it on a logarithmic scale. The realized errors against BSM are of the order of those standard errors: *Fig. — Monte Carlo standard error of the ORX call by number of paths, with three estimators.* On log scales, the standard error of the three estimators falls along parallel lines of slope −1/2: plain, from 6.44 points with 1,000 paths to 0.20 with a million; antithetic, from 3.29 to 0.10; and with the index as control variate (highlighted), from 3.07 to 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 antithetic 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, with 200,000 antithetic paths and its own generator, gives 132.37 ± 0.23: within one standard error of BSM. ### An Asian option: where Monte Carlo is the tool An **Asian** call pays $\max(\overline{S} - K, 0)$, with $\overline{S}$ the arithmetic mean of the closes over a period. It has no closed formula, but its sibling with a **geometric** mean does (the logarithm of a geometric mean of lognormals is normal). It is the ideal control variate: almost identical to the arithmetic one and with an exact value, which QuantLib computes analytically. An Asian call on ORX is priced with the mean of the 20 NYSE sessions left until the July expiry: ```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), "fixings, from", ses[0].date(), "to", 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=["geometric MC", "geometric exact (QL)", "arithmetic MC simple", "arithmetic MC control", "arithmetic QL MC control"]).round(4) ``` ```output 20 fixings, from 2025-06-20 to 2025-07-18 precio error_tipico geometric MC 81.8746 0.3860 geometric exact (QL) 81.6093 0.0000 arithmetic MC simple 82.7353 0.3893 arithmetic MC control 82.4677 0.0032 arithmetic QL MC control 82.4743 0.0046 ``` The simulated geometric one matches the exact one within its standard error, which validates the simulation. With it as a control, the standard error of the arithmetic one drops from 0.39 to 0.003 points: about 120 times less, the equivalent of multiplying the paths by almost 15,000. QuantLib, with the same control, gives the same price. The Asian option is worth much less than the European call (82.47 against 132.59) because the mean over a month varies less than the close of the last day. ## Common errors - **Units.** Business days / 252 for time, volatility in percent, a simple rate as if it were continuous or the forgotten multiplier. Each one, on its own, moves the price more than the difference between two models. - **Spot and an assumed dividend for an index.** The chain's implied forward already discounts the dividend the market expects; with Black-76 on that forward there is no need to guess it. - **European formulas for American options.** With discrete dividends or deep in-the-money puts, BSM undervalues; a tree or finite differences must be used. - **Trees with few steps.** CRR oscillates with the parity of $n$; with a library, the engine is checked against a known case (like LR with even $n$ in the convergence cell). - **Monte Carlo without a standard error.** An estimate without its error cannot be compared with anything. The seed is fixed so that it is reproducible and the error is always reported. - **Volatility.** Everything above assumes a given volatility. Which one to use (implied per strike, realized, from a surface) is the subject of lessons 3.4.4 to 3.4.7; here the data model's volatility was used to isolate the method. ## Exercises 1. Price with `crr` the ORX put from the example (strike 5,150, July expiry, on June 18, 2025) as if it were American, with 100, 500 and 1,000 steps. How much is early exercise worth against the European one? **Solution.** ```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"American {am:.4f} | European BSM {put_eu:.4f} | " f"premium {am - put_eu:.4f}") ``` ```output 100 American 123.1397 | European BSM 122.5030 | premium 0.6367 500 American 123.2939 | European BSM 122.5030 | premium 0.7908 1000 American 123.3085 | European BSM 122.5030 | premium 0.8055 ``` Early exercise of an at-the-money put at 30 days is worth 0.81 points with 1,000 steps: the interest from collecting the strike a few days earlier, minus the index dividend that is forgone. It is little because the at-the-money put still has a lot of time value. Index options such as SPX are European, so this premium does not exist in them; in American options on ETFs on the same index, it does. 2. Estimate by Monte Carlo a **digital** ORX call (it pays 100 points if the index ends above 5,150) with 100,000 paths and compare it with the formula $100\,e^{-rT}N(d_2)$. **Solution.** ```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"formula {100 * DF * norm.cdf(d2):.3f}") ``` ```output MC 49.688 ± 0.158 | formula 49.846 ``` The digital is worth the risk-neutral probability of ending in the money, discounted and times 100. Monte Carlo recovers it within its standard error. Its relative error is larger than that of the call because the payoff is discontinuous; a quoted digital call is replicated in practice with a narrow call spread (lesson 3.4.1). 3. A colleague prices KRN's calls of November 3, 2025 with BSM and a continuous yield $q = 4 \times 1.04 / S$. Explain which two errors they make and at which strikes they are most noticeable. **Solution.** The first is spreading over the whole year a dividend that falls on a specific day, 10 days away: the 18-day forward carries the whole dividend (1.04 USD), not 18/365 of four dividends, so the continuous yield underestimates the drop the market discounts and overvalues the European calls. The second is treating American options as European: it ignores early exercise the day before the ex-date, which for in-the-money calls is worth almost the whole dividend (0.65 USD on the 165 one). The two errors have opposite signs and do not offset equally at all strikes: the first affects the whole chain; the second, mostly the in-the-money calls. ## Self-assessment 1. In a one-step tree with u = 1.10, d = 0.90 and a 5% rate, the risk-neutral probability of an up move is 0.756. What does it mean? - A. That the market expects the underlying to rise with a 75.6% probability - B. That it is the weight that makes the discounted underlying a martingale; the option price comes from replication and does not depend on the real probability of an up move - C. That the option is worth 75.6% of its intrinsic value - D. That there is arbitrage, because it should be 0.5 *Answer B: That it is the weight that makes the discounted underlying a martingale; the option price comes from replication and does not depend on the real probability of an up move.* The price is set by replication: the portfolio of 0.5 units of the underlying and a loan of 42.81 pays the same as the call in both states and costs 7.19. The risk-neutral probability is just a way of writing that price as a discounted expectation. 2. The CRR tree of an at-the-money call gives errors of +0.62 with 51 steps and −0.53 with 50. What should be done? - A. Nothing: the tree is coded wrong - B. Always use 50 steps, because the error is smaller - C. Change the interest rate - D. Use many more steps, average two consecutive trees or use a tree that converges monotonically, such as Leisen-Reimer with an odd number of steps *Answer D: Use many more steps, average two consecutive trees or use a tree that converges monotonically, such as Leisen-Reimer with an odd number of steps.* The even-odd oscillation comes from where the strike falls relative to the final nodes. The CRR error shrinks as 1/n but oscillates; QuantLib's Leisen-Reimer tree with 51 steps is already within 0.002 points of Black-Scholes-Merton. 3. When do Black-Scholes-Merton and Black-76 coincide? - A. Whenever the forward passed to Black-76 is F = S·e^((r − q)T) for the same expiry: then they are the same formula - B. Never: one is for stocks and the other for futures - C. Only with zero interest rates - D. Only for American options *Answer A: Whenever the forward passed to Black-76 is F = S·e^((r − q)T) for the same expiry: then they are the same formula.* Black-Scholes-Merton with continuous yield q is Black-76 applied to the forward of that expiry. On the ORX option in the lesson, the six implementations agree to ten decimals. That is why with indices it is more robust to price on the chain's implied forward than on spot with an assumed dividend. 4. A stock pays a dividend of 1.04 USD in 10 days. What error is made by pricing its American calls with European Black-Scholes-Merton? - A. None: the dividend is already in the stock price - B. Only a rounding error - C. The value of early exercise just before the ex-date is ignored, which for KRN's in-the-money calls reaches 0.59 USD per share - D. The call is overvalued *Answer C: The value of early exercise just before the ex-date is ignored, which for KRN's in-the-money calls reaches 0.59 USD per share.* The holder of an American call can exercise the day before the ex-date and collect the dividend. The European one cannot. In KRN's chain of November 3, 2025 the American with strike 167.5 is worth 17.31 and the European 16.71: 0.59 USD per share, 59 USD per contract. 5. A Monte Carlo estimate with 100,000 paths has a standard error of 0.62 points. How many are needed to bring it down to 0.062 without changing method? - A. 1,000,000 - B. 200,000 - C. It is impossible - D. 10,000,000 *Answer D: 10,000,000.* The standard error falls with the square root of the number of paths: dividing it by 10 requires 100 times more. That is why variance reduction is used: with the underlying as a control variate, the error with 100,000 paths drops from 0.62 to 0.30: the same as quadrupling the paths. ## Key points - The price of an option comes from replication (underlying plus bond); the risk-neutral probability is a way of writing it, not a forecast. - The CRR tree converges to Black-Scholes-Merton with an error that falls as 1/n and oscillates with the parity of n; Leisen-Reimer with odd n converges much sooner. - Black-Scholes-Merton with continuous yield is Black-76 on the forward of the same expiry: with indices and futures it pays to price on the forward. - Explicit conventions: time in calendar days / 365, continuous rate from the curve at that maturity, annual volatility as a decimal, price per point and multiplier separately. - American options with discrete dividends call for a tree with dividends (“escrowed”): it recovers KRN's prices with an error of thousandths and QuantLib confirms it. - The Monte Carlo error falls with the square root of the number of paths; antithetic and control variates speed it up, and it is the natural method for path-dependent options, such as Asian ones. ## Resources ### Documentation - [QuantLib: library documentation (option pricing engines)](https://www.quantlib.org/docs.shtml) · quantlib.org - [py_vollib (vollib): documentation](https://vollib.org/) · vollib.org ### Books - [Cox, J. C., Ross, S. A. and 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. and 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 ### Data - [Level 3 data: American options on the KRN stock with model price (Parquet)](https://pythonparatrading.com/curso/datos/nivel-3/accion3_opciones.parquet) · Download · accion3_opciones.parquet Previous lesson: [3.4.1 · Payoffs, put-call parity and arbitrage bounds](https://pythonparatrading.com/en/python-trading-course/level-3/option-payoffs-put-call-parity-arbitrage/) Next lesson: [3.4.3 · Greeks: delta, gamma, vega, theta, rho, vanna, volga and charm](https://pythonparatrading.com/en/python-trading-course/level-3/first-and-second-order-option-greeks/)