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

nautilus_model/data/
black_scholes.rs
