Option Pricing Models

Option pricing models

../../../fin-core/skills/derivatives-pricing/SKILL.md tells you which library and what the Greeks mean. This skill is about the models themselves — what happens inside, and the four failure modes that produce a number rather than an exception.

Every number below is printed by scripts/option_models.py (runs in 6.1 s, QuantLib optional). ✅ Measured means this file produced it on 2026-09-09 with QuantLib 1.43, numpy 2.2.6, scipy 1.13.0, Python 3.11.3.

The four rules, in one line each: price Heston with the Albrecher/Gatheral form; average CRR odd and even step counts; report Monte Carlo bias and standard error separately; and check which CRR probability your reference library uses before you call your tree wrong.

1. 🚨 The little Heston trap — the 1993 formula is discontinuous, and it does not raise

Heston (1993) eq. (17)–(18) writes the log-characteristic function with

g = (b - rho*sigma*i*phi + d) / (b - rho*sigma*i*phi - d)     and     exp(+d*tau)

Albrecher, Mayer, Schoutens & Tistaert (2007), The Little Heston Trap, showed the algebraically identical form with 1/g and exp(-d*tau) is continuous under the principal branch of the complex log, while the original is not. Under the original, the inversion integral integrates across a jump, and you get a finite, plausible, wrong price.

✅ Measured (v0=0.04, kappa=1.5, theta=0.04, sigma=0.3, rho=-0.7, r=3%, q=0, ATM call S=K=100, fixed 8-point composite Gauss-Legendre quadrature so a discontinuity shows up as a wrong number rather than an adaptive-quadrature warning):

• First maturity that breaks: T = 2.0 years. Not an exotic long-dated case — a two-year option.
• The Albrecher form has 0 branch jumps at T = 1, 10 and 30.
• ✅ The NaN, located: at T=20 the 1993 C_1 stops being finite at phi = 165.3, where exp(+d*tau) overflows to inf; inf - inf in the log term is the NaN. The Albrecher form's exp(-d*tau) decays instead and is finite everywhere on the grid.
• The error is not monotone in T and it flips sign. At T=3 the 1993 form is too high by 1.39, at T=5 too low by 2.65. You cannot calibrate your way past it or spot it from a smooth residual plot.

1b. ✅ Maturity is not the only axis — kappa*theta breaks it too

The offending log is multiplied by a/sigma^2 = kappa*theta/sigma^2. ✅ Measured at T held at 2 years, varying kappa only:

Exact up to kappa = 1.00, wrong from kappa = 1.50 on, at an ordinary 2-year maturity. A fast mean-reverting calibration breaks the original formulation on short-dated options.

1c. ✅ QuantLib names both formulations, and agrees with the safe one to 1e-13

ql.AnalyticHestonEngine takes a ComplexLogFormula. ✅ Verified against QuantLib 1.43:

🚨 BranchCorrection is not the naive 1993 form. QuantLib repairs it by tracking how many times the log wraps. The unrepaired textbook transcription — the one in most blog posts and most first implementations — is what produces the table in §1. If you write your own Heston, write the Albrecher form. If you must use the 1993 form, you owe it a rotation counter.

2. 🚨 A CRR tree does not converge — it oscillates, and "more steps" buys almost nothing

✅ Measured, European put, S=K=100, T=1, r=5%, q=2%, sigma=20%, Black-Scholes limit 6.330081:

Even step counts sit below the limit, odd counts above it — the strike falls between two terminal nodes on one parity and near one on the other. Going 100 → 800 steps is 64× the work for an 8.0× smaller error, and the sign of your error is decided by the parity of a number you picked arbitrarily.

✅ The two fixes, measured against the same limit:

🚨 Richardson extrapolation only works if both legs share the parity. 2*C(2n) - C(n) doubles an even n to an even 2n, so the two errors are on the same side and the O(1/n) term cancels. Start from an odd n and 2n is even: the legs straddle the limit, and doubling the leg on the wrong side amplifies the oscillation. This is a one-character bug that makes a "converged" tree worse than the raw one.

✅ American early exercise, same tree, 800 steps: put 6.659528 vs European 6.327651 — an early-exercise premium of 0.331877. ✅ QuantLib's QdFpAmericanEngine prices the same American put at 6.660684, a difference of −0.001156 (the tree's own discretisation error, consistent with §2's magnitudes). derivatives-pricing §3 ranks the QuantLib American engines; this is only the tree.

3. 🚨 Monte Carlo: paths shrink the error bar, steps shrink the bias, and the error bar cannot see the bias

The single most expensive Monte Carlo mistake is reading a small standard error as "converged". ✅ Measured, European call, exact price 9.227006, seeded (numpy.random.default_rng):

• ✅ 4× the paths halves the error bar: 0.1371 → 0.0688 → 0.0345, ratios 1.99, 1.99 against the 2.00 that 1/sqrt(N) predicts.
• ✅ 52× the time steps does nothing to it: 0.0688 → 0.0686, ratio 1.00.
• ✅ Antithetic variates at the same 40,000 paths: 0.0688 → 0.0511, 1.35× smaller. Modest — antithetics help a monotone payoff, and a call's kink caps the gain.
• 🚨 The trap, in one row: with an Euler scheme and 1 step, 400,000 antithetic paths report a standard error of 0.0126 around a price that is 0.1238 too low — 10 standard errors of bias. Nothing in the output says so. Four steps cut the bias by 9×; the error bar barely moved.

Report both, always. The standard error measures how much your estimator wobbles; it says nothing about whether the thing it converges to is the price. For plain Black-Scholes use exact lognormal increments (scheme="log") and the bias is zero at any step count — Euler is only needed when the SDE has no exact solution.

⚠️ Compute the standard error over antithetic PAIRS, not over 2N correlated draws. The pair average is the i.i.d. unit; treating both legs as independent understates the error.

4. ✅ SABR (Hagan et al. 2002) — reproduced exactly against QuantLib

Hagan, Kumar, Lesniewski & Woodward (2002) eq. (2.17a)–(2.17c), lognormal implied vol, F=100, T=1, alpha=0.25, beta=0.6, nu=0.4, rho=−0.25:

✅ Worst |mine − ql.sabrVolatility| over all seven strikes: 0.0e+00 — bit-identical. That is the documented reproduction; there is no ambiguity left about the formula.

• 🚨 The formula is 0/0 at the money. z/x(z) with z -> 0 returns NaN if you transcribe it literally. Below |z| ~ 1e-7 use the series 1 - rho*z/2 + (2 - 3*rho^2)*z^2/12; that is the branch QuantLib takes, which is why the ATM strike also matches to 0.0e+00. The threshold is a choice and it differs between libraries — a smile that is smooth in one library and has a spike at K=F in another is this.
• ✅ The general formula's ATM value and eq. (2.18) agree to 0.0e+00 (0.040077965390 both ways). If your ATM special case disagrees with your general formula, one of the two is mistyped.
• 🚨 alpha is not the ATM vol. alpha=0.25 here gives an ATM vol of 0.0401, because the scaling is alpha / F^(1-beta) = 0.25/100^0.4. Seeding a calibration with "alpha ≈ ATM vol" diverges for any beta < 1. derivatives-pricing §5 measures the same trap from the calibration side; this is the formula side.

5. ✅ QuantLib's "crr" is not the textbook CRR probability

Found while benchmarking §2 against ql.BinomialVanillaEngine(process, "crr", n). Same u and d, different p:

✅ Reproduced to 5.8e-13 by using pu = 1/2 + (r−q−sigma²/2)·dt / (2·sigma·sqrt(dt)). Both probabilities converge to Black-Scholes, so neither is wrong — but the textbook one misses QuantLib by 1.2e-04 at n=100, six times the tolerance a careful person would pick for a tree test. ⚠️ This is verified by reproduction against the installed 1.43 wheel, not by reading the C++ (QuantLib ships no sdist — derivatives-pricing §3). crr_price(..., prob="quantlib") switches.

Rule: when a tree disagrees with a library at O(1/n), suspect the probability convention before the tree.

6. ✅ Also measured while cross-checking

• 🚨 MCEuropeanEngine(..., antitheticVariate=True, requiredSamples=100_000) counts PAIRS. ✅ Inverting QuantLib's own errorEstimate() against the measured payoff standard deviation (13.826 plain, 7.287 per antithetic pair) implies 100,241 payoffs in the plain run and 100,171 pairs in the antithetic one — about 200,342 paths, not 100,000. A timing comparison at equal requiredSamples is a 2× unfair comparison.
• ✅ QuantLib's MCEuropeanEngine reproduces §3's two facts in its own error estimate: 0.0437 at 1 step vs 0.0436 at 52 steps (steps do not move the error bar), 0.0230 antithetic.
• ✅ Black-Scholes-Merton put-call parity residual: +0.00e+00 (exactly zero in double precision) at S=K=100, T=1, r=5%, q=2%, sigma=20%; call 9.2270055082, put 6.3300806275.

7. What the script gives you

scripts/option_models.py is importable and library-free at import time (numpy + scipy only; QuantLib is imported inside quantlib_cross_checks, so the module works on a bare install and the demo degrades to a one-line notice):

Where this sits

• ../../../fin-core/skills/derivatives-pricing/SKILL.md — which library, Greek units (vega 100×, theta 365×), licences, and the ranked QuantLib American engines. Read that first.
• ../implied-vol-surface/SKILL.md — fitting a whole smile, SVI, and the static no-arbitrage checks. This skill prices one option; that one builds the surface.
• ../../../fin-libraries/skills/lib-quantlib/SKILL.md — 🚨 Settings.instance().evaluationDate is a global and a stale one returns NPV exactly 0.0 with no warning. Set it before you construct anything, including everything in §1c.
• ../../../fin-libraries/skills/lib-vollib/SKILL.md — vollib is the live package; py_vollib is a deprecated shim.
• ../../../fin-core/skills/options-backtesting/SKILL.md — assignment, expiry, pin risk, and historical chains. A model price is not a fill.