Skip to content

Commit

Permalink
chore: Upgrade to 1.0.0
Browse files Browse the repository at this point in the history
chore: Apply formatter (as many as possible)
  • Loading branch information
nakashima-hikaru committed Aug 11, 2024
1 parent 66504fe commit 5cd21dd
Show file tree
Hide file tree
Showing 8 changed files with 331 additions and 71 deletions.
2 changes: 1 addition & 1 deletion Cargo.toml
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[package]
name = "implied-vol"
version = "0.2.2"
version = "1.0.0"
authors = ["Hikaru Nakashima <nakashima.alg57@gmail.com>"]
description = "A pure rust implementation of Peter Jäckel's implied volatility calculation"
edition = "2021"
Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,7 @@
[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)

More information about this crate can be found in
the [crate documentation](https://docs.rs/implied-vol/0.2/implied_vol/).
the [crate documentation](https://docs.rs/implied-vol/1.0/implied_vol/).

## About

Expand Down
90 changes: 71 additions & 19 deletions src/bachelier.rs
Original file line number Diff line number Diff line change
@@ -1,16 +1,27 @@
use std::cmp::Ordering;
use crate::constants::{ONE_OVER_SQRT_TWO_PI, SQRT_TWO_PI};
use crate::normal_distribution::norm_pdf;
use std::cmp::Ordering;

#[inline]
fn intrinsic_value(forward: f64, strike: f64, q: bool) -> f64 {
(if !q { strike - forward } else { forward - strike }).max(0.0).abs()
(if !q {
strike - forward
} else {
forward - strike
})
.max(0.0)
.abs()
}

fn phi_tilde_times_x(x: f64) -> f64 {
if x.abs() <= 0.612_003_180_962_480_7 {
let h = (x * x - 1.872_739_467_540_974_8E-1) * 5.339_771_053_755_08;
let g = (1.964_154_984_377_470_3E-1 + h * (2.944_481_222_626_891_4E-3 + 3.095_828_855_856_471E-5 * h)) / (1.0 + h * (3.026_101_684_659_232_6E-2 + h * (3.373_546_191_189_62E-4 + h * (1.290_112_376_540_573_2E-6 - 1.671_197_583_524_420_5E-9 * h))));
let g = (1.964_154_984_377_470_3E-1
+ h * (2.944_481_222_626_891_4E-3 + 3.095_828_855_856_471E-5 * h))
/ (1.0
+ h * (3.026_101_684_659_232_6E-2
+ h * (3.373_546_191_189_62E-4
+ h * (1.290_112_376_540_573_2E-6 - 1.671_197_583_524_420_5E-9 * h))));
return ONE_OVER_SQRT_TWO_PI + x * (0.5 + x * g);
}

Expand All @@ -19,12 +30,42 @@ fn phi_tilde_times_x(x: f64) -> f64 {
}

if x >= -3.5 {
let g = (3.989_422_804_009_617_5E-1 + x * ((-2.882_725_012_271_64E-1) + x * (1.174_893_477_005_507_4E-1 + x * ((-2.920_893_049_832_423_4E-2) + x * (4.670_481_708_734_893E-3 + x * ((-4.444_840_548_247_636E-4) + x * (1.986_526_744_238_593_6E-5 + x * (7.638_739_347_414_361E-10 + 1.329_152_522_013_758_3E-11 * x)))))))) / (1.0 + x * ((-1.975_906_139_672_860_6) + x * (1.770_933_219_893_362_5 + x * ((-9.435_025_002_644_624E-1) + x * (3.281_611_814_538_859_5E-1 + x * ((-7.669_740_808_821_474E-2) + x * (1.184_322_430_309_622_3E-2 + x * ((-1.115_141_636_552_486_1E-3) + 4.974_100_533_375_869E-5 * x))))))));
let g = (3.989_422_804_009_617_5E-1
+ x * ((-2.882_725_012_271_64E-1)
+ x * (1.174_893_477_005_507_4E-1
+ x * ((-2.920_893_049_832_423_4E-2)
+ x * (4.670_481_708_734_893E-3
+ x * ((-4.444_840_548_247_636E-4)
+ x * (1.986_526_744_238_593_6E-5
+ x * (7.638_739_347_414_361E-10
+ 1.329_152_522_013_758_3E-11 * x))))))))
/ (1.0
+ x * ((-1.975_906_139_672_860_6)
+ x * (1.770_933_219_893_362_5
+ x * ((-9.435_025_002_644_624E-1)
+ x * (3.281_611_814_538_859_5E-1
+ x * ((-7.669_740_808_821_474E-2)
+ x * (1.184_322_430_309_622_3E-2
+ x * ((-1.115_141_636_552_486_1E-3)
+ 4.974_100_533_375_869E-5 * x))))))));
return (-0.5 * (x * x)).exp() * g;
}

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)))))));
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 @@ -44,23 +85,26 @@ fn inv_phi_tilde(phi_tilde_star: f64) -> f64 {
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))) /
(1.0 - g2 * (0.6635646938 - g2 * (0.14528712196 - 0.010472855461 * g2)));
let xi_bar = (0.032114372355
- g2 * (0.016969777977 - g2 * (0.002620733246 - 0.000096066952861 * g2)))
/ (1.0 - g2 * (0.6635646938 - g2 * (0.14528712196 - 0.010472855461 * g2)));
// Equation (2.3)
g * (ONE_OVER_SQRT_TWO_PI + xi_bar * g2)
} else {
// Equation (2.4)
let h = (-(-phi_tilde_star).ln()).sqrt();
// Equation (2.5)
(9.4883409779 - h * (9.6320903635 - h * (0.58556997323 + 2.1464093351 * h))) /
(1.0 - h * (0.65174820867 + h * (1.5120247828 + 0.000066437847132 * h)))
(9.4883409779 - h * (9.6320903635 - h * (0.58556997323 + 2.1464093351 * h)))
/ (1.0 - h * (0.65174820867 + h * (1.5120247828 + 0.000066437847132 * h)))
};
// Equation (2.7)
let q = (phi_tilde(x_bar) - phi_tilde_star) / norm_pdf(x_bar);
let x2 = x_bar * x_bar;
// Equation (2.6)
x_bar + 3.0 * q * x2 * (2.0 - q * x_bar * (2.0 + x2)) /
(6.0 + q * x_bar * (-12.0 + x_bar * (6.0 * q + x_bar * (-6.0 + q * x_bar * (3.0 + x2)))))
x_bar
+ 3.0 * q * x2 * (2.0 - q * x_bar * (2.0 + x2))
/ (6.0
+ q * x_bar * (-12.0 + x_bar * (6.0 * q + x_bar * (-6.0 + q * x_bar * (3.0 + x2)))))
}

/// Calculates the price of an option using Bachelier's model.
Expand All @@ -81,19 +125,29 @@ 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 moneyness = if q { forward - strike } else { strike - forward };
let moneyness = if q {
forward - strike
} else {
strike - forward
};
let x = moneyness / s;
s * phi_tilde_times_x(x)
}

pub(crate) fn implied_normal_volatility(price: f64, forward: f64, strike: f64, t: f64, q: bool) -> f64 {
pub(crate) fn implied_normal_volatility(
price: f64,
forward: f64,
strike: f64,
t: f64,
q: bool,
) -> f64 {
if forward == strike {
return price * SQRT_TWO_PI / t.sqrt();
}
let intrinsic = intrinsic_value(forward, strike, q);
match price.total_cmp(&intrinsic) {
Ordering::Less => { f64::NEG_INFINITY }
Ordering::Equal => { 0.0 }
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;
Expand All @@ -103,11 +157,10 @@ pub(crate) fn implied_normal_volatility(price: f64, forward: f64, strike: f64, t
}
}


#[cfg(test)]
mod tests {
use rand::Rng;
use super::*;
use rand::Rng;

#[test]
fn reconstruction_call_atm() {
Expand Down Expand Up @@ -191,7 +244,6 @@ mod tests {
}
}


#[test]
fn reconstruction_random_put_itm() {
let n = 100_000;
Expand Down Expand Up @@ -227,4 +279,4 @@ mod tests {
assert!((price - reprice).abs() <= 2.0 * f64::EPSILON);
}
}
}
}
2 changes: 0 additions & 2 deletions src/constants.rs
Original file line number Diff line number Diff line change
Expand Up @@ -17,5 +17,3 @@ 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;


4 changes: 1 addition & 3 deletions src/erf_cody.rs
Original file line number Diff line number Diff line change
Expand Up @@ -216,9 +216,7 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 {

#[cfg(test)]
mod tests {
use crate::erf_cody::{
erfc_cody, erfcx_cody, THRESH, XBIG, XHUGE, XMAX, XNEG,
};
use crate::erf_cody::{erfc_cody, erfcx_cody, THRESH, XBIG, XHUGE, XMAX, XNEG};

#[test]
fn calerf_1() {
Expand Down
45 changes: 34 additions & 11 deletions src/lib.rs
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,7 @@
//!
//! ```toml
//! [dependencies]
//! implied-vol = "0.2.2"
//! implied-vol = "1.0.0"
//! ```
//!
//! Then, in your code, bring the functions you need into scope with:
Expand All @@ -44,18 +44,17 @@
//!
//! ```toml
//! [dependencies]
//! implied-vol = { versions = "0.2.2", features = ["normal-distribution", "error-function"] }
//! implied-vol = { versions = "1.0.0", features = ["normal-distribution", "error-function"] }
//! ```
//!
//! For detailed explanations of each feature, please refer to the README.md file.

mod bachelier;
mod constants;
mod erf_cody;
mod lets_be_rational;
mod normal_distribution;
mod rational_cubic;
mod lets_be_rational;
mod bachelier;
mod constants;

/// Calculates the implied black volatility using a transformed rational guess with limited iterations.
///
Expand All @@ -78,7 +77,13 @@ mod constants;
/// assert_eq!(black_vol, 0.07011701801482094);
/// ```
#[inline]
pub fn implied_black_volatility(option_price: f64, forward: f64, strike: f64, expiry: f64, is_call: bool) -> f64 {
pub fn implied_black_volatility(
option_price: f64,
forward: f64,
strike: f64,
expiry: f64,
is_call: bool,
) -> f64 {
lets_be_rational::implied_black_volatility(option_price, forward, strike, expiry, is_call)
}

Expand All @@ -103,7 +108,13 @@ pub fn implied_black_volatility(option_price: f64, forward: f64, strike: f64, ex
/// assert!((price - 20.0).abs()<= 2.0 * f64::EPSILON * 20.0);
/// ```
#[inline]
pub fn calculate_european_option_price_by_black_scholes(forward: f64, strike: f64, volatility: f64, expiry: f64, is_call: bool) -> f64 {
pub fn calculate_european_option_price_by_black_scholes(
forward: f64,
strike: f64,
volatility: f64,
expiry: f64,
is_call: bool,
) -> f64 {
lets_be_rational::black(forward, strike, volatility, expiry, is_call)
}

Expand All @@ -127,7 +138,13 @@ pub fn calculate_european_option_price_by_black_scholes(forward: f64, strike: f6
/// let normal_vol = implied_vol::implied_normal_volatility(20.0, 100.0, 90.0, 30.0, true);
/// assert_eq!(normal_vol, 6.614292466299764);
/// ```
pub fn implied_normal_volatility(option_price: f64, forward: f64, strike: f64, expiry: f64, is_call: bool) -> f64 {
pub fn implied_normal_volatility(
option_price: f64,
forward: f64,
strike: f64,
expiry: f64,
is_call: bool,
) -> f64 {
bachelier::implied_normal_volatility(option_price, forward, strike, expiry, is_call)
}

Expand All @@ -152,7 +169,13 @@ pub fn implied_normal_volatility(option_price: f64, forward: f64, strike: f64, e
/// assert!((price - 20.0).abs()<= 2.0 * f64::EPSILON * 20.0);
/// ```
#[inline]
pub fn calculate_european_option_price_by_bachelier(forward: f64, strike: f64, volatility: f64, expiry: f64, is_call: bool) -> f64 {
pub fn calculate_european_option_price_by_bachelier(
forward: f64,
strike: f64,
volatility: f64,
expiry: f64,
is_call: bool,
) -> f64 {
bachelier::bachelier(forward, strike, volatility, expiry, is_call)
}

Expand Down Expand Up @@ -270,4 +293,4 @@ pub fn norm_cdf(x: f64) -> f64 {
#[inline]
pub fn inverse_norm_cdf(x: f64) -> f64 {
normal_distribution::inverse_norm_cdf(x)
}
}
34 changes: 18 additions & 16 deletions src/normal_distribution.rs
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
use std::f64::consts::FRAC_1_SQRT_2;
use crate::erf_cody::erfc_cody;
use std::f64::consts::FRAC_1_SQRT_2;

const NORM_CDF_ASYMPTOTIC_EXPANSION_FIRST_THRESHOLD: f64 = -10.0;
const NORM_CDF_ASYMPTOTIC_EXPANSION_SECOND_THRESHOLD: f64 = -67108864.0;
Expand Down Expand Up @@ -41,7 +41,6 @@ pub(crate) fn norm_cdf(z: f64) -> f64 {
0.5 * erfc_cody(-z * FRAC_1_SQRT_2)
}


pub(crate) fn inverse_norm_cdf(u: f64) -> f64 {
//
// ALGORITHM AS241 APPL. STATIST. (1988) VOL. 37, NO. 3
Expand Down Expand Up @@ -113,21 +112,24 @@ pub(crate) fn inverse_norm_cdf(u: f64) -> f64 {
let q = u - 0.5;
if q.abs() <= SPLIT1 {
let r = CONST1 - q * q;
q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0) /
(((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0)
q * (((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0)
/ (((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + 1.0)
} else {
let mut r = if q.is_sign_negative() { u } else { 1.0 - u };
r = (-r.ln()).sqrt();
let ret =
if r < SPLIT2 {
r -= CONST2;
(((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) * r + C0) /
(((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) * r + 1.0)
} else {
r -= SPLIT2;
(((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0) /
(((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + 1.0)
};
if q.is_sign_negative() { -ret } else { ret }
let ret = if r < SPLIT2 {
r -= CONST2;
(((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) * r + C0)
/ (((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) * r + 1.0)
} else {
r -= SPLIT2;
(((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0)
/ (((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + 1.0)
};
if q.is_sign_negative() {
-ret
} else {
ret
}
}
}
}
Loading

0 comments on commit 5cd21dd

Please sign in to comment.