Skip to content

Commit

Permalink
fix: Fix to use f64::MIN_POSITIVE
Browse files Browse the repository at this point in the history 
refactor: Change to use more f64 methods
  • Loading branch information
@nakashima-hikaru
nakashima-hikaru committed Jan 21, 2024
1 parent db619a3 commit 5016341
Show file tree
Hide file tree
Showing 6 changed files with 109 additions and 107 deletions.
30 changes: 15 additions & 15 deletions src/bachelier.rs
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
use std::cmp::Ordering;
use crate::constants::{ONE_OVER_SQRT_TWO_PI, SQRT_TWO_PI};
use crate::normal_distribution::norm_pdf;

Expand All @@ -13,7 +14,7 @@ fn phi_tilde_times_x(x: f64) -> f64 {
return ONE_OVER_SQRT_TWO_PI + x * (0.5 + x * g);
}

if x > 0.0 {
if x.is_sign_positive() {
return phi_tilde_times_x(-x) + x;
}

Expand All @@ -22,7 +23,7 @@ fn phi_tilde_times_x(x: f64) -> f64 {
return (-0.5 * (x * x)).exp() * g;
}

let w = 1.0 / (x * x);
let w = (x * x).recip();
let g = (2.999_999_999_999_991 + w * (2.365_455_662_782_315E2 + w * (6.812_677_344_935_879E3 + w * (8.969_794_159_836_079E4 + w * (5.516_392_059_126_862E5 + w * (1.434_506_112_333_566_2E6 + w * (1.150_498_824_634_488_2E6 + 1.186_760_040_099_769_1E4 * w))))))) / (1.0 + w * (8.384_852_209_273_714E1 + w * (2.655_135_058_780_958E3 + w * (4.055_529_088_467_379E4 + w * (3.166_737_476_299_376_6E5 + w * (1.232_979_595_802_432_2E6 + w * (2.140_981_054_061_905E6 + 1.214_566_780_409_316E6 * w)))))));
ONE_OVER_SQRT_TWO_PI * (-0.5 * (x * x)).exp() * w * (1.0 - g * w)
}
Expand All @@ -35,12 +36,12 @@ fn inv_phi_tilde(phi_tilde_star: f64) -> f64 {
if phi_tilde_star > 1.0 {
return -inv_phi_tilde(1.0 - phi_tilde_star);
}
if phi_tilde_star >= 0.0 {
if !phi_tilde_star.is_sign_negative() {
return f64::NAN;
}
let x_bar = if phi_tilde_star < -0.00188203927 {
// Equation (2.1)
let g = 1.0 / (phi_tilde_star - 0.5);
let g = (phi_tilde_star - 0.5).recip();
let g2 = g * g;
// Equation (2.2)
let xi_bar = (0.032114372355 - g2 * (0.016969777977 - g2 * (0.002620733246 - 0.000096066952861 * g2))) /
Expand Down Expand Up @@ -80,8 +81,7 @@ pub(crate) fn bachelier(forward: f64, strike: f64, sigma: f64, t: f64, q: bool)
if s < f64::MIN_POSITIVE {
return intrinsic_value(forward, strike, q);
}
let theta = if q { 1.0 } else { -1.0 };
let moneyness = theta * (forward - strike);
let moneyness = if q { forward - strike } else { strike - forward };
let x = moneyness / s;
s * phi_tilde_times_x(x)
}
Expand All @@ -91,16 +91,16 @@ pub(crate) fn implied_normal_volatility(price: f64, forward: f64, strike: f64, t
return price * SQRT_TWO_PI / t.sqrt();
}
let intrinsic = intrinsic_value(forward, strike, q);
if price == intrinsic {
return 0.0;
}
if price < intrinsic {
return f64::NEG_INFINITY;
match price.total_cmp(&intrinsic) {
Ordering::Less => { f64::NEG_INFINITY }
Ordering::Equal => { 0.0 }
Ordering::Greater => {
let absolute_moneyness = (forward - strike).abs();
let phi_tilde_star = (intrinsic - price) / absolute_moneyness;
let x_star = inv_phi_tilde(phi_tilde_star);
absolute_moneyness / (x_star * t.sqrt()).abs()
}
}
let absolute_moneyness = (forward - strike).abs();
let phi_tilde_star = (intrinsic - price) / absolute_moneyness;
let x_star = inv_phi_tilde(phi_tilde_star);
absolute_moneyness / (x_star * t.sqrt()).abs()
}


Expand Down
7 changes: 3 additions & 4 deletions src/constants.rs
View file
Original file line number Diff line number Diff line change
@@ -1,19 +1,18 @@
pub(crate) const TWO_PI: f64 = std::f64::consts::TAU;
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 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 PI_OVER_SIX: f64 = std::f64::consts::FRAC_PI_6;
pub(crate) const FOURTH_ROOT_DBL_EPSILON: f64 = 0.0001220703125;
pub(crate) const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433;
pub(crate) const SQRT_DBL_MIN: f64 = 1.4916681462400413e-154;
pub(crate) const SQRT_MIN_POSITIVE: f64 = 1.4916681462400413e-154;
pub(crate) const SQRT_DBL_MAX: f64 = 1.3407807929942596e154;
// pub(crate) const SQRT_DBL_MIN: f64 = 1.3407807929942596e154;
// Set this to 0 if you want positive results for (positive) denormalised inputs, else to DBL_MIN.
// Note that you cannot achieve full machine accuracy from denormalised inputs!
pub(crate) const DENORMALISATION_CUTOFF: f64 = 0.0;
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 LN_TWO_PI: f64 = 1.8378770664093453;
pub(crate) const HALF_OF_LN_TWO_PI: f64 = 0.918_938_533_204_672_8;
76 changes: 35 additions & 41 deletions src/erf_cody.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,14 +48,8 @@ const Q: [f64; 5] = [
0.002_335_204_976_268_691_8,
];

const ZERO: f64 = 0.0;
// const HALF: f64 = 0.5;
const ONE: f64 = 1.0;
const TWO: f64 = 2.0;
const FOUR: f64 = 4.0;
const SQRPI: f64 = 0.564_189_583_547_756_3;
const THRESH: f64 = 0.46875;
const SIXTEN: f64 = 16.0;
const XINF: f64 = f64::MAX;
const XNEG: f64 = -26.628;
const XSMALL: f64 = 1.11e-16;
Expand All @@ -70,10 +64,10 @@ pub(crate) fn erfc_cody(x: f64) -> f64 {
/* Author/date: W. J. Cody, January 8, 1985 */
/* -------------------------------------------------------------------- */
let y = x.abs();
let mut ysq = ZERO;
let mut ysq = 0.0;
let mut xden;
let mut xnum;
let mut result = ZERO;
let mut result = 0.0;

if y <= THRESH {
if y > XSMALL {
Expand All @@ -88,9 +82,9 @@ pub(crate) fn erfc_cody(x: f64) -> f64 {
}
result = x * (xnum + A[3]) / (xden + B[3]);

result = ONE - result;
result = 1.0 - result;
return result;
} else if y <= FOUR {
} else if y <= 4.0 {
xnum = C[8] * y;
xden = y;

Expand All @@ -100,16 +94,16 @@ pub(crate) fn erfc_cody(x: f64) -> f64 {
}
result = (xnum + C[7]) / (xden + D[7]);

ysq = (y * SIXTEN).trunc() / SIXTEN;
ysq = (y * 16.0).trunc() / 16.0;
let del = (y - ysq) * (y + ysq);
result *= (-ysq * ysq).exp() * (-del).exp();
} else if y >= XBIG {
if x < ZERO {
result = TWO - result;
if x.is_sign_negative() {
result = 2.0 - result;
}
return result;
} else {
ysq = ONE / (y * y);
ysq = (y * y).recip();
xnum = P[5] * ysq;
xden = ysq;

Expand All @@ -120,12 +114,12 @@ pub(crate) fn erfc_cody(x: f64) -> f64 {
result = ysq * (xnum + P[4]) / (xden + Q[4]);
result = (SQRPI - result) / y;

ysq = (y * SIXTEN).trunc() / SIXTEN;
ysq = (y * 16.0).trunc() / 16.0;
let del = (y - ysq) * (y + ysq);
result *= (-ysq * ysq).exp() * (-del).exp();
}
if x < ZERO {
result = TWO - result;
if x.is_sign_negative() {
result = 2.0 - result;
}
result
}
Expand All @@ -137,10 +131,10 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
/* Author/date: W. J. Cody, March 30, 1987 */
/* ------------------------------------------------------------------ */
let y = x.abs();
let mut ysq = ZERO;
let mut ysq = 0.0;
let mut xden;
let mut xnum;
let mut result = ZERO;
let mut result = 0.0;

if y <= THRESH {
if y > XSMALL {
Expand All @@ -155,11 +149,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
}
result = x * (xnum + A[3]) / (xden + B[3]);

result = ONE - result;
result = 1.0 - result;

result *= ysq.exp();
return result;
} else if y <= FOUR {
} else if y <= 4.0 {
xnum = C[8] * y;
xden = y;

Expand All @@ -170,11 +164,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
result = (xnum + C[7]) / (xden + D[7]);
} else if y >= XBIG {
if y >= XMAX {
if x < ZERO {
if x.is_sign_negative() {
if x < XNEG {
result = XINF;
} else {
let ysq = (x * SIXTEN).trunc() / SIXTEN;
let ysq = (x * 16.0).trunc() / 16.0;
let del = (x - ysq) * (x + ysq);
let y = (ysq * ysq).exp() * del.exp();
result = (y + y) - result;
Expand All @@ -183,11 +177,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
return result;
} else if y >= XHUGE {
result = SQRPI / y;
if x < ZERO {
if x.is_sign_negative() {
if x < XNEG {
result = XINF;
} else {
let ysq = (x * SIXTEN).trunc() / SIXTEN;
let ysq = (x * 16.0).trunc() / 16.0;
let del = (x - ysq) * (x + ysq);
let y = (ysq * ysq).exp() * del.exp();
result = (y + y) - result;
Expand All @@ -196,7 +190,7 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
return result;
}
} else {
ysq = ONE / (y * y);
ysq = (y * y).recip();
xnum = P[5] * ysq;
xden = ysq;

Expand All @@ -207,11 +201,11 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
result = ysq * (xnum + P[4]) / (xden + Q[4]);
result = (SQRPI - result) / y;
}
if x < ZERO {
if x.is_sign_negative() {
if x < XNEG {
result = XINF;
} else {
let ysq = (x * SIXTEN).trunc() / SIXTEN;
let ysq = (x * 16.0).trunc() / 16.0;
let del = (x - ysq) * (x + ysq);
let y = (ysq * ysq).exp() * del.exp();
result = (y + y) - result;
Expand All @@ -223,7 +217,7 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {
#[cfg(test)]
mod tests {
use crate::erf_cody::{
erfc_cody, erfcx_cody, FOUR, THRESH, XBIG, XHUGE, XMAX, XNEG, ZERO,
erfc_cody, erfcx_cody, THRESH, XBIG, XHUGE, XMAX, XNEG,
};

#[test]
Expand All @@ -237,13 +231,13 @@ mod tests {
let x = erfc_cody(-THRESH + f64::EPSILON);
assert_eq!(x, 1.4926134732179377);

let x = erfc_cody(FOUR + f64::EPSILON);
let x = erfc_cody(4.0 + f64::EPSILON);
assert_eq!(x, 1.541725790028002e-8);
let x = erfc_cody(FOUR - f64::EPSILON);
let x = erfc_cody(4.0 - f64::EPSILON);
assert_eq!(x, 1.541725790028002e-8);
let x = erfc_cody(-FOUR - f64::EPSILON);
let x = erfc_cody(-4.0 - f64::EPSILON);
assert_eq!(x, 1.999999984582742);
let x = erfc_cody(-FOUR + f64::EPSILON);
let x = erfc_cody(-4.0 + f64::EPSILON);
assert_eq!(x, 1.999999984582742);

let x = erfc_cody(XBIG + f64::EPSILON);
Expand Down Expand Up @@ -273,9 +267,9 @@ mod tests {
let x = erfc_cody(-XHUGE + f64::EPSILON);
assert_eq!(x, 2.0);

let x = erfc_cody(ZERO + f64::EPSILON);
let x = erfc_cody(0.0 + f64::EPSILON);
assert_eq!(x, 0.9999999999999998);
let x = erfc_cody(ZERO - f64::EPSILON);
let x = erfc_cody(0.0 - f64::EPSILON);
assert_eq!(x, 1.0000000000000002);

let x = erfc_cody(XNEG + f64::EPSILON);
Expand All @@ -295,13 +289,13 @@ mod tests {
let x = erfcx_cody(-THRESH + f64::EPSILON);
assert_eq!(x, 1.8594024168714214);

let x = erfcx_cody(FOUR + f64::EPSILON);
let x = erfcx_cody(4.0 + f64::EPSILON);
assert_eq!(x, 0.1369994576250614);
let x = erfcx_cody(FOUR - f64::EPSILON);
let x = erfcx_cody(4.0 - f64::EPSILON);
assert_eq!(x, 0.1369994576250614);
let x = erfcx_cody(-FOUR - f64::EPSILON);
let x = erfcx_cody(-4.0 - f64::EPSILON);
assert_eq!(x, 17772220.904016286);
let x = erfcx_cody(-FOUR + f64::EPSILON);
let x = erfcx_cody(-4.0 + f64::EPSILON);
assert_eq!(x, 17772220.904016286);

let x = erfcx_cody(XBIG + f64::EPSILON);
Expand Down Expand Up @@ -331,9 +325,9 @@ mod tests {
let x = erfcx_cody(-XHUGE + f64::EPSILON);
assert_eq!(x, 1.7976931348623157e308);

let x = erfcx_cody(ZERO + f64::EPSILON);
let x = erfcx_cody(0.0 + f64::EPSILON);
assert_eq!(x, 0.9999999999999998);
let x = erfcx_cody(ZERO - f64::EPSILON);
let x = erfcx_cody(0.0 - f64::EPSILON);
assert_eq!(x, 1.0000000000000002);

let x = erfcx_cody(XNEG + f64::EPSILON);
Expand Down
Loading

0 comments on commit 5016341

Please sign in to comment.