From 85da785571826fd89a9df45613ac10c6c5976613 Mon Sep 17 00:00:00 2001 From: hnakashima Date: Sun, 21 Jan 2024 19:47:17 +0900 Subject: [PATCH] Update usage of f64 methods and constants The calculation in the rational model have been optimized by updating f64 methods and constants across different functions in src/lets_be_rational.rs, src/normal_distribution.rs and src/constants.rs. This includes switching from SQRT_2.recip() to FRAC_1_SQRT_2, and from SQRT_2 * x.abs().sqrt() to SQRT_TWO_OVER_PI for clearer and more efficient calculations. --- src/constants.rs | 5 ++++- src/lets_be_rational.rs | 16 ++++++++-------- src/normal_distribution.rs | 3 ++- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/constants.rs b/src/constants.rs index b79a7aa..5eaab90 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,7 +1,8 @@ pub(crate) const SQRT_PI_OVER_TWO: f64 = 1.253_314_137_315_500_3; pub(crate) const SQRT_TWO_PI: f64 = 2.506_628_274_631_000_7; pub(crate) const SQRT_THREE: f64 = 1.732_050_807_568_877_2; -pub(crate) const SQRT_ONE_OVER_THREE: f64 = 0.577_350_269_189_625_7; +pub(crate) const SQRT_TWO_OVER_PI: f64 = 0.797_884_560_802_865_4; +pub(crate) const ONE_OVER_SQRT_THREE: f64 = 0.577_350_269_189_625_7; pub(crate) const TWO_PI_OVER_SQRT_TWENTY_SEVEN: f64 = 1.209_199_576_156_145_2; pub(crate) const SQRT_THREE_OVER_THIRD_ROOT_TWO_PI: f64 = 0.938_643_487_427_383_6; pub(crate) const FOURTH_ROOT_DBL_EPSILON: f64 = 0.0001220703125; @@ -16,3 +17,5 @@ pub(crate) const ONE_OVER_SQRT_TWO_PI: f64 = 0.3989422804014327; pub(crate) const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC: f64 = f64::NEG_INFINITY; pub(crate) const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM: f64 = f64::INFINITY; pub(crate) const HALF_OF_LN_TWO_PI: f64 = 0.918_938_533_204_672_8; + + diff --git a/src/lets_be_rational.rs b/src/lets_be_rational.rs index ddad641..3f55b4f 100644 --- a/src/lets_be_rational.rs +++ b/src/lets_be_rational.rs @@ -1,5 +1,5 @@ -use std::f64::consts::SQRT_2; -use crate::constants::{DENORMALISATION_CUTOFF, FOURTH_ROOT_DBL_EPSILON, HALF_OF_LN_TWO_PI, SIXTEENTH_ROOT_DBL_EPSILON, SQRT_DBL_MAX, SQRT_MIN_POSITIVE, SQRT_ONE_OVER_THREE, SQRT_PI_OVER_TWO, SQRT_THREE, SQRT_THREE_OVER_THIRD_ROOT_TWO_PI, SQRT_TWO_PI, TWO_PI_OVER_SQRT_TWENTY_SEVEN, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC}; +use std::f64::consts::{FRAC_1_SQRT_2, SQRT_2}; +use crate::constants::{DENORMALISATION_CUTOFF, FOURTH_ROOT_DBL_EPSILON, HALF_OF_LN_TWO_PI, SIXTEENTH_ROOT_DBL_EPSILON, SQRT_DBL_MAX, SQRT_MIN_POSITIVE, ONE_OVER_SQRT_THREE, SQRT_PI_OVER_TWO, SQRT_THREE, SQRT_THREE_OVER_THIRD_ROOT_TWO_PI, SQRT_TWO_PI, TWO_PI_OVER_SQRT_TWENTY_SEVEN, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM, VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC, SQRT_TWO_OVER_PI}; use crate::erf_cody::{erfc_cody, erfcx_cody}; use crate::normal_distribution::{inverse_norm_cdf, norm_cdf, norm_pdf}; use crate::rational_cubic::{convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side, convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side, rational_cubic_interpolation}; @@ -62,7 +62,7 @@ fn asymptotic_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f6 fn small_t_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f64 { let w = t * t; let h2 = h * h; - let a = 1_f64 + h * SQRT_PI_OVER_TWO * erfcx_cody(-SQRT_2.recip() * h); + let a = 1_f64 + h * SQRT_PI_OVER_TWO * erfcx_cody(-FRAC_1_SQRT_2 * h); let b_over_vega = 2.0 * t * (a + w * ((-1.0 + 3.0 * a + a * h2) / 6.0 + w * ((-7.0 + 15.0 * a + h2 * (-1.0 + 10.0 * a + a * h2)) / 120.0 + w * ((-57.0 + 105.0 * a + h2 * (-18.0 + 105.0 * a + h2 * (-1.0 + 21.0 * a + a * h2))) / 5040.0 + w * ((-561.0 + 945.0 * a + h2 * (-285.0 + 1260.0 * a + h2 * (-33.0 + 378.0 * a + h2 * (-1.0 + 36.0 * a + a * h2)))) / 362880.0 + w * ((-6555.0 + 10395.0 * a + h2 * (-4680.0 + 17325.0 * a + h2 * (-840.0 + 6930.0 * a + h2 * (-52.0 + 990.0 * a + h2 * (-1.0 + 55.0 * a + a * h2))))) / 39916800.0 + ((-89055.0 + 135135.0 * a + h2 * (-82845.0 + 270270.0 * a + h2 * (-20370.0 + 135135.0 * a + h2 * (-1926.0 + 25740.0 * a + h2 * (-75.0 + 2145.0 * a + h2 * (-1.0 + 78.0 * a + a * h2)))))) * w) / 6227020800.0)))))); b_over_vega.max(0.0).abs() } @@ -71,8 +71,8 @@ fn normalised_black_call_with_optimal_use_of_codys_functions(x: f64, s: f64) -> const CODYS_THRESHOLD: f64 = 0.46875; let h = x / s; let t = 0.5 * s; - let q1 = -SQRT_2.recip() * (h + t); - let q2 = -SQRT_2.recip() * (h - t); + let q1 = -FRAC_1_SQRT_2 * (h + t); + let q2 = -FRAC_1_SQRT_2 * (h - t); let two_b: f64 = if q1 < CODYS_THRESHOLD { if q2 < CODYS_THRESHOLD { @@ -160,7 +160,7 @@ pub(crate) fn black(f: f64, k: f64, sigma: f64, t: f64, q: bool) -> f64 { fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f64) { let ax = x.abs(); - let z = SQRT_ONE_OVER_THREE * ax / s; + let z = ONE_OVER_SQRT_THREE * ax / s; let y = z * z; let s2 = s * s; let phi_m = norm_cdf(-z); @@ -238,7 +238,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let mut ds = f64::MIN; let mut s_left = f64::MIN_POSITIVE; let mut s_right = f64::MAX; - let s_c = SQRT_2 * x.abs().sqrt(); + let s_c = SQRT_2 * x.abs().sqrt(); let b_c = normalised_black_call(x, s_c); let v_c = normalised_vega(x, s_c); if beta < b_c { @@ -316,7 +316,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit while iterations < n && ds.abs() > f64::EPSILON * s { let h = x / s; let t = s / 2.0; - let gp = (2.0 / SQRT_TWO_PI) / (erfcx_cody((t + h) * SQRT_2.recip()) + erfcx_cody((t - h) * SQRT_2.recip())); + let gp = SQRT_TWO_OVER_PI / (erfcx_cody((t + h) * FRAC_1_SQRT_2) + erfcx_cody((t - h) * FRAC_1_SQRT_2)); let b_bar = normalised_vega(x, s) / gp; let g = (beta_bar / b_bar).ln(); let x_over_s_square = (h * h) / s; diff --git a/src/normal_distribution.rs b/src/normal_distribution.rs index 81e0b66..8545179 100644 --- a/src/normal_distribution.rs +++ b/src/normal_distribution.rs @@ -1,3 +1,4 @@ +use std::f64::consts::FRAC_1_SQRT_2; use crate::erf_cody::erfc_cody; const NORM_CDF_ASYMPTOTIC_EXPANSION_FIRST_THRESHOLD: f64 = -10.0; @@ -37,7 +38,7 @@ pub(crate) fn norm_cdf(z: f64) -> f64 { } return -norm_pdf(z) * sum / z; } - 0.5 * erfc_cody(-z * std::f64::consts::SQRT_2.recip()) + 0.5 * erfc_cody(-z * FRAC_1_SQRT_2) }