jetblack-options 0.4.0

use crate::distributions::{cdf, pdf};
use crate::fdm::FdmWithoutCarry;
use crate::implied_volatility::solve_ivol;

/// # Black (1976)
///
/// Options on futures/forwards.
///
/// * The discounted futures price $ F $,
/// * Strike price $ K $,
/// * Risk-free rate $ r $,
/// * Annual dividend yield $ q $,
/// * Time to maturity $ \tau = t - t $
/// * Volatility $ \sigma $.
///
/// Most of the formula use one or both of the following terms.
///
/// $$
/// d_1 = \frac{\ln(F/K) + (\sigma^2/2)t}{\sigma\sqrt{t}}
/// $$
///
/// $$
/// d_2 = \frac{\ln(F/K) - (\sigma^2/2)t}{\sigma\sqrt{t}} = d_1 - \sigma\sqrt{t}
/// $$
///
/// $$
/// \varphi(x) &= \frac{1}{\sqrt{2\pi}} e^{-\frac{1}{2} x^2}
/// $$
///
/// $$
/// \Phi(x) &= \frac{1}{\sqrt{2\pi}} \int_{-\infty}^x e^{-\frac{1}{2} y^2} \,dy = 1 - \frac{1}{\sqrt{2\pi}} \int_x^\infty e^{-\frac{1}{2} y^2} \,dy
/// $$
pub struct Black76 {}

/// The following arguments are used in this module.
///
/// * is_call (bool): True for a call, false for a put.
/// * F (f64): The price of the future.
/// * K (f64): The strike price.
/// * t (f64): The time to expiry in years.
/// * r (f64): The risk free rate.
/// * v (f64): The asset volatility.
/// * max_iterations (usize): The maximum number of iterations before a price is returned.
/// * epsilon (f64): The largest acceptable error.
impl Black76 {
    /// The fair value of an option on a futures/forward using Black 76.
    ///
    /// For a call:
    ///
    /// $$
    /// C = e^{-r \tau}[F\Phi(d_1) - K\Phi(d_2)]
    /// $$
    ///
    /// For a put:
    ///
    /// $$
    /// P = e^{-r \tau} [K\Phi(-d_2) -  F\Phi(-d_1)]
    /// $$
    #[allow(non_snake_case)]
    pub fn price(is_call: bool, F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + ((v * v) / 2.0) * t) / (v * t.sqrt());
        let d2 = d1 - v * t.sqrt();

