Lesson 3.4.2 · 75 min

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.

Fig. — Convergence of the binomial tree to Black-Scholes-Merton: ORX 5,150 call, 30 days, 18 June 2025

02/07

Module: 3.4 · Option theory and pricing
Duration: 75 min
Prerequisites: 3.4.1 · Payoffs, put-call parity and arbitrage bounds, 0.7.2 · Probability and distributions, correlation versus causation and simple regression, 1.6.1 · Returns, volatility, rolling correlations and drawdown
Author: Jesús Cuesta
Updated: 04 Nov 2026

Level 3 lessons

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 source

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.

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"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}")

OUT [01]One step: replication versus risk-neutral expectation

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:

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))  # 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}")

OUT [02]European and American CRR tree

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):

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} 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")

OUT [03]BSM on an ORX call

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:

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,
                  "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)

OUT [04]Tree convergence: our own and QuantLib

     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

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:

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({
    "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)

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

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:

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)} 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")

OUT [06]One day's chain with Black-76

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):

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}, 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)

OUT [07]Black-76 on the PECG26 future, European and American

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 source

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]

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} 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")

OUT [08]KRN: the tree with discrete dividends

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):

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 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)

OUT [09]KRN with QuantLib finite differences

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

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.

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 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)

OUT [10]Monte Carlo of the ORX call: three estimators

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

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]Realized errors and QuantLib

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:

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), "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)

OUT [12]Arithmetic Asian with a geometric control

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

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

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"American {am:.4f} | European BSM {put_eu:.4f} | "
          f"premium {am - put_eu:.4f}")

OUT [13]American put on ORX

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.

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

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"formula {100 * DF * norm.cdf(d2):.3f}")

OUT [14]Digital by Monte Carlo

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).

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.

Each answer is folded under its question

Self-assessment

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?
- AThat the market expects the underlying to rise with a 75.6% probability
- BThat 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
- CThat the option is worth 75.6% of its intrinsic value
- DThat there is arbitrage, because it should be 0.5
Show answer
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.
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?
- ANothing: the tree is coded wrong
- BAlways use 50 steps, because the error is smaller
- CChange the interest rate
- DUse many more steps, average two consecutive trees or use a tree that converges monotonically, such as Leisen-Reimer with an odd number of steps
Show answer
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.
When do Black-Scholes-Merton and Black-76 coincide?
- AWhenever the forward passed to Black-76 is F = S·e^((r − q)T) for the same expiry: then they are the same formula
- BNever: one is for stocks and the other for futures
- COnly with zero interest rates
- DOnly for American options
Show answer
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.
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?
- ANone: the dividend is already in the stock price
- BOnly a rounding error
- CThe 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
- DThe call is overvalued
Show answer
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.
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?
- A1,000,000
- B200,000
- CIt is impossible
- D10,000,000
Show answer
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)quantlib.org
py_vollib (vollib): documentationvollib.org

Books

Data

Level 3 data: American options on the KRN stock with model price (Parquet)Download · accion3_opciones.parquet

Back to level 3

Conventions before pricing anything

One step of the tree: replication and risk-neutral probability

The Cox-Ross-Rubinstein tree

Black-Scholes-Merton as the limit of the tree

Value before expiry

py_vollib and Black-76

An option on a crude oil future

American options with discrete dividends

Monte Carlo

An Asian option: where Monte Carlo is the tool

Common errors

Exercises

Self-assessment.

Key points.

Resources.

Documentation

Books

Data

Self-assessment

Key points

Resources