regit_svi/

density.rs

// Copyright 2026 Regit.io — Nicolas Koenig
// SPDX-License-Identifier: Apache-2.0

//! Risk-neutral density implied by a raw SVI slice.
//!
//! The butterfly function `g(k)` (see [`crate::no_arb::butterfly`]) is the
//! sign-controlling factor in the risk-neutral density, not the density
//! itself. With
//!
//! ```text
//! d_minus(k) = -k/sqrt(w(k)) - sqrt(w(k))/2
//! d_plus(k)  = -k/sqrt(w(k)) + sqrt(w(k))/2
//! ```
//!
//! the risk-neutral density of the log-strike is
//!
//! ```text
//! p(k) = g(k) / sqrt(2*pi*w(k)) * exp(-d_minus(k)^2 / 2)
//! ```
//!
//! On a regular positive slice, `p(k) >= 0` for all `k` is equivalent to
//! `g(k) >= 0`. Unit continuous mass additionally requires the left-tail
//! boundary that excludes an atom at zero. [`integral`] is a bounded numerical
//! diagnostic for that continuous mass.
//!
//! # References
//!
//! - Breeden, D. & Litzenberger, R., "Prices of state-contingent claims
//!   implicit in option prices", *Journal of Business* 51(4):621-651 (1978).
//! - Gatheral, J. & Jacquier, A., "Arbitrage-free SVI volatility surfaces",
//!   *Quantitative Finance* 14(1):59-71 (2014), eq. (2.2).

use crate::no_arb::butterfly::{g, g_with_scale};
use crate::numerics::index_to_f64;
use crate::smile::raw::RawSvi;
use core::fmt;

/// `2*pi` — normalising constant for the density.
const TWO_PI: f64 = std::f64::consts::TAU;

/// The Black `d_minus` quantity: `d_minus(k) = -k/sqrt(w) - sqrt(w)/2`.
///
/// This infallible kernel assumes finite `k` and positive finite `w(k)` and
/// propagates IEEE non-finite results when those preconditions fail.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::density::d_minus;
///
/// // At k = 0 with total variance w, d_minus = -sqrt(w)/2.
/// let svi = RawSvi::new(0.04, 0.0, 0.0, 0.0, 0.1)?;
/// assert!((d_minus(&svi, 0.0) + 0.1).abs() < 1e-12);
/// # Ok(())
/// # }
/// ```
#[must_use]
#[inline]
pub fn d_minus(svi: &RawSvi, k: f64) -> f64 {
    let w = svi.total_variance(k);
    let sqrt_w = w.sqrt();
    -k / sqrt_w - sqrt_w / 2.0
}

/// The Black `d_plus` quantity: `d_plus(k) = -k/sqrt(w) + sqrt(w)/2`.
///
/// This infallible kernel assumes finite `k` and positive finite `w(k)` and
/// propagates IEEE non-finite results when those preconditions fail.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::density::d_plus;
///
/// // At k = 0 with total variance w, d_plus = +sqrt(w)/2.
/// let svi = RawSvi::new(0.04, 0.0, 0.0, 0.0, 0.1)?;
/// assert!((d_plus(&svi, 0.0) - 0.1).abs() < 1e-12);
/// # Ok(())
/// # }
/// ```
#[must_use]
#[inline]
pub fn d_plus(svi: &RawSvi, k: f64) -> f64 {
    let w = svi.total_variance(k);
    let sqrt_w = w.sqrt();
    -k / sqrt_w + sqrt_w / 2.0
}