        if is_call {
            (-r * t).exp() * (F * cdf(d1) - K * cdf(d2))
        } else {
            (-r * t).exp() * (K * cdf(-d2) - F * cdf(-d1))
        }
    }

    /// Calculate the volatility of a Black 76 option that is implied by the price.
    #[allow(non_snake_case)]
    pub fn ivol(
        is_call: bool,
        F: f64,
        K: f64,
        t: f64,
        r: f64,
        p: f64,
        max_iterations: usize,
        epsilon: f64,
    ) -> f64 {
        solve_ivol(
            p,
            |v| Black76::price(is_call, F, K, t, r, v),
            max_iterations,
            epsilon,
        )
    }

    /// Return a struct to calculate greeks numerically using finite difference methods.
    pub fn fdm_greeks<'a>(is_call: bool) -> FdmWithoutCarry<'a> {
        #[allow(non_snake_case)]
        FdmWithoutCarry::new(move |S: f64, K: f64, t: f64, r: f64, b: f64| {
            Black76::price(is_call, S, K, t, r, b)
        })
    }

    /// The sensitivity of the option to a change in the asset price
    /// using Black 76.
    ///
    /// For the call.
    ///
    /// $$
    /// \frac{\partial C}{\partial S} = e^{-r \tau} \Phi(d_1)
    /// $$
    ///
    /// For the put.
    ///
    /// $$
    /// \frac{\partial P}{\partial S} = -e^{-r \tau} \Phi(-d_1)
    /// $$
    #[allow(non_snake_case)]
    pub fn delta(is_call: bool, F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + t * ((v * v) / 2.0)) / (v * t.sqrt());
        if is_call {
            return (-r * t).exp() * cdf(d1);
        } else {
            return -(-r * t).exp() * cdf(-d1);
        }
    }

    /// The second derivative to the change in asset price using Black 76.
    ///
    /// The gamma for both calls and puts.
    ///
    /// $$
    /// \frac{\partial^2 V}{\partial S^2} = e^{-r \tau} \frac{\varphi(d_1)}{F\sigma\sqrt{\tau}} = K e^{-r \tau} \frac{\varphi(d_2)}{F^2\sigma\sqrt{\tau}}
    /// $$
    #[allow(non_snake_case)]
    pub fn gamma(F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + t * ((v * v) / 2.0)) / (v * t.sqrt());

        (-r * t).exp() * pdf(d1) / (F * v * t.sqrt())
    }

    /// The change in the value of the option with respect to time to expiry
    /// using Black 76.
    ///
    /// For the call.
    ///
    /// $$
    /// \frac{\partial C}{\partial t} - \frac{F e^{-r \tau} \varphi(d_1) \sigma}{2 \sqrt{\tau}} - rKe^{-r \tau}\Phi(d_2) + rFe^{-r \tau}\Phi(d_1)
    /// $$
    ///
    /// For the put.
    ///
    /// $$
    /// \frac{\partial P}{\partial t} - \frac{F e^{-r \tau} \varphi(d_1) \sigma}{2 \sqrt{\tau}} + rKe^{-r \tau}\Phi(-d_2) - rFe^{-r \tau}\Phi(-d_1)
    /// $$
    #[allow(non_snake_case)]
    pub fn theta(is_call: bool, F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + ((v * v) / 2.0) * t) / (v * t.sqrt());
        let d2 = d1 - v * t.sqrt();

        if is_call {
            -F * (-r * t).exp() * pdf(d1) * v / (2.0 * t.sqrt()) + r * F * (-r * t).exp() * cdf(d1)
                - r * K * (-r * t).exp() * cdf(d2)
        } else {
            -F * (-r * t).exp() * pdf(d1) * v / (2.0 * t.sqrt()) - r * F * (-r * t).exp() * cdf(-d1)
                + r * K * (-r * t).exp() * cdf(-d2)
        }
    }

    /// The sensitivity of the options price or a change in the asset volatility
    /// using Black 76.
    ///
    /// For both calls and puts.
    ///
    /// $$
    /// \frac{\partial V}{\partial \sigma} = F e^{-r \tau} \varphi(d_1) \sqrt{\tau} = K e^{-r \tau} \varphi(d_2) \sqrt{\tau}
    /// $$
    #[allow(non_snake_case)]
    pub fn vega(F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + ((v * v) / 2.0) * t) / (v * t.sqrt());

        F * (-r * t).exp() * pdf(d1) * t.sqrt()
    }

    /// The sensitivity of the option price to a change in the risk free rate
    /// using Black 76.
    ///
    /// For a call:
    ///
    /// $$
    /// \frac{\partial C}{\partial r} = -\tau e^{-r \tau}[F\Phi(d_1) - K\Phi(d_2)]
    /// $$
    ///
    /// For a put:
    ///
    /// $$
    /// \frac{\partial P}{\partial r} = -\tau e^{-r \tau} [K\Phi(-d_2) -  F\Phi(-d_1)]
    /// $$
    #[allow(non_snake_case)]
    pub fn rho(is_call: bool, F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + ((v * v) / 2.0) * t) / (v * t.sqrt());
        let d2 = d1 - v * t.sqrt();
        if is_call {
            -t * (-r * t).exp() * (F * cdf(d1) - K * cdf(d2))
        } else {
            -t * (-r * t).exp() * (K * cdf(-d2) - F * cdf(-d1))
        }
    }

    /// The sensitivity of the option value to the underlying
    /// asset price and the volatility.
    ///
    /// For both calls and puts.
    ///
    /// $$
    /// \frac{\partial^2 V}{\partial F \partial \sigma} = -e^{-r \tau} \varphi(d_1) \frac{d_2}{\sigma} \, = \frac{\mathcal{V}}{F}\left[1 - \frac{d_1}{\sigma\sqrt{\tau}} \right]
    /// $$
    #[allow(non_snake_case)]
    pub fn vanna(F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + t * ((v * v) / 2.0)) / (v * t.sqrt());
        let d2 = d1 - v * t.sqrt();

        -(-r * t).exp() * pdf(d1) * d2 / v
    }

    /// The second order sensitivity to volatility.
    ///
    /// For both puts and calls.
    ///
    /// $$
    /// \frac{\partial^2 V}{\partial \sigma^2} = F e^{-r \tau} \varphi(d_1) \sqrt{\tau} \frac{d_1 d_2}{\sigma} = \mathcal{V}  \frac{d_1 d_2}{\sigma}
    /// $$
    #[allow(non_snake_case)]
    pub fn vomma(F: f64, K: f64, t: f64, r: f64, v: f64) -> f64 {
        let d1 = ((F / K).ln() + t * ((v * v) / 2.0)) / (v * t.sqrt());
        let d2 = d1 - v * t.sqrt();

        F * (-r * t).exp() * pdf(d1) * t.sqrt() * d1 * d2 / v
    }
}

