RadonLab C++

Monte Carlo Greeks for discontinuous payoffs, in header-only C++20, validated against closed-form references and cross-checked to 1e-9 against the companion Python library.

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 estimator for its sensitivities breaks: the payoff's derivative is zero almost everywhere, so the estimate collapses toward zero. This library implements the established remedies (likelihood-ratio, smoothing, Brownian-bridge conditional Monte Carlo), tests each against a known analytic answer, and benchmarks them so you can see which method to use and what it costs.

Header-only. C++20. Standard library only, no third-party dependencies. Apache-2.0.

Why this exists (and why C++)

This is the native counterpart to the Python radonlab package (github.com/quantsingularity/radonlab, on PyPI as radonlab). The two share the same methods and the same closed-form references, and the C++ analytics are checked against the Python library's machine-precision-validated values to 1e-9 in the test suite.

	Python (radonlab)	C++ (radonlab-cpp)
Role	Reference implementation; carries the research extensions	Header-only production port
Build	pip install, pure NumPy/SciPy	Header-only; C++20 concepts constrain payoff types
Extras	Merton jump-diffusion, path-dependent Malliavin delta	Adjoint algorithmic differentiation, exact second-order Greeks, threaded estimator

The C++ build keeps the hot loops in plain vectorizable form.

What's inside

Estimators (each returns a value with its Monte Carlo standard error):

Estimator	Identity / mechanism	Notes
Likelihood-ratio	$\dfrac{\partial P}{\partial \theta} = e^{-rT},\mathbb{E}!\left[f(S_T)\cdot \dfrac{\partial}{\partial \theta}\log p(S_T;\theta)\right]$	Differentiates the density, never the payoff, so it is unbiased for digitals. Delta, vega, gamma.
Smoothing	differentiate a $C^1$ approximation of bandwidth $h$	A tunable bias/variance trade-off.
Conditional Monte Carlo (barriers)	multiplies the terminal payoff by the Brownian-bridge non-crossing probability $\displaystyle\prod_i\left(1-\exp!\left(\dfrac{-2(X_i-b)(X_{i+1}-b)}{\sigma^2\Delta t}\right)\right)$	Unbiased for the continuously-monitored value at any number of steps, and Lipschitz in spot, so its delta comes out pathwise.
Pathwise / finite-difference	baselines	Kept so their failure modes are visible.

Closed-form Black-Scholes analytics are provided for European call/put, cash-or-nothing and asset-or-nothing digital calls, and continuously-monitored down-and-out / down-and-in calls.

Adjoint algorithmic differentiation (AAD)

A compact reverse-mode automatic differentiation tape (autodiff.hpp) computes a full gradient from a single reverse sweep. bs_greeks_aad differentiates the closed-form price and reproduces the analytic delta, vega and rho to 1e-9 (this also makes the analytics test suite self-contained: the Greeks are checked two independent ways, explicit formula and reverse-mode differentiation, with no reliance on external constants). mc_greeks_aad differentiates a Monte Carlo pricer path by path for pathwise or smoothed Greeks.

bs_second_order adds exact second-order Greeks (gamma, vanna, volga) via hyper-dual numbers, forward-mode second-order AD with no truncation error, validated against the closed form. A hyper-dual number carries a value and three infinitesimal parts:

$$ x = a + b,\varepsilon_1 + c,\varepsilon_2 + d,\varepsilon_1\varepsilon_2, \qquad \varepsilon_1^2=\varepsilon_2^2=0, \quad \varepsilon_1\varepsilon_2\neq 0. $$

Evaluating $f$ at $x = a+\varepsilon_1+\varepsilon_2$ gives

$$ f(x) = f(a) + f'(a)(\varepsilon_1+\varepsilon_2) + f''(a),\varepsilon_1\varepsilon_2, $$

so the $\varepsilon_1\varepsilon_2$ component is the exact second derivative, with no truncation error, unlike a finite difference of a first derivative. For a large Hessian a forward-over-reverse scheme scales better and is the natural next step.

AAD's value is the whole gradient at a cost independent of the number of inputs. Pricing a book of $N$ vanilla calls and taking every delta (examples/perf_aad.cpp):

