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

nautilus_model/data/
black_scholes.rs