/// The risk-neutral density `p(k)` implied by a raw SVI slice (MATH.md §9).
///
/// ```text
/// p(k) = g(k) / sqrt(2*pi*w(k)) * exp(-d_minus(k)^2 / 2)
/// ```
///
/// On a regular positive slice, a negative `p(k)` is a butterfly-arbitrage
/// witness. Non-negativity of the factor must be combined with the applicable
/// call-price tail condition for a full slice conclusion.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::density::risk_neutral_density;
///
/// // A benign slice has a positive density at the money.
/// let svi = RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3)?;
/// assert!(risk_neutral_density(&svi, 0.0)? > 0.0);
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`DensityError::InvalidDomain`] for non-finite `k`,
/// [`DensityError::NonPositiveVariance`] at a zero/singular variance point,
/// [`DensityError::IllConditionedVariance`] when a positive variance is inside
/// the scale-aware cancellation band,
/// or [`DensityError::NonFiniteEvaluation`] when binary64 evaluation overflows.
pub fn risk_neutral_density(svi: &RawSvi, k: f64) -> Result<f64, DensityError> {
    if !k.is_finite() {
        return Err(DensityError::InvalidDomain);
    }
    let w = svi.total_variance(k);
    if !w.is_finite() {
        return Err(DensityError::NonFiniteEvaluation { k });
    }
    if w <= 0.0 {
        return Err(DensityError::NonPositiveVariance { k, w });
    }
    let variance_scale = 1.0 + svi.a().abs() + svi.b().abs() * svi.sigma().abs();
    let near_minimum = (k - svi.k_min()).abs() <= 16.0 * f64::EPSILON * (1.0 + k.abs());
    let uncertainty = 16.0 * f64::EPSILON * variance_scale;
    if svi.b() > 0.0 && near_minimum && w <= uncertainty {
        return Err(DensityError::IllConditionedVariance { k, w, uncertainty });
    }
    let dm = d_minus(svi, k);
    let value = g(svi, k) / (TWO_PI * w).sqrt() * (-0.5 * dm * dm).exp();
    if value.is_finite() {
        Ok(value)
    } else {
        Err(DensityError::NonFiniteEvaluation { k })
    }
}

/// Error returned when density evaluation or integration has no finite domain.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum DensityError {
    /// The requested range, panel count, or point is invalid.
    InvalidDomain,
    /// Total variance is non-positive at the requested log-moneyness.
    NonPositiveVariance {
        /// Log-moneyness of the singularity.
        k: f64,
        /// Evaluated total variance.
        w: f64,
    },
    /// Positive variance is too close to floating-point cancellation for a
    /// reliable regular-density evaluation.
    IllConditionedVariance {
        /// Log-moneyness of the ill-conditioned evaluation.
        k: f64,
        /// Positive total variance obtained in binary64.
        w: f64,
        /// Scale-aware uncertainty threshold used for classification.
        uncertainty: f64,
    },
    /// Floating-point evaluation overflowed or produced `NaN`.
    NonFiniteEvaluation {
        /// Log-moneyness of the failed evaluation.
        k: f64,
    },
}

impl fmt::Display for DensityError {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::InvalidDomain => {
                write!(f, "density domain must be finite, ordered, and non-empty")
            }
            Self::NonPositiveVariance { k, w } => {
                write!(f, "density is singular at k={k}: total variance is {w}")
            }
            Self::IllConditionedVariance { k, w, uncertainty } => write!(
                f,
                "density is ill-conditioned at k={k}: total variance {w} is within uncertainty {uncertainty}"
            ),
            Self::NonFiniteEvaluation { k } => {
                write!(f, "density evaluation is non-finite at k={k}")
            }
        }
    }
}

impl std::error::Error for DensityError {}

/// A diagnostic report on the risk-neutral density of a raw SVI slice.
///
/// When the continuous positive-strike mass is one, [`integral`] over a wide
/// window should approach `1`. A value far from `1` can reflect truncation, a
/// zero-strike atom, arbitrage, or numerical ill-conditioning.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct DensityReport {
    /// Requested lower log-moneyness bound.
    lower: f64,
    /// Requested upper log-moneyness bound.
    upper: f64,
    /// Number of panels in the finer composite-Simpson estimate.
    panels: usize,
    /// Numerical integral of `p` over the diagnostic window.
    integral: f64,
    /// Richardson error estimate from nested composite-Simpson grids.
    integration_error: f64,
    /// Smallest density value observed on the integration grid.
    min_density: f64,
    /// Whether a negative density was observed on the bounded grid.
    violation_observed: bool,
}

impl DensityReport {
    /// Returns the requested finite integration bounds.
    #[must_use]
    pub const fn domain(self) -> (f64, f64) {
        (self.lower, self.upper)
    }
    /// Returns the panel count used by the reported fine-grid integral.
    #[must_use]
    pub const fn panels(self) -> usize {
        self.panels
    }
    /// Returns the bounded numerical integral.
    #[must_use]
    pub const fn integral(self) -> f64 {
        self.integral
    }
    /// Returns the nested-grid Simpson error estimate.
    ///
    /// This estimates quadrature discretization error on the declared finite
    /// window; it does not include probability mass outside that window.
    #[must_use]
    pub const fn integration_error(self) -> f64 {
        self.integration_error
    }
    /// Returns the minimum sampled density.
    #[must_use]
    pub const fn min_density(self) -> f64 {
        self.min_density
    }
    /// Returns whether the bounded grid contained a negative density.
    #[must_use]
    pub const fn violation_observed(self) -> bool {
        self.violation_observed
    }
}

