feat: Add implied normal volatility and Bachelier model

refactor: Change `q` f64 value to `is_call` flag refactor: Collect mathematical constants to constants.rs refactor: Change module visibility refactor: Add `inline` to thin-wrapper functions refactor: Remove `is_below_horizon` and use `is_subnormal` method instead refactor: Use f64::MIN instead of -f64::MAX refactor: Use `square` function instead of `powi(2)` method refactor: Remove comment outed code fix: Fix wrong constant value refactor: Remove mutable references from functions chore: Update LICENSE and README.md chore: Add doctest to CI
nakashima-hikaru · Jan 16, 2024 · ede561a · ede561a
1 parent f19d2c3
commit ede561a
Show file tree

Hide file tree

Showing 10 changed files with 529 additions and 335 deletions.
diff --git 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
diff --git 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.
+    ==========================================================================================
diff --git 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.
+## 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).
diff --git a/src/bachelier.rs 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);
+        }
+    }
+}
diff --git a/src/constants.rs 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;