#[cfg(test)]
mod tests {
    use std::{collections::HashMap, f64};

    use crate::fdm::DifferenceMethod;

    use super::*;

    fn is_close_to(actual: f64, expected: f64, threshold: f64) -> bool {
        (actual - expected).abs() < threshold
    }

    #[test]
    fn it_should_calc_price() {
        #[allow(non_snake_case)]
        for (is_call, F, K, t, r, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                6.0 / 12.0,
                0.1,
                0.125,
                10.143390791460092,
                1e-12,
            ),
            (
                false,
                110.0,
                100.0,
                6.0 / 12.0,
                0.1,
                0.125,
                0.63109654645295865,
                1e-12,
            ),
            (
                true,
                100.0,
                100.0,
                6.0 / 12.0,
                0.1,
                0.125,
                3.3531192847248605,
                1e-12,
            ),
            (
                false,
                100.0,
                100.0,
                6.0 / 12.0,
                0.1,
                0.125,
                3.3531192847248534,
                1e-12,
            ),
            (
                true,
                100.0,
                110.0,
                6.0 / 12.0,
                0.1,
                0.125,
                0.63109654645296376,
                1e-12,
            ),
            (
                false,
                100.0,
                110.0,
                6.0 / 12.0,
                0.1,
                0.125,
                10.143390791460105,
                1e-12,
            ),
        ] {
            let actual = Black76::price(is_call, F, K, t, r, v);
            assert!(
                is_close_to(actual, expected, f64::EPSILON),
                "price({}, {}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                actual,
                expected,
                (expected - actual).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_ivol() {
        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, p, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                10.143390791460092,
                0.125,
                1e-12,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.6310965464529535,
                0.125,
                1e-12,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                3.3531192847248605,
                0.125,
                1e-12,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                3.3531192847248534,
                0.125,
                1e-12,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.6310965464529654,
                0.125,
                1e-12,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                10.143390791460092,
                0.125,
                1e-12,
            ),
        ] {
            let actual = Black76::ivol(is_call, F, K, t, r, p, 100, 1e-12);
            assert!(
                is_close_to(actual, expected, threshold),
                "ivol({}, {}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                p,
                actual,
                expected,
                (expected - actual).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_delta() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.82678604419568191,
                1e-11,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.12444338030503208,
                1e-10,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.4923803086739813,
                1e-10,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.45884911582673277,
                1e-10,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.14319868380006498,
                1e-10,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.80803074070064895,
                1e-9,
            ),
        ] {
            let analytic = Black76::delta(is_call, F, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, f64::EPSILON),
                "delta({}, {}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                f64::EPSILON
            );

            let numeric = ng[&is_call].delta(F, K, t, r, v, 0.0001, DifferenceMethod::Central);
            assert!(
                is_close_to(numeric, analytic, threshold),
                "[{}].delta({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                numeric,
                analytic,
                (analytic - numeric).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_gamma() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.020787293603316447,
                1e-8,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.020787293603316447,
                1e-8,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.042891991459829734,
                1e-8,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.042891991459829734,
                1e-8,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.025152625260012922,
                1e-9,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                0.025152625260012922,
                1e-9,
            ),
        ] {
            let analytic = Black76::gamma(F, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, f64::EPSILON),
                "gamma({}, {}, {}, {}, {}) -> {} (expected: {}, diff {:e}, threshold: {:e})",
                F,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                f64::EPSILON
            );