/// Numerically integrates the risk-neutral density over `[k_lo, k_hi]` by the
/// composite Simpson rule with `2n` panels.
///
/// When the continuous positive-strike mass is one, a sufficiently wide
/// window approaches `1`. The caller controls the window and panel count;
/// this routine makes no truncation-error guarantee.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::density::integral;
///
/// // A benign, low-variance slice integrates to near 1 over a wide window.
/// let svi = RawSvi::new(0.04, 0.05, -0.1, 0.0, 0.4)?;
/// let mass = integral(&svi, -6.0, 6.0, 2000)?;
/// assert!((mass - 1.0).abs() < 1e-2, "mass = {mass}");
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`DensityError::InvalidDomain`] for invalid bounds/counts and
/// propagates pointwise density errors from [`risk_neutral_density`].
pub fn integral(svi: &RawSvi, k_lo: f64, k_hi: f64, n: usize) -> Result<f64, DensityError> {
    if !k_lo.is_finite() || !k_hi.is_finite() || k_lo >= k_hi || n == 0 || n > usize::MAX / 2 {
        return Err(DensityError::InvalidDomain);
    }
    let panels = 2 * n;
    let h = (k_hi - k_lo) / index_to_f64(panels);
    let mut sum = risk_neutral_density(svi, k_lo)? + risk_neutral_density(svi, k_hi)?;
    if !sum.is_finite() {
        return Err(DensityError::NonFiniteEvaluation { k: k_lo });
    }
    for i in 1..panels {
        let k = h.mul_add(index_to_f64(i), k_lo);
        let weight = if i % 2 == 1 { 4.0 } else { 2.0 };
        sum += weight * risk_neutral_density(svi, k)?;
        if !sum.is_finite() {
            return Err(DensityError::NonFiniteEvaluation { k });
        }
    }
    let value = sum * h / 3.0;
    if value.is_finite() {
        Ok(value)
    } else {
        Err(DensityError::NonFiniteEvaluation { k: k_hi })
    }
}

/// Builds a [`DensityReport`] for a slice over `[k_lo, k_hi]`.
///
/// Integrates the density and samples the same bounded grid for its minimum.
/// A clean report is numerical diagnostic evidence, not a global conclusion.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::smile::raw::RawSvi;
/// use regit_svi::density::density_report;
///
/// let svi = RawSvi::new(0.04, 0.05, -0.1, 0.0, 0.4)?;
/// let report = density_report(&svi, -6.0, 6.0, 2000)?;
/// assert!(!report.violation_observed());
/// assert!((report.integral() - 1.0).abs() < 1e-2);
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`DensityError::InvalidDomain`] for invalid bounds/counts and
/// propagates pointwise density errors from [`risk_neutral_density`].
pub fn density_report(
    svi: &RawSvi,
    k_lo: f64,
    k_hi: f64,
    n: usize,
) -> Result<DensityReport, DensityError> {
    if !k_lo.is_finite() || !k_hi.is_finite() || k_lo >= k_hi || n == 0 || n > usize::MAX / 4 {
        return Err(DensityError::InvalidDomain);
    }
    let panels = 4 * n;
    let h = (k_hi - k_lo) / index_to_f64(panels);

    let mut min_density = f64::INFINITY;
    let mut violation_observed = false;
    for i in 0..=panels {
        let k = h.mul_add(index_to_f64(i), k_lo);
        let p = risk_neutral_density(svi, k)?;
        if p < min_density {
            min_density = p;
        }
        let (density_factor, evaluation_scale) = g_with_scale(svi, k);
        violation_observed |= density_factor < -128.0 * f64::EPSILON * evaluation_scale;
    }

    let coarse_integral = integral(svi, k_lo, k_hi, n)?;
    let fine_integral = integral(svi, k_lo, k_hi, 2 * n)?;
    let integration_error = (fine_integral - coarse_integral).abs() / 15.0;
    if !integration_error.is_finite() {
        return Err(DensityError::NonFiniteEvaluation { k: k_hi });
    }
    Ok(DensityReport {
        lower: k_lo,
        upper: k_hi,
        panels,
        integral: fine_integral,
        integration_error,
        min_density,
        violation_observed,
    })
}

