libitofin 0.52.1

//! Black 1976 formula family.
//!
//! Port of the value and undiscounted-form subset of
//! `ql/pricingengines/blackformula.{hpp,cpp}`: [`black_formula`], its forward
//! derivative, the cash/asset in-the-money probabilities and the standard
//! deviation first and second derivatives. Every function takes the *standard
//! deviation* over the option life, `volatility * sqrt(time_to_maturity)`,
//! not the volatility itself, and an optional lognormal `displacement`
//! shifting both forward and strike.
//!
//! The Bachelier (normal-model) pricing pair, [`bachelier_black_formula`] and
//! its forward derivative [`bachelier_black_formula_forward_derivative`], sit
//! beside the Black family; the optionlet volatility surface selects between
//! the two through its volatility type. Unlike the lognormal family the normal
//! model admits negative forwards and strikes and applies no displacement, so
//! these two functions deliberately skip [`check_parameters`]: they validate
//! only the standard deviation and discount, exactly as the C++ reference does.
//!
//! Black approximation-seeded and normal implied-volatility inversions support
//! caplet stripping; other specialized approximation families remain separate.
//!
//! At `std_dev == 0` the Black cash and asset and the Bachelier asset
//! in-the-money probabilities all return QuantLib v1.44's indicator
//! (`blackformula.cpp:580,611,970`): one when `sign * (forward - strike) > 0`,
//! zero otherwise. That is the `std_dev -> 0` limit off the money; exactly at
//! the money it is zero, where the limit would be one half.

use crate::errors::QlResult;
use crate::fail;
use crate::math::comparison::close_enough;
use crate::math::distributions::normal::{CumulativeNormalDistribution, NormalDistribution};
use crate::math::solver1d::{DerivativeSolver, func1d};
use crate::math::solvers1d::newtonsafe::NewtonSafe;
use crate::option::OptionType;
use crate::types::{Natural, Real};

/// QuantLib's `checkParameters` (`blackformula.cpp:44,47,50`): the
/// `displacement >= 0`, `strike + displacement >= 0` and
/// `forward + displacement > 0` requirements, kept intact.
///
/// Divergence: the standalone finiteness checks on `strike` and `forward`, and
/// the `!is_finite()` clauses replacing C++'s implicit NaN handling. In C++ a
/// NaN argument fails every comparison, so `QL_REQUIRE(x >= 0.0)` already
/// throws; an infinite one does not, and `+inf - inf` in the shifted sums then
/// yields NaN downstream. Rejecting both here keeps the failure at the boundary.
fn check_parameters(strike: Real, forward: Real, displacement: Real) -> QlResult<()> {
    if !displacement.is_finite() || displacement < 0.0 {
        fail!("displacement ({displacement}) must be non-negative");
    }
    if !strike.is_finite() {
        fail!("strike ({strike}) must be finite");
    }
    if !forward.is_finite() {
        fail!("forward ({forward}) must be finite");
    }
    let shifted_strike = strike + displacement;
    if !shifted_strike.is_finite() || shifted_strike < 0.0 {
        fail!("strike + displacement ({strike} + {displacement}) must be non-negative");
    }
    let shifted_forward = forward + displacement;
    if !shifted_forward.is_finite() || shifted_forward <= 0.0 {
        fail!("forward + displacement ({forward} + {displacement}) must be positive");
    }
    Ok(())
}

fn check_std_dev_and_discount(std_dev: Real, discount: Real) -> QlResult<()> {
    check_std_dev(std_dev)?;
    if !discount.is_finite() || discount <= 0.0 {
        fail!("discount ({discount}) must be positive");
    }
    Ok(())
}

/// QuantLib's `QL_REQUIRE(stdDev >= 0.0)` (`blackformula.cpp:67`), extended to
/// reject `+inf`.
///
/// Divergence: `blackFormulaCashItmProbability` and
/// `blackFormulaAssetItmProbability` call only `checkParameters` and never
/// validate `stdDev`, so a negative one silently flips the sign of `d2`. This
/// port applies the same check there as in `black_formula`.
fn check_std_dev(std_dev: Real) -> QlResult<()> {
    if !std_dev.is_finite() || std_dev < 0.0 {
        fail!("stdDev ({std_dev}) must be non-negative");
    }
    Ok(())
}

fn sign_of(option_type: OptionType) -> Real {
    Real::from(option_type as i32)
}

/// Black 1976 value of a European option on the given forward.
pub fn black_formula(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev_and_discount(std_dev, discount)?;

    let sign = sign_of(option_type);

    if std_dev == 0.0 {
        let intrinsic = (forward - strike) * sign;
        let intrinsic = if intrinsic < 0.0 { 0.0 } else { intrinsic };
        return Ok(intrinsic * discount);
    }

    let forward = forward + displacement;
    let strike = strike + displacement;

    if strike == 0.0 {
        return Ok(match option_type {
            OptionType::Call => forward * discount,
            OptionType::Put => 0.0,
        });
    }

    let d1 = (forward / strike).ln() / std_dev + 0.5 * std_dev;
    let d2 = d1 - std_dev;
    let phi = CumulativeNormalDistribution::standard();
    let nd1 = phi.value(sign * d1);
    let nd2 = phi.value(sign * d2);
    let result = discount * sign * (forward * nd1 - strike * nd2);
    if result.is_nan() || result < 0.0 {
        fail!(
            "negative value ({result}) for {std_dev} stdDev, {option_type} option, \
             {strike} strike, {forward} forward"
        );
    }
    Ok(result)
}

/// Derivative of [`black_formula`] with respect to the forward.
pub fn black_formula_forward_derivative(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev_and_discount(std_dev, discount)?;

    let sign = sign_of(option_type);

    if std_dev == 0.0 {
        let moneyness = (forward - strike) * sign;
        return Ok(if moneyness > 0.0 {
            sign * discount
        } else {
            0.0
        });
    }

    let forward = forward + displacement;
    let strike = strike + displacement;

    if strike == 0.0 {
        return Ok(match option_type {
            OptionType::Call => discount,
            OptionType::Put => 0.0,
        });
    }

    let d1 = (forward / strike).ln() / std_dev + 0.5 * std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(sign * phi.value(sign * d1) * discount)
}

/// Risk-neutral probability of exercise in the bond martingale measure, `N(d2)`.
pub fn black_formula_cash_itm_probability(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev(std_dev)?;

    let sign = sign_of(option_type);

    if std_dev == 0.0 {
        return Ok(if forward * sign > strike * sign {
            1.0
        } else {
            0.0
        });
    }

    let forward = forward + displacement;
    let strike = strike + displacement;
    if strike == 0.0 {
        return Ok(match option_type {
            OptionType::Call => 1.0,
            OptionType::Put => 0.0,
        });
    }
    let d2 = (forward / strike).ln() / std_dev - 0.5 * std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(phi.value(sign * d2))
}

/// Risk-neutral probability of exercise in the asset martingale measure, `N(d1)`.
pub fn black_formula_asset_itm_probability(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev(std_dev)?;

    let sign = sign_of(option_type);

    if std_dev == 0.0 {
        return Ok(if forward * sign > strike * sign {
            1.0
        } else {
            0.0
        });
    }

    let forward = forward + displacement;
    let strike = strike + displacement;
    if strike == 0.0 {
        return Ok(match option_type {
            OptionType::Call => 1.0,
            OptionType::Put => 0.0,
        });
    }
    let d1 = (forward / strike).ln() / std_dev + 0.5 * std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(phi.value(sign * d1))
}

/// Derivative of [`black_formula`] with respect to the standard deviation.
///
/// Multiplying by `sqrt(time_to_maturity)` turns this into the Black vega.
pub fn black_formula_std_dev_derivative(
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev_and_discount(std_dev, discount)?;

    let forward = forward + displacement;
    let strike = strike + displacement;

    if std_dev == 0.0 || strike == 0.0 {
        return Ok(0.0);
    }

    let d1 = (forward / strike).ln() / std_dev + 0.5 * std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(discount * forward * phi.derivative(d1))
}

/// Derivative of [`black_formula`] with respect to the implied volatility.
///
/// This is the Black vega: [`black_formula_std_dev_derivative`] times
/// `sqrt(expiry)`.
pub fn black_formula_vol_derivative(
    strike: Real,
    forward: Real,
    std_dev: Real,
    expiry: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    let derivative =
        black_formula_std_dev_derivative(strike, forward, std_dev, discount, displacement)?;
    Ok(derivative * expiry.sqrt())
}

/// Second derivative of [`black_formula`] with respect to the standard deviation.
pub fn black_formula_std_dev_second_derivative(
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    check_std_dev_and_discount(std_dev, discount)?;

    let forward = forward + displacement;
    let strike = strike + displacement;

    if std_dev == 0.0 || strike == 0.0 {
        return Ok(0.0);
    }

    let d1 = (forward / strike).ln() / std_dev + 0.5 * std_dev;
    let d1_prime = -(forward / strike).ln() / (std_dev * std_dev) + 0.5;
    let density = NormalDistribution::standard();
    Ok(discount * forward * density.derivative(d1) * d1_prime)
}

/// Implied Black standard deviation that reprices `black_price` under
/// [`black_formula`].
///
/// Port of `blackFormulaImpliedStdDev` (`blackformula.cpp:382-440`): checks the
/// price against put-call parity, converts an in-the-money option to its
/// out-of-the-money complement (greater vega/price ratio, hence numerically more
/// robust), then inverts [`black_formula`] for the standard deviation with a
/// [`NewtonSafe`] solve bracketed on `[0, 24]` (`24 = 300% * sqrt(60)`).
///
/// The undiscounted objective and its derivative mirror C++'s
/// `BlackImpliedStdDevHelper` (`blackformula.cpp:330-379`): the value is
/// `black_formula(type, strike, forward, std_dev, 1, 0) - black_price/discount`
/// on the displacement-shifted strike and forward, and the derivative is the
/// helper's signed `signedForward * phi(signed d1)`, reproduced here as
/// `sign * black_formula_std_dev_derivative(...)`. The sign is negative for the
/// out-of-the-money put branch, exactly as the helper computes it; the solve
/// still converges because [`NewtonSafe`] bisects whenever a Newton step leaves
/// the bracket.
///
/// Pass `None` for QuantLib's approximation-seeded null-guess branch; a finite
/// scalar guess remains accepted for backwards compatibility.
///
/// # Errors
///
/// Rejects a non-positive `discount`, a negative `black_price`, a price implying
/// a negative complementary option under put-call parity, and a negative or
/// non-finite `guess`; propagates a solver that fails to converge.
#[allow(clippy::too_many_arguments)]
pub fn black_formula_implied_std_dev(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    black_price: Real,
    discount: Real,
    displacement: Real,
    guess: impl Into<Option<Real>>,
    accuracy: Real,
    max_iterations: Natural,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    if !discount.is_finite() || discount <= 0.0 {
        fail!("discount ({discount}) must be positive");
    }
    if !black_price.is_finite() || black_price < 0.0 {
        fail!("option price ({black_price}) must be non-negative");
    }

    let mut option_type = option_type;
    let mut black_price = black_price;
    let other_option_price = black_price - sign_of(option_type) * (forward - strike) * discount;
    if other_option_price < 0.0 {
        fail!(
            "negative complementary-option price ({other_option_price}) implied by put-call \
             parity: no solution exists for {option_type} strike {strike}, forward {forward}, \
             price {black_price}, deflator {discount}"
        );
    }

    if matches!(option_type, OptionType::Put) && strike > forward {
        option_type = OptionType::Call;
        black_price = other_option_price;
    }
    if matches!(option_type, OptionType::Call) && strike < forward {
        option_type = OptionType::Put;
        black_price = other_option_price;
    }

    let guess = match guess.into() {
        Some(guess) => guess,
        None => black_formula_implied_std_dev_approximation(
            option_type,
            strike,
            forward,
            black_price,
            discount,
            displacement,
        )?,
    };
    if !guess.is_finite() || guess < 0.0 {
        fail!("stdDev guess ({guess}) must be non-negative");
    }

    let shifted_strike = strike + displacement;
    let shifted_forward = forward + displacement;
    let target = black_price / discount;
    let otm_sign = sign_of(option_type);

    let value = |std_dev: Real| {
        black_formula(
            option_type,
            shifted_strike,
            shifted_forward,
            std_dev,
            1.0,
            0.0,
        )
        .map(|price| price - target)
        .unwrap_or(Real::NAN)
    };
    let derivative = |std_dev: Real| {
        black_formula_std_dev_derivative(shifted_strike, shifted_forward, std_dev, 1.0, 0.0)
            .map(|d| otm_sign * d)
            .unwrap_or(Real::NAN)
    };

    let std_dev = NewtonSafe::new()
        .with_max_evaluations(max_iterations as usize)
        .solve_bracketed(func1d(value, derivative), accuracy, guess, 0.0, 24.0)?;
    if std_dev < 0.0 {
        fail!("stdDev ({std_dev}) must be non-negative");
    }
    Ok(std_dev)
}

/// Corrado-Miller approximation, with the Brenner-Subrahmanyan ATM limit.
///
/// This is QuantLib's seed for an unspecified Black implied-standard-deviation guess.
pub fn black_formula_implied_std_dev_approximation(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    price: Real,
    discount: Real,
    displacement: Real,
) -> QlResult<Real> {
    check_parameters(strike, forward, displacement)?;
    if !price.is_finite() || price < 0.0 || !discount.is_finite() || discount <= 0.0 {
        fail!("invalid option price or discount for implied deviation approximation");
    }
    let forward = forward + displacement;
    let strike = strike + displacement;
    let scale = (2.0 * std::f64::consts::PI).sqrt();
    if strike == forward {
        return Ok(price / discount * scale / forward);
    }
    let delta = sign_of(option_type) * (forward - strike);
    let value = price / discount - delta / 2.0;
    Ok((value
        + (value * value - delta * delta / std::f64::consts::PI)
            .max(0.0)
            .sqrt())
        * scale
        / (forward + strike))
}

/// Normal implied volatility, solved on out-of-the-money time value.
///
/// A bracketed inversion avoids subtracting intrinsic value inside the solver.
/// As in QuantLib v1.44 (`blackformula.cpp:908-926`), a strike `close_enough` to
/// the forward takes the at-the-money closed form, and a premium `close_enough`
/// to the intrinsic value returns zero before the negative-time-value check, so
/// a price that rounds a few ulps below intrinsic still inverts.
/// Invalid prices and nonpositive expiry/discount return an error.
pub fn bachelier_black_formula_implied_vol(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    expiry: Real,
    price: Real,
    discount: Real,
) -> QlResult<Real> {
    use crate::math::solver1d::Solver1D;
    use crate::math::solvers1d::brent::Brent;
    if ![strike, forward, expiry, price, discount]
        .iter()
        .all(|v| v.is_finite())
        || expiry <= 0.0
        || price < 0.0
        || discount <= 0.0
    {
        fail!("invalid normal implied-volatility inputs");
    }
    let premium = price / discount;
    let intrinsic = (sign_of(option_type) * (forward - strike)).max(0.0);
    let distance = (forward - strike).abs();
    let scale = (2.0 * std::f64::consts::PI).sqrt();
    if close_enough(strike, forward) {
        return Ok(premium * scale / expiry.sqrt());
    }
    if close_enough(premium, intrinsic) {
        return Ok(0.0);
    }
    let time_value = premium - intrinsic;
    if time_value < 0.0 {
        fail!("option price implies negative time value ({time_value})");
    }
    let upper = (time_value * scale + distance) * 2.0;
    let objective = |std_dev| {
        bachelier_black_formula(OptionType::Call, distance, 0.0, std_dev, 1.0).unwrap_or(Real::NAN)
            - time_value
    };
    let std_dev = Brent::new().with_max_evaluations(200).solve_bracketed(
        objective,
        1e-14,
        upper / 2.0,
        0.0,
        upper,
    )?;
    Ok(std_dev / expiry.sqrt())
}

/// Bachelier (normal-model) value of a European option on the given forward.
///
/// Port of `bachelierBlackFormula` (`blackformula.cpp:705`). With
/// `d = (forward - strike) * sign` and `h = d / std_dev`, the premium is
/// `discount * (std_dev * phi(h) + d * Phi(h))`, where `phi`/`Phi` are the
/// standard normal density and cumulative distribution. The normal model
/// prices negative forwards and strikes, so only `std_dev` and `discount` are
/// validated; the terminal non-negativity check mirrors C++'s `QL_ENSURE`.
pub fn bachelier_black_formula(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
) -> QlResult<Real> {
    check_std_dev_and_discount(std_dev, discount)?;

    let sign = sign_of(option_type);
    let d = (forward - strike) * sign;

    if std_dev == 0.0 {
        return Ok(discount * d.max(0.0));
    }

    let h = d / std_dev;
    let phi = CumulativeNormalDistribution::standard();
    let result = discount * (std_dev * phi.derivative(h) + d * phi.value(h));
    if result.is_nan() || result < 0.0 {
        fail!(
            "negative value ({result}) for {std_dev} stdDev, {option_type} option, \
             {strike} strike, {forward} forward"
        );
    }
    Ok(result)
}

/// Derivative of [`bachelier_black_formula`] with respect to the forward.
///
/// Port of `bachelierBlackFormulaForwardDerivative` (`blackformula.cpp:738`),
/// which equals `sign * Phi(h) * discount`. At `std_dev == 0` the reference
/// collapses this to `sign * max(boost_sign((forward - strike) * sign), 0) *
/// discount`; `boost::math::sign` yields `0` at the money, so the derivative
/// is `0` exactly when `forward == strike` (distinct from `f64::signum`, which
/// never returns zero).
pub fn bachelier_black_formula_forward_derivative(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
) -> QlResult<Real> {
    check_std_dev_and_discount(std_dev, discount)?;

    let sign = sign_of(option_type);
    let moneyness = (forward - strike) * sign;

    if std_dev == 0.0 {
        let boost_sign: Real = if moneyness > 0.0 {
            1.0
        } else if moneyness < 0.0 {
            -1.0
        } else {
            0.0
        };
        return Ok(sign * boost_sign.max(0.0) * discount);
    }

    let h = moneyness / std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(sign * phi.value(h) * discount)
}

/// Derivative of [`bachelier_black_formula`] with respect to the standard
/// deviation.
///
/// Port of `bachelierBlackFormulaStdDevDerivative` (`blackformula.cpp:923`),
/// which equals `discount * phi((forward - strike) / std_dev)` where `phi` is
/// the standard normal density. The normal model prices negative forwards and
/// strikes, so only `std_dev` and `discount` are validated. Multiplying by
/// `sqrt(exercise_time)` turns this into the Bachelier vega.
pub fn bachelier_black_formula_std_dev_derivative(
    strike: Real,
    forward: Real,
    std_dev: Real,
    discount: Real,
) -> QlResult<Real> {
    check_std_dev_and_discount(std_dev, discount)?;

    if std_dev == 0.0 {
        return Ok(0.0);
    }

    let d1 = (forward - strike) / std_dev;
    let phi = CumulativeNormalDistribution::standard();
    Ok(discount * phi.derivative(d1))
}

/// QuantLib's Bachelier asset ITM helper (`blackformula.cpp:962`).
///
/// For positive `std_dev`, returns `Phi(sign * (forward - strike) / std_dev)`.
/// At exactly zero, returns QuantLib v1.44's indicator (`blackformula.cpp:970`):
/// one when `sign * (forward - strike) > 0`, zero otherwise. Its at-the-money
/// value is zero, not the positive-volatility limit of one half. Negative
/// forwards and strikes are supported.
///
/// Divergence: rejects non-finite inputs. For positive volatility, an
/// overflowing difference is divided before subtraction, retaining a
/// representable standardized distance. An infinite standardized distance
/// legitimately yields a normal-CDF endpoint.
///
/// # Errors
///
/// Returns an error for non-finite inputs or negative `std_dev`.
///
/// # Examples
///
/// ```
/// use libitofin::option::OptionType;
/// use libitofin::pricingengines::blackformula::bachelier_black_formula_asset_itm_probability;
/// let value = bachelier_black_formula_asset_itm_probability(
///     OptionType::Call, -0.01, -0.01, 0.005,
/// )?;
/// assert_eq!(value, 0.5);
/// # Ok::<(), libitofin::errors::QlError>(())
/// ```
pub fn bachelier_black_formula_asset_itm_probability(
    option_type: OptionType,
    strike: Real,
    forward: Real,
    std_dev: Real,
) -> QlResult<Real> {
    check_std_dev(std_dev)?;
    if !strike.is_finite() || !forward.is_finite() {
        fail!("strike ({strike}) and forward ({forward}) must be finite");
    }

    let sign = sign_of(option_type);
    let difference = forward - strike;
    if std_dev == 0.0 {
        return Ok(if difference * sign > 0.0 { 1.0 } else { 0.0 });
    }

    let h = if difference.is_finite() {
        difference / std_dev
    } else {
        forward / std_dev - strike / std_dev
    };
    Ok(CumulativeNormalDistribution::standard().value(sign * h))
}

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

    use crate::pricingengines::hull_fixture::{
        DISCOUNT as HULL_DISCOUNT, FORWARD as HULL_FORWARD, STD_DEV as HULL_STD_DEV,
    };

    fn assert_close(actual: Real, expected: Real, tolerance: Real) {
        assert!(
            (actual - expected).abs() <= tolerance,
            "actual {actual} vs expected {expected} (tolerance {tolerance})"
        );
    }

    #[test]
    fn known_values_match_black_scholes() {
        let call = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let put = black_formula(
            OptionType::Put,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        assert_close(call, 4.759422392871536, 1e-10);
        assert_close(put, 0.8085993729000926, 1e-10);
    }

    #[test]
    fn put_call_parity_holds() {
        for strike in [20.0, 40.0, 44.15338604779301, 60.0] {
            let call = black_formula(
                OptionType::Call,
                strike,
                HULL_FORWARD,
                HULL_STD_DEV,
                HULL_DISCOUNT,
                0.0,
            )
            .expect("valid inputs");
            let put = black_formula(
                OptionType::Put,
                strike,
                HULL_FORWARD,
                HULL_STD_DEV,
                HULL_DISCOUNT,
                0.0,
            )
            .expect("valid inputs");
            assert_close(call - put, HULL_DISCOUNT * (HULL_FORWARD - strike), 1e-12);
        }
    }

    #[test]
    fn zero_std_dev_returns_discounted_intrinsic() {
        let call = black_formula(OptionType::Call, 40.0, 44.0, 0.0, 0.95, 0.0).expect("valid");
        assert_close(call, 0.95 * 4.0, 1e-15);
        let put = black_formula(OptionType::Put, 40.0, 44.0, 0.0, 0.95, 0.0).expect("valid");
        assert_close(put, 0.0, 0.0);
    }

    #[test]
    fn zero_strike_prices_the_forward() {
        let call = black_formula(OptionType::Call, 0.0, 44.0, 0.2, 0.95, 0.0).expect("valid");
        assert_close(call, 44.0 * 0.95, 1e-15);
        let put = black_formula(OptionType::Put, 0.0, 44.0, 0.2, 0.95, 0.0).expect("valid");
        assert_close(put, 0.0, 0.0);
    }

    #[test]
    fn invalid_inputs_are_rejected() {
        assert!(black_formula(OptionType::Call, 40.0, 44.0, -0.1, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, 44.0, 0.1, 0.0, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, 44.0, Real::NAN, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, Real::NAN, 0.1, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, -1.0, 44.0, 0.1, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, -44.0, 0.1, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, 44.0, 0.1, 0.95, -0.01).is_err());
    }

    fn assert_forward_derivative_consistency(option_type: OptionType, strikes: &[Real], vol: Real) {
        let forward = 1.0;
        let tte: Real = 10.0;
        let std_dev = vol * tte.sqrt();
        let discount = 0.95;
        let displacement = 0.01;
        let bump = 0.0001;
        let epsilon = 1.0e-10;

        for &strike in strikes {
            let delta = black_formula_forward_derivative(
                option_type,
                strike,
                forward,
                std_dev,
                discount,
                displacement,
            )
            .expect("valid inputs");
            let bumped_delta = black_formula_forward_derivative(
                option_type,
                strike,
                forward + bump,
                std_dev,
                discount,
                displacement,
            )
            .expect("valid inputs");

            let base_premium = black_formula(
                option_type,
                strike,
                forward,
                std_dev,
                discount,
                displacement,
            )
            .expect("valid inputs");
            let bumped_premium = black_formula(
                option_type,
                strike,
                forward + bump,
                std_dev,
                discount,
                displacement,
            )
            .expect("valid inputs");
            let delta_approx = (bumped_premium - base_premium) / bump;

            assert!(
                delta.max(bumped_delta) + epsilon > delta_approx
                    && delta_approx > delta.min(bumped_delta) - epsilon,
                "forward derivative inconsistent with bump for {option_type} at strike {strike}: \
                 analytical {delta}, approximated {delta_approx}"
            );
        }
    }

    #[test]
    fn forward_derivative_is_consistent_with_bumping() {
        let strikes = [0.1, 0.5, 1.0, 2.0, 3.0];
        assert_forward_derivative_consistency(OptionType::Call, &strikes, 0.1);
        assert_forward_derivative_consistency(OptionType::Put, &strikes, 0.1);
    }

    #[test]
    fn forward_derivative_is_consistent_with_bumping_at_zero_strike() {
        assert_forward_derivative_consistency(OptionType::Call, &[0.0], 0.1);
        assert_forward_derivative_consistency(OptionType::Put, &[0.0], 0.1);
    }

    #[test]
    fn forward_derivative_is_consistent_with_bumping_at_zero_volatility() {
        let strikes = [0.1, 0.5, 1.0, 2.0, 3.0];
        assert_forward_derivative_consistency(OptionType::Call, &strikes, 0.0);
        assert_forward_derivative_consistency(OptionType::Put, &strikes, 0.0);
    }

    #[test]
    fn std_dev_derivative_matches_bumped_value() {
        let bump = 1.0e-6;
        let up = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV + bump,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let down = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV - bump,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let analytical =
            black_formula_std_dev_derivative(40.0, HULL_FORWARD, HULL_STD_DEV, HULL_DISCOUNT, 0.0)
                .expect("valid inputs");
        assert_close((up - down) / (2.0 * bump), analytical, 1e-6);

        let wide = 1.0e-4;
        let up_wide = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV + wide,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let down_wide = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV - wide,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let second = black_formula_std_dev_second_derivative(
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        let base = black_formula(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV,
            HULL_DISCOUNT,
            0.0,
        )
        .expect("valid inputs");
        assert_close(
            (up_wide - 2.0 * base + down_wide) / (wide * wide),
            second,
            1e-4,
        );
    }

    #[test]
    fn vol_derivative_scales_by_sqrt_expiry() {
        let std_dev_derivative =
            black_formula_std_dev_derivative(40.0, HULL_FORWARD, HULL_STD_DEV, HULL_DISCOUNT, 0.0)
                .expect("valid inputs");
        let vol_derivative =
            black_formula_vol_derivative(40.0, HULL_FORWARD, HULL_STD_DEV, 0.5, HULL_DISCOUNT, 0.0)
                .expect("valid inputs");
        assert_close(vol_derivative, std_dev_derivative * 0.5_f64.sqrt(), 1e-12);
        assert_close(vol_derivative, 8.813415059602862, 1e-10);
    }

    #[test]
    fn itm_probabilities_match_normal_quantiles() {
        let cash = black_formula_cash_itm_probability(
            OptionType::Call,
            40.0,
            HULL_FORWARD,
            HULL_STD_DEV,
            0.0,
        )
        .expect("valid inputs");
        assert_close(cash, 0.7349460368459086, 1e-10);
    }

    #[test]
    fn asset_itm_probability_is_continuous_at_zero_std_dev() {
        for (option_type, forward) in [
            (OptionType::Call, 44.0),
            (OptionType::Call, 36.0),
            (OptionType::Put, 44.0),
            (OptionType::Put, 36.0),
        ] {
            let limit = black_formula_asset_itm_probability(option_type, 40.0, forward, 1e-12, 0.0)
                .expect("valid inputs");
            let at_zero = black_formula_asset_itm_probability(option_type, 40.0, forward, 0.0, 0.0)
                .expect("valid inputs");
            assert_close(at_zero, limit, 1e-9);
        }
    }

    #[test]
    fn itm_probabilities_at_the_money_keep_the_zero_convention() {
        for option_type in [OptionType::Call, OptionType::Put] {
            let asset = black_formula_asset_itm_probability(option_type, 40.0, 40.0, 0.0, 0.0)
                .expect("valid inputs");
            assert_close(asset, 0.0, 0.0);
            let cash = black_formula_cash_itm_probability(option_type, 40.0, 40.0, 0.0, 0.0)
                .expect("valid inputs");
            assert_close(cash, 0.0, 0.0);
        }
    }

    #[test]
    fn itm_probabilities_reject_invalid_standard_deviation() {
        assert!(
            black_formula_cash_itm_probability(OptionType::Call, 40.0, 44.0, -0.1, 0.0).is_err()
        );
        assert!(
            black_formula_asset_itm_probability(OptionType::Call, 40.0, 44.0, -0.1, 0.0).is_err()
        );
        assert!(
            black_formula_cash_itm_probability(OptionType::Call, 40.0, 44.0, Real::NAN, 0.0)
                .is_err()
        );
        assert!(
            black_formula_asset_itm_probability(OptionType::Call, 40.0, 44.0, Real::NAN, 0.0)
                .is_err()
        );
        assert!(
            black_formula_cash_itm_probability(OptionType::Call, 40.0, 44.0, Real::INFINITY, 0.0)
                .is_err()
        );
    }

    #[test]
    fn black_formula_rejects_non_finite_inputs() {
        assert!(black_formula(OptionType::Call, Real::INFINITY, 44.0, 0.2, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, Real::INFINITY, 0.2, 0.95, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, 44.0, 0.2, Real::INFINITY, 0.0).is_err());
        assert!(black_formula(OptionType::Call, 40.0, 44.0, 0.2, 0.95, Real::INFINITY).is_err());
    }

    #[test]
    fn bachelier_put_call_parity_is_exact() {
        // C - P == discount * (forward - strike) for the normal model, derived from
        // C - P = discount * [d*Phi(h) + d*(1 - Phi(h))] with d = forward - strike.
        for (forward, strike) in [(0.03, 0.02), (0.02, 0.02), (-0.01, 0.01), (0.05, -0.02)] {
            for std_dev in [0.001, 0.02, 0.1] {
                let call =
                    bachelier_black_formula(OptionType::Call, strike, forward, std_dev, 0.95)
                        .expect("valid inputs");
                let put = bachelier_black_formula(OptionType::Put, strike, forward, std_dev, 0.95)
                    .expect("valid inputs");
                assert_close(call - put, 0.95 * (forward - strike), 1e-15);
            }
        }
    }

    #[test]
    fn bachelier_zero_volatility_is_discounted_intrinsic() {
        for (forward, strike) in [(0.03, 0.02), (0.01, 0.02), (-0.02, -0.01)] {
            let call = bachelier_black_formula(OptionType::Call, strike, forward, 0.0, 0.95)
                .expect("valid inputs");
            let put = bachelier_black_formula(OptionType::Put, strike, forward, 0.0, 0.95)
                .expect("valid inputs");
            assert_close(call, 0.95 * (forward - strike).max(0.0), 0.0);
            assert_close(put, 0.95 * (strike - forward).max(0.0), 0.0);
        }
    }

    #[test]
    fn bachelier_at_the_money_matches_the_closed_form() {
        // At forward == strike, d = 0 and h = 0, so the premium collapses to
        // discount * std_dev * phi(0) = discount * std_dev / sqrt(2*pi).
        for std_dev in [0.001, 0.05, 0.2] {
            let expected = 0.95 * std_dev / (2.0 * std::f64::consts::PI).sqrt();
            for option_type in [OptionType::Call, OptionType::Put] {
                let premium = bachelier_black_formula(option_type, 1.0, 1.0, std_dev, 0.95)
                    .expect("valid inputs");
                assert_close(premium, expected, 1e-16);
            }
        }
    }

    #[test]
    fn bachelier_matches_the_normal_pricing_definition() {
        // premium = discount * (std_dev * phi(h) + d * Phi(h)), d = (forward-strike)*sign,
        // h = d / std_dev, reconstructed from the normal primitives directly.
        let phi = CumulativeNormalDistribution::standard();
        for (forward, strike) in [(0.03, 0.025), (-0.01, 0.005), (0.02, 0.02)] {
            for std_dev in [0.001, 0.01, 0.05] {
                for option_type in [OptionType::Call, OptionType::Put] {
                    let sign = Real::from(option_type as i32);
                    let d = (forward - strike) * sign;
                    let h = d / std_dev;
                    let expected = 0.95 * (std_dev * phi.derivative(h) + d * phi.value(h));
                    let actual =
                        bachelier_black_formula(option_type, strike, forward, std_dev, 0.95)
                            .expect("valid inputs");
                    assert_close(actual, expected, 1e-16);
                }
            }
        }
    }

    #[test]
    fn bachelier_prices_negative_forward_and_strike() {
        let premium = bachelier_black_formula(OptionType::Call, -0.005, -0.01, 0.01, 1.0)
            .expect("normal model admits negative rates");
        assert!(premium > 0.0);
    }

    #[test]
    fn bachelier_rejects_invalid_std_dev_and_discount() {
        assert!(bachelier_black_formula(OptionType::Call, 1.0, 1.0, -0.1, 0.95).is_err());
        assert!(bachelier_black_formula(OptionType::Call, 1.0, 1.0, Real::INFINITY, 0.95).is_err());
        assert!(bachelier_black_formula(OptionType::Call, 1.0, 1.0, 0.1, 0.0).is_err());
        assert!(bachelier_black_formula(OptionType::Call, 1.0, 1.0, 0.1, -1.0).is_err());
        assert!(
            bachelier_black_formula_forward_derivative(OptionType::Call, 1.0, 1.0, -0.1, 0.95)
                .is_err()
        );
        assert!(
            bachelier_black_formula_forward_derivative(OptionType::Call, 1.0, 1.0, 0.1, 0.0)
                .is_err()
        );
    }

    fn assert_bachelier_forward_derivative(option_type: OptionType, strikes: &[Real], bpvol: Real) {
        // Mean-value-theorem check ported from blackformula.cpp:357: the analytical
        // forward derivative must bracket the bumped finite difference of the premium.
        let forward = 1.0;
        let tte: Real = 10.0;
        let std_dev = bpvol * tte.sqrt();
        let discount = 0.95;
        let bump = 0.0001;
        let epsilon = 1.0e-10;
        for &strike in strikes {
            let delta = bachelier_black_formula_forward_derivative(
                option_type,
                strike,
                forward,
                std_dev,
                discount,
            )
            .expect("valid inputs");
            let bumped_delta = bachelier_black_formula_forward_derivative(
                option_type,
                strike,
                forward + bump,
                std_dev,
                discount,
            )
            .expect("valid inputs");
            let base = bachelier_black_formula(option_type, strike, forward, std_dev, discount)
                .expect("valid inputs");
            let bumped =
                bachelier_black_formula(option_type, strike, forward + bump, std_dev, discount)
                    .expect("valid inputs");
            let delta_approx = (bumped - base) / bump;
            assert!(delta.max(bumped_delta) + epsilon > delta_approx);
            assert!(delta_approx > delta.min(bumped_delta) - epsilon);
        }
    }

    #[test]
    fn bachelier_forward_derivative_brackets_the_bumped_premium() {
        let strikes = [-3.0, -2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 3.0];
        assert_bachelier_forward_derivative(OptionType::Call, &strikes, 0.001);
        assert_bachelier_forward_derivative(OptionType::Put, &strikes, 0.001);
    }

    #[test]
    fn bachelier_forward_derivative_brackets_the_bumped_premium_zero_vol() {
        let strikes = [-3.0, -2.0, -1.0, -0.5, 0.0, 0.5, 1.0, 2.0, 3.0];
        assert_bachelier_forward_derivative(OptionType::Call, &strikes, 0.0);
        assert_bachelier_forward_derivative(OptionType::Put, &strikes, 0.0);
    }

    #[test]
    fn bachelier_std_dev_derivative_matches_bumped_value() {
        let forward = 1.0;
        let std_dev = 0.001 * 10.0_f64.sqrt();
        let discount = 0.95;
        let bump = 1.0e-6;
        for strike in [-1.0, 0.0, 0.9, 1.0, 1.1, 2.0] {
            let up = bachelier_black_formula(
                OptionType::Call,
                strike,
                forward,
                std_dev + bump,
                discount,
            )
            .expect("valid inputs");
            let down = bachelier_black_formula(
                OptionType::Call,
                strike,
                forward,
                std_dev - bump,
                discount,
            )
            .expect("valid inputs");
            let analytical =
                bachelier_black_formula_std_dev_derivative(strike, forward, std_dev, discount)
                    .expect("valid inputs");
            assert_close((up - down) / (2.0 * bump), analytical, 1e-6);
        }
    }

    #[test]
    fn bachelier_std_dev_derivative_edge_cases() {
        assert_eq!(
            bachelier_black_formula_std_dev_derivative(1.0, 1.0, 0.0, 0.95).expect("zero std dev"),
            0.0
        );
        assert!(bachelier_black_formula_std_dev_derivative(1.0, 1.0, -0.1, 0.95).is_err());
        assert!(bachelier_black_formula_std_dev_derivative(1.0, 1.0, 0.1, 0.0).is_err());
    }

    /// Round-trips vol -> price -> implied stddev across Call/Put and ITM/OTM,
    /// with and without a displacement, so both the put-call-parity conversion
    /// (the ITM cases flip to the OTM complement) and the shifted-lognormal
    /// branch are exercised. This is the supporting oracle for the enabler; the
    /// discriminating stripper round-trip is #576.
    #[test]
    fn implied_std_dev_round_trips_black_formula() {
        let discount = 0.95;
        let true_std_dev = 0.15;
        let accuracy = 1e-10;
        for displacement in [0.0, 0.01] {
            for (option_type, strike, forward) in [
                (OptionType::Call, 0.03, 0.025),
                (OptionType::Call, 0.02, 0.025),
                (OptionType::Put, 0.02, 0.025),
                (OptionType::Put, 0.03, 0.025),
            ] {
                let price = black_formula(
                    option_type,
                    strike,
                    forward,
                    true_std_dev,
                    discount,
                    displacement,
                )
                .expect("valid inputs");
                let recovered = black_formula_implied_std_dev(
                    option_type,
                    strike,
                    forward,
                    price,
                    discount,
                    displacement,
                    0.14,
                    accuracy,
                    100,
                )
                .expect("inverts black_formula");
                assert_close(recovered, true_std_dev, 1e-6);
            }
        }
    }

    #[test]
    fn implied_std_dev_rejects_a_negative_guess() {
        let price = black_formula(OptionType::Call, 0.03, 0.025, 0.15, 0.95, 0.0).expect("valid");
        assert!(
            black_formula_implied_std_dev(
                OptionType::Call,
                0.03,
                0.025,
                price,
                0.95,
                0.0,
                -1.0,
                1e-10,
                100
            )
            .is_err()
        );
    }

    #[test]
    fn implied_std_dev_rejects_invalid_price_and_discount() {
        assert!(
            black_formula_implied_std_dev(
                OptionType::Call,
                0.03,
                0.025,
                -0.1,
                0.95,
                0.0,
                0.14,
                1e-10,
                100
            )
            .is_err()
        );
        assert!(
            black_formula_implied_std_dev(
                OptionType::Call,
                0.03,
                0.025,
                0.01,
                0.0,
                0.0,
                0.14,
                1e-10,
                100
            )
            .is_err()
        );
    }

    #[test]
    fn itm_probabilities_with_zero_volatility_are_the_limit() {
        let cases = [
            (OptionType::Call, 0.03, 0.05),
            (OptionType::Call, 0.05, 0.03),
            (OptionType::Put, 0.03, 0.05),
            (OptionType::Put, 0.05, 0.03),
            (OptionType::Call, 90.0, 100.0),
            (OptionType::Put, 90.0, 100.0),
        ];
        let tiny_std_dev = 1.0e-8;
        for (option_type, strike, forward) in cases {
            let expected = if sign_of(option_type) * (forward - strike) > 0.0 {
                1.0
            } else {
                0.0
            };
            let probabilities = |std_dev| {
                [
                    black_formula_cash_itm_probability(option_type, strike, forward, std_dev, 0.0),
                    black_formula_asset_itm_probability(option_type, strike, forward, std_dev, 0.0),
                    bachelier_black_formula_asset_itm_probability(
                        option_type,
                        strike,
                        forward,
                        std_dev,
                    ),
                ]
            };
            for value in probabilities(tiny_std_dev)
                .into_iter()
                .chain(probabilities(0.0))
            {
                assert_close(value.expect("valid inputs"), expected, 1.0e-12);
            }
        }
    }

    #[test]
    fn bachelier_implied_vol_at_intrinsic_inverts() {
        let tte: Real = 0.01;
        let vol = 1e-5;
        for forward in [-0.01, 0.005, 0.05, 1.0] {
            for offset in [0.01, 0.03, 0.07] {
                for option_type in [OptionType::Call, OptionType::Put] {
                    let strike = match option_type {
                        OptionType::Call => forward - offset,
                        OptionType::Put => forward + offset,
                    };
                    let price = bachelier_black_formula(
                        option_type,
                        strike,
                        forward,
                        vol * tte.sqrt(),
                        1.0,
                    )
                    .expect("valid inputs");
                    let implied = bachelier_black_formula_implied_vol(
                        option_type,
                        strike,
                        forward,
                        tte,
                        price,
                        1.0,
                    )
                    .expect("a price at intrinsic must invert");
                    assert!(implied >= 0.0, "implied vol {implied}");
                }
            }
        }
    }

    #[test]
    fn bachelier_implied_vol_one_ulp_below_intrinsic_is_zero() {
        for (option_type, strike, forward) in [
            (OptionType::Call, 0.03, 0.05),
            (OptionType::Put, 0.05, 0.03),
            (OptionType::Call, -0.02, 0.01),
        ] {
            let intrinsic = sign_of(option_type) * (forward - strike);
            let price = intrinsic.next_down();
            let implied =
                bachelier_black_formula_implied_vol(option_type, strike, forward, 0.5, price, 1.0)
                    .expect("a price within close_enough of intrinsic must invert");
            assert_eq!(implied, 0.0);
        }
    }

    #[test]
    fn bachelier_implied_vol_one_ulp_above_intrinsic_is_zero() {
        for (option_type, strike, forward) in [
            (OptionType::Call, 0.03, 0.05),
            (OptionType::Put, 0.05, 0.03),
        ] {
            let intrinsic = sign_of(option_type) * (forward - strike);
            let price = intrinsic.next_up();
            let implied =
                bachelier_black_formula_implied_vol(option_type, strike, forward, 0.5, price, 1.0)
                    .expect("a price within close_enough of intrinsic must invert");
            assert_eq!(implied, 0.0);
        }
    }

    #[test]
    fn bachelier_implied_vol_of_a_negligible_out_of_the_money_premium_is_zero() {
        let price = 1e-30;
        for (option_type, strike, forward) in [
            (OptionType::Call, 0.05, 0.03),
            (OptionType::Put, 0.03, 0.05),
        ] {
            let implied =
                bachelier_black_formula_implied_vol(option_type, strike, forward, 0.5, price, 1.0)
                    .expect("a premium within close_enough of zero intrinsic must invert");
            assert_eq!(implied, 0.0);
        }
    }

    #[test]
    fn bachelier_implied_vol_near_the_money_takes_the_atm_closed_form() {
        let forward: Real = 0.05;
        let strike = forward.next_up();
        let (tte, price): (Real, Real) = (0.5, 0.004);
        let expected = price * (2.0 * std::f64::consts::PI).sqrt() / tte.sqrt();
        for option_type in [OptionType::Call, OptionType::Put] {
            let implied =
                bachelier_black_formula_implied_vol(option_type, strike, forward, tte, price, 1.0)
                    .expect("valid inputs");
            assert_eq!(implied.to_bits(), expected.to_bits());
        }
    }

    #[test]
    fn deep_in_the_money_prices_are_not_below_intrinsic() {
        for (option_type, strike, forward) in [
            (OptionType::Call, 50.0, 100.0),
            (OptionType::Put, 100.0, 50.0),
            (OptionType::Call, 0.03, 0.05),
            (OptionType::Put, 0.05, 0.03),
            (OptionType::Call, 1000.0, 1234.5),
            (OptionType::Put, 1234.5, 1000.0),
        ] {
            let moneyness = sign_of(option_type) * (forward - strike);
            for discount in [1.0, 0.97] {
                let intrinsic = discount * moneyness.max(0.0);
                let black = black_formula(option_type, strike, forward, 0.01, discount, 0.0)
                    .expect("valid inputs");
                assert!(
                    black >= intrinsic,
                    "Black {black} below intrinsic {intrinsic}"
                );
                let bachelier = bachelier_black_formula(
                    option_type,
                    strike,
                    forward,
                    moneyness / 20.0,
                    discount,
                )
                .expect("valid inputs");
                assert!(
                    bachelier >= intrinsic,
                    "Bachelier {bachelier} below intrinsic {intrinsic}"
                );
            }
        }
    }
}