            let numeric = ng[&is_call].gamma(F, K, t, r, v, 0.01, DifferenceMethod::Central);
            assert!(
                is_close_to(numeric, analytic, threshold),
                "[{}].gamma({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                numeric,
                analytic,
                (analytic - numeric).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_theta() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.95070976929249973,
                1e-8,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -1.901939193793212,
                1e-8,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -3.0156249043267125,
                1e-8,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -3.0156249043267138,
                1e-8,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -1.9019391937932129,
                1e-9,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.95070976929249973,
                1e-9,
            ),
        ] {
            let analytic = Black76::theta(is_call, F, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, f64::EPSILON),
                "theta({}, {}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                f64::EPSILON
            );

            let numerical = ng[&is_call].theta(
                F,
                K,
                t,
                r,
                v,
                1.0 / 365.0 / 24.0 / 60.0,
                DifferenceMethod::Central,
            );
            assert!(
                is_close_to(numerical, analytic, threshold),
                "[{}].theta({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                numerical,
                analytic,
                (analytic - numerical).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_vega() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                15.720390787508064,
                1e-7,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                15.720390787508064,
                1e-7,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                26.807494662393587,
                1e-9,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                26.807494662393587,
                1e-9,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                15.720390787508077,
                1e-7,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                15.720390787508077,
                1e-7,
            ),
        ] {
            let analytic = Black76::vega(F, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, f64::EPSILON),
                "vega({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e}",
                F,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                threshold
            );

            let numerical = ng[&is_call].vega(F, K, t, r, v, 0.00001, DifferenceMethod::Central);
            assert!(
                is_close_to(numerical, analytic, threshold),
                "[{}].vega({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                numerical,
                analytic,
                (analytic - numerical).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_rho() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, F, K, r, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -5.0716953957300461,
                1e-10,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.31554827322647933,
                1e-11,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -1.6765596423624303,
                1e-10,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -1.6765596423624267,
                1e-11,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -0.31554827322648188,
                1e-11,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                6.0 / 12.0,
                0.125,
                -5.0716953957300523,
                1e-10,
            ),
        ] {
            let analytic = Black76::rho(is_call, F, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, f64::EPSILON),
                "rho({}, {}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                threshold
            );

            let numeric = ng[&is_call].rho(F, K, t, r, v, 0.00001, DifferenceMethod::Central);
            assert!(
                is_close_to(numeric, analytic, threshold),
                "[{}].rho({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                F,
                K,
                t,
                r,
                v,
                numeric,
                analytic,
                (analytic - numeric).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_vanna() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, S, K, r, _q, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                -1.6720354862986977,
                1e-4,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                -1.6720354862986977,
                1e-4,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                0.13403747331196794,
                1e-4,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                0.13403747331196794,
                1e-4,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                1.996442942803649,
                1e-4,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                1.996442942803649,
                1e-4,
            ),
        ] {
            let analytic = Black76::vanna(S, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, 1e-12),
                "vanna({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                S,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                1e-12
            );

            let numeric = ng[&is_call].vanna(S, K, t, r, v, 0.01, 0.001, DifferenceMethod::Central);
            assert!(
                is_close_to(numeric, analytic, threshold),
                "[{}].vanna({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                S,
                K,
                t,
                r,
                v,
                numeric,
                analytic,
                (expected - analytic).abs(),
                threshold
            );
        }
    }

    #[test]
    fn it_should_calc_vomma() {
        let ng = HashMap::from([
            (true, Black76::fdm_greeks(true)),
            (false, Black76::fdm_greeks(false)),
        ]);

        #[allow(non_snake_case)]
        for (is_call, S, K, r, _q, t, v, expected, threshold) in [
            (
                true,
                110.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                145.98618448189166,
                1e-2,
            ),
            (
                false,
                110.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                145.98618448189166,
                1e-2,
            ),
            (
                true,
                100.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                -0.41886710409989986,
                1e-2,
            ),
            (
                false,
                100.0,
                100.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                -0.41886710409989986,
                1e-2,
            ),
            (
                true,
                100.0,
                110.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                145.98618448189166,
                1e-2,
            ),
            (
                false,
                100.0,
                110.0,
                0.1,
                0.08,
                6.0 / 12.0,
                0.125,
                145.98618448189166,
                1e-2,
            ),
        ] {
            let analytic = Black76::vomma(S, K, t, r, v);
            assert!(
                is_close_to(analytic, expected, 1e-12),
                "vomma({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                S,
                K,
                t,
                r,
                v,
                analytic,
                expected,
                (expected - analytic).abs(),
                1e-12
            );

            let numeric = ng[&is_call].vomma(S, K, t, r, v, 0.001, DifferenceMethod::Central);
            assert!(
                is_close_to(numeric, analytic, threshold),
                "[{}].vomma({}, {}, {}, {}, {}) -> {} (expected: {}, diff: {:e}, threshold: {:e})",
                is_call,
                S,
                K,
                t,
                r,
                v,
                numeric,
                analytic,
                (analytic - numeric).abs(),
                threshold
            );
        }
    }
}