#[cfg(test)]
#[allow(clippy::expect_used)] // Validated fixtures use contextual expectations.
mod tests {
    use super::*;

    #[test]
    fn d_plus_d_minus_differ_by_sqrt_w() {
        let svi =
            RawSvi::new(0.04, 0.2, -0.3, 0.05, 0.12).expect("valid test or documentation fixture");
        for &k in &[-0.5, 0.0, 0.3] {
            let w = svi.total_variance(k);
            assert!((d_plus(&svi, k) - d_minus(&svi, k) - w.sqrt()).abs() < 1e-12);
        }
    }

    #[test]
    fn density_positive_for_benign_slice() {
        let svi =
            RawSvi::new(0.04, 0.1, -0.2, 0.0, 0.3).expect("valid test or documentation fixture");
        for &k in &[-1.0, -0.3, 0.0, 0.3, 1.0] {
            assert!(
                risk_neutral_density(&svi, k).expect("valid test or documentation fixture") > 0.0,
                "p({k})"
            );
        }
    }

    #[test]
    fn density_integrates_to_one() {
        let svi =
            RawSvi::new(0.04, 0.05, -0.1, 0.0, 0.4).expect("valid test or documentation fixture");
        let mass = integral(&svi, -8.0, 8.0, 4000).expect("valid test or documentation fixture");
        assert!((mass - 1.0).abs() < 1e-3, "mass = {mass}");
    }

    #[test]
    fn density_integrates_to_one_low_vol() {
        let svi =
            RawSvi::new(0.02, 0.04, -0.15, 0.0, 0.3).expect("valid test or documentation fixture");
        let mass = integral(&svi, -6.0, 6.0, 4000).expect("valid test or documentation fixture");
        assert!((mass - 1.0).abs() < 1e-3, "mass = {mass}");
    }

    #[test]
    fn density_report_benign_slice() {
        let svi =
            RawSvi::new(0.04, 0.05, -0.1, 0.0, 0.4).expect("valid test or documentation fixture");
        let report =
            density_report(&svi, -8.0, 8.0, 4000).expect("valid test or documentation fixture");
        assert!(!report.violation_observed());
        assert!((report.integral() - 1.0).abs() < 1e-3);
        assert!(report.integration_error().is_finite());
        assert!(report.integration_error() >= 0.0);
        let (lower, upper) = report.domain();
        assert!((lower + 8.0).abs() < f64::EPSILON);
        assert!((upper - 8.0).abs() < f64::EPSILON);
        assert_eq!(report.panels(), 16_000);
        assert!(report.min_density() >= 0.0);
    }

    #[test]
    fn density_report_flags_vogt_slice() {
        // The Vogt slice has butterfly arbitrage -> negative density region.
        let vogt = RawSvi::new(-0.0410, 0.1331, 0.3060, 0.3586, 0.4153)
            .expect("valid test or documentation fixture");
        let report =
            density_report(&vogt, -2.0, 2.0, 2000).expect("valid test or documentation fixture");
        assert!(report.violation_observed());
        assert!(report.min_density() < 0.0);
    }

    #[test]
    fn density_handles_zero_variance_gracefully() {
        // A slice whose w_min is exactly 0 should not produce NaN.
        let svi = RawSvi::new(-0.125, 0.5, 0.0, 0.0, 0.25).expect("valid exact-zero fixture");
        assert!(matches!(
            risk_neutral_density(&svi, svi.k_min()),
            Err(DensityError::NonPositiveVariance { .. })
        ));
    }

    #[test]
    fn tiny_positive_variance_is_not_classified_as_non_positive() {
        let slice =
            RawSvi::new(-0.02 + 1e-16, 0.1, 0.0, 0.0, 0.2).expect("valid tiny-positive fixture");
        assert!(slice.w_min() > 0.0);
        assert!(matches!(
            risk_neutral_density(&slice, slice.k_min()),
            Err(DensityError::IllConditionedVariance { w, .. }) if w > 0.0
        ));
    }

    #[test]
    fn non_finite_variance_is_not_classified_as_non_positive() {
        let slice = RawSvi::new_unchecked(0.0, f64::MAX, 0.0, -f64::MAX, f64::MIN_POSITIVE);
        assert!(matches!(
            risk_neutral_density(&slice, 0.0),
            Err(DensityError::NonFiniteEvaluation { k: 0.0 })
        ));
    }
}