diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index c61c508..f0ab0fc 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -19,4 +19,5 @@ jobs: # - run: cargo fmt --all --check - run: cargo check --workspace - run: cargo clippy --all-targets --all-features -- -D warnings -Aclippy::multiple_crate_versions + - run: cargo test --doc -- --show-output - run: cargo nextest run \ No newline at end of file diff --git a/LICENSE b/LICENSE index 7b04c08..4e88e47 100644 --- a/LICENSE +++ b/LICENSE @@ -20,7 +20,7 @@ LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -The externally maintained library from which parts of the Software is derived is: +The externally maintained libraries from which parts of the Software is derived are: - LetsBeRational, licensed as follows: @@ -35,3 +35,17 @@ The externally maintained library from which parts of the Software is derived is including without limitation any implied warranties of condition, uninterrupted use, merchantability, fitness for a particular purpose, or non-infringement. ========================================================================================== + + -　ImpliedNormalVolatility, licensed as follows: + + ========================================================================================= + Copyright © 2022 Peter Jäckel. + + Permission to use, copy, modify, and distribute this software is freely granted, + provided that this notice is preserved. + + WARRANTY DISCLAIMER + The Software is provided "as is" without warranty of any kind, either express or implied, + including without limitation any implied warranties of condition, uninterrupted use, + merchantability, fitness for a particular purpose, or non-infringement. + ========================================================================================== \ No newline at end of file diff --git a/README.md b/README.md index 263b79e..95319ec 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,42 @@ # Implied Vol -Welcome to `implied-vol`, a lightning-fast, pure Rust implementation of Peter Jäckel's renowned implied volatility -calculations. +[![Crates.io](https://img.shields.io/crates/v/implied-vol)](https://crates.io/crates/implied-vol) +[![Actions status](https://github.com/nakashima-hikaru/implied-vol/actions/workflows/ci.yaml/badge.svg)](https://github.com/nakashima-hikaru/implied-vol/actions) +[![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) -Originally found at [Peter Jäckel's homepage](http://www.jaeckel.org/), this library is a robust Rust reimplementation -of the pioneering work laid out in Jäckel's article, [Let's Be Rational](http://www.jaeckel.org/LetsBeRational.pdf). +More information about this crate can be found in the [crate documentation](https://docs.rs/implied-vol/0.1.0/implied_vol/). -Our aim is to offer a tool that provides blazingly fast and precise calculations of implied Black volatility. -With `implied-vol`, speed and accuracy in your quantitative finance projects need not be mutually exclusive. \ No newline at end of file +## About + +`implied-vol` is a high-performance, pure Rust implementation of Peter Jäckel's implied volatility calculations. This +library serves as a robust Rust reimplementation of the methodologies presented in Jäckel's works. + +## Source Works + +Our library follows the methods presented in two pivotal papers by Peter Jäckel: + +1. [Let's Be Rational](http://www.jaeckel.org/LetsBeRational.pdf): This work presents an approach to deduce Black’s + volatility from option prices with high precision. + +2. [Implied Normal Volatility](http://www.jaeckel.org/ImpliedNormalVolatility.pdf): Here, Jäckel provides an analytical + formula to calculate implied normal volatility (also known as Bachelier volatility) from vanilla option prices. + +Both resources can be accessed at [Peter Jäckel's homepage](http://www.jaeckel.org/). + +## Features + +`implied-vol` gains the benefits of being implemented in Rust, such as cross-compilation with Cargo and powerful static +analysis using Clippy lint. This library stands out by offering: + +- Rapid, precise calculations of both implied Black volatility and normal (Bachelier) implied volatility. +- Exceptional testability and maintainability due to its implementation in Rust. +- Unit tests aiding error checking. + +And most importantly, `implied-vol` follows the original C++ implementations closely, maintaining the same performance +characteristics and output precision. + +Community contributions are always welcome! + +## License + +This project is licensed under the [MIT license](https://github.com/nakashima-hikaru/implied-vol/blob/main/LICENSE). \ No newline at end of file diff --git a/src/bachelier.rs b/src/bachelier.rs new file mode 100644 index 0000000..e31c93a --- /dev/null +++ b/src/bachelier.rs @@ -0,0 +1,230 @@ +use crate::constants::{ONE_OVER_SQRT_TWO_PI, SQRT_TWO_PI}; +use crate::normal_distribution::norm_pdf; + +#[inline] +fn intrinsic_value(forward: f64, strike: f64, q: bool) -> f64 { + (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)))); + return ONE_OVER_SQRT_TWO_PI + x * (0.5 + x * g); + } + + if x > 0.0 { + return phi_tilde_times_x(-x) + x; + } + + 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)))))))); + return (-0.5 * (x * x)).exp() * g; + } + + let w = 1.0 / (x * x); + 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) +} + +fn phi_tilde(x: f64) -> f64 { + phi_tilde_times_x(x) / x +} + +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 { + 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 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))); + // 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))) + }; + // 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))))) +} + +/// Calculates the price of an option using Bachelier's model. +/// +/// # Arguments +/// +/// * `forward` - The forward price of the underlying asset. +/// * `strike` - The strike price of the option. +/// * `sigma` - The volatility of the underlying asset. +/// * `t` - The time to expiration of the option. +/// * `q` - A boolean flag indicating whether the option is a put (true) or a call (false). +/// +/// # Returns +/// +/// The price of the option. +pub(crate) fn bachelier(forward: f64, strike: f64, sigma: f64, t: f64, q: bool) -> f64 { + let s = sigma.abs() * t.sqrt(); + 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 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 { + if forward == strike { + 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; + } + 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() +} + + +#[cfg(test)] +mod tests { + use rand::Rng; + use super::*; + + #[test] + fn reconstruction_call_atm() { + for i in 1..100 { + let price = 0.01 * i as f64; + let f = 100.0; + let k = f; + let t = 1.0; + let q = true; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-14); + } + } + + #[test] + fn reconstruction_put_atm() { + for i in 1..100 { + let price = 0.01 * i as f64; + let f = 100.0; + let k = f; + let t = 1.0; + let q = false; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-14); + } + } + + #[test] + fn reconstruction_random_call_intrinsic() { + let n = 100_000; + let seed: [u8; 32] = [13; 32]; + let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); + for _ in 0..n { + let (r, r2, r3): (f64, f64, f64) = rng.gen(); + let price = 1e5 * r2; + let f = r + 1e5 * r2; + let k = f - price; + let t = 1e5 * r3; + let q = true; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-16); + } + } + + #[test] + fn reconstruction_random_call_itm() { + let n = 100_000; + let seed: [u8; 32] = [13; 32]; + let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); + for _ in 0..n { + let (r, r2, r3): (f64, f64, f64) = rng.gen(); + let price = 1.0 * (1.0 - r) + 1.0 * r * r2; + let f = 1.0; + let k = 1.0 * r; + let t = 1e5 * r3; + let q = true; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-16); + } + } + + #[test] + fn reconstruction_random_call_otm() { + let n = 100_000; + let seed: [u8; 32] = [13; 32]; + let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); + for _ in 0..n { + let (r, r2, r3): (f64, f64, f64) = rng.gen(); + let price = 1.0 * r * r2; + let f = 1.0 * r; + let k = 1.0; + let t = 1e5 * r3; + let q = true; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-16); + } + } + + + #[test] + fn reconstruction_random_put_itm() { + let n = 100_000; + let seed: [u8; 32] = [13; 32]; + let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); + for _ in 0..n { + let (r, r2, r3): (f64, f64, f64) = rng.gen(); + let price = 1.0 * r * r2; + let f = 1.0; + let k = 1.0 * r; + let t = 1e5 * r3; + let q = false; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-16); + } + } + + #[test] + fn reconstruction_random_put_otm() { + let n = 100_000; + let seed: [u8; 32] = [13; 32]; + let mut rng: rand::rngs::StdRng = rand::SeedableRng::from_seed(seed); + for _ in 0..n { + let (r, r2, r3): (f64, f64, f64) = rng.gen(); + let price = 1.0 * (1.0 - r) + 1.0 * r * r2; + let f = 1.0 * r; + let k = 1.0; + let t = 1e5 * r3; + let q = false; + let sigma = implied_normal_volatility(price, f, k, t, q); + let reprice = bachelier(f, k, sigma, t, q); + assert!((price - reprice).abs() < 5e-16); + } + } +} \ No newline at end of file diff --git a/src/constants.rs b/src/constants.rs new file mode 100644 index 0000000..8dd886d --- /dev/null +++ b/src/constants.rs @@ -0,0 +1,19 @@ +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_DBL_MAX: 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; \ No newline at end of file diff --git a/src/erf_cody.rs b/src/erf_cody.rs index c8f79ea..db75b74 100644 --- a/src/erf_cody.rs +++ b/src/erf_cody.rs @@ -63,73 +63,6 @@ const XBIG: f64 = 26.543; const XHUGE: f64 = 6.71e7; const XMAX: f64 = 2.53e307; -// pub(crate) fn erf_cody(x: f64) -> f64 { -// /* -------------------------------------------------------------------- */ -// /* This subprogram computes approximate values for erf(x). */ -// /* (see comments heading CALERF). */ -// /* Author/date: W. J. Cody, January 8, 1985 */ -// /* -------------------------------------------------------------------- */ -// let y = x.abs(); -// let mut ysq = ZERO; -// let mut xden; -// let mut xnum; -// let mut result = ZERO; -// -// if y <= THRESH { -// if y > XSMALL { -// ysq = y * y; -// } -// xnum = A[4] * ysq; -// xden = ysq; -// -// for i in 0..3 { -// xnum = (xnum + A[i]) * ysq; -// xden = (xden + B[i]) * ysq; -// } -// result = x * (xnum + A[3]) / (xden + B[3]); -// return result; -// } else if y <= FOUR { -// xnum = C[8] * y; -// xden = y; -// -// for i in 0..7 { -// xnum = (xnum + C[i]) * y; -// xden = (xden + D[i]) * y; -// } -// result = (xnum + C[7]) / (xden + D[7]); -// -// ysq = (y * SIXTEN).trunc() / SIXTEN; -// let del = (y - ysq) * (y + ysq); -// result *= (-ysq * ysq).exp() * (-del).exp(); -// } else if y >= XBIG { -// result = (HALF - result) + HALF; -// if x < ZERO { -// result = -result; -// } -// return result; -// } else { -// ysq = ONE / (y * y); -// xnum = P[5] * ysq; -// xden = ysq; -// -// for i in 0..4 { -// xnum = (xnum + P[i]) * ysq; -// xden = (xden + Q[i]) * ysq; -// } -// result = ysq * (xnum + P[4]) / (xden + Q[4]); -// result = (SQRPI - result) / y; -// -// ysq = (y * SIXTEN).trunc() / SIXTEN; -// let del = (y - ysq) * (y + ysq); -// result *= (-ysq * ysq).exp() * (-del).exp(); -// } -// result = (HALF - result) + HALF; -// if x < ZERO { -// result = -result; -// } -// result -// } - pub(crate) fn erfc_cody(x: f64) -> f64 { /* -------------------------------------------------------------------- */ /* This subprogram computes approximate values for erfc(x). */ @@ -197,7 +130,6 @@ pub(crate) fn erfc_cody(x: f64) -> f64 { result } - pub(crate) fn erfcx_cody(x: f64) -> f64 { /* ------------------------------------------------------------------ */ /* This subprogram computes approximate values for exp(x*x) * erfc(x). */ @@ -288,71 +220,12 @@ pub(crate) fn erfcx_cody(x: f64) -> f64 { result } - #[cfg(test)] mod tests { use crate::erf_cody::{ erfc_cody, erfcx_cody, FOUR, THRESH, XBIG, XHUGE, XMAX, XNEG, ZERO, }; - // #[test] - // fn calerf_0() { - // let x = erf_cody(THRESH + f64::EPSILON); - // assert_eq!(x, 0.49261347321793825); - // let x = erf_cody(THRESH - f64::EPSILON); - // assert_eq!(x, 0.49261347321793775); - // let x = erf_cody(-THRESH - f64::EPSILON); - // assert_eq!(x, -0.49261347321793825); - // let x = erf_cody(-THRESH + f64::EPSILON); - // assert_eq!(x, -0.49261347321793775); - // - // let x = erf_cody(FOUR + f64::EPSILON); - // assert_eq!(x, 0.9999999845827421); - // let x = erf_cody(FOUR - f64::EPSILON); - // assert_eq!(x, 0.9999999845827421); - // let x = erf_cody(-FOUR - f64::EPSILON); - // assert_eq!(x, -0.9999999845827421); - // let x = erf_cody(-FOUR + f64::EPSILON); - // assert_eq!(x, -0.9999999845827421); - // - // let x = erf_cody(XBIG + f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(XBIG - f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(-XBIG - f64::EPSILON); - // assert_eq!(x, -1.0); - // let x = erf_cody(-XBIG + f64::EPSILON); - // assert_eq!(x, -1.0); - // - // let x = erf_cody(XMAX + f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(XMAX - f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(-XMAX - f64::EPSILON); - // assert_eq!(x, -1.0); - // let x = erf_cody(-XMAX + f64::EPSILON); - // assert_eq!(x, -1.0); - // - // let x = erf_cody(XHUGE + f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(XHUGE - f64::EPSILON); - // assert_eq!(x, 1.0); - // let x = erf_cody(-XHUGE - f64::EPSILON); - // assert_eq!(x, -1.0); - // let x = erf_cody(-XHUGE + f64::EPSILON); - // assert_eq!(x, -1.0); - // - // let x = erf_cody(ZERO + f64::EPSILON); - // assert_eq!(x, 2.5055050636335897e-16); - // let x = erf_cody(ZERO - f64::EPSILON); - // assert_eq!(x, -2.5055050636335897e-16); - // - // let x = erf_cody(XNEG + f64::EPSILON); - // assert_eq!(x, -1.0); - // let x = erf_cody(XNEG - f64::EPSILON); - // assert_eq!(x, -1.0); - // } - #[test] fn calerf_1() { let x = erfc_cody(THRESH + f64::EPSILON); diff --git a/src/lets_be_rational.rs b/src/lets_be_rational.rs index 78d7428..18451ca 100644 --- a/src/lets_be_rational.rs +++ b/src/lets_be_rational.rs @@ -1,52 +1,35 @@ use std::f64::consts::SQRT_2; +use crate::constants::{DENORMALISATION_CUTOFF, FOURTH_ROOT_DBL_EPSILON, LN_TWO_PI, PI_OVER_SIX, SIXTEENTH_ROOT_DBL_EPSILON, SQRT_DBL_MAX, SQRT_DBL_MIN, SQRT_ONE_OVER_THREE, SQRT_PI_OVER_TWO, SQRT_THREE, SQRT_THREE_OVER_THIRD_ROOT_TWO_PI, SQRT_TWO_PI, 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 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}; -const TWO_PI: f64 = std::f64::consts::TAU; -const SQRT_PI_OVER_TWO: f64 = 1.253_314_137_315_500_3; -const SQRT_THREE: f64 = 1.732_050_807_568_877_2; -const SQRT_ONE_OVER_THREE: f64 = 0.577_350_269_189_625_7; -const TWO_PI_OVER_SQRT_TWENTY_SEVEN: f64 = 1.209_199_576_156_145_2; -const SQRT_THREE_OVER_THIRD_ROOT_TWO_PI: f64 = 0.938_643_487_427_383_6; -const PI_OVER_SIX: f64 = std::f64::consts::FRAC_PI_6; - -// const SQRT_DBL_EPSILON: f64 = 1.4901161193847656e-8; -// const FOURTH_ROOT_DBL_EPSILON: f64 = 0.0001220703125; -// const EIGHTH_ROOT_DBL_EPSILON: f64 = 0.011048543456039806; -const SIXTEENTH_ROOT_DBL_EPSILON: f64 = 0.10511205190671433; -const SQRT_DBL_MIN: f64 = 1.4916681462400413e-154; -const SQRT_DBL_MAX: f64 = 1.3407807929942596e154; -const DENORMALISATION_CUTOFF: f64 = 0.0; - -const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC: f64 = f64::NEG_INFINITY; -const VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM: f64 = f64::INFINITY; - - -fn is_below_horizon(x: f64) -> bool { x.abs() < DENORMALISATION_CUTOFF } +#[inline] fn householder3_factor(v: f64, h2: f64, h3: f64) -> f64 { (1.0 + 0.5 * h2 * v) / (1.0 + v * (h2 + h3 * v / 6.0)) } +#[inline] fn householder4_factor(v: f64, h2: f64, h3: f64, h4: f64) -> f64 { (1.0 + v * (h2 + v * h3 / 6.0)) / (1.0 + v * (1.5 * h2 + v * (h2 * h2 / 4.0 + h3 / 3.0 + v * h4 / 24.0))) } -fn normalised_intrinsic(x: f64, q: f64 /* q=±1 */) -> f64 { - if q * x <= 0.0 { +fn normalised_intrinsic(x: f64, q: bool) -> f64 { + if (q && (x <= 0.0)) || (!q && (x >= 0.0)) { return 0.0; } let x2 = x * x; - const FOURTH_ROOT_DBL_EPSILON: f64 = 0.0001; // Just assuming a value here. Please replace with actual constant. if x2 < 98.0 * FOURTH_ROOT_DBL_EPSILON { - return ((q < 0.0) as i32 as f64 * -2.0 + 1.0) * x * (1.0 + x2 * ((1.0 / 24.0) + x2 * ((1.0 / 1920.0) + x2 * ((1.0 / 322560.0) + (1.0 / 92897280.0) * x2)))).abs().max(0.0); + return (if q { 1.0 } else { -1.0 }) * x * (1.0 + x2 * ((1.0 / 24.0) + x2 * ((1.0 / 1920.0) + x2 * ((1.0 / 322560.0) + (1.0 / 92897280.0) * x2)))).abs().max(0.0); } let b_max = (0.5 * x).exp(); - let one_over_b_max = 1.0 / b_max; - ((q < 0.0) as i32 as f64 * -2.0 + 1.0) * (b_max - one_over_b_max).abs().max(0.0) + let one_over_b_max = b_max.recip(); + (if q { 1.0 } else { -1.0 }) * (b_max - one_over_b_max).abs().max(0.0) } +#[inline] fn normalised_intrinsic_call(x: f64) -> f64 { - normalised_intrinsic(x, 1.0) + normalised_intrinsic(x, true) } +#[inline] fn square(x: f64) -> f64 { x * x } @@ -56,9 +39,9 @@ const SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD: f64 = 2.0 * SIXTEENTH_ROO fn asymptotic_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f64 { assert!((h < -f64::abs(ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD)) && (h + t < -f64::abs(SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD + ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD))); - let e = (t / h).powi(2); + let e = square(t / h); let r = (h + t) * (h - t); - let q = (h / r).powi(2); + let q = square(h / r); let asymptotic_expansion_sum = 2.0 + q * (-6.0E0 - 2.0 * e + 3.0 * q * (1.0E1 + e * (2.0E1 + 2.0 * e) + 5.0 * q * (-1.4E1 + e * (-7.0E1 + e * (-4.2E1 - 2.0 * e)) + 7.0 * q * (1.8E1 + e * (1.68E2 + e * (2.52E2 + e * (7.2E1 + 2.0 * e))) + 9.0 * q * (-2.2E1 + e * (-3.3E2 + e * (-9.24E2 + e * (-6.6E2 + e * (-1.1E2 - 2.0 * e)))) + 1.1E1 * q * (2.6E1 + e * (5.72E2 + e * (2.574E3 + e * (3.432E3 + e * (1.43E3 + e * (1.56E2 + 2.0 * e))))) + 1.3E1 * q * (-3.0E1 + e * (-9.1E2 + e * (-6.006E3 + e * (-1.287E4 + e * (-1.001E4 + e * (-2.73E3 + e * (-2.1E2 - 2.0 * e)))))) + 1.5E1 * q * (3.4E1 + e * (1.36E3 + e * (1.2376E4 + e * (3.8896E4 + e * (4.862E4 + e * (2.4752E4 + e * (4.76E3 + e * (2.72E2 + 2.0 * e))))))) + 1.7E1 * q * (-3.8E1 + e * (-1.938E3 + e * (-2.3256E4 + e * (-1.00776E5 + e * (-1.84756E5 + e * (-1.51164E5 + e * (-5.4264E4 + e * (-7.752E3 + e * (-3.42E2 - 2.0 * e)))))))) + 1.9E1 * q * (4.2E1 + e * (2.66E3 + e * (4.0698E4 + e * (2.3256E5 + e * (5.8786E5 + e * (7.05432E5 + e * (4.0698E5 + e * (1.08528E5 + e * (1.197E4 + e * (4.2E2 + 2.0 * e))))))))) + 2.1E1 * q * (-4.6E1 + e * (-3.542E3 + e * (-6.7298E4 + e * (-4.90314E5 + e * (-1.63438E6 + e * (-2.704156E6 + e * (-2.288132E6 + e * (-9.80628E5 + e * (-2.01894E5 + e * (-1.771E4 + e * (-5.06E2 - 2.0 * e)))))))))) + 2.3E1 * q * (5.0E1 + e * (4.6E3 + e * (1.0626E5 + e * (9.614E5 + e * (4.08595E6 + e * (8.9148E6 + e * (1.04006E7 + e * (6.53752E6 + e * (2.16315E6 + e * (3.542E5 + e * (2.53E4 + e * (6.0E2 + 2.0 * e))))))))))) + 2.5E1 * q * (-5.4E1 + e * (-5.85E3 + e * (-1.6146E5 + e * (-1.77606E6 + e * (-9.37365E6 + e * (-2.607579E7 + e * (-4.01166E7 + e * (-3.476772E7 + e * (-1.687257E7 + e * (-4.44015E6 + e * (-5.9202E5 + e * (-3.51E4 + e * (-7.02E2 - 2.0 * e)))))))))))) + 2.7E1 * q * (5.8E1 + e * (7.308E3 + e * (2.3751E5 + e * (3.12156E6 + e * (2.003001E7 + e * (6.919458E7 + e * (1.3572783E8 + e * (1.5511752E8 + e * (1.0379187E8 + e * (4.006002E7 + e * (8.58429E6 + e * (9.5004E5 + e * (4.7502E4 + e * (8.12E2 + 2.0 * e))))))))))))) + 2.9E1 * q * (-6.2E1 + e * (-8.99E3 + e * (-3.39822E5 + e * (-5.25915E6 + e * (-4.032015E7 + e * (-1.6934463E8 + e * (-4.1250615E8 + e * (-6.0108039E8 + e * (-5.3036505E8 + e * (-2.8224105E8 + e * (-8.870433E7 + e * (-1.577745E7 + e * (-1.472562E6 + e * (-6.293E4 + e * (-9.3E2 - 2.0 * e)))))))))))))) + 3.1E1 * q * (6.6E1 + e * (1.0912E4 + e * (4.74672E5 + e * (8.544096E6 + e * (7.71342E7 + e * (3.8707344E8 + e * (1.14633288E9 + e * (2.07431664E9 + e * (2.33360622E9 + e * (1.6376184E9 + e * (7.0963464E8 + e * (1.8512208E8 + e * (2.7768312E7 + e * (2.215136E6 + e * (8.184E4 + e * (1.056E3 + 2.0 * e))))))))))))))) + 3.3E1 * (-7.0E1 + e * (-1.309E4 + e * (-6.49264E5 + e * (-1.344904E7 + e * (-1.4121492E8 + e * (-8.344518E8 + e * (-2.9526756E9 + e * (-6.49588632E9 + e * (-9.0751353E9 + e * (-8.1198579E9 + e * (-4.6399188E9 + e * (-1.6689036E9 + e * (-3.67158792E8 + e * (-4.707164E7 + e * (-3.24632E6 + e * (-1.0472E5 + e * (-1.19E3 - 2.0 * e))))))))))))))))) * q)))))))))))))))); @@ -66,29 +49,14 @@ fn asymptotic_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f6 f64::abs(f64::max(b_over_vega, 0.)) } - -// fn normalised_black_call_using_erfcx(h: f64, t: f64) -> f64 { -// let b = 0.5 * (-0.5 * (h * h + t * t)).exp() * (erfcx_cody(-(1.0 / SQRT_2) * (h + t)) - erfcx_cody(-(1.0 / SQRT_2) * (h - t))); -// b.abs().max(0.0) -// } - fn small_t_expansion_of_normalised_black_call_over_vega(h: f64, t: f64) -> f64 { - let w = t.powi(2); - let h2 = h.powi(2); + let w = t * t; + let h2 = h * h; let a = 1_f64 + h * SQRT_PI_OVER_TWO * erfcx_cody(-(1.0 / 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)))))); f64::abs(f64::max(b_over_vega, 0_f64)) } -// fn normalised_black_call_using_norm_cdf(x: f64, s: f64) -> f64 { -// let h = x / s; -// let t = 0.5 * s; -// let b_max = (0.5 * x).exp(); -// let b = norm_cdf(h + t) * b_max - norm_cdf(h - t) / b_max; -// -// b.max(0.0).abs() -// } - fn normalised_black_call_with_optimal_use_of_codys_functions(x: f64, s: f64) -> f64 { const CODYS_THRESHOLD: f64 = 0.46875; let h = x / s; @@ -110,9 +78,6 @@ fn normalised_black_call_with_optimal_use_of_codys_functions(x: f64, s: f64) -> two_b.abs().max(0.0) } -const SQRT_TWO_PI: f64 = 2.5066282746310002; -const LN_TWO_PI: f64 = 1.8378770664093453; - fn normalised_vega(x: f64, s: f64) -> f64 { let ax = x.abs(); if ax <= 0.0 { @@ -120,7 +85,7 @@ fn normalised_vega(x: f64, s: f64) -> f64 { } else if s <= 0.0 || s <= ax * SQRT_DBL_MIN { 0.0 } else { - (1.0 / SQRT_TWO_PI) * (-0.5 * ((x / s).powi(2) + (0.5 * s).powi(2))).exp() + (1.0 / SQRT_TWO_PI) * (-0.5 * (square(x / s) + square(0.5 * s))).exp() } } @@ -129,9 +94,9 @@ fn ln_normalised_vega(x: f64, s: f64) -> f64 { if ax <= 0.0 { -(LN_TWO_PI / 2.0) - 0.125 * s * s } else if s <= 0.0 { - -f64::MAX + f64::MIN } else { - -(LN_TWO_PI / 2.0) - 0.5 * ((x / s).powi(2) + (0.5 * s).powi(2)) + -(LN_TWO_PI / 2.0) - 0.5 * (square(x / s) + square(0.5 * s)) } } @@ -142,7 +107,7 @@ fn normalised_black_call(x: f64, s: f64) -> f64 { if s <= x.abs() * DENORMALISATION_CUTOFF { return normalised_intrinsic_call(x); } - if x < s * ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD && (0.5 * s).powi(2) + x < s * (SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD + ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD) { + if x < s * ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD && square(0.5 * s) + x < s * (SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD + ASYMPTOTIC_EXPANSION_ACCURACY_THRESHOLD) { return asymptotic_expansion_of_normalised_black_call_over_vega(x / s, 0.5 * s) * normalised_vega(x, s); } if 0.5 * s < SMALL_T_EXPANSION_OF_NORMALISED_BLACK_THRESHOLD { @@ -170,65 +135,45 @@ fn normalised_black_call_over_vega_and_ln_vega(x: f64, s: f64) -> (f64, f64) { (normalised_black_call_with_optimal_use_of_codys_functions(x, s) * (-ln_vega).exp(), ln_vega) } -// fn normalised_volga(x: f64, s: f64) -> f64 { -// let ax = x.abs(); -// if ax <= 0f64 { return (1f64 / SQRT_TWO_PI) * (-0.125 * s * s).exp(); } -// if s <= 0f64 || s <= ax * SQRT_DBL_MIN { return 0f64; } -// let h2 = x.powi(2) / s.powi(2); -// let t2 = (0.5 * s).powi(2); -// return (1f64 / SQRT_TWO_PI) * (-0.5 * (h2 + t2)).exp() * (h2 - t2) / s; -// } - -fn normalised_black(x: f64, s: f64, theta: f64) -> f64 { - normalised_black_call(if theta < 0f64 { -x } else { x }, s) +#[inline] +fn normalised_black(x: f64, s: f64, theta: bool) -> f64 { + normalised_black_call(if !theta { -x } else { x }, s) } -/// Calculates the price of a European option using the Black-Scholes formula. -/// -/// # Arguments -/// -/// * `f` - The current value of the underlying asset. -/// * `k` - The strike price of the option. -/// * `sigma` - The volatility of the underlying asset. -/// * `t` - The time to expiration of the option. -/// * `q` - Positive if call else put. -/// -/// # Returns -/// -/// The price of the European option. -pub fn black(f: f64, k: f64, sigma: f64, t: f64, q: f64) -> f64 { - let intrinsic = if q < 0f64 { k - f } else { f - k }.max(0f64).abs(); - if q * (f - k) > 0f64 { - return intrinsic + black(f, k, sigma, t, -q); +pub(crate) fn black(f: f64, k: f64, sigma: f64, t: f64, q: bool) -> f64 { + let intrinsic = if !q { k - f } else { f - k }.max(0f64).abs(); + if (q && ((f - k) > 0.0)) || (!q && ((f - k) < 0.0)) { + return intrinsic + black(f, k, sigma, t, !q); } intrinsic.max((f.sqrt() * k.sqrt()) * normalised_black((f / k).ln(), sigma * t.sqrt(), q)) } -fn compute_f_lower_map_and_first_two_derivatives(x: f64, s: f64, f: &mut f64, fp: &mut f64, fpp: &mut 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 y = z * z; let s2 = s * s; let phi_m = norm_cdf(-z); let phi = norm_pdf(z); - *fpp = PI_OVER_SIX * y / (s2 * s) * phi_m * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_m / phi) * (2.0 * y + 0.25 * s2).exp(); - (*fp, *f) = if is_below_horizon(s) { + let fpp = PI_OVER_SIX * y / (s2 * s) * phi_m * (8.0 * SQRT_THREE * s * ax + (3.0 * s2 * (s2 - 8.0) - 8.0 * x * x) * phi_m / phi) * (2.0 * y + 0.25 * s2).exp(); + let (fp, f) = if s.is_subnormal() { (1.0, 0.0) } else { let phi2 = phi_m * phi_m; let fp_val = TWO_PI * y * phi2 * (y + 0.125 * s * s).exp(); - let f_val = if is_below_horizon(x) { + let f_val = if x.is_subnormal() { 0.0 } else { TWO_PI_OVER_SQRT_TWENTY_SEVEN * ax * (phi2 * phi_m) }; (fp_val, f_val) }; + (f, fp, fpp) } fn inverse_f_lower_map(x: f64, f: f64) -> f64 { - if is_below_horizon(f) { + if f.is_subnormal() { 0.0 } else { (x / (SQRT_THREE * inverse_norm_cdf(SQRT_THREE_OVER_THIRD_ROOT_TWO_PI * f.cbrt() / x.abs().cbrt()))).abs() @@ -239,11 +184,11 @@ fn compute_f_upper_map_and_first_two_derivatives(x: f64, s: f64) -> (f64, f64, f let f = norm_cdf(-0.5 * s); let (fp, fpp); - if is_below_horizon(x) { + if x.is_subnormal() { fp = -0.5; fpp = 0.0; } else { - let w = (x / s).powi(2); + let w = square(x / s); fp = -0.5 * (0.5 * w).exp(); fpp = SQRT_PI_OVER_TWO * ((w + 0.125 * s * s).exp()) * w / s; } @@ -256,26 +201,19 @@ fn inverse_f_upper_map(f: f64) -> f64 { -2.0 * inverse_norm_cdf(f) } -fn take_step(x_min: f64, x_max: f64, x: f64, dx: &mut f64) -> f64 { - let new_x = x_min.max(x_max.min(x + *dx)); - *dx = new_x - x; - new_x +fn take_step(x_min: f64, x_max: f64, x: f64, dx: f64) -> (f64, f64) { + let new_x = x_min.max(x_max.min(x + dx)); + (new_x, new_x - x) } -// pub fn complementary_normalised_black_inner(h: f64, t: f64) -> f64 { -// 0.5 * (erfcx_cody((t + h) * (1.0 / SQRT_2)) + erfcx_cody((t - h) * (1.0 / SQRT_2))) * (-0.5 * (t * t + h * h)).exp() -// } - fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( - mut beta: f64, mut x: f64, q: f64, n: i32, + mut beta: f64, mut x: f64, q: bool, n: i32, ) -> f64 { - if q * x > 0. { + if (q && (x > 0.0)) || (!q && (x < 0.0)) { beta = (beta - normalised_intrinsic(x, q)).max(0.).abs(); - // q = -q; } - if q < 0. { + if !q { x = -x; - // q = -q; } assert!(x <= 0.); if beta <= 0. || beta < DENORMALISATION_CUTOFF { @@ -286,9 +224,9 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM; } let mut iterations = 0; - let mut f = -f64::MAX; + let mut f = f64::MIN; let mut s; - let mut ds = -f64::MAX; + let mut ds = f64::MIN; let mut s_left = f64::MIN; let mut s_right = f64::MAX; let s_c = (2. * x.abs()).sqrt(); @@ -298,10 +236,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let s1 = s_c - b_c / v_c; let b1 = normalised_black_call(x, s1); if beta < b1 { - let mut f_lower_map_l: f64 = 0.0; - let mut d_f_lower_map_l_d_beta: f64 = 0.0; - let mut d2_f_lower_map_l_d_beta2: f64 = 0.0; - compute_f_lower_map_and_first_two_derivatives(x, s1, &mut f_lower_map_l, &mut d_f_lower_map_l_d_beta, &mut d2_f_lower_map_l_d_beta2); + let (f_lower_map_l, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2) = compute_f_lower_map_and_first_two_derivatives(x, s1); let r2 = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_right_side(0.0, b1, 0.0, f_lower_map_l, 1.0, d_f_lower_map_l_d_beta, d2_f_lower_map_l_d_beta2, true); f = rational_cubic_interpolation(beta, 0.0, b1, 0.0, f_lower_map_l, 1.0, d_f_lower_map_l_d_beta, r2); if f <= 0.0 { @@ -323,7 +258,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let lambda = 1.0 / ln_b; let otlambda = 1.0 + 2.0 * lambda; let h2 = b_h2 - bpob * otlambda; - let c = 3.0 * (h / s).powi(2); + let c = 3.0 * square(h / s); let b_h3 = b_h2 * b_h2 - c - 0.25; let sq_bpob = bpob * bpob; let mu = 6.0 * lambda * (1.0 + lambda); @@ -333,7 +268,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit } else { nu * householder3_factor(nu, h2, h3) }; - s = take_step(s_left, s_right, s, &mut ds); + (s, ds) = take_step(s_left, s_right, s, ds); iterations += 1; } return s; @@ -355,11 +290,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit s_left = s_c; s_right = s_u; } else { - let f_upper_map_h: f64; - let d_f_upper_map_h_d_beta: f64; - let d2_f_upper_map_h_d_beta2: f64; - - (f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2) = compute_f_upper_map_and_first_two_derivatives(x, s_u); + let (f_upper_map_h, d_f_upper_map_h_d_beta, d2_f_upper_map_h_d_beta2) = compute_f_upper_map_and_first_two_derivatives(x, s_u); if d2_f_upper_map_h_d_beta2 > -SQRT_DBL_MAX && d2_f_upper_map_h_d_beta2 < SQRT_DBL_MAX { let r_uu = convex_rational_cubic_control_parameter_to_fit_second_derivative_at_left_side(b_u, b_max, f_upper_map_h, 0.0, d_f_upper_map_h_d_beta, -0.5, d2_f_upper_map_h_d_beta2, true); @@ -391,7 +322,7 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit } else { nu * householder3_factor(nu, h2, h3) }; - s = take_step(s_left, s_right, s, &mut ds); + (s, ds) = take_step(s_left, s_right, s, ds); iterations += 1; } return s; @@ -408,10 +339,10 @@ fn unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_wit let nu = (beta - b) / bp; let h = x / s; let h2 = (h * h) / s - s / 4.0; - let h3 = h2 * h2 - 3.0 * (h / s).powi(2) - 0.25; + let h3 = h2 * h2 - 3.0 * square(h / s) - 0.25; ds = nu * householder3_factor(nu, h2, h3); // Never leave the branch (or bracket) - s = take_step(s_left, s_right, s, &mut ds); + (s, ds) = take_step(s_left, s_right, s, ds); } s } @@ -421,24 +352,24 @@ fn implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( f: f64, k: f64, t: f64, - mut q: f64, // q=±1 + mut q: bool, n: i32, ) -> f64 { - let intrinsic = f64::abs(f64::max(if q < 0.0 { k - f } else { f - k }, 0.0)); + let intrinsic = f64::abs(f64::max(if !q { k - f } else { f - k }, 0.0)); if price < intrinsic { return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC; } - let max_price = if q < 0.0 { k } else { f }; + let max_price = if !q { k } else { f }; if price >= max_price { return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_ABOVE_MAXIMUM; } let x = (f / k).ln(); // Map in-the-money to out-of-the-money - if q * x > 0.0 { + if (q && (x > 0.0)) || (!q && (x < 0.0)) { // q && x > 0.0 or !q && x < 0.0 price = (price - intrinsic).max(0.0).abs(); - q = -q; + q = !q; } unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( @@ -449,67 +380,11 @@ fn implied_volatility_from_a_transformed_rational_guess_with_limited_iterations( ) / t.sqrt() } -// fn complementary_normalised_black(x: f64, s: f64) -> f64 { -// complementary_normalised_black_inner(x / s, s / 2.0) -// } - -// fn normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(mut beta: f64, x: f64, mut q: f64 /* q=±1 */, n: i32) -> f64 { -// // Map in-the-money to out-of-the-money -// if q * x > 0.0 { -// beta -= normalised_intrinsic(x, q); -// q = -q; -// } -// if beta < 0.0 { -// return VOLATILITY_VALUE_TO_SIGNAL_PRICE_IS_BELOW_INTRINSIC; -// } -// return unchecked_normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(beta, x, q, n); -// } - -/// Calculates the implied black volatility using a transformed rational guess with limited iterations. -/// -/// # Arguments -/// -/// * `price` - The current price of the option. -/// * `f` - The current forward price of the underlying asset. -/// * `k` - The strike price of the option. -/// * `t` - The time to expiration in years. -/// * `q` - Positive if call else put. -/// -/// # Returns -/// -/// The implied black volatility. -pub fn implied_black_volatility(price: f64, f: f64, k: f64, t: f64, q: f64) -> f64 { +pub(crate) fn implied_black_volatility(price: f64, f: f64, k: f64, t: f64, q: bool) -> f64 { implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(price, f, k, t, q, 2) } -// fn normalised_implied_black_volatility(beta: f64, x: f64, q: f64 /* q=±1 */) -> f64 { -// normalised_implied_volatility_from_a_transformed_rational_guess_with_limited_iterations(beta, x, q, 2) -// } - -// fn vega(f: f64, k: f64, sigma: f64, t: f64) -> f64 { -// (f.sqrt() * k.sqrt()) * normalised_vega((f / k).ln(), sigma * t.sqrt()) * t.sqrt() -// } - - -// fn volga(f: f64, k: f64, sigma: f64, t: f64) -> f64 { -// (f.sqrt() * k.sqrt()) * normalised_volga((f / k).ln(), sigma * t.sqrt()) * t -// } - - -// fn black_accuracy_factor(x: f64, s: f64, theta: f64) -> f64 { -// s / normalised_black_call_over_vega_and_ln_vega(if theta < 0.0 { -x } else { x }, s).0 -// } -// fn implied_volatility_attainable_accuracy(x: f64, s: f64, theta: f64) -> f64 { -// let (bx, ln_vega) = normalised_black_call_over_vega_and_ln_vega(if theta < 0.0 { -x } else { x }, s); -// let b = bx * ln_vega.exp(); -// let b_max = (if theta < 0.0 { -x } else { x } / 2.0).exp(); -// if b > f64::MIN && b < b_max { -// f64::EPSILON * (1.0 + (bx / s).abs()) -// } else { -// 1.0 -// } -// } #[cfg(test)] mod tests { use rand::Rng; @@ -522,7 +397,7 @@ mod tests { let f = 100.0; let k = f; let t = 1.0; - let q = 1.0; + let q = true; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-14); @@ -536,7 +411,7 @@ mod tests { let f = 100.0; let k = f; let t = 1.0; - let q = -1.0; + let q = false; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-14); @@ -554,7 +429,7 @@ mod tests { let f = r + 1e5 * r2; let k = f - price; let t = 1e5 * r3; - let q = 1.0; + let q = true; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-16); @@ -572,7 +447,7 @@ mod tests { let f = 1.0; let k = 1.0 * r; let t = 1e5 * r3; - let q = 1.0; + let q = true; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-16); @@ -590,7 +465,7 @@ mod tests { let f = 1.0 * r; let k = 1.0; let t = 1e5 * r3; - let q = 1.0; + let q = true; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-16); @@ -609,7 +484,7 @@ mod tests { let f = 1.0; let k = 1.0 * r; let t = 1e5 * r3; - let q = -1.0; + let q = false; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-16); @@ -627,7 +502,7 @@ mod tests { let f = 1.0 * r; let k = 1.0; let t = 1e5 * r3; - let q = -1.0; + let q = false; let sigma = implied_black_volatility(price, f, k, t, q); let reprice = black(f, k, sigma, t, q); assert!((price - reprice).abs() < 5e-16); diff --git a/src/lib.rs b/src/lib.rs index f4043ef..bea5be0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,4 +1,152 @@ +//! This module provides a Rust implementation of Peter Jäckel's original C++ code for calculating +//! implied volatilities in financial derivatives. +//! To learn more about the algorithms, please refer to Peter Jäckel's papers [Let's Be Rational](http://www.jaeckel.org/LetsBeRational.pdf) and [Implied Normal Volatility](http://www.jaeckel.org/ImpliedNormalVolatility.pdf). +//! +//! # Features +//! +//! - Calculation of implied Black volatility +//! - Calculation of the price of a European option using the Black-Scholes model +//! - Calculation of implied normal volatility +//! - Calculation of the price of an option using Bachelier's model +//! +//! All models support both call and put options. +//! +//! # Examples +//! +//! Check out the documentation for each function for practical examples. +//! +//! # Usage +//! +//! Import the crate in your Rust project by adding the following to your `Cargo.toml` +//! +//! ```toml +//! [dependencies] +//! implied-vol = "0.1.0" +//! ``` +//! +//! Then, in your code, bring the functions you need into scope with: +//! +//! ```rust +//! use implied_vol::{implied_black_volatility, calculate_european_option_price_by_black_scholes, implied_normal_volatility, calculate_european_option_price_by_bachelier}; +//! +//! let black_vol = implied_black_volatility(20.0, 100.0, 90.0, 30.0, true); +//! assert_eq!(black_vol, 0.07011701801482094); +//! +//! let price = calculate_european_option_price_by_black_scholes(100.0, 90.0, 0.07011701801482094, 30.0, true); +//! assert!((price - 20.0).abs() < 5e-16 * 20.0); +//! +//! let normal_vol = implied_normal_volatility(20.0, 100.0, 90.0, 30.0, true); +//! assert_eq!(normal_vol, 6.614292466299764); +//! +//! let price = calculate_european_option_price_by_bachelier(100.0, 90.0, 6.614292466299764, 30.0, true); +//! assert!((price - 20.0).abs() < 5e-16 * 20.0); +//! ``` mod erf_cody; mod normal_distribution; mod rational_cubic; -pub mod lets_be_rational; +pub(crate) mod lets_be_rational; +pub(crate) mod bachelier; +pub mod constants; + +/// Calculates the implied black volatility using a transformed rational guess with limited iterations. +/// +/// # Arguments +/// +/// * `option_price` - The current price of the option. +/// * `forward` - The current forward price of the underlying asset. +/// * `strike` - The strike price of the option. +/// * `expiry` - The time to expiration in years. +/// * `is_call` - A boolean flag indicating whether the option is a call (true) or put (false). +/// +/// # Returns +/// +/// The implied black volatility. +/// +/// # Examples +/// +/// ``` +/// use implied_vol::implied_black_volatility; +/// let black_vol = implied_black_volatility(20.0, 100.0, 90.0, 30.0, true); +/// 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 { + lets_be_rational::implied_black_volatility(option_price, forward, strike, expiry, is_call) +} + +/// Calculates the price of a European option using the Black-Scholes formula. +/// +/// # Arguments +/// +/// * `forward` - The current value of the underlying asset. +/// * `strike` - The strike price of the option. +/// * `volatility` - The volatility of the underlying asset. +/// * `expiry` - The time to expiration of the option. +/// * `is_call` - A boolean flag indicating whether the option is a call (true) or put (false). +/// +/// # Returns +/// +/// The price of the European option. +/// +/// # Examples +/// +/// ``` +/// use implied_vol::calculate_european_option_price_by_black_scholes; +/// let price = calculate_european_option_price_by_black_scholes(100.0, 90.0, 0.07011701801482094, 30.0, true); +/// assert!((price - 20.0).abs() < 5e-16 * 20.0); +/// ``` +#[inline] +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) +} + +/// Calculates the implied normal volatility. +/// +/// # Arguments +/// +/// * `price` - The market price of the option. +/// * `forward` - The forward price of the underlying asset. +/// * `strike` - The strike price of the option. +/// * `expiry` - The time to expiration in years. +/// * `is_call` - A boolean flag indicating whether the option is a call (true) or put (false). +/// +/// # Returns +/// +/// The implied normal volatility as a `f64` value. +/// +/// # Examples +/// +/// ``` +/// use implied_vol::implied_normal_volatility; +/// let normal_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 { + bachelier::implied_normal_volatility(option_price, forward, strike, expiry, is_call) +} + +/// Calculates the price of an option using Bachelier's model. +/// +/// # Arguments +/// +/// * `forward` - The forward price of the underlying asset. +/// * `strike` - The strike price of the option. +/// * `volatility` - The volatility of the underlying asset. +/// * `expiry` - The time to expiration in years. +/// * `is_call` - A boolean flag indicating whether the option is a call (true) or a put (false). +/// +/// # Returns +/// +/// The price of the European option. +/// +/// # Examples +/// +/// ``` +/// use implied_vol::calculate_european_option_price_by_bachelier; +/// let price = calculate_european_option_price_by_bachelier(100.0, 90.0, 6.614292466299764, 30.0, true); +/// assert!((price - 20.0).abs() < 5e-16 * 20.0); +/// ``` +#[inline] +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) +} \ No newline at end of file diff --git a/src/normal_distribution.rs b/src/normal_distribution.rs index 87d7581..6e3299d 100644 --- a/src/normal_distribution.rs +++ b/src/normal_distribution.rs @@ -5,6 +5,7 @@ const NORM_CDF_ASYMPTOTIC_EXPANSION_SECOND_THRESHOLD: f64 = -67108864.0; // 1.0 / f64::sqrt(f64::EPSILON); const FRAC_SQRT_2_PI: f64 = 0.398_942_280_401_432_7; +#[inline] pub(crate) fn norm_pdf(x: f64) -> f64 { FRAC_SQRT_2_PI * (-0.5 * x * x).exp() } @@ -36,7 +37,7 @@ pub(crate) fn norm_cdf(z: f64) -> f64 { } return -norm_pdf(z) * sum / z; } - 0.5 * erfc_cody(-z * (1.0 / std::f64::consts::SQRT_2)) + 0.5 * erfc_cody(-z * std::f64::consts::SQRT_2.recip()) } @@ -117,15 +118,15 @@ pub(crate) fn inverse_norm_cdf(u: f64) -> f64 { let mut r = if q < 0.0 { 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 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 < 0.0 { -ret } else { ret } } } \ No newline at end of file diff --git a/src/rational_cubic.rs b/src/rational_cubic.rs index 9e56dae..15e1f8a 100644 --- a/src/rational_cubic.rs +++ b/src/rational_cubic.rs @@ -1,6 +1,7 @@ const MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE: f64 = -(1f64 - 0.000000014901161193847656); const MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE: f64 = 2f64 / (f64::EPSILON * f64::EPSILON); +#[inline] fn is_zero(x: f64) -> bool { x.abs() < f64::MIN_POSITIVE } @@ -106,7 +107,7 @@ pub(crate) fn minimum_rational_cubic_control_parameter( } else if monotonic && s == 0.0 && prefer_shape_preservation_over_smoothness { r1 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE; } else { - r1 = -f64::MAX; + r1 = f64::MIN; } let r2; if convex || concave { @@ -115,12 +116,12 @@ pub(crate) fn minimum_rational_cubic_control_parameter( } else if prefer_shape_preservation_over_smoothness { r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE; } else { - r2 = -f64::MAX; + r2 = f64::MIN; } } else if monotonic && prefer_shape_preservation_over_smoothness { r2 = MAXIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE; } else { - r2 = -f64::MAX; + r2 = f64::MIN; } r1.max(r2) .max(MINIMUM_RATIONAL_CUBIC_CONTROL_PARAMETER_VALUE)