nautilus_model/data/
black_scholes.rs

1// -------------------------------------------------------------------------------------------------
2//  Copyright (C) 2015-2026 Nautech Systems Pty Ltd. All rights reserved.
3//  https://nautechsystems.io
4//
5//  Licensed under the GNU Lesser General Public License Version 3.0 (the "License");
6//  You may not use this file except in compliance with the License.
7//  You may obtain a copy of the License at https://www.gnu.org/licenses/lgpl-3.0.en.html
8//
9//  Unless required by applicable law or agreed to in writing, software
10//  distributed under the License is distributed on an "AS IS" BASIS,
11//  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12//  See the License for the specific language governing permissions and
13//  limitations under the License.
14// -------------------------------------------------------------------------------------------------
15
16pub trait BlackScholesReal:
17    Sized
18    + Copy
19    + Send
20    + Sync
21    + Default
22    + std::ops::Add<Output = Self>
23    + std::ops::Sub<Output = Self>
24    + std::ops::Mul<Output = Self>
25    + std::ops::Div<Output = Self>
26    + std::ops::Neg<Output = Self>
27{
28    type Mask: Copy;
29    fn splat(val: f64) -> Self;
30    #[must_use]
31    fn abs(self) -> Self;
32    #[must_use]
33    fn sqrt(self) -> Self;
34    #[must_use]
35    fn ln(self) -> Self;
36    #[must_use]
37    fn exp(self) -> Self;
38    #[must_use]
39    fn cdf(self) -> Self;
40    fn cdf_with_pdf(self) -> (Self, Self);
41    #[must_use]
42    fn mul_add(self, a: Self, b: Self) -> Self;
43    #[must_use]
44    fn recip_precise(self) -> Self;
45    fn select(mask: Self::Mask, t: Self, f: Self) -> Self;
46    fn cmp_gt(self, other: Self) -> Self::Mask;
47    #[must_use]
48    fn max(self, other: Self) -> Self;
49    #[must_use]
50    fn min(self, other: Self) -> Self;
51    #[must_use]
52    fn signum(self) -> Self;
53}
54
55impl BlackScholesReal for f32 {
56    type Mask = bool;
57    #[inline(always)]
58    fn splat(val: f64) -> Self {
59        val as Self
60    }
61    #[inline(always)]
62    fn abs(self) -> Self {
63        self.abs()
64    }
65    #[inline(always)]
66    fn sqrt(self) -> Self {
67        self.sqrt()
68    }
69    #[inline(always)]
70    fn select(mask: bool, t: Self, f: Self) -> Self {
71        if mask { t } else { f }
72    }
73    #[inline(always)]
74    fn cmp_gt(self, other: Self) -> bool {
75        self > other
76    }
77    #[inline(always)]
78    fn recip_precise(self) -> Self {
79        1.0 / self
80    }
81    #[inline(always)]
82    fn max(self, other: Self) -> Self {
83        self.max(other)
84    }
85    #[inline(always)]
86    fn min(self, other: Self) -> Self {
87        self.min(other)
88    }
89    #[inline(always)]
90    fn signum(self) -> Self {
91        self.signum()
92    }
93    #[inline(always)]
94    fn mul_add(self, a: Self, b: Self) -> Self {
95        self.mul_add(a, b)
96    }
97
98    #[inline(always)]
99    fn ln(self) -> Self {
100        // Minimax polynomial approximation for ln(x) on [1, 2)
101        // Optimized for f32 precision with max error ~1e-7
102        // Uses range reduction: ln(mantissa) = ln(1 + x) where x = (mantissa - 1) / (mantissa + 1)
103        // See: J.-M. Muller et al., "Handbook of Floating-Point Arithmetic", 2018, Section 10.2
104        //      A. J. Salgado & S. M. Wise, "Classical Numerical Analysis", 2023, Chapter 10
105        let bits = self.to_bits();
106
107        // Positive normal bit patterns occupy one contiguous interval, so the ordinary
108        // path reaches the polynomial through a single range test. Everything else -
109        // zeros, negatives, subnormals, infinity, NaN - is handled out of line.
110        if !(0x0080_0000..0x7f80_0000).contains(&bits) {
111            return ln_f32_outside_normal_range(self);
112        }
113
114        let exponent = ((bits >> 23) as i32 - 127) as Self;
115        let mantissa = Self::from_bits((bits & 0x007F_FFFF) | 0x3f80_0000);
116        let x = (mantissa - 1.0) / (mantissa + 1.0);
117        let x2 = x * x;
118        let mut res = 0.239_282_85_f32;
119        res = x2.mul_add(res, 0.285_182_11);
120        res = x2.mul_add(res, 0.400_005_83);
121        res = x2.mul_add(res, 0.666_666_7);
122        res = x2.mul_add(res, 2.0);
123        x.mul_add(res, exponent * std::f32::consts::LN_2)
124    }
125
126    #[inline(always)]
127    fn exp(self) -> Self {
128        // Minimax polynomial approximation for exp(x) on [-0.5*ln(2), 0.5*ln(2))
129        // Optimized for f32 precision with max error ~1e-7
130        // Uses range reduction: exp(x) = 2^k * exp(r) where k = round(x / ln(2)) and r = x - k*ln(2)
131        // See: J.-M. Muller et al., "Handbook of Floating-Point Arithmetic", 2018, Section 10.3
132        //      A. J. Salgado & S. M. Wise, "Classical Numerical Analysis", 2023, Chapter 10
133        let k = (self.mul_add(
134            std::f32::consts::LOG2_E,
135            if self > 0.0 { 0.5 } else { -0.5 },
136        )) as i32;
137
138        if k <= -151 {
139            return 0.0;
140        }
141
142        if k >= 129 {
143            return Self::INFINITY;
144        }
145
146        let r = self - (k as Self * 0.693_145_75) - (k as Self * 1.428_606_8e-6);
147        let mut res = 0.001_388_89_f32;
148        res = r.mul_add(res, 0.008_333_33);
149        res = r.mul_add(res, 0.041_666_67);
150        res = r.mul_add(res, 0.166_666_67);
151        res = r.mul_add(res, 0.5);
152        res = r.mul_add(res, 1.0);
153        if (-126..=127).contains(&k) {
154            r.mul_add(res, 1.0) * Self::from_bits(((k + 127) as u32) << 23)
155        } else {
156            // Split 2^k across two valid biased exponents; a single one would be
157            // out of range here. Apply them in sequence: pre-multiplying the two
158            // factors underflows to zero at k = -150 before the polynomial can
159            // round the result up into the subnormal range.
160            let ka = k >> 1;
161            let kb = k - ka;
162            r.mul_add(res, 1.0)
163                * Self::from_bits(((ka + 127) as u32) << 23)
164                * Self::from_bits(((kb + 127) as u32) << 23)
165        }
166    }
167
168    #[inline(always)]
169    fn cdf(self) -> Self {
170        self.cdf_with_pdf().0
171    }
172
173    #[inline(always)]
174    fn cdf_with_pdf(self) -> (Self, Self) {
175        // Minimax rational approximation for normal CDF
176        // Optimized for f32 precision with max error ~1e-7
177        // Uses transformation t = 1 / (1 + 0.2316419 * |x|) for numerical stability
178        // See: M. Abramowitz & I. A. Stegun (eds.), "Handbook of Mathematical Functions
179        //      with Formulas, Graphs, and Mathematical Tables", 1972, Section 26.2.17
180        let abs_x = self.abs();
181        let t = 1.0 / (1.0 + 0.231_641_9 * abs_x);
182        let mut poly = 1.330_274_5_f32.mul_add(t, -1.821_255_9);
183        poly = t.mul_add(poly, 1.781_477_9);
184        poly = t.mul_add(poly, -0.356_563_78);
185        poly = t.mul_add(poly, 0.319_381_54);
186        let pdf = 0.398_942_3 * (-0.5 * self * self).exp();
187        let res = 1.0 - pdf * (poly * t);
188        // Use >= to ensure CDF(0) = 0.5 exactly (maintains symmetry)
189        (if self >= 0.0 { res } else { 1.0 - res }, pdf)
190    }
191}
192
193#[derive(Debug, Clone, Copy, Default, PartialEq)]
194pub struct Greeks<T> {
195    pub price: T,
196    pub vol: T,
197    pub delta: T,
198    pub gamma: T,
199    pub vega: T,
200    pub theta: T,
201    pub itm_prob: T,
202}
203
204/// Lightweight kernel for IV search - only computes price and vega.
205/// `phi` is +1 for call, -1 for put (caller does the select once).
206#[inline(always)]
207fn pricing_kernel_price_vega<T: BlackScholesReal>(
208    s_forward: T,
209    k: T,
210    df_r: T,
211    d1: T,
212    d2: T,
213    sqrt_t: T,
214    phi: T,
215) -> (T, T) {
216    let (cdf_phi_d1, pdf_d1) = (phi * d1).cdf_with_pdf();
217    let cdf_phi_d2 = (phi * d2).cdf();
218
219    let price = phi * (s_forward * cdf_phi_d1 - k * df_r * cdf_phi_d2);
220    let vega = s_forward * sqrt_t * pdf_d1;
221
222    (price, vega)
223}
224
225#[expect(clippy::too_many_arguments)]
226#[inline(always)]
227fn pricing_kernel<T: BlackScholesReal>(
228    s_forward: T,
229    k: T,
230    df_r: T,
231    d1: T,
232    d2: T,
233    inv_scaled_vol: T,
234    vol: T,
235    sqrt_t: T,
236    t: T,
237    r: T,
238    b: T,
239    s: T,
240    phi: T,
241) -> Greeks<T> {
242    let (cdf_phi_d1, pdf_d1) = (phi * d1).cdf_with_pdf();
243    let cdf_phi_d2 = (phi * d2).cdf();
244
245    let df_b = ((b - r) * t).exp();
246    let price = phi * (s_forward * cdf_phi_d1 - k * df_r * cdf_phi_d2);
247    let delta = phi * df_b * cdf_phi_d1;
248    let vega = s_forward * sqrt_t * pdf_d1;
249    let gamma = df_b * pdf_d1 * inv_scaled_vol / s;
250
251    let theta_v = -(s_forward * pdf_d1 * vol) * (T::splat(2.0) * sqrt_t).recip_precise();
252    let theta_b = -phi * (b - r) * s_forward * cdf_phi_d1;
253    let theta_r = -phi * r * k * df_r * cdf_phi_d2;
254    let theta = theta_v + theta_b + theta_r;
255
256    Greeks {
257        price,
258        vol,
259        delta,
260        gamma,
261        vega,
262        theta,
263        itm_prob: cdf_phi_d2,
264    }
265}
266
267#[inline(always)]
268pub fn compute_greeks<T: BlackScholesReal>(
269    s: T,
270    k: T,
271    t: T,
272    r: T,
273    b: T,
274    vol: T,
275    is_call: T::Mask,
276) -> Greeks<T> {
277    let sqrt_t = t.sqrt();
278    let scaled_vol = vol * sqrt_t;
279    let inv_scaled_vol = scaled_vol.recip_precise();
280    let df_r = (-r * t).exp();
281    let df_b = ((b - r) * t).exp();
282    let d1 = ((s / k).ln() + (b + T::splat(0.5) * vol * vol) * t) * inv_scaled_vol;
283    let d2 = d1 - scaled_vol;
284    let s_forward = s * df_b;
285    let phi = T::select(is_call, T::splat(1.0), T::splat(-1.0));
286
287    pricing_kernel(
288        s_forward,
289        k,
290        df_r,
291        d1,
292        d2,
293        inv_scaled_vol,
294        vol,
295        sqrt_t,
296        t,
297        r,
298        b,
299        s,
300        phi,
301    )
302}
303
304/// Performs a single Halley iteration to refine an implied volatility estimate and compute greeks.
305///
306/// # Important Notes
307///
308/// This function is intended as a **refinement step** when a good initial guess for volatility
309/// is available (e.g., from a previous calculation or a fast approximation). It performs only
310/// a single Halley iteration and does not implement a full convergence loop.
311///
312/// **This is NOT a standalone implied volatility solver.** For production use, prefer
313/// `imply_vol_and_greeks` which uses the robust `implied_vol` crate for full convergence.
314///
315/// # Parameters
316///
317/// - `initial_guess`: Must be a reasonable estimate of the true volatility. Poor initial guesses
318///   (especially for deep ITM/OTM options) may result in significant errors.
319///
320/// # Accuracy
321///
322/// With a good initial guess (within ~25% of true vol), one Halley step typically achieves
323/// ~1% relative error. For deep ITM/OTM options or poor initial guesses, multiple iterations
324/// or a better initial estimate may be required.
325#[expect(clippy::too_many_arguments)]
326#[inline(always)]
327pub fn compute_iv_and_greeks<T: BlackScholesReal>(
328    mkt_price: T,
329    s: T,
330    k: T,
331    t: T,
332    r: T,
333    b: T,
334    is_call: T::Mask,
335    initial_guess: T,
336) -> Greeks<T> {
337    let sqrt_t = t.sqrt();
338    let inv_sqrt_t = sqrt_t.recip_precise();
339    let ln_sk_bt = (s.ln() - k.ln()) + (b * t);
340    let half_t = T::splat(0.5) * t;
341    let df_r = (-r * t).exp();
342    let vol = initial_guess;
343
344    let inv_vol = vol.recip_precise();
345    let inv_scaled_vol = inv_vol * inv_sqrt_t;
346    let d1 = (ln_sk_bt + half_t * vol * vol) * inv_scaled_vol;
347    let d2 = d1 - vol * sqrt_t;
348    let s_forward = s * ((b - r) * t).exp();
349    let phi = T::select(is_call, T::splat(1.0), T::splat(-1.0));
350    let (price, vega_raw) = pricing_kernel_price_vega(s_forward, k, df_r, d1, d2, sqrt_t, phi);
351
352    let diff = price - mkt_price;
353    let vega = vega_raw.abs().max(T::splat(1e-9));
354    let volga = (vega * d1 * d2) * inv_vol;
355    let num = T::splat(2.0) * diff * vega;
356    let den = T::splat(2.0) * vega * vega - diff * volga;
357
358    // Clamp denominator magnitude while preserving sign
359    let den_safe = den.signum() * den.abs().max(T::splat(1e-9));
360
361    // Clamp volatility between 0.0001% and 1000% annualized to keep outputs finite
362    let vol = (vol - (num * den_safe.recip_precise()))
363        .max(T::splat(1e-6))
364        .min(T::splat(10.0));
365
366    let inv_vol_f = vol.recip_precise();
367    let inv_scaled_vol_f = inv_vol_f * inv_sqrt_t;
368    let scaled_vol_f = vol * sqrt_t;
369    let d1_f = (ln_sk_bt + half_t * vol * vol) * inv_scaled_vol_f;
370    let d2_f = d1_f - scaled_vol_f;
371    pricing_kernel(
372        s_forward,
373        k,
374        df_r,
375        d1_f,
376        d2_f,
377        inv_scaled_vol_f,
378        vol,
379        sqrt_t,
380        t,
381        r,
382        b,
383        s,
384        phi,
385    )
386}
387
388/// Returns `ln(value)` for the inputs outside the positive normal range: zeros, negatives,
389/// subnormals, infinity, and NaN.
390///
391/// Kept out of line because `<f32 as BlackScholesReal>::ln` is `#[inline(always)]`, so any
392/// classification left in its body is duplicated at every call site on the pricing path.
393#[cold]
394#[inline(never)]
395fn ln_f32_outside_normal_range(value: f32) -> f32 {
396    let bits = value.to_bits();
397    let magnitude = bits & 0x7fff_ffff;
398    if magnitude == 0 {
399        return f32::NEG_INFINITY;
400    }
401
402    let exponent_bits = magnitude & 0x7f80_0000;
403    let fraction_bits = magnitude & 0x007f_ffff;
404    if exponent_bits == 0x7f80_0000 && fraction_bits != 0 {
405        return f32::NAN;
406    }
407
408    if bits & 0x8000_0000 != 0 {
409        return f32::NAN;
410    }
411
412    if exponent_bits == 0x7f80_0000 {
413        return f32::INFINITY;
414    }
415
416    // Only positive subnormals remain. Scaling by 2^23 is exact and carries every one of
417    // them into the normal range, so this re-enters the ordinary path exactly once.
418    <f32 as BlackScholesReal>::ln(value * 8_388_608.0) - 23.0 * std::f32::consts::LN_2
419}
420
421#[cfg(test)]
422mod tests {
423    use rstest::*;
424
425    use super::*;
426    use crate::data::greeks::black_scholes_greeks_exact;
427
428    /// The positive-normal `ln` formula exactly as it stood before the domain guard.
429    /// Frozen deliberately: it is the oracle for the byte-identity characterization
430    /// below, so it must not be updated alongside the implementation.
431    fn old_ln_formula(input: f32) -> f32 {
432        let bits = input.to_bits();
433        let exponent = ((bits >> 23) as i32 - 127) as f32;
434        let mantissa = f32::from_bits((bits & 0x007F_FFFF) | 0x3f80_0000);
435        let x = (mantissa - 1.0) / (mantissa + 1.0);
436        let x2 = x * x;
437        let mut res = 0.239_282_85_f32;
438        res = x2.mul_add(res, 0.285_182_11);
439        res = x2.mul_add(res, 0.400_005_83);
440        res = x2.mul_add(res, 0.666_666_7);
441        res = x2.mul_add(res, 2.0);
442        x.mul_add(res, exponent * std::f32::consts::LN_2)
443    }
444
445    fn assert_ln_close(actual: f32, expected: f32, max_ulps: u32) {
446        assert_eq!(actual.is_nan(), expected.is_nan());
447        assert_eq!(actual.is_infinite(), expected.is_infinite());
448        assert_eq!(actual.is_sign_negative(), expected.is_sign_negative());
449        if actual.is_finite() {
450            assert!(
451                actual.to_bits().abs_diff(expected.to_bits()) <= max_ulps,
452                "ln mismatch: actual={actual:e} ({:#010x}), expected={expected:e} ({:#010x})",
453                actual.to_bits(),
454                expected.to_bits(),
455            );
456        }
457    }
458
459    fn assert_exp_close(actual: f32, expected: f32, max_ulps: u32) {
460        assert_eq!(actual.is_nan(), expected.is_nan());
461        assert_eq!(actual.is_infinite(), expected.is_infinite());
462        assert_eq!(actual.is_sign_negative(), expected.is_sign_negative());
463        if actual.is_finite() {
464            assert!(
465                actual.to_bits().abs_diff(expected.to_bits()) <= max_ulps,
466                "exp mismatch: actual={actual:e} ({:#010x}), expected={expected:e} ({:#010x})",
467                actual.to_bits(),
468                expected.to_bits(),
469            );
470        }
471    }
472
473    #[rstest]
474    fn test_ln_special_values() {
475        for input in [0.0_f32, -0.0] {
476            let actual = <f32 as BlackScholesReal>::ln(input);
477            assert_eq!(actual, f32::NEG_INFINITY, "input={input:?}");
478        }
479
480        for input in [-1.0_f32, -2.0, -f32::MIN_POSITIVE, f32::NEG_INFINITY] {
481            assert!(
482                <f32 as BlackScholesReal>::ln(input).is_nan(),
483                "input={input:?}"
484            );
485        }
486
487        assert_eq!(<f32 as BlackScholesReal>::ln(f32::INFINITY), f32::INFINITY);
488
489        for input in [f32::from_bits(0x7fc0_1234), f32::from_bits(0xffc0_5678)] {
490            assert!(
491                <f32 as BlackScholesReal>::ln(input).is_nan(),
492                "input_bits={:#010x}",
493                input.to_bits()
494            );
495        }
496    }
497
498    #[rstest]
499    fn test_ln_positive_subnormals() {
500        for input in [
501            f32::from_bits(1),
502            f32::from_bits(0x0040_0000),
503            f32::from_bits(0x007f_ffff),
504        ] {
505            assert_ln_close(<f32 as BlackScholesReal>::ln(input), input.ln(), 3);
506        }
507    }
508
509    #[rstest]
510    fn test_ln_positive_normal_path_is_unchanged() {
511        // Every normal exponent field crossed with a spread of mantissa patterns.
512        // Exhausting the exponent matters because the final fused addition combines
513        // the polynomial with `exponent * LN_2`, so rounding can differ per exponent
514        // even when mantissa handling is untouched.
515        let fractions = [
516            0x0000_0000,
517            0x0000_0001,
518            0x001f_ffff,
519            0x003f_ffff,
520            0x0040_0000,
521            0x0055_5555,
522            0x007f_fffe,
523            0x007f_ffff,
524        ];
525
526        for exponent in 1_u32..=254 {
527            for fraction in fractions {
528                let input = f32::from_bits((exponent << 23) | fraction);
529                assert_eq!(
530                    <f32 as BlackScholesReal>::ln(input).to_bits(),
531                    old_ln_formula(input).to_bits(),
532                    "input={input:e} ({:#010x})",
533                    input.to_bits()
534                );
535            }
536        }
537    }
538
539    #[rstest]
540    fn test_compute_greeks_negative_strike_returns_nan_price() {
541        let greeks = compute_greeks::<f32>(100.0, -100.0, 1.0, 0.05, 0.05, 0.2, true);
542
543        assert!(greeks.price.is_nan());
544    }
545
546    #[rstest]
547    fn test_compute_iv_and_greeks_negative_strike_returns_nan_price() {
548        let greeks = compute_iv_and_greeks::<f32>(10.0, 100.0, -100.0, 1.0, 0.05, 0.05, true, 0.2);
549
550        assert!(greeks.price.is_nan());
551    }
552
553    #[rstest]
554    fn test_exp_exponent_boundaries() {
555        let ln_2 = std::f32::consts::LN_2;
556        let inputs = [
557            -126.5 * ln_2 - 0.000_1,
558            -126.5 * ln_2 + 0.000_1,
559            -127.5 * ln_2 - 0.000_1,
560            -127.5 * ln_2 + 0.000_1,
561            -150.0 * ln_2 + 0.1,
562            88.7,
563            88.8,
564            129.0 * ln_2,
565        ];
566
567        for input in inputs {
568            assert_exp_close(<f32 as BlackScholesReal>::exp(input), input.exp(), 3);
569        }
570    }
571
572    #[rstest]
573    fn test_exp_tail_sweep() {
574        let mut previous = 0.0;
575        let mut input = -105.0_f32;
576        while input <= 89.0 {
577            let actual = <f32 as BlackScholesReal>::exp(input);
578            assert!(!actual.is_nan(), "exp({input}) returned NaN");
579            assert!(actual >= 0.0, "exp({input}) returned {actual}");
580            assert!(
581                actual >= previous,
582                "exp is not monotonic at {input}: {actual} < {previous}"
583            );
584            assert_exp_close(actual, input.exp(), 3);
585            previous = actual;
586            input += 0.031_25;
587        }
588    }
589
590    #[rstest]
591    fn test_exp_extreme_inputs() {
592        for input in [-1e10_f32, -1e30, f32::NEG_INFINITY] {
593            assert_eq!(<f32 as BlackScholesReal>::exp(input), 0.0, "input={input}");
594        }
595
596        for input in [1e10_f32, 1e30, f32::INFINITY] {
597            assert_eq!(
598                <f32 as BlackScholesReal>::exp(input),
599                f32::INFINITY,
600                "input={input}"
601            );
602        }
603        assert!(<f32 as BlackScholesReal>::exp(f32::NAN).is_nan());
604    }
605
606    #[rstest]
607    fn test_accuracy_1e7() {
608        let s = 100.0;
609        let k = 100.0;
610        let t = 1.0;
611        let r = 0.05;
612        let vol = 0.2;
613        let g = compute_greeks::<f32>(s, k, t, r, r, vol, true); // Use r as b
614        assert!((g.price - 10.45058).abs() < 1e-5);
615        let solved = compute_iv_and_greeks::<f32>(g.price, s, k, t, r, r, true, vol); // Use r as b
616        assert!((solved.vol - vol).abs() < 1e-6);
617    }
618
619    #[rstest]
620    fn test_compute_greeks_accuracy_vs_exact() {
621        let s = 100.0f64;
622        let k = 100.0f64;
623        let t = 1.0f64;
624        let r = 0.05f64;
625        let b = 0.05f64; // cost of carry
626        let vol = 0.2f64;
627        let multiplier = 1.0f64;
628
629        // Compute using fast f32 method
630        let g_fast = compute_greeks::<f32>(
631            s as f32, k as f32, t as f32, r as f32, b as f32, vol as f32, true,
632        );
633
634        // Compute using exact f64 method
635        let g_exact = black_scholes_greeks_exact(s, r, b, vol, true, k, t);
636
637        // Compare with tolerance for f32 precision
638        let price_tol = 1e-4;
639        let greeks_tol = 1e-3;
640
641        assert!(
642            (f64::from(g_fast.price) - g_exact.price).abs() < price_tol,
643            "Price mismatch: fast={}, exact={}",
644            g_fast.price,
645            g_exact.price
646        );
647        assert!(
648            (f64::from(g_fast.delta) - g_exact.delta).abs() < greeks_tol,
649            "Delta mismatch: fast={}, exact={}",
650            g_fast.delta,
651            g_exact.delta
652        );
653        assert!(
654            (f64::from(g_fast.gamma) - g_exact.gamma).abs() < greeks_tol,
655            "Gamma mismatch: fast={}, exact={}",
656            g_fast.gamma,
657            g_exact.gamma
658        );
659        // Vega units differ: exact uses multiplier * 0.01, fast uses raw units
660        let vega_exact_raw = g_exact.vega / (multiplier * 0.01);
661        assert!(
662            (f64::from(g_fast.vega) - vega_exact_raw).abs() < greeks_tol,
663            "Vega mismatch: fast={}, exact_raw={}, exact_scaled={}",
664            g_fast.vega,
665            vega_exact_raw,
666            g_exact.vega
667        );
668        // Theta units differ: exact uses multiplier * daily_factor (0.0027378507871321013), fast uses raw units
669        let theta_daily_factor = 0.002_737_850_787_132_101_3;
670        let theta_exact_raw = g_exact.theta / (multiplier * theta_daily_factor);
671        assert!(
672            (f64::from(g_fast.theta) - theta_exact_raw).abs() < greeks_tol,
673            "Theta mismatch: fast={}, exact_raw={}, exact_scaled={}",
674            g_fast.theta,
675            theta_exact_raw,
676            g_exact.theta
677        );
678    }
679
680    #[rstest]
681    fn test_put_theta_with_cost_of_carry_not_equal_to_rate() {
682        let s = 100.0f64;
683        let k = 100.0f64;
684        let t = 1.0f64;
685        let r = 0.05f64;
686        let b = 0.0f64; // cost of carry != r (e.g. futures option)
687        let vol = 0.2f64;
688        let multiplier = 1.0f64;
689
690        let g_fast = compute_greeks::<f32>(
691            s as f32, k as f32, t as f32, r as f32, b as f32, vol as f32, false,
692        );
693
694        let g_exact = black_scholes_greeks_exact(s, r, b, vol, false, k, t);
695
696        let theta_daily_factor = 0.002_737_850_787_132_101_3;
697        let theta_exact_raw = g_exact.theta / (multiplier * theta_daily_factor);
698        assert!(
699            (f64::from(g_fast.theta) - theta_exact_raw).abs() < 1e-3,
700            "Put theta mismatch with b!=r: fast={}, exact_raw={}",
701            g_fast.theta,
702            theta_exact_raw
703        );
704    }
705
706    #[rstest]
707    fn test_compute_iv_and_greeks_halley_accuracy() {
708        let s = 100.0f64;
709        let k = 100.0f64;
710        let t = 1.0f64;
711        let r = 0.05f64;
712        let b = 0.05f64; // cost of carry
713        let vol_true = 0.2f64; // True volatility
714        let initial_guess = 0.25f64; // Initial guess (25% higher than true)
715        let multiplier = 1.0f64;
716
717        // Compute the exact price using the true volatility
718        let g_exact = black_scholes_greeks_exact(s, r, b, vol_true, true, k, t);
719        let mkt_price = g_exact.price;
720
721        // Compute implied vol using one Halley step with initial guess
722        let g_halley = compute_iv_and_greeks::<f32>(
723            mkt_price as f32,
724            s as f32,
725            k as f32,
726            t as f32,
727            r as f32,
728            b as f32,
729            true,
730            initial_guess as f32,
731        );
732
733        // Check that one Halley step gets close to the true volatility
734        let vol_error = (f64::from(g_halley.vol) - vol_true).abs();
735
736        // One Halley step should get within ~1% of true vol for a 25% initial error
737        assert!(
738            vol_error < 0.01,
739            "Halley step accuracy: vol_error={}, initial_guess={}, vol_true={}, computed_vol={}",
740            vol_error,
741            initial_guess,
742            vol_true,
743            g_halley.vol
744        );
745
746        // Check that the computed greeks are close to exact
747        let price_tol = 5e-3; // Relaxed for one Halley step
748        let greeks_tol = 5e-3; // Relaxed for one-step approximation
749
750        assert!(
751            (f64::from(g_halley.price) - g_exact.price).abs() < price_tol,
752            "Price mismatch after Halley: halley={}, exact={}, diff={}",
753            g_halley.price,
754            g_exact.price,
755            (f64::from(g_halley.price) - g_exact.price).abs()
756        );
757        assert!(
758            (f64::from(g_halley.delta) - g_exact.delta).abs() < greeks_tol,
759            "Delta mismatch after Halley: halley={}, exact={}",
760            g_halley.delta,
761            g_exact.delta
762        );
763        assert!(
764            (f64::from(g_halley.gamma) - g_exact.gamma).abs() < greeks_tol,
765            "Gamma mismatch after Halley: halley={}, exact={}",
766            g_halley.gamma,
767            g_exact.gamma
768        );
769        // Vega units differ: exact uses multiplier * 0.01, fast uses raw units
770        let vega_exact_raw = g_exact.vega / (multiplier * 0.01);
771        assert!(
772            (f64::from(g_halley.vega) - vega_exact_raw).abs() < greeks_tol,
773            "Vega mismatch after Halley: halley={}, exact_raw={}",
774            g_halley.vega,
775            vega_exact_raw
776        );
777        // Theta units differ: exact uses multiplier * daily_factor (0.0027378507871321013), fast uses raw units
778        let theta_daily_factor = 0.002_737_850_787_132_101_3;
779        let theta_exact_raw = g_exact.theta / (multiplier * theta_daily_factor);
780        assert!(
781            (f64::from(g_halley.theta) - theta_exact_raw).abs() < greeks_tol,
782            "Theta mismatch after Halley: halley={}, exact_raw={}",
783            g_halley.theta,
784            theta_exact_raw
785        );
786    }
787
788    #[rstest]
789    fn test_print_halley_iv() {
790        let s = 100.0f64;
791        let k = 100.0f64;
792        let t = 1.0f64;
793        let r = 0.05f64;
794        let b = 0.05f64;
795        let vol_true = 0.2f64;
796
797        let g_exact = black_scholes_greeks_exact(s, r, b, vol_true, true, k, t);
798        let mkt_price = g_exact.price;
799
800        println!("\n=== Halley Step IV Test (Using True Vol as Initial Guess) ===");
801        println!("True volatility: {vol_true}");
802        println!("Market price: {mkt_price:.8}");
803        println!("Initial guess: {vol_true} (using true vol)");
804
805        let g_halley = compute_iv_and_greeks::<f32>(
806            mkt_price as f32,
807            s as f32,
808            k as f32,
809            t as f32,
810            r as f32,
811            b as f32,
812            true,
813            vol_true as f32, // Using true vol as initial guess
814        );
815
816        println!("\nAfter one Halley step:");
817        println!("Computed volatility: {:.8}", g_halley.vol);
818        println!("True volatility: {vol_true:.8}");
819        println!(
820            "Absolute error: {:.8}",
821            (f64::from(g_halley.vol) - vol_true).abs()
822        );
823        println!(
824            "Relative error: {:.4}%",
825            (f64::from(g_halley.vol) - vol_true).abs() / vol_true * 100.0
826        );
827    }
828
829    #[rstest]
830    fn test_compute_iv_and_greeks_deep_itm_otm() {
831        let t = 1.0f64;
832        let r = 0.05f64;
833        let b = 0.05f64;
834        let vol_true = 0.2f64;
835
836        // Deep ITM: s=150, k=100 (spot is 50% above strike)
837        let s_itm = 150.0f64;
838        let k_itm = 100.0f64;
839        let g_exact_itm = black_scholes_greeks_exact(s_itm, r, b, vol_true, true, k_itm, t);
840        let mkt_price_itm = g_exact_itm.price;
841
842        println!("\n=== Deep ITM Test ===");
843        println!("Spot: {s_itm}, Strike: {k_itm}, True vol: {vol_true}");
844        println!("Market price: {mkt_price_itm:.8}");
845
846        let g_recovered_itm = compute_iv_and_greeks::<f32>(
847            mkt_price_itm as f32,
848            s_itm as f32,
849            k_itm as f32,
850            t as f32,
851            r as f32,
852            b as f32,
853            true,
854            vol_true as f32, // Using true vol as initial guess
855        );
856
857        let vol_error_itm = (f64::from(g_recovered_itm.vol) - vol_true).abs();
858        let rel_error_itm = vol_error_itm / vol_true * 100.0;
859
860        println!("Recovered volatility: {:.8}", g_recovered_itm.vol);
861        println!("Absolute error: {vol_error_itm:.8}");
862        println!("Relative error: {rel_error_itm:.4}%");
863
864        // Deep OTM: s=50, k=100 (spot is 50% below strike)
865        let s_otm = 50.0f64;
866        let k_otm = 100.0f64;
867        let g_exact_otm = black_scholes_greeks_exact(s_otm, r, b, vol_true, true, k_otm, t);
868        let mkt_price_otm = g_exact_otm.price;
869
870        println!("\n=== Deep OTM Test ===");
871        println!("Spot: {s_otm}, Strike: {k_otm}, True vol: {vol_true}");
872        println!("Market price: {mkt_price_otm:.8}");
873
874        let g_recovered_otm = compute_iv_and_greeks::<f32>(
875            mkt_price_otm as f32,
876            s_otm as f32,
877            k_otm as f32,
878            t as f32,
879            r as f32,
880            b as f32,
881            false,
882            vol_true as f32, // Using true vol as initial guess
883        );
884
885        let vol_error_otm = (f64::from(g_recovered_otm.vol) - vol_true).abs();
886        let rel_error_otm = vol_error_otm / vol_true * 100.0;
887
888        println!("Recovered volatility: {:.8}", g_recovered_otm.vol);
889        println!("Absolute error: {vol_error_otm:.8}");
890        println!("Relative error: {rel_error_otm:.4}%");
891
892        // Assertions: Deep ITM and OTM are challenging cases
893        // One Halley step with Corrado-Miller initial guess may not be sufficient
894        // We use a more relaxed tolerance to verify the method still converges in the right direction
895        // For production use, multiple iterations or better initial guesses would be needed
896        let vol_tol_itm = 50.0; // 50% relative error tolerance for deep ITM
897        let vol_tol_otm = 150.0; // 150% relative error tolerance for deep OTM (very challenging)
898
899        // Check that we at least get a reasonable volatility (not NaN or extreme values)
900        assert!(
901            g_recovered_itm.vol.is_finite()
902                && g_recovered_itm.vol > 0.0
903                && g_recovered_itm.vol < 2.0,
904            "Deep ITM vol recovery: invalid result={}",
905            g_recovered_itm.vol
906        );
907
908        assert!(
909            g_recovered_otm.vol.is_finite()
910                && g_recovered_otm.vol > 0.0
911                && g_recovered_otm.vol < 2.0,
912            "Deep OTM vol recovery: invalid result={}",
913            g_recovered_otm.vol
914        );
915
916        // Verify the error is within acceptable bounds (one step may not be enough)
917        assert!(
918            rel_error_itm < vol_tol_itm,
919            "Deep ITM vol recovery error too large: recovered={}, true={}, error={:.4}%",
920            g_recovered_itm.vol,
921            vol_true,
922            rel_error_itm
923        );
924
925        assert!(
926            rel_error_otm < vol_tol_otm,
927            "Deep OTM vol recovery error too large: recovered={}, true={}, error={:.4}%",
928            g_recovered_otm.vol,
929            vol_true,
930            rel_error_otm
931        );
932
933        println!("\n=== Summary ===");
934        println!("Deep ITM: One Halley iteration error = {rel_error_itm:.2}%");
935        println!(
936            "Deep OTM: One Halley iteration error = {rel_error_otm:.2}% (still challenging, deep OTM is difficult)"
937        );
938    }
939}
nautilus_model/data/black_scholes.rs

nautilus_model/data/
black_scholes.rs
