rustyqlib/core/aad/
mod.rs

1//! Adjoint Algorithmic Differentiation (AAD): tape-based reverse-mode
2//! differentiation for pricing code.
3//!
4//! Write a pricer over [`Var`] instead of `f64` — the operator
5//! overloading records every operation on a [`Tape`] — and one backward
6//! sweep ([`Var::grad`]) returns the sensitivity of the output to
7//! **every** input at once, at a fixed small multiple of the pricing
8//! cost. That is the AAD trade against bump-and-reprice: bumping costs
9//! one full reprice *per input*, the adjoint sweep costs ~one reprice
10//! *total*, however many inputs there are — the difference between
11//! seconds and hours on a book with per-pillar curve and surface
12//! sensitivities.
13//!
14//! Two worked and tested quant applications live in this module's tests:
15//!
16//! - [`black_scholes`]: the closed form written over `Var`; a single
17//!   sweep produces delta, dual delta, rho, carry rho, vega and the
18//!   maturity sensitivity simultaneously, matching the library's
19//!   closed-form Greeks to near machine precision;
20//! - **pathwise Monte Carlo Greeks**: differentiate through the
21//!   simulation itself (one small tape per path), giving delta, vega
22//!   and rho from the same paths that price — the standard pathwise
23//!   estimator, validated against the closed forms.
24//!
25//! The `max(x, 0)` payoff kink is handled by the almost-everywhere
26//! derivative ([`Var::maxf`]), which is exactly the classical pathwise
27//! estimator's requirement (fine for vanillas and smooth-density
28//! payoffs; digitals need smoothing or likelihood-ratio methods).
29
30pub mod tape;
31pub mod var;
32
33pub use tape::{Gradients, Tape};
34pub use var::Var;
35
36use crate::core::trade::PutOrCall;
37
38/// Black-Scholes price recorded on the tape: differentiate to get every
39/// first-order Greek from one backward sweep.
40pub fn black_scholes<'a>(
41    s: Var<'a>,
42    k: Var<'a>,
43    r: Var<'a>,
44    q: Var<'a>,
45    sigma: Var<'a>,
46    t: Var<'a>,
47    put_or_call: PutOrCall,
48) -> Var<'a> {
49    let sqrt_t = t.sqrt();
50    let st = sigma * sqrt_t;
51    let d1 = ((s / k).ln() + (r - q + sigma * sigma * 0.5) * t) / st;
52    let d2 = d1 - st;
53    let df_q = (-(q * t)).exp();
54    let df_r = (-(r * t)).exp();
55    match put_or_call {
56        PutOrCall::Call => s * df_q * d1.norm_cdf() - k * df_r * d2.norm_cdf(),
57        PutOrCall::Put => k * df_r * (-d2).norm_cdf() - s * df_q * (-d1).norm_cdf(),
58    }
59}
60
61#[cfg(test)]
62mod tests {
63    use super::*;
64    use crate::equity::blackscholes::{bs_price, bs_vega};
65
66    const S: f64 = 100.0;
67    const K: f64 = 105.0;
68    const R: f64 = 0.05;
69    const Q: f64 = 0.02;
70    const SIG: f64 = 0.3;
71    const T: f64 = 1.0;
72
73    #[test]
74    fn one_sweep_reproduces_every_black_scholes_greek() {
75        let tape = Tape::new();
76        let (s, k, r, q, sigma, t) =
77            (tape.var(S), tape.var(K), tape.var(R), tape.var(Q), tape.var(SIG), tape.var(T));
78        let price = black_scholes(s, k, r, q, sigma, t, PutOrCall::Call);
79        assert!((price.value() - bs_price(S, K, R, Q, SIG, T, PutOrCall::Call)).abs() < 1e-12);
80
81        // ONE backward pass: six sensitivities
82        let g = price.grad();
83
84        // vega against the library closed form
85        assert!((g.wrt(sigma) - bs_vega(S, K, R, Q, SIG, T)).abs() < 1e-10, "vega");
86        // the rest against tight central differences of bs_price
87        let h = 1e-6;
88        let fd = |f: &dyn Fn(f64) -> f64| (f(h) - f(-h)) / (2.0 * h);
89        let cases: [(f64, Box<dyn Fn(f64) -> f64>); 5] = [
90            (g.wrt(s), Box::new(|e| bs_price(S + e, K, R, Q, SIG, T, PutOrCall::Call))),
91            (g.wrt(k), Box::new(|e| bs_price(S, K + e, R, Q, SIG, T, PutOrCall::Call))),
92            (g.wrt(r), Box::new(|e| bs_price(S, K, R + e, Q, SIG, T, PutOrCall::Call))),
93            (g.wrt(q), Box::new(|e| bs_price(S, K, R, Q + e, SIG, T, PutOrCall::Call))),
94            (g.wrt(t), Box::new(|e| bs_price(S, K, R, Q, SIG, T + e, PutOrCall::Call))),
95        ];
96        for (i, (aad, f)) in cases.iter().enumerate() {
97            let numeric = fd(f);
98            assert!((aad - numeric).abs() < 1e-7, "greek {i}: aad {aad} vs fd {numeric}");
99        }
100        // put side too
101        let put = black_scholes(s, k, r, q, sigma, t, PutOrCall::Put);
102        let gp = put.grad();
103        let put_delta_fd =
104            (bs_price(S + h, K, R, Q, SIG, T, PutOrCall::Put)
105                - bs_price(S - h, K, R, Q, SIG, T, PutOrCall::Put))
106                / (2.0 * h);
107        assert!((gp.wrt(s) - put_delta_fd).abs() < 1e-7);
108    }
109
110    #[test]
111    fn pathwise_monte_carlo_greeks_from_one_sweep_per_path() {
112        // differentiate straight through a GBM simulation: delta, vega
113        // and rho of a European call from the same paths that price it
114        use crate::core::montecarlo::path_rng;
115        use rand::Rng;
116
117        let n_paths = 40_000;
118        // accumulate mean and variance of each estimator so the
119        // assertions are proper statistical bands, not magic numbers
120        let mut acc = [[0.0f64; 2]; 4]; // [sum, sum_sq] x {price, delta, vega, rho}
121        for i in 0..n_paths {
122            let mut rng = path_rng(2026, i);
123            let z: f64 = rng.sample(rand_distr::StandardNormal);
124            let tape = Tape::new();
125            let s0 = tape.var(S);
126            let sigma = tape.var(SIG);
127            let r = tape.var(R);
128            let drift = (r - Q - sigma * sigma * 0.5) * T;
129            let s_t = s0 * (drift + sigma * (T.sqrt() * z)).exp();
130            let payoff = (s_t - K).maxf(0.0) * (-(r * T)).exp();
131            let g = payoff.grad();
132            for (slot, x) in
133                [payoff.value(), g.wrt(s0), g.wrt(sigma), g.wrt(r)].into_iter().enumerate()
134            {
135                acc[slot][0] += x;
136                acc[slot][1] += x * x;
137            }
138        }
139        let n = n_paths as f64;
140        let stats = |slot: usize| -> (f64, f64) {
141            let mean = acc[slot][0] / n;
142            let var = (acc[slot][1] / n - mean * mean).max(0.0);
143            (mean, (var / n).sqrt())
144        };
145        let h = 1e-5;
146        let bs = |s: f64, sig: f64, r: f64| bs_price(s, K, r, Q, sig, T, PutOrCall::Call);
147        let truths = [
148            bs(S, SIG, R),
149            (bs(S + h, SIG, R) - bs(S - h, SIG, R)) / (2.0 * h),
150            bs_vega(S, K, R, Q, SIG, T),
151            (bs(S, SIG, R + h) - bs(S, SIG, R - h)) / (2.0 * h),
152        ];
153        for (slot, name) in ["price", "delta", "vega", "rho"].iter().enumerate() {
154            let (mean, se) = stats(slot);
155            assert!(
156                (mean - truths[slot]).abs() < 4.0 * se + 1e-10,
157                "{name}: {mean} vs {} (se {se})",
158                truths[slot]
159            );
160        }
161    }
162
163    #[test]
164    fn tape_cost_is_a_small_constant_multiple_of_pricing() {
165        // the adjoint promise: node count (a proxy for work) does not
166        // grow with the number of sensitivities requested
167        let tape = Tape::new();
168        let (s, k, r, q, sigma, t) =
169            (tape.var(S), tape.var(K), tape.var(R), tape.var(Q), tape.var(SIG), tape.var(T));
170        let price = black_scholes(s, k, r, q, sigma, t, PutOrCall::Call);
171        let nodes = tape.len();
172        assert!(nodes < 60, "tape has {nodes} nodes");
173        // one sweep serves all six inputs
174        let g = price.grad();
175        let six = [g.wrt(s), g.wrt(k), g.wrt(r), g.wrt(q), g.wrt(sigma), g.wrt(t)];
176        assert!(six.iter().all(|x| x.is_finite() && *x != 0.0));
177    }
178}
rustyqlib/core/aad/mod.rs

rustyqlib/core/aad/
mod.rs
