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        let exponent = ((bits >> 23) as i32 - 127) as Self;
109        let mantissa = Self::from_bits((bits & 0x007FFFFF) | 0x3f800000);
110        let x = (mantissa - 1.0) / (mantissa + 1.0);
111        let x2 = x * x;
112        let mut res = 0.23928285_f32;
113        res = x2.mul_add(res, 0.28518211);
114        res = x2.mul_add(res, 0.40000583);
115        res = x2.mul_add(res, 0.666_666_7);
116        res = x2.mul_add(res, 2.0);
117        x.mul_add(res, exponent * std::f32::consts::LN_2)
118    }
119
120    #[inline(always)]
121    fn exp(self) -> Self {
122        // Minimax polynomial approximation for exp(x) on [-0.5*ln(2), 0.5*ln(2))
123        // Optimized for f32 precision with max error ~1e-7
124        // Uses range reduction: exp(x) = 2^k * exp(r) where k = round(x / ln(2)) and r = x - k*ln(2)
125        // See: J.-M. Muller et al., "Handbook of Floating-Point Arithmetic", 2018, Section 10.3
126        //      A. J. Salgado & S. M. Wise, "Classical Numerical Analysis", 2023, Chapter 10
127        let k = (self.mul_add(
128            std::f32::consts::LOG2_E,
129            if self > 0.0 { 0.5 } else { -0.5 },
130        )) as i32;
131        let r = self - (k as Self * 0.69314575) - (k as Self * 1.4286068e-6);
132        let mut res = 0.00138889_f32;
133        res = r.mul_add(res, 0.00833333);
134        res = r.mul_add(res, 0.04166667);
135        res = r.mul_add(res, 0.16666667);
136        res = r.mul_add(res, 0.5);
137        res = r.mul_add(res, 1.0);
138        r.mul_add(res, 1.0) * Self::from_bits(((k + 127) as u32) << 23)
139    }
140
141    #[inline(always)]
142    fn cdf(self) -> Self {
143        self.cdf_with_pdf().0
144    }
145
146    #[inline(always)]
147    fn cdf_with_pdf(self) -> (Self, Self) {
148        // Minimax rational approximation for normal CDF
149        // Optimized for f32 precision with max error ~1e-7
150        // Uses transformation t = 1 / (1 + 0.2316419 * |x|) for numerical stability
151        // See: M. Abramowitz & I. A. Stegun (eds.), "Handbook of Mathematical Functions
152        //      with Formulas, Graphs, and Mathematical Tables", 1972, Section 26.2.17
153        let abs_x = self.abs();
154        let t = 1.0 / (1.0 + 0.2316419 * abs_x);
155        let mut poly = 1.330_274_5_f32.mul_add(t, -1.821_255_9);
156        poly = t.mul_add(poly, 1.781_477_9);
157        poly = t.mul_add(poly, -0.356_563_78);
158        poly = t.mul_add(poly, 0.319_381_54);
159        let pdf = 0.398_942_3 * (-0.5 * self * self).exp();
160        let res = 1.0 - pdf * (poly * t);
161        // Use >= to ensure CDF(0) = 0.5 exactly (maintains symmetry)
162        (if self >= 0.0 { res } else { 1.0 - res }, pdf)
163    }
164}
165
166// 3. DATA STRUCTURES & CORE KERNEL
167#[derive(Debug, Clone, Copy, Default, PartialEq)]
168pub struct Greeks<T> {
169    pub price: T,
170    pub vol: T,
171    pub delta: T,
172    pub gamma: T,
173    pub vega: T,
174    pub theta: T,
175}
176
177/// Lightweight kernel for IV search - only computes price and vega
178#[inline(always)]
179fn pricing_kernel_price_vega<T: BlackScholesReal>(
180    s: T,
181    k: T,
182    disc: T,
183    d1: T,
184    d2: T,
185    sqrt_t: T,
186    is_call: T::Mask,
187) -> (T, T) {
188    let (n_d1, pdf_d1) = d1.cdf_with_pdf();
189    let n_d2 = d2.cdf();
190
191    let c_price = s * n_d1 - k * disc * n_d2;
192    let p_price = c_price - s + k * disc;
193    let price = T::select(is_call, c_price, p_price);
194
195    let vega = s * sqrt_t * pdf_d1;
196
197    (price, vega)
198}
199
200#[allow(clippy::too_many_arguments)]
201#[inline(always)]
202fn pricing_kernel<T: BlackScholesReal>(
203    s: T,
204    k: T,
205    disc: T,
206    d1: T,
207    d2: T,
208    inv_vol_sqrt_t: T,
209    vol: T,
210    sqrt_t: T,
211    t: T,
212    r: T,
213    b: T,
214    s_orig: T,
215    is_call: T::Mask,
216) -> Greeks<T> {
217    let (n_d1, pdf_d1) = d1.cdf_with_pdf();
218    let n_d2 = d2.cdf();
219
220    let c_price = s * n_d1 - k * disc * n_d2;
221    let p_price = c_price - s + k * disc;
222
223    let vega = s * sqrt_t * pdf_d1;
224    let df = ((b - r) * t).exp();
225    let gamma = df * pdf_d1 * inv_vol_sqrt_t / s_orig;
226
227    let theta_base = -(vega * vol) * (T::splat(2.0) * t).recip_precise();
228    let phi_theta = T::select(is_call, T::splat(1.0), -T::splat(1.0));
229    let c_theta = theta_base - r * k * disc * n_d2 - phi_theta * (b - r) * s * n_d1;
230    let p_theta =
231        theta_base + r * k * disc * (T::splat(1.0) - n_d2) - phi_theta * (b - r) * s * n_d1;
232
233    Greeks {
234        price: T::select(is_call, c_price, p_price),
235        vol,
236        delta: T::select(is_call, n_d1, n_d1 - T::splat(1.0)),
237        gamma,
238        vega,
239        theta: T::select(is_call, c_theta, p_theta),
240    }
241}
242
243// 5. SOLVERS: STANDALONE GREEKS & IV SEARCH
244#[inline(always)]
245pub fn compute_greeks<T: BlackScholesReal>(
246    s: T,
247    k: T,
248    t: T,
249    r: T,
250    b: T,
251    vol: T,
252    is_call: T::Mask,
253) -> Greeks<T> {
254    let sqrt_t = t.sqrt();
255    let vol_sqrt_t = vol * sqrt_t;
256    let inv_vol_sqrt_t = vol_sqrt_t.recip_precise();
257    let disc = (-r * t).exp();
258    let d1 = ((s / k).ln() + (b + T::splat(0.5) * vol * vol) * t) * inv_vol_sqrt_t;
259    let s_forward = s * ((b - r) * t).exp();
260
261    pricing_kernel(
262        s_forward,
263        k,
264        disc,
265        d1,
266        d1 - vol_sqrt_t,
267        inv_vol_sqrt_t,
268        vol,
269        sqrt_t,
270        t,
271        r,
272        b,
273        s,
274        is_call,
275    )
276}
277
278/// Performs a single Halley iteration to refine an implied volatility estimate and compute greeks.
279///
280/// # Important Notes
281///
282/// This function is intended as a **refinement step** when a good initial guess for volatility
283/// is available (e.g., from a previous calculation or a fast approximation). It performs only
284/// a single Halley iteration and does not implement a full convergence loop.
285///
286/// **This is NOT a standalone implied volatility solver.** For production use, prefer
287/// `imply_vol_and_greeks` which uses the robust `implied_vol` crate for full convergence.
288///
289/// # Parameters
290///
291/// * `initial_guess`: Must be a reasonable estimate of the true volatility. Poor initial guesses
292///   (especially for deep ITM/OTM options) may result in significant errors.
293///
294/// # Accuracy
295///
296/// With a good initial guess (within ~25% of true vol), one Halley step typically achieves
297/// ~1% relative error. For deep ITM/OTM options or poor initial guesses, multiple iterations
298/// or a better initial estimate may be required.
299#[allow(clippy::too_many_arguments)]
300#[inline(always)]
301pub fn compute_iv_and_greeks<T: BlackScholesReal>(
302    mkt_price: T,
303    s: T,
304    k: T,
305    t: T,
306    r: T,
307    b: T,
308    is_call: T::Mask,
309    initial_guess: T,
310) -> Greeks<T> {
311    // PRE-CALCULATION (Hoisted outside iteration)
312    let sqrt_t = t.sqrt();
313    let inv_sqrt_t = sqrt_t.recip_precise();
314    let ln_sk_bt = (s.ln() - k.ln()) + (b * t); // Numerical Idea 1: Merged constant with b
315    let half_t = T::splat(0.5) * t; // Numerical Idea 2: Hoisted half-time
316    let disc = (-r * t).exp();
317    let mut vol = initial_guess;
318
319    // SINGLE HALLEY PASS
320    let inv_vol = vol.recip_precise();
321    let inv_vst = inv_vol * inv_sqrt_t;
322    let d1 = (ln_sk_bt + half_t * vol * vol) * inv_vst;
323    let d2 = d1 - vol * sqrt_t;
324    let s_forward = s * ((b - r) * t).exp();
325    let (price, vega_raw) = pricing_kernel_price_vega(s_forward, k, disc, d1, d2, sqrt_t, is_call);
326
327    let diff = price - mkt_price;
328    let vega = vega_raw.abs().max(T::splat(1e-9));
329    let volga = (vega * d1 * d2) * inv_vol;
330    let num = T::splat(2.0) * diff * vega;
331    let den = T::splat(2.0) * vega * vega - diff * volga;
332    // Clamp denominator magnitude while preserving sig
333    let den_safe = den.signum() * den.abs().max(T::splat(1e-9));
334    vol = vol - (num * den_safe.recip_precise());
335
336    // Clamp volatility to reasonable bounds to prevent negative or infinite values
337    // Lower bound: 1e-6 (0.0001% annualized), Upper bound: 10.0 (1000% annualized)
338    // Using max/min compiles to single instructions for f32
339    vol = vol.max(T::splat(1e-6)).min(T::splat(10.0));
340
341    // FINAL RE-SYNC
342    let inv_vol_f = vol.recip_precise();
343    let inv_vst_f = inv_vol_f * inv_sqrt_t;
344    let d1_f = (ln_sk_bt + half_t * vol * vol) * inv_vst_f;
345    let mut g_final = pricing_kernel(
346        s_forward,
347        k,
348        disc,
349        d1_f,
350        d1_f - vol * sqrt_t,
351        inv_vst_f,
352        vol,
353        sqrt_t,
354        t,
355        r,
356        b,
357        s,
358        is_call,
359    );
360    g_final.vol = vol;
361
362    g_final
363}
364
365// 4. UNIT TESTS
366#[cfg(test)]
367mod tests {
368    use rstest::*;
369
370    use super::*;
371    use crate::data::greeks::black_scholes_greeks_exact;
372
373    #[rstest]
374    fn test_accuracy_1e7() {
375        let s = 100.0;
376        let k = 100.0;
377        let t = 1.0;
378        let r = 0.05;
379        let vol = 0.2;
380        let g = compute_greeks::<f32>(s, k, t, r, r, vol, true); // Use r as b
381        assert!((g.price - 10.45058).abs() < 1e-5);
382        let solved = compute_iv_and_greeks::<f32>(g.price, s, k, t, r, r, true, vol); // Use r as b
383        assert!((solved.vol - vol).abs() < 1e-6);
384    }
385
386    #[rstest]
387    fn test_compute_greeks_accuracy_vs_exact() {
388        let s = 100.0f64;
389        let k = 100.0f64;
390        let t = 1.0f64;
391        let r = 0.05f64;
392        let b = 0.05f64; // cost of carry
393        let vol = 0.2f64;
394        let multiplier = 1.0f64;
395
396        // Compute using fast f32 method
397        let g_fast = compute_greeks::<f32>(
398            s as f32, k as f32, t as f32, r as f32, b as f32, vol as f32, true,
399        );
400
401        // Compute using exact f64 method
402        let g_exact = black_scholes_greeks_exact(s, r, b, vol, true, k, t, multiplier);
403
404        // Compare with tolerance for f32 precision
405        let price_tol = 1e-4;
406        let greeks_tol = 1e-3;
407
408        assert!(
409            (g_fast.price as f64 - g_exact.price).abs() < price_tol,
410            "Price mismatch: fast={}, exact={}",
411            g_fast.price,
412            g_exact.price
413        );
414        assert!(
415            (g_fast.delta as f64 - g_exact.delta).abs() < greeks_tol,
416            "Delta mismatch: fast={}, exact={}",
417            g_fast.delta,
418            g_exact.delta
419        );
420        assert!(
421            (g_fast.gamma as f64 - g_exact.gamma).abs() < greeks_tol,
422            "Gamma mismatch: fast={}, exact={}",
423            g_fast.gamma,
424            g_exact.gamma
425        );
426        // Vega units differ: exact uses multiplier * 0.01, fast uses raw units
427        let vega_exact_raw = g_exact.vega / (multiplier * 0.01);
428        assert!(
429            (g_fast.vega as f64 - vega_exact_raw).abs() < greeks_tol,
430            "Vega mismatch: fast={}, exact_raw={}, exact_scaled={}",
431            g_fast.vega,
432            vega_exact_raw,
433            g_exact.vega
434        );
435        // Theta units differ: exact uses multiplier * daily_factor (0.0027378507871321013), fast uses raw units
436        let theta_daily_factor = 0.0027378507871321013;
437        let theta_exact_raw = g_exact.theta / (multiplier * theta_daily_factor);
438        assert!(
439            (g_fast.theta as f64 - theta_exact_raw).abs() < greeks_tol,
440            "Theta mismatch: fast={}, exact_raw={}, exact_scaled={}",
441            g_fast.theta,
442            theta_exact_raw,
443            g_exact.theta
444        );
445    }
446
447    #[rstest]
448    fn test_compute_iv_and_greeks_halley_accuracy() {
449        let s = 100.0f64;
450        let k = 100.0f64;
451        let t = 1.0f64;
452        let r = 0.05f64;
453        let b = 0.05f64; // cost of carry
454        let vol_true = 0.2f64; // True volatility
455        let initial_guess = 0.25f64; // Initial guess (25% higher than true)
456        let multiplier = 1.0f64;
457
458        // Compute the exact price using the true volatility
459        let g_exact = black_scholes_greeks_exact(s, r, b, vol_true, true, k, t, multiplier);
460        let mkt_price = g_exact.price;
461
462        // Compute implied vol using one Halley step with initial guess
463        let g_halley = compute_iv_and_greeks::<f32>(
464            mkt_price as f32,
465            s as f32,
466            k as f32,
467            t as f32,
468            r as f32,
469            b as f32,
470            true,
471            initial_guess as f32,
472        );
473
474        // Check that one Halley step gets close to the true volatility
475        let vol_error = (g_halley.vol as f64 - vol_true).abs();
476
477        // One Halley step should get within ~1% of true vol for a 25% initial error
478        assert!(
479            vol_error < 0.01,
480            "Halley step accuracy: vol_error={}, initial_guess={}, vol_true={}, computed_vol={}",
481            vol_error,
482            initial_guess,
483            vol_true,
484            g_halley.vol
485        );
486
487        // Check that the computed greeks are close to exact
488        let price_tol = 5e-3; // Relaxed for one Halley step
489        let greeks_tol = 5e-3; // Relaxed for one-step approximation
490
491        assert!(
492            (g_halley.price as f64 - g_exact.price).abs() < price_tol,
493            "Price mismatch after Halley: halley={}, exact={}, diff={}",
494            g_halley.price,
495            g_exact.price,
496            (g_halley.price as f64 - g_exact.price).abs()
497        );
498        assert!(
499            (g_halley.delta as f64 - g_exact.delta).abs() < greeks_tol,
500            "Delta mismatch after Halley: halley={}, exact={}",
501            g_halley.delta,
502            g_exact.delta
503        );
504        assert!(
505            (g_halley.gamma as f64 - g_exact.gamma).abs() < greeks_tol,
506            "Gamma mismatch after Halley: halley={}, exact={}",
507            g_halley.gamma,
508            g_exact.gamma
509        );
510        // Vega units differ: exact uses multiplier * 0.01, fast uses raw units
511        let vega_exact_raw = g_exact.vega / (multiplier * 0.01);
512        assert!(
513            (g_halley.vega as f64 - vega_exact_raw).abs() < greeks_tol,
514            "Vega mismatch after Halley: halley={}, exact_raw={}",
515            g_halley.vega,
516            vega_exact_raw
517        );
518        // Theta units differ: exact uses multiplier * daily_factor (0.0027378507871321013), fast uses raw units
519        let theta_daily_factor = 0.0027378507871321013;
520        let theta_exact_raw = g_exact.theta / (multiplier * theta_daily_factor);
521        assert!(
522            (g_halley.theta as f64 - theta_exact_raw).abs() < greeks_tol,
523            "Theta mismatch after Halley: halley={}, exact_raw={}",
524            g_halley.theta,
525            theta_exact_raw
526        );
527    }
528
529    #[rstest]
530    fn test_print_halley_iv() {
531        let s = 100.0f64;
532        let k = 100.0f64;
533        let t = 1.0f64;
534        let r = 0.05f64;
535        let b = 0.05f64;
536        let vol_true = 0.2f64;
537        let multiplier = 1.0f64;
538
539        let g_exact = black_scholes_greeks_exact(s, r, b, vol_true, true, k, t, multiplier);
540        let mkt_price = g_exact.price;
541
542        println!("\n=== Halley Step IV Test (Using True Vol as Initial Guess) ===");
543        println!("True volatility: {vol_true}");
544        println!("Market price: {mkt_price:.8}");
545        println!("Initial guess: {vol_true} (using true vol)");
546
547        let g_halley = compute_iv_and_greeks::<f32>(
548            mkt_price as f32,
549            s as f32,
550            k as f32,
551            t as f32,
552            r as f32,
553            b as f32,
554            true,
555            vol_true as f32, // Using true vol as initial guess
556        );
557
558        println!("\nAfter one Halley step:");
559        println!("Computed volatility: {:.8}", g_halley.vol);
560        println!("True volatility: {vol_true:.8}");
561        println!(
562            "Absolute error: {:.8}",
563            (g_halley.vol as f64 - vol_true).abs()
564        );
565        println!(
566            "Relative error: {:.4}%",
567            (g_halley.vol as f64 - vol_true).abs() / vol_true * 100.0
568        );
569    }
570
571    #[rstest]
572    fn test_compute_iv_and_greeks_deep_itm_otm() {
573        let t = 1.0f64;
574        let r = 0.05f64;
575        let b = 0.05f64;
576        let vol_true = 0.2f64;
577        let multiplier = 1.0f64;
578
579        // Deep ITM: s=150, k=100 (spot is 50% above strike)
580        let s_itm = 150.0f64;
581        let k_itm = 100.0f64;
582        let g_exact_itm =
583            black_scholes_greeks_exact(s_itm, r, b, vol_true, true, k_itm, t, multiplier);
584        let mkt_price_itm = g_exact_itm.price;
585
586        println!("\n=== Deep ITM Test ===");
587        println!("Spot: {s_itm}, Strike: {k_itm}, True vol: {vol_true}");
588        println!("Market price: {mkt_price_itm:.8}");
589
590        let g_recovered_itm = compute_iv_and_greeks::<f32>(
591            mkt_price_itm as f32,
592            s_itm as f32,
593            k_itm as f32,
594            t as f32,
595            r as f32,
596            b as f32,
597            true,
598            vol_true as f32, // Using true vol as initial guess
599        );
600
601        let vol_error_itm = (g_recovered_itm.vol as f64 - vol_true).abs();
602        let rel_error_itm = vol_error_itm / vol_true * 100.0;
603
604        println!("Recovered volatility: {:.8}", g_recovered_itm.vol);
605        println!("Absolute error: {vol_error_itm:.8}");
606        println!("Relative error: {rel_error_itm:.4}%");
607
608        // Deep OTM: s=50, k=100 (spot is 50% below strike)
609        let s_otm = 50.0f64;
610        let k_otm = 100.0f64;
611        let g_exact_otm =
612            black_scholes_greeks_exact(s_otm, r, b, vol_true, true, k_otm, t, multiplier);
613        let mkt_price_otm = g_exact_otm.price;
614
615        println!("\n=== Deep OTM Test ===");
616        println!("Spot: {s_otm}, Strike: {k_otm}, True vol: {vol_true}");
617        println!("Market price: {mkt_price_otm:.8}");
618
619        let g_recovered_otm = compute_iv_and_greeks::<f32>(
620            mkt_price_otm as f32,
621            s_otm as f32,
622            k_otm as f32,
623            t as f32,
624            r as f32,
625            b as f32,
626            false,
627            vol_true as f32, // Using true vol as initial guess
628        );
629
630        let vol_error_otm = (g_recovered_otm.vol as f64 - vol_true).abs();
631        let rel_error_otm = vol_error_otm / vol_true * 100.0;
632
633        println!("Recovered volatility: {:.8}", g_recovered_otm.vol);
634        println!("Absolute error: {vol_error_otm:.8}");
635        println!("Relative error: {rel_error_otm:.4}%");
636
637        // Assertions: Deep ITM and OTM are challenging cases
638        // One Halley step with Corrado-Miller initial guess may not be sufficient
639        // We use a more relaxed tolerance to verify the method still converges in the right direction
640        // For production use, multiple iterations or better initial guesses would be needed
641        let vol_tol_itm = 50.0; // 50% relative error tolerance for deep ITM
642        let vol_tol_otm = 150.0; // 150% relative error tolerance for deep OTM (very challenging)
643
644        // Check that we at least get a reasonable volatility (not NaN or extreme values)
645        assert!(
646            g_recovered_itm.vol.is_finite()
647                && g_recovered_itm.vol > 0.0
648                && g_recovered_itm.vol < 2.0,
649            "Deep ITM vol recovery: invalid result={}",
650            g_recovered_itm.vol
651        );
652
653        assert!(
654            g_recovered_otm.vol.is_finite()
655                && g_recovered_otm.vol > 0.0
656                && g_recovered_otm.vol < 2.0,
657            "Deep OTM vol recovery: invalid result={}",
658            g_recovered_otm.vol
659        );
660
661        // Verify the error is within acceptable bounds (one step may not be enough)
662        assert!(
663            rel_error_itm < vol_tol_itm,
664            "Deep ITM vol recovery error too large: recovered={}, true={}, error={:.4}%",
665            g_recovered_itm.vol,
666            vol_true,
667            rel_error_itm
668        );
669
670        assert!(
671            rel_error_otm < vol_tol_otm,
672            "Deep OTM vol recovery error too large: recovered={}, true={}, error={:.4}%",
673            g_recovered_otm.vol,
674            vol_true,
675            rel_error_otm
676        );
677
678        println!("\n=== Summary ===");
679        println!("Deep ITM: One Halley iteration error = {rel_error_itm:.2}%");
680        println!(
681            "Deep OTM: One Halley iteration error = {rel_error_otm:.2}% (still challenging, deep OTM is difficult)"
682        );
683    }
684}
nautilus_model/data/black_scholes.rs

nautilus_model/data/
black_scholes.rs
