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

nautilus_model/data/
black_scholes.rs
