sanos 0.2.1

// src/backbone/bs.rs
//
// Minimal Black–Scholes utilities for SANOS.
//
// Conventions:
// - Forward-normalized calls: F = 1
// - Strike is moneyness k > 0
// - Total variance var >= 0
//
// Formula:
//   C(k, var) = N(d1) - k N(d2)
//   d1 = (-ln k + 0.5 var) / sqrt(var)
//   d2 = d1 - sqrt(var)
//
// Edge case:
// - If var == 0, C = max(1 - k, 0)

use crate::error::{SanosError, SanosResult};
use statrs::distribution::{ContinuousCDF, Normal};

/// Standard normal CDF N(x).
#[inline]
pub fn norm_cdf(x: f64) -> SanosResult<f64> {
    if !x.is_finite() {
        return Err(SanosError::NonFinite {
            field: "x",
            value: x,
        });
    }
    // Safe unwrap: std normal params are valid.
    let n = Normal::new(0.0, 1.0).expect("Normal(0,1) must be constructible");
    Ok(n.cdf(x))
}

/// Forward-normalized Black–Scholes call price with forward F = 1.
///
/// Inputs:
/// - `k`: strike in forward moneyness (k > 0)
/// - `var`: total variance (var >= 0)
///
/// Returns:
/// - `C` in [0, 1] for k >= 0 (numerically, may have tiny eps violations)
pub fn bs_call_forward_norm(k: f64, var: f64) -> SanosResult<f64> {
    if !k.is_finite() {
        return Err(SanosError::NonFinite {
            field: "k",
            value: k,
        });
    }
    if !var.is_finite() {
        return Err(SanosError::NonFinite {
            field: "var",
            value: var,
        });
    }
    if k <= 0.0 {
        return Err(SanosError::InvalidBound {
            field: "k",
            value: k,
            min: f64::MIN_POSITIVE,
            max: f64::INFINITY,
        });
    }
    if var < 0.0 {
        return Err(SanosError::InvalidBound {
            field: "var",
            value: var,
            min: 0.0,
            max: f64::INFINITY,
        });
    }

    // Deterministic limit as var -> 0.
    if var == 0.0 {
        return Ok((1.0 - k).max(0.0));
    }

    let sqrt_var = var.sqrt();
    // d1, d2 with forward-normalized convention.
    let ln_k = k.ln();
    let d1 = (-ln_k + 0.5 * var) / sqrt_var;
    let d2 = d1 - sqrt_var;

    let n = Normal::new(0.0, 1.0).expect("Normal(0,1) must be constructible");
    let nd1 = n.cdf(d1);
    let nd2 = n.cdf(d2);

    Ok(nd1 - k * nd2)
}

/// Black–Scholes ATM implied total variance W from forward-normalized ATM call price C(1, W).
///
/// For k = 1 (ATM forward), we have:
///   C = N(d1) - N(d2) with d1 = 0.5*sqrt(W), d2 = -0.5*sqrt(W)
/// => C = 2*N(0.5*sqrt(W)) - 1
/// => N(0.5*sqrt(W)) = (1 + C)/2
/// => W = 4 * (Phi^{-1}((1+C)/2))^2
pub fn bs_implied_atm_var_from_call(call_atm: f64) -> SanosResult<f64> {
    if !call_atm.is_finite() {
        return Err(SanosError::NonFinite {
            field: "call_atm",
            value: call_atm,
        });
    }
    if !(0.0..=1.0).contains(&call_atm) {
        return Err(SanosError::InvalidBound {
            field: "call_atm",
            value: call_atm,
            min: 0.0,
            max: 1.0,
        });
    }

    let p = 0.5 * (1.0 + call_atm);

    // clamp away from {0,1} to avoid +/- infinity inverse_cdf
    let eps = 1e-15;
    let p = p.clamp(eps, 1.0 - eps);

    let n = Normal::new(0.0, 1.0).expect("Normal(0,1) must be constructible");
    let x = n.inverse_cdf(p); // x = 0.5*sqrt(W)
    Ok(4.0 * x * x)
}

/// Black implied volatility from a call price using Jaeckel's solver.
///
/// Inputs follow Black's forward convention:
/// - `forward` > 0
/// - `strike` > 0
/// - `maturity` > 0
/// - `call` in [intrinsic, forward]
#[cfg(feature = "iv-jaeckel")]
pub fn bs_implied_vol_from_call(
    call: f64,
    forward: f64,
    strike: f64,
    maturity: f64,
) -> SanosResult<f64> {
    if !call.is_finite() {
        return Err(SanosError::NonFinite {
            field: "call",
            value: call,
        });
    }
    if !forward.is_finite() {
        return Err(SanosError::NonFinite {
            field: "forward",
            value: forward,
        });
    }
    if !strike.is_finite() {
        return Err(SanosError::NonFinite {
            field: "strike",
            value: strike,
        });
    }
    if !maturity.is_finite() {
        return Err(SanosError::NonFinite {
            field: "maturity",
            value: maturity,
        });
    }
    if forward <= 0.0 {
        return Err(SanosError::InvalidBound {
            field: "forward",
            value: forward,
            min: f64::MIN_POSITIVE,
            max: f64::INFINITY,
        });
    }
    if strike <= 0.0 {
        return Err(SanosError::InvalidBound {
            field: "strike",
            value: strike,
            min: f64::MIN_POSITIVE,
            max: f64::INFINITY,
        });
    }
    if maturity <= 0.0 {
        return Err(SanosError::InvalidBound {
            field: "maturity",
            value: maturity,
            min: f64::MIN_POSITIVE,
            max: f64::INFINITY,
        });
    }

    let intrinsic = (forward - strike).max(0.0);
    if call < intrinsic - 1e-14 || call > forward + 1e-14 {
        return Err(SanosError::InvalidBound {
            field: "call",
            value: call,
            min: intrinsic,
            max: forward,
        });
    }

    let iv = jaeckel::implied_black_volatility(call, forward, strike, maturity, 1.0);
    if !iv.is_finite() || iv == f64::MAX || iv == -f64::MAX || iv < 0.0 {
        return Err(SanosError::External {
            msg: format!(
                "implied vol inversion failed (call={call}, forward={forward}, strike={strike}, maturity={maturity}, iv={iv})"
            ),
        });
    }

    Ok(iv)
}

#[cfg(not(feature = "iv-jaeckel"))]
pub fn bs_implied_vol_from_call(
    _call: f64,
    _forward: f64,
    _strike: f64,
    _maturity: f64,
) -> SanosResult<f64> {
    Err(SanosError::NotImplemented {
        what: "Enable feature `iv-jaeckel` to use implied volatility inversion.",
    })
}

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

    #[test]
    fn bs_call_var_zero_limit() {
        let c = bs_call_forward_norm(1.2, 0.0).unwrap();
        assert!((c - 0.0).abs() < 1e-15);

        let c = bs_call_forward_norm(0.8, 0.0).unwrap();
        assert!((c - 0.2).abs() < 1e-15);
    }

    #[test]
    fn bs_implied_atm_var_roundtrip() {
        let w = 0.09;
        let c = bs_call_forward_norm(1.0, w).unwrap();
        let w2 = bs_implied_atm_var_from_call(c).unwrap();
        assert!((w - w2).abs() < 1e-10);
    }

    #[cfg(feature = "iv-jaeckel")]
    #[test]
    fn bs_implied_vol_roundtrip() {
        let k = 1.15;
        let t = 0.7;
        let sigma = 0.32;
        let var = sigma * sigma * t;
        let price = bs_call_forward_norm(k, var).unwrap();
        let implied = bs_implied_vol_from_call(price, 1.0, k, t).unwrap();
        assert!((implied - sigma).abs() < 1e-10);
    }
}