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

nautilus_model/data/
black_scholes.rs