N	AAD (1 sweep)	Bump (2N revaluations)	Speedup	Max grad diff
1	0.17 s	0.12 s	0.72x	7.7e-10
4	0.23 s	0.50 s	2.17x	6.7e-09
16	0.27 s	2.14 s	8.05x	7.5e-09
64	0.30 s	8.83 s	29.9x	7.6e-09
256	0.30 s	35.6 s	120x	7.7e-09
1024	0.30 s	141 s	468x	7.7e-09

Read honestly: for a single sensitivity the tape overhead makes AAD slower than a bump (0.72x at N = 1). The win is asymptotic in the number of sensitivities, which is exactly the regime of a real risk book. The gradient agrees with bump to the bump's own truncation error; AAD itself is machine-precision.

Threading

parallel_likelihood_ratio splits paths across threads, each with an independent spawned RNG, and matches the single-threaded estimate within Monte Carlo noise (verified in the test suite). Scaling numbers are not reported here because the build host has a single core; the code is written to scale on multicore hardware.

Build and test

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j
ctest --test-dir build --output-on-failure

Requires CMake 3.20+ and a C++20 compiler (tested with GCC 13). To consume it from another CMake project after cmake --install build:

find_package(radonlab CONFIG REQUIRED)
target_link_libraries(your_target PRIVATE radonlab::radonlab)

Because it is header-only you can also just add include/ to your include path.

Package managers

Via Conan:

conan create . --build=missing

# your conanfile.py
def requirements(self):
    self.requires("radonlab/0.1.1")

Via vcpkg as an overlay port (until the port is published in the curated registry):

vcpkg install radonlab --overlay-ports=<path-to-this-repo>/ports

Both install the headers plus a CMake package config, so find_package(radonlab CONFIG REQUIRED) works the same way regardless of how the library was installed.

Quick start

#include "radonlab/radonlab.hpp"
#include <cstdio>

int main() {
    using namespace radonlab;
    const GeometricBrownianMotion model(100.0, 0.03, 0.20, 1.0);
    const CashOrNothingCall payoff{100.0};

    NormalSampler rng(0);
    const GreekResult d = likelihood_ratio(model, payoff, Greek::Delta, 1'000'000, rng);
    const auto [lo, hi] = d.confidence_interval(0.95);
    std::printf("digital delta = %.6f  95%% CI [%.6f, %.6f]\n", d.value, lo, hi);
}

Adding a payoff is just a struct satisfying the TerminalPayoff concept (operator(), pathwise_derivative, smoothed, smoothed_derivative).

Benchmark

From examples/benchmark_digital_delta.cpp (1,000,000 paths). Delta of an at-the-money cash-or-nothing digital call, $S_0=K=100$, $r=3%$, $\sigma=20%$, $T=1$. Analytic reference delta = 0.019333. Ranked by standard error, the metric that reflects accuracy:

Method	Value	Std. error	Biased?
Likelihood-ratio	0.019373	2.838e-05	No
Smoothing	0.019263	4.568e-05	No
Finite-difference (1% bump)	0.019331	9.490e-05	No
Pathwise	0.000000	0.000e+00	Yes

Naive pathwise returns exactly zero, a 100% error. Among the unbiased methods the likelihood-ratio estimator has the lowest standard error (about 3.3 times lower than a 1% finite difference) and needs no bump to tune. Sweeping the bump size shows the finite-difference standard error diverging as the bump shrinks:

Rel. bump	Value	Std. error	Abs. error
1e-01	0.018608	2.359e-05	7.25e-04 (bias)
1e-02	0.019296	9.482e-05	3.76e-05
1e-03	0.019409	3.063e-04	7.55e-05
1e-04	0.018778	9.544e-04	5.55e-04

These numbers track the Python library's to within Monte Carlo noise, as they should: same estimators, same math.

Scope and honesty

This release covers the Black-Scholes / geometric Brownian motion world and the payoffs above, now with reverse-mode AAD and a threaded estimator. It is a rigorous, tested core, not a kitchen sink. Remaining natural next steps: stochastic volatility and jump models, Asian and lookback payoffs, forward-over- reverse AAD for large Hessians (exact second-order Greeks for a single option are already provided via hyper-dual numbers), and a SIMD-vectorized kernel. Where a method is expected to be biased (pathwise on a digital) the library shows it rather than hiding it.

License

Apache-2.0. See LICENSE.