options-pricing-engine-rs

A Rust options pricing library covering Black-Scholes-Merton, Black-76, Heston (1993), Bates (1996), and Dupire local volatility, with full analytic Greeks where closed forms exist, a Halley-iteration implied vol solver, Levenberg-Marquardt calibration (single-start, multistart, and differential-evolution global search) for both Heston and Bates, no-arbitrage surface repair, and a Monte Carlo engine (full truncation Euler or Andersen QE) for path-dependent payoffs. Built for a vol surface update cycle, not a scripting exercise.

License: MIT. See LICENSE.

Contents

• Models
• Design
• Build
• Usage
• Testing
• Performance
• Known limitations and roadmap
• Dependencies
• References

Models

| Model | Pricing method | Greeks | |---|---|---| | Black-Scholes-Merton | Closed form | Full analytic: Δ, Γ, ν, Θ, ρ, vanna, volga | | Black-76 | Closed form | Full analytic | | Heston (1993) | Albrecher et al. (2007) stable characteristic function, adaptive Gauss-Kronrod-15 quadrature | Bump-and-reprice (heston_price_and_greeks), or forward-mode automatic differentiation (heston_greeks_ad) | | Bates (1996) | Heston CF × Merton (1976) log-normal jump CF | Bump-and-reprice (bates_price_and_greeks), or forward-mode AD (bates_greeks_ad) | | Local Vol (Dupire 1994) | Fritsch-Butland monotone cubic spline, differentiated through the spline, not the raw grid | Numerical (local vol surface) | | Monte Carlo (Heston/Bates) | Full truncation Euler (default) or Andersen (2008) QE, exact per-step Poisson jump counts, antithetic variates | N/A, path-dependent payoffs only (European, Asian, up-and-out barrier) |

All five analytic models share the same OptionContract/PricingResult conventions where applicable, so switching models in a caller doesn't mean rewriting the call site.

Design

Dispatch. No Box<dyn Model> anywhere in a pricing path. Every model is a free function, statically dispatched and monomorphized. Model selection happens at the call site, not through a trait object indirection that shows up in a profiler.

Memory layout. LocalVolSurface stores local_vols as a flat Vec<f64> indexed i_strike * n_expiry + j_expiry, not Vec<Vec<f64>>. One allocation, contiguous, cache-friendly for the row/column sweeps the Dupire finite differences need.

Characteristic function stability. The Heston and Bates CF use the Albrecher et al. (2007) formulation, which removes the branch-cut discontinuities of the original 1993 formula. stable_cf and the adaptive Gauss-Kronrod integrator (gk_integrate) live in heston.rs and are shared, unmodified, by bates.rs and ad.rs. The quadrature is adaptive (substitution to a finite interval, panel subdivision until the Kronrod/Gauss error estimate is below tolerance), which matters most in the wings and at short expiries, where a fixed panel under-resolves the integrand and silently produces prices that violate static arbitrage bounds. stable_cf's one sqrt() call per evaluation uses fast_csqrt, a closed-form algebraic complex square root, instead of num_complex::Complex64::sqrt()'s general branch, which goes through to_polar()/from_polar() (hypot + atan2 + sqrt + cos + sin, verified by reading the crate source). fast_csqrt needs hypot + 2 sqrt + a sign, no trig. In isolation this is ~5x faster (measured), but stable_cf also calls exp() and ln() and does several complex multiply/divides, so the end-to-end win on a full stable_cf call is a real but modest ~4%, not 5x, isolated micro-benchmarks of one operation don't linearly predict aggregate impact when a CPU can overlap independent work. Kept anyway: verified correct (680+ point sweep against the builtin, concentrated near the axes where a naive version of this formula catastrophically cancels, see Testing), strictly no downside, and it's the same fix ad.rs's dual csqrt needed for the same reason.

Local vol. dvar_dk, dvar_dt, and d2var_dk2 differentiate the Fritsch-Butland monotone cubic spline (Fritsch & Butland, 1984), not the raw IV grid. Central differences on raw quotes amplify quote noise into negative local variances; sampling the spline at a symmetric offset around each node means the F-B overshoot limiter is actually load-bearing in the derivative, not just present for interpolation queries between nodes. Boundary nodes use a one-sided offset within the grid range, a symmetric step at the edge would sample outside [strikes[0], strikes[n-1]] and hit the spline's flat clamp, which silently halves the estimated slope.

Surface repair. check_and_repair_surface is multi-pass: it loops fixing calendar-spread and butterfly violations until the surface is clean or a pass cap is hit. Fixing one violation can create another next to it (bumping an IV to kill a calendar violation can turn a previously-fine butterfly into a violation), so a single pass is not sufficient on a surface with more than one problem.

Monte Carlo. mc_heston/mc_bates default to full truncation Euler (Lord, Koekkoek, van Dijk 2010) for the variance process, correlated via the standard two-normal construction. McConfig::scheme = VarianceScheme::QuadraticExponential switches to Andersen's (2008) QE scheme: samples v(t+Δ) from a moment-matched distribution (quadratic-in-normal when ψ≤1.5, an exponential/point-mass mixture above) instead of discretizing and truncating the CIR SDE, and prices using his equation (33) for the log-price update (K0-K4 coefficients, central discretization γ1=γ2=0.5), verified against the primary source PDF, not a secondary writeup. The price innovation in (33) uses an independent normal, correlation with V is already analytic in K1/K2, reusing Euler's rho-correlated shock there was the bug in the first version of this, it gave a worse price than plain Euler (off by 6.2 vs analytic, Euler was off by 1.8) until fixed against the actual paper. Measured, not assumed: in a badly Feller-violating case (2κθ=0.4 ≪ σ²=1.44) with a coarse 8-step/year grid, QE cuts the bias against the analytic price roughly 20x versus Euler (qe_reduces_bias_in_feller_violating_regime). Bates jumps use an exact per-step Poisson draw (Knuth's algorithm), not the "coin flip with probability λdt" shortcut that silently drops the probability of two or more jumps landing in the same step. Parallelized over path chunks via rayon, with a splitmix64-hashed seed per chunk so parallel runs are reproducible and statistically independent. Use this for path-dependent payoffs the CF-inversion pricers can't touch (Asian averages, barriers), not for vanillas, heston_price/bates_price are exact and far cheaper for those.

Automatic differentiation. heston_greeks_ad and bates_greeks_ad propagate dual numbers (Dual<f64>, Complex<Dual>) through the characteristic function and integrate the derivative alongside the value via the Leibniz rule, giving exact vega/vanna without a finite-difference bump size to tune. The Bates version composes the Heston CF and the Merton jump CF before differentiating (same order bates_call uses), so price and the five Heston-driven Greeks are exact through the jump-adjusted CF. Jump-parameter sensitivities (d/dλ, d/dμⱼ, d/dσⱼ) are a separate function, ad::bates_jump_sensitivities_ad, same forward_pass machinery with the Heston side pinned constant and a jump parameter carrying the active derivative instead, it's what calibrate_bates's Jacobian uses for its jump columns now instead of FD.

heston_greeks_ad5 is a second, experimental AD path: Dual5 carries all 5 Heston-parameter tangent directions at once (dot: [f64; 5] instead of f64) instead of running 5 separate scalar-Dual passes, so the CF's value gets computed once instead of five times redundantly. Measured (profile_dual5_vs_five_scalar_passes, #[ignore]d, same methodology as the other profile tests): the joint pass is a real ~3x faster than 5 scalar passes at the integration level (211µs vs 623µs on this box). But heston_greeks_ad/ad5 spend most of their wall-clock time in the FD-bumped delta/gamma/theta/rho/vanna/volga, identical between both versions, so the win at the full-function level is a much more modest ~6%, and heston_greeks_ad5 is still slower than bump-and-reprice overall. The 3x is real and would matter if delta/gamma/rho/vanna/volga also moved onto AD (a further, larger Dual5-style tangent space covering spot and rate too), that's the natural next step this result points to, not implemented here.

Profiled properly (see ad::tests::profile_*, #[ignore]d, run with cargo test --release -- --ignored --nocapture --test-threads=1 ad::tests::profile), not guessed at: Complex<Dual> multiply is ~5.5x a plain Complex64 multiply, divide is ~11x, both measured with varying inputs cycled through the benchmark so LLVM can't hoist a closure-captured constant out of the loop and report a fake sub-nanosecond number (the first version of this benchmark did exactly that for sqrt/ln, caught by the results being physically impossible, not by inspection). But exp is only ~1.3x and ln ~1.1x, dual bookkeeping for those reuses the already-computed value via the standard AD reuse trick instead of redoing the whole computation twice. A full CF evaluation is dominated by the cheap transcendentals, not by raw multiply/divide count, so the aggregate overhead lands around 1.4-1.5x per GK panel and per full pricing pass, nowhere near the 11x the division number alone would suggest. That 1.4-1.5x is why heston_greeks_ad/bates_greeks_ad are exact but not currently faster wall-clock than bump-and-re