radonlab 0.1.0

Monte Carlo Greeks for discontinuous payoffs: likelihood-ratio, smoothing and conditional Monte Carlo, benchmarked against closed-form references.

RadonLab

Monte Carlo Greeks for discontinuous payoffs, done correctly and benchmarked against closed-form references.

When an option's payoff is discontinuous (a digital that pays a fixed amount if the underlying finishes above a strike, a barrier that knocks out on touch), the textbook pathwise (infinitesimal perturbation) estimator for its sensitivities quietly breaks: the derivative of the payoff is zero almost everywhere, so the estimate collapses toward zero. radonlab implements the established remedies, tests each one against a known analytic answer, and benchmarks them side by side so you can see which method to use and what it costs.

Pure NumPy/SciPy. Apache-2.0. No copyleft dependencies.

Companion C++ implementation

A header-only C++20 port lives at github.com/quantsingularity/radonlab-cpp. The two libraries share the same methods and the same closed-form references, and the C++ analytics are cross-checked against this library's values (which are verified to machine precision by complex-step differentiation), so the two agree to 1e-9 across languages. The C++ side additionally provides adjoint algorithmic differentiation (a full Greek vector from one reverse sweep) and exact second-order Greeks. This Python library is the reference and carries the research extensions (Merton jump-diffusion and a path-dependent Malliavin delta).

What's inside

Estimators (all with per-path standard errors):

• Likelihood-ratio method (LRM): differentiates the density, never the payoff, so it is unbiased for digitals and other discontinuous payoffs. Delta, vega and gamma.
• Smoothing: replaces the payoff with a C¹ approximation of bandwidth h and applies the pathwise method to it. A tunable bias/variance trade-off.
• Conditional Monte Carlo: for barriers, integrates out the crossing event with the Brownian-bridge non-crossing probability, giving an estimator that is unbiased for the continuously-monitored price at any number of steps and Lipschitz in spot (so its delta comes out by the pathwise method).
• Pathwise and finite-difference (bump) baselines, kept so their failure modes on discontinuous payoffs are visible and measurable.

Closed-form Black-Scholes analytics (the ground truth) for European calls/puts, cash-or-nothing and asset-or-nothing digitals, and continuously-monitored down-and-out / down-and-in calls. Every explicit Greek formula is checked against complex-step differentiation of the price to machine precision.

Install

pip install -e ".[benchmark,test]"

Runtime dependencies are just numpy and scipy. pandas/matplotlib are only needed for the benchmark table and plots.

Quick start

from radonlab import GeometricBrownianMotion, CashOrNothingCall, Greek
from radonlab import likelihood_ratio, pathwise, analytics

model = GeometricBrownianMotion(S0=100, r=0.03, sigma=0.20, T=1.0)
payoff = CashOrNothingCall(strike=100)

# Ground truth
ref = analytics.cash_or_nothing_call(100, 100, 0.03, 0.20, 1.0, greek="delta")

# Likelihood ratio: unbiased for the digital
lrm = likelihood_ratio(model, payoff, Greek.DELTA, n_paths=1_000_000)
print(lrm)                    # delta=... +/- ...  [likelihood-ratio]
print(lrm.confidence_interval(0.95))

# Naive pathwise: collapses to ~0, badly biased
pw = pathwise(model, payoff, Greek.DELTA, n_paths=1_000_000)
print(pw.value, "vs analytic", ref)

The math, briefly

Under geometric Brownian motion S_T = S₀·exp((r − q − ½σ²)T + σ√T·Z), with Z ~ N(0,1), a price is P = e^{−rT} E[f(S_T)].

• LRM writes ∂P/∂θ = e^{−rT} E[f(S_T)·∂_θ log p(S_T;θ)]. The score functions are Z/(S₀σ√T) for delta, (Z²−1)/σ − √T·Z for vega, and (Z² − σ√T·Z − 1)/(S₀²σ²T) for gamma. f is never differentiated, so a jump in f is no problem.
• Pathwise writes ∂P/∂θ = e^{−rT} E[f′(S_T)·∂_θ S_T]. Fine when f is Lipschitz (vanilla), degenerate when f jumps (digital).
• Conditional MC for a down-and-out call multiplies the terminal payoff by ∏ᵢ(1 − exp(−2(Xᵢ−b)(Xᵢ₊₁−b)/(σ²Δt))), the exact bridge probability of not breaching b = ln B between simulated points.

Full references are in each module's docstring (Broadie-Glasserman 1996, Glasserman 2003, Reiner-Rubinstein 1991, and the smoothed-perturbation and Malliavin literature).

Benchmark

Reproduce with python examples/benchmark_digital_delta.py (1,000,000 paths, seed 0). Delta of an at-the-money cash-or-nothing call, S₀=K=100, r=3%, σ=20%, T=1. Analytic reference delta = 0.019333. Ranked by standard error, the metric that actually reflects an estimator's accuracy:

method	value	std error	|error| vs analytic	biased?
likelihood-ratio	0.01939	2.84e-05	5.9e-05	no
smoothing (h=default)	0.01921	4.56e-05	1.2e-04	no
finite-difference (1% bump)	0.01934	9.49e-05	8.5e-06	no
pathwise	0.00000	0.00e+00	1.9e-02	yes

Two things to read off it. First, pathwise is 100% wrong: it returns exactly zero, because the digital's derivative is zero almost everywhere. Second, among the unbiased methods the likelihood-ratio estimator has the lowest variance (about 3.3× smaller standard error than the 1% finite difference) and needs no bump to tune.

And finite difference forces a bad choice: the same run sweeping the bump size shows the standard error exploding as the bump shrinks, while a large bump introduces its own bias:

rel. bump	value	std error	|error|
1e-01	0.018620	2.36e-05	7.1e-04 (bias)
1e-02	0.019311	9.49e-05	2.2e-05
1e-03	0.019161	3.04e-04	1.7e-04
1e-04	0.017274	9.15e-04	2.1e-03

The likelihood-ratio method sidesteps that trade-off entirely. This is the whole reason the library exists.

Testing

pytest

The suite validates the analytics by complex-step differentiation, confirms each Monte Carlo estimator agrees with the analytic Greek within a few standard errors, and documents the pathwise failure on digitals as an explicit test.

Scope and honesty

This release covers the Black-Scholes / GBM world and, in radonlab.merton, the Merton jump-diffusion model with a closed-form digital reference: it shows the likelihood-ratio delta weight is unchanged by jumps and validates the estimator against the closed form. It is a rigorous, tested core rather than a kitchen sink: the estimators and their validation are real, and where a method is expected to be biased (pathwise on a digital), the library shows it rather than hiding it. It also includes a genuinely path-dependent frontier estimator in radonlab.asian: the likelihood-ratio / Malliavin delta of a geometric Asian digital, validated against the closed form, with a weight built from the whole path. Stochastic volatility, arithmetic-Asian and lookback payoffs, and a JIT/vectorized backend are natural next steps.

License

Apache-2.0. See LICENSE.