regit-svi 2.0.0

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

//! Surface SVI (SSVI) parametrisation.
//!
//! SSVI parametrises the entire surface at once, as a function of
//! log-moneyness `k` and the ATM total variance `theta = w(0, t)`:
//!
//! ```text
//! w(k, theta) = (theta / 2)
//!             * ( 1 + rho*phi(theta)*k
//!                   + sqrt( (phi(theta)*k + rho)^2 + (1 - rho^2) ) )
//! ```
//!
//! The free objects are a single global correlation `rho in (-1, 1)`, the ATM
//! total-variance term structure `theta_t`, and a smoothing function
//! `phi : R+ -> R+`. Two standard choices of `phi` are supported (MATH.md §6):
//!
//! ```text
//! Heston-like:  phi(theta) = (1/(lambda*theta))
//!                          * (1 - (1 - exp(-lambda*theta))/(lambda*theta))
//! Power-law:    phi(theta) = eta / (theta^gamma * (1 + theta)^(1 - gamma))
//! ```
//!
//! For fixed `theta`, an SSVI slice equals a raw SVI slice — the closed-form
//! map is [`Ssvi::slice_at`].
//!
//! SSVI's value is that static-arbitrage conditions admit closed-form
//! inequalities (Gatheral & Jacquier 2014, Theorems 4.1 and 4.2), checked here by
//! [`Ssvi::butterfly_assessment`] and [`Ssvi::calendar_assessment`].
//!
//! # References
//!
//! - Gatheral, J. & Jacquier, A., "Arbitrage-free SVI volatility surfaces",
//!   *Quantitative Finance* 14(1):59-71 (2014), Section 4.

use crate::error::ParamError;
use crate::no_arb::evidence::{ArbitrageAssessment, ArbitrageEvidence, ArbitrageStatus};
use crate::smile::raw::{RawSvi, stable_shape};

/// SSVI smoothing function `phi : R+ -> R+`.
///
/// `phi` controls how the smile curvature evolves with maturity. Both
/// supported forms are positive and decreasing in `theta`; arbitrage evidence
/// is assessed separately because positivity alone is not sufficient.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::surface::ssvi::Phi;
///
/// let heston = Phi::heston(1.0)?;
/// let power = Phi::modified_power_law(0.5, 0.5)?;
/// assert!(heston.eval(0.04) > 0.0);
/// assert!(power.eval(0.04) > 0.0);
/// # Ok(())
/// # }
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Phi(PhiKind);

#[derive(Debug, Clone, Copy, PartialEq)]
enum PhiKind {
    Heston { lambda: f64 },
    ModifiedPowerLaw { eta: f64, gamma: f64 },
}

impl Phi {
    /// Creates a validated Heston-like smoothing function.
    ///
    /// # Errors
    ///
    /// Returns [`ParamError::InvalidPhiParameter`] if `lambda <= 0`, or
    /// [`ParamError::NonFinite`] if `lambda` is not finite.
    ///
    /// # Examples
    ///
    /// ```
    /// use regit_svi::surface::ssvi::Phi;
    ///
    /// assert!(Phi::heston(1.5).is_ok());
    /// assert!(Phi::heston(0.0).is_err());
    /// ```
    pub fn heston(lambda: f64) -> Result<Self, ParamError> {
        if !lambda.is_finite() {
            return Err(ParamError::NonFinite { name: "lambda" });
        }
        if lambda <= 0.0 {
            return Err(ParamError::InvalidPhiParameter {
                name: "lambda",
                value: lambda,
            });
        }
        Ok(Self(PhiKind::Heston { lambda }))
    }

    /// Creates a validated power-law smoothing function.
    ///
    /// # Errors
    ///
    /// Returns [`ParamError::InvalidPhiParameter`] if `eta <= 0` or
    /// `gamma` is outside `(0, 1)`, or [`ParamError::NonFinite`] if either
    /// parameter is not finite.
    ///
    /// # Examples
    ///
    /// ```
    /// use regit_svi::surface::ssvi::Phi;
    ///
    /// assert!(Phi::modified_power_law(0.5, 0.5).is_ok());
    /// assert!(Phi::modified_power_law(0.5, 1.0).is_err());
    /// ```
    pub fn modified_power_law(eta: f64, gamma: f64) -> Result<Self, ParamError> {
        if !eta.is_finite() {
            return Err(ParamError::NonFinite { name: "eta" });
        }
        if !gamma.is_finite() {
            return Err(ParamError::NonFinite { name: "gamma" });
        }
        if eta <= 0.0 {
            return Err(ParamError::InvalidPhiParameter {
                name: "eta",
                value: eta,
            });
        }
        if gamma <= 0.0 || gamma >= 1.0 {
            return Err(ParamError::InvalidPhiParameter {
                name: "gamma",
                value: gamma,
            });
        }
        Ok(Self(PhiKind::ModifiedPowerLaw { eta, gamma }))
    }

    /// Evaluates `phi(theta)` for `theta > 0`.
    ///
    /// The Heston-like form uses `x = lambda*theta` and
    /// `phi = (1/x)*(1 - (1 - exp(-x))/x)`; the power-law form uses
    /// `phi = eta / (theta^gamma * (1 + theta)^(1-gamma))`.
    /// This infallible numerical kernel assumes finite, positive `theta` and
    /// propagates IEEE non-finite results outside that domain or range.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::Phi;
    ///
    /// // Power-law at theta = 1: phi = eta / (1 * 2^(1-gamma)).
    /// let p = Phi::modified_power_law(1.0, 0.5)?;
    /// assert!((p.eval(1.0) - 1.0 / 2.0_f64.powf(0.5)).abs() < 1e-12);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    pub fn eval(&self, theta: f64) -> f64 {
        match self.0 {
            PhiKind::Heston { lambda } => {
                let x = lambda * theta;
                if x.abs() <= 1.0e-4 {
                    0.5 + x * (-1.0 / 6.0 + x * (1.0 / 24.0 + x * (-1.0 / 120.0 + x / 720.0)))
                } else {
                    (x + (-x).exp_m1()) / (x * x)
                }
            }
            PhiKind::ModifiedPowerLaw { eta, gamma } => {
                eta / (theta.powf(gamma) * (1.0 + theta).powf(1.0 - gamma))
            }
        }
    }

    /// Returns the Heston-like `lambda`, or `None` for another family.
    #[must_use]
    pub const fn heston_lambda(self) -> Option<f64> {
        match self.0 {
            PhiKind::Heston { lambda } => Some(lambda),
            PhiKind::ModifiedPowerLaw { .. } => None,
        }
    }

    /// Returns modified-power-law `(eta, gamma)`, or `None` for another family.
    #[must_use]
    pub const fn modified_power_law_parameters(self) -> Option<(f64, f64)> {
        match self.0 {
            PhiKind::ModifiedPowerLaw { eta, gamma } => Some((eta, gamma)),
            PhiKind::Heston { .. } => None,
        }
    }
}

/// A Surface SVI parametrisation: a global correlation plus a smoothing
/// function.
///
/// The ATM total-variance term structure `theta_t` is supplied per evaluation
/// (it is interpolated from market quotes by the surface layer), so an
/// [`Ssvi`] value is the maturity-independent part of the surface.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::surface::ssvi::{Phi, Ssvi};
///
/// let ssvi = Ssvi::new(-0.4, Phi::modified_power_law(0.5, 0.5)?)?;
/// let w = ssvi.total_variance(0.1, 0.04);
/// assert!(w > 0.0);
/// # Ok(())
/// # }
/// ```
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Ssvi {
    /// Global correlation `rho in (-1, 1)`.
    pub(crate) rho: f64,
    /// The smoothing function `phi`.
    pub(crate) phi: Phi,
}

impl Ssvi {
    /// Creates a validated SSVI surface.
    ///
    /// # Errors
    ///
    /// - [`ParamError::NonFinite`] if `rho` is not finite.
    /// - [`ParamError::CorrelationOutOfRange`] if `|rho| >= 1`.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// assert!(Ssvi::new(-0.4, Phi::heston(1.0)?).is_ok());
    /// assert!(Ssvi::new(1.0, Phi::heston(1.0)?).is_err());
    /// # Ok(())
    /// # }
    /// ```
    pub fn new(rho: f64, phi: Phi) -> Result<Self, ParamError> {
        if !rho.is_finite() {
            return Err(ParamError::NonFinite { name: "rho" });
        }
        if rho.abs() >= 1.0 {
            return Err(ParamError::CorrelationOutOfRange { rho });
        }
        Ok(Self { rho, phi })
    }

    /// Returns the global correlation parameter.
    #[must_use]
    pub const fn rho(self) -> f64 {
        self.rho
    }

    /// Returns the validated smoothing function.
    #[must_use]
    pub const fn phi(self) -> Phi {
        self.phi
    }

    /// Total implied variance `w(k, theta)` from the SSVI form.
    ///
    /// This infallible numerical kernel assumes finite `k` and finite,
    /// positive `theta`; it propagates IEEE non-finite results outside that
    /// domain or range. [`Self::slice_at`] provides validated slice mapping.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(0.0, Phi::modified_power_law(0.5, 0.5)?)?;
    /// // With rho = 0 and k = 0, w = theta/2 * (1 + sqrt(1)) = theta.
    /// assert!((ssvi.total_variance(0.0, 0.04) - 0.04).abs() < 1e-15);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    pub fn total_variance(&self, k: f64, theta: f64) -> f64 {
        let phi = self.phi.eval(theta);
        let pk = phi * k;
        let one_minus_rho2 = 1.0 - self.rho * self.rho;
        let shifted = pk + self.rho;
        let rho_scale = one_minus_rho2.sqrt();
        let radial = shifted.hypot(rho_scale);
        let bracket = one_minus_rho2 + stable_shape(shifted, rho_scale, self.rho, radial);
        (theta / 2.0) * bracket
    }

    /// The raw SVI slice equal to this SSVI surface at fixed `theta`.
    ///
    /// Closed-form map (MATH.md §6):
    ///
    /// ```text
    /// a     = (theta/2) * (1 - rho^2)
    /// b     = theta * phi / 2
    /// rho   = rho
    /// m     = -rho / phi
    /// sigma = sqrt(1 - rho^2) / phi
    /// ```
    ///
    /// # Errors
    ///
    /// Returns [`ParamError::NonPositiveTheta`] if `theta` is non-finite or
    /// non-positive, [`ParamError::InvalidPhiParameter`] if `phi(theta)` is
    /// non-finite or non-positive, and [`ParamError::NonFinite`],
    /// [`ParamError::NonPositiveSigma`], or
    /// [`ParamError::NegativeMinVariance`] if the mapped Raw parameters or
    /// their validated derived invariants are not representable.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
    /// let raw = ssvi.slice_at(0.04)?;
    /// // The raw slice reproduces the SSVI total variance.
    /// let direct = ssvi.total_variance(0.1, 0.04);
    /// assert!((raw.total_variance(0.1) - direct).abs() < 1e-12);
    /// # Ok(())
    /// # }
    /// ```
    pub fn slice_at(&self, theta: f64) -> Result<RawSvi, ParamError> {
        if theta <= 0.0 || !theta.is_finite() {
            return Err(ParamError::NonPositiveTheta { theta });
        }
        let phi = self.phi.eval(theta);
        if phi <= 0.0 || !phi.is_finite() {
            return Err(ParamError::InvalidPhiParameter {
                name: "phi(theta)",
                value: phi,
            });
        }
        let one_minus_rho2 = 1.0 - self.rho * self.rho;
        let a = (theta / 2.0) * one_minus_rho2;
        let b = theta * phi / 2.0;
        let m = -self.rho / phi;
        let sigma = one_minus_rho2.sqrt() / phi;
        RawSvi::new(a, b, self.rho, m, sigma)
    }

    /// Tests the SSVI sufficient no-butterfly-arbitrage condition at one
    /// `theta` (Gatheral & Jacquier 2014, Theorem 4.2).
    ///
    /// The slice at `theta` is free of butterfly arbitrage if both hold:
    ///
    /// ```text
    /// theta * phi      * (1 + |rho|) < 4
    /// theta * phi^2    * (1 + |rho|) <= 4
    /// ```
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
    /// assert_eq!(ssvi.butterfly_assessment_at(0.04).status(), regit_svi::ArbitrageStatus::NoViolationDetected);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    pub fn butterfly_assessment_at(&self, theta: f64) -> ArbitrageAssessment {
        let evidence = ArbitrageEvidence::AnalyticSufficient {
            theorem: "Gatheral–Jacquier Theorem 4.2",
        };
        if theta <= 0.0 || !theta.is_finite() {
            return ArbitrageAssessment::new(
                ArbitrageStatus::Indeterminate,
                evidence,
                f64::NAN,
                None,
            );
        }
        let phi = self.phi.eval(theta);
        if phi <= 0.0 || !phi.is_finite() {
            return ArbitrageAssessment::new(
                ArbitrageStatus::Indeterminate,
                evidence,
                f64::NAN,
                None,
            );
        }
        let factor = 1.0 + self.rho.abs();
        let tp = theta * phi;
        let first = tp * factor;
        let second = tp * phi * factor;
        let margin = (4.0 - first).min(4.0 - second);
        let boundary_band = 128.0 * f64::EPSILON * (1.0 + first.abs().max(second.abs()).max(4.0));
        let status = if 4.0 - first > boundary_band && 4.0 - second > boundary_band {
            ArbitrageStatus::NoViolationDetected
        } else {
            ArbitrageStatus::Indeterminate
        };
        ArbitrageAssessment::new(status, evidence, margin, None)
    }

    /// Tests the SSVI sufficient no-butterfly condition across a set of
    /// `theta` values (Gatheral & Jacquier 2014, Theorem 4.2).
    ///
    /// Every supplied value is checked. This is a conclusion only on the
    /// declared set of `theta` values, not between them or beyond them.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
    /// assert_eq!(ssvi.butterfly_assessment(&[0.01, 0.04, 0.09]).status(), regit_svi::ArbitrageStatus::NoViolationDetected);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    pub fn butterfly_assessment(&self, thetas: &[f64]) -> ArbitrageAssessment {
        if thetas.is_empty() {
            return ArbitrageAssessment::new(
                ArbitrageStatus::Indeterminate,
                ArbitrageEvidence::AnalyticSufficient {
                    theorem: "Gatheral–Jacquier Theorem 4.2 on declared theta support",
                },
                f64::NAN,
                None,
            );
        }
        let mut margin = f64::INFINITY;
        for &theta in thetas {
            let assessment = self.butterfly_assessment_at(theta);
            margin = margin.min(assessment.margin());
            if assessment.status() != ArbitrageStatus::NoViolationDetected {
                return ArbitrageAssessment::new(
                    ArbitrageStatus::Indeterminate,
                    assessment.evidence(),
                    margin,
                    None,
                );
            }
        }
        ArbitrageAssessment::new(
            ArbitrageStatus::NoViolationDetected,
            ArbitrageEvidence::AnalyticSufficient {
                theorem: "Gatheral–Jacquier Theorem 4.2 on declared theta support",
            },
            margin,
            None,
        )
    }

    /// Assesses the global positive-theta Theorem 4.2 envelope for the
    /// supported smoothing family.
    ///
    /// Passing is analytic sufficient evidence. Failing a sufficient envelope
    /// is [`ArbitrageStatus::Indeterminate`], never a detected violation.
    #[must_use]
    pub fn global_butterfly_assessment(&self) -> ArbitrageAssessment {
        let evidence = ArbitrageEvidence::AnalyticSufficient {
            theorem: "Gatheral–Jacquier Theorem 4.2 global family envelope",
        };
        let c = 1.0 + self.rho.abs();
        let margin = if let Some(lambda) = self.phi.heston_lambda() {
            lambda - c / 4.0
        } else if let Some((eta, gamma)) = self.phi.modified_power_law_parameters() {
            if gamma > 0.5 {
                return ArbitrageAssessment::new(
                    ArbitrageStatus::Indeterminate,
                    evidence,
                    f64::NEG_INFINITY,
                    None,
                );
            }
            let h_gamma = if (gamma - 0.5).abs() <= 64.0 * f64::EPSILON {
                1.0
            } else {
                let x = 1.0 - 2.0 * gamma;
                x.powf(x) / (2.0 - 2.0 * gamma).powf(2.0 - 2.0 * gamma)
            };
            (4.0 / c).min(2.0 / (c * h_gamma).sqrt()) - eta
        } else {
            f64::NAN
        };
        let boundary_band = 128.0
            * f64::EPSILON
            * (1.0
                + margin.abs()
                + self
                    .phi
                    .heston_lambda()
                    .unwrap_or_else(|| {
                        self.phi
                            .modified_power_law_parameters()
                            .map_or(0.0, |(eta, _)| eta)
                    })
                    .abs());
        let status = if margin > boundary_band {
            ArbitrageStatus::NoViolationDetected
        } else {
            ArbitrageStatus::Indeterminate
        };
        ArbitrageAssessment::new(status, evidence, margin, None)
    }

    /// Tests the SSVI no-calendar-spread-arbitrage condition at one `theta`
    /// (Gatheral & Jacquier 2014, Theorem 4.1, condition (ii)).
    ///
    /// Given a non-decreasing `theta_t` curve, the surface is calendar-free
    /// if, for every `theta`,
    ///
    /// ```text
    /// 0 <= d(theta * phi)/d(theta)
    ///        <= (1/rho^2) * (1 + sqrt(1 - rho^2)) * phi
    /// ```
    ///
    /// The upper bound is `+infinity` when `rho = 0`. For both supported phi
    /// families the derivative inequalities hold globally on `theta > 0` by
    /// direct differentiation; no finite theta sampling is used as a proof.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
    /// assert_eq!(ssvi.calendar_assessment(&[0.04]).status(), regit_svi::ArbitrageStatus::NoViolationDetected);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    fn calendar_condition_at(&self, theta: f64) -> bool {
        if theta <= 0.0 || !theta.is_finite() {
            return false;
        }
        match self.phi.0 {
            PhiKind::ModifiedPowerLaw { gamma, .. } => {
                // d(theta*phi)/dtheta = phi*(1-gamma)/(1+theta), whose
                // ratio to phi lies in (0,1). The theorem's upper multiplier
                // is at least one for every |rho|<1.
                gamma > 0.0 && gamma < 1.0
            }
            PhiKind::Heston { lambda } => {
                // The Heston-like form is the admissible family constructed
                // in Gatheral–Jacquier §4.1; validated lambda>0 satisfies the
                // derivative condition globally.
                lambda > 0.0
            }
        }
    }

    /// Tests the SSVI no-calendar condition across a set of `theta` values
    /// (Gatheral & Jacquier 2014, Theorem 4.1).
    ///
    /// Verifies condition (i) — `theta_t` non-decreasing — on the supplied
    /// sequence, then condition (ii) at each `theta`.
    ///
    /// # Examples
    ///
    /// ```
    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
    /// use regit_svi::surface::ssvi::{Phi, Ssvi};
    ///
    /// let ssvi = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
    /// assert_eq!(ssvi.calendar_assessment(&[0.01, 0.04, 0.09]).status(), regit_svi::ArbitrageStatus::NoViolationDetected);
    /// # Ok(())
    /// # }
    /// ```
    #[must_use]
    pub fn calendar_assessment(&self, thetas: &[f64]) -> ArbitrageAssessment {
        const BOUNDARY_TOLERANCE: f64 = 1e-12;
        let evidence = ArbitrageEvidence::AnalyticNecessaryAndSufficient {
            theorem: "Gatheral–Jacquier Theorem 4.1 for supported phi family",
            boundary_tolerance: BOUNDARY_TOLERANCE,
        };
        if thetas.is_empty()
            || thetas
                .iter()
                .any(|theta| !theta.is_finite() || *theta <= 0.0)
        {
            return ArbitrageAssessment::new(
                ArbitrageStatus::Indeterminate,
                evidence,
                f64::NAN,
                None,
            );
        }
        // Condition (i): theta_t non-decreasing in maturity order.
        let mut calendar_margin = f64::INFINITY;
        for pair in thetas.windows(2) {
            let difference = pair[1] - pair[0];
            calendar_margin = calendar_margin.min(difference);
            if difference < -BOUNDARY_TOLERANCE {
                return ArbitrageAssessment::new(
                    ArbitrageStatus::ViolationDetected,
                    evidence,
                    difference,
                    None,
                );
            }
            if difference < 0.0 {
                return ArbitrageAssessment::new(
                    ArbitrageStatus::Indeterminate,
                    evidence,
                    difference,
                    None,
                );
            }
        }
        // Condition (ii): the slope bound at every theta.
        if thetas
            .iter()
            .all(|&theta| self.calendar_condition_at(theta))
        {
            ArbitrageAssessment::new(
                ArbitrageStatus::NoViolationDetected,
                evidence,
                calendar_margin,
                None,
            )
        } else {
            ArbitrageAssessment::new(ArbitrageStatus::ViolationDetected, evidence, -0.0, None)
        }
    }
}

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

    #[test]
    fn phi_heston_validation() {
        assert!(Phi::heston(1.0).is_ok());
        assert!(Phi::heston(0.0).is_err());
        assert!(Phi::heston(-1.0).is_err());
        assert!(Phi::heston(f64::NAN).is_err());
    }

    #[test]
    fn heston_phi_small_theta_uses_cancellation_free_series() {
        let phi = Phi::heston(1.0).expect("valid test or documentation fixture");
        let value = phi.eval(1e-12);
        assert!(value.is_finite());
        assert!((value - 0.5).abs() < 1e-12);
    }

    #[test]
    fn total_variance_avoids_intermediate_square_overflow() {
        let ssvi = Ssvi::new(0.0, Phi::heston(1.0).expect("valid extreme-value fixture"))
            .expect("valid extreme-value fixture");
        let value = ssvi.total_variance(1e308, 0.04);
        assert!(value.is_finite());
        assert!((value - 9.867_989_404_040_118e305).abs() / value < 1e-14);
    }

    #[test]
    fn total_variance_preserves_near_monotone_wing_correction() {
        let ssvi = Ssvi::new(
            -1.0 + f64::EPSILON,
            Phi::heston(1.0).expect("valid near-wing fixture"),
        )
        .expect("valid near-wing fixture");
        let theta = 0.04;
        let k = 1e8 / ssvi.phi().eval(theta);
        let expected_bracket = 2.220_446_071_454_774e-8;
        let expected = (theta / 2.0) * expected_bracket;
        assert!((ssvi.total_variance(k, theta) - expected).abs() < 1e-23);
    }

    #[test]
    fn empty_butterfly_support_is_indeterminate() {
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test fixture"),
        )
        .expect("valid test fixture");
        assert_eq!(
            ssvi.butterfly_assessment(&[]).status(),
            ArbitrageStatus::Indeterminate
        );
    }

    #[test]
    fn boundary_near_calendar_decrease_is_indeterminate() {
        let ssvi = Ssvi::new(-0.3, Phi::heston(1.0).expect("valid test fixture"))
            .expect("valid test fixture");
        let theta = 0.04_f64;
        let smaller = f64::from_bits(theta.to_bits() - 1);
        let assessment = ssvi.calendar_assessment(&[theta, smaller]);
        assert_eq!(assessment.status(), ArbitrageStatus::Indeterminate);
        assert!(assessment.margin() < 0.0);
    }

    #[test]
    fn material_calendar_decrease_is_a_violation() {
        let ssvi = Ssvi::new(-0.3, Phi::heston(1.0).expect("valid test fixture"))
            .expect("valid test fixture");
        let assessment = ssvi.calendar_assessment(&[0.04, 0.039]);
        assert_eq!(assessment.status(), ArbitrageStatus::ViolationDetected);
    }

    #[test]
    fn phi_power_law_validation() {
        assert!(Phi::modified_power_law(0.5, 0.5).is_ok());
        assert!(Phi::modified_power_law(0.0, 0.5).is_err());
        assert!(Phi::modified_power_law(0.5, 0.0).is_err());
        assert!(Phi::modified_power_law(0.5, 1.0).is_err());
        assert!(Phi::modified_power_law(f64::INFINITY, 0.5).is_err());
    }

    #[test]
    fn phi_eval_positive_and_decreasing() {
        for phi in [
            Phi::heston(1.0).expect("valid test or documentation fixture"),
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        ] {
            let a = phi.eval(0.01);
            let b = phi.eval(0.04);
            let c = phi.eval(0.16);
            assert!(a > 0.0 && b > 0.0 && c > 0.0);
            assert!(a > b && b > c, "phi should be decreasing in theta");
        }
    }

    #[test]
    fn phi_power_law_golden() {
        // phi(1) = eta / (1^gamma * 2^(1-gamma)).
        let p = Phi::modified_power_law(1.0, 0.5).expect("valid test or documentation fixture");
        assert!((p.eval(1.0) - 1.0 / 2.0_f64.sqrt()).abs() < 1e-12);
    }

    #[test]
    fn ssvi_new_validation() {
        let phi = Phi::heston(1.0).expect("valid test or documentation fixture");
        assert!(Ssvi::new(-0.4, phi).is_ok());
        assert!(Ssvi::new(1.0, phi).is_err());
        assert!(Ssvi::new(f64::NAN, phi).is_err());
    }

    #[test]
    fn ssvi_total_variance_atm_zero_rho() {
        // rho = 0, k = 0: w = theta/2 * (1 + 1) = theta.
        let ssvi = Ssvi::new(
            0.0,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert!((ssvi.total_variance(0.0, 0.04) - 0.04).abs() < 1e-15);
    }

    #[test]
    fn slice_at_reproduces_total_variance() {
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        let raw = ssvi
            .slice_at(0.04)
            .expect("valid test or documentation fixture");
        for &k in &[-0.5, -0.1, 0.0, 0.1, 0.5] {
            let direct = ssvi.total_variance(k, 0.04);
            assert!((raw.total_variance(k) - direct).abs() < 1e-12, "k = {k}");
        }
    }

    #[test]
    fn slice_at_rejects_non_positive_theta() {
        let ssvi = Ssvi::new(
            -0.3,
            Phi::heston(1.0).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert!(matches!(
            ssvi.slice_at(0.0),
            Err(ParamError::NonPositiveTheta { .. })
        ));
    }

    #[test]
    fn slice_at_heston_reproduces_total_variance() {
        let ssvi = Ssvi::new(
            0.2,
            Phi::heston(2.0).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        let raw = ssvi
            .slice_at(0.09)
            .expect("valid test or documentation fixture");
        for &k in &[-0.4, 0.0, 0.4] {
            let direct = ssvi.total_variance(k, 0.09);
            assert!((raw.total_variance(k) - direct).abs() < 1e-12, "k = {k}");
        }
    }

    #[test]
    fn butterfly_free_holds_for_small_eta() {
        // Power-law corollary: eta*(1+|rho|) <= 2 implies butterfly-free.
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.butterfly_assessment(&[0.01, 0.04, 0.09, 0.25])
                .status(),
            ArbitrageStatus::NoViolationDetected
        );
    }

    #[test]
    fn butterfly_violation_for_large_phi() {
        // A power-law phi with large eta gives huge phi near theta -> 0,
        // which violates Theorem 4.2 (theta*phi^2*(1+|rho|) > 4).
        let ssvi = Ssvi::new(
            0.5,
            Phi::modified_power_law(20.0, 0.9).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.butterfly_assessment_at(1e-3).status(),
            ArbitrageStatus::Indeterminate
        );
    }

    #[test]
    fn theorem_42_binary64_boundary_is_indeterminate() {
        let theta = 4.0_f64;
        let gamma = 0.5_f64;
        let boundary_eta = theta.powf(gamma) * (1.0 + theta).powf(1.0 - gamma);
        let boundary = Ssvi::new(
            0.0,
            Phi::modified_power_law(boundary_eta, gamma).expect("positive boundary phi"),
        )
        .expect("valid SSVI");
        assert_eq!(
            boundary.butterfly_assessment_at(theta).status(),
            ArbitrageStatus::Indeterminate
        );

        let interior = Ssvi::new(
            0.0,
            Phi::modified_power_law(boundary_eta * (1.0 - 1024.0 * f64::EPSILON), gamma)
                .expect("strictly interior phi"),
        )
        .expect("valid SSVI");
        assert_eq!(
            interior.butterfly_assessment_at(theta).status(),
            ArbitrageStatus::NoViolationDetected
        );
    }

    #[test]
    fn global_butterfly_envelope_boundary_is_indeterminate() {
        let boundary = Ssvi::new(0.0, Phi::heston(0.25).expect("positive boundary lambda"))
            .expect("valid SSVI");
        assert_eq!(
            boundary.global_butterfly_assessment().status(),
            ArbitrageStatus::Indeterminate
        );
    }

    #[test]
    fn butterfly_free_rejects_bad_theta() {
        let ssvi = Ssvi::new(
            0.0,
            Phi::heston(1.0).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.butterfly_assessment_at(0.0).status(),
            ArbitrageStatus::Indeterminate
        );
        assert_eq!(
            ssvi.butterfly_assessment_at(-1.0).status(),
            ArbitrageStatus::Indeterminate
        );
    }

    #[test]
    fn calendar_free_for_monotone_thetas() {
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.calendar_assessment(&[0.01, 0.04, 0.09]).status(),
            ArbitrageStatus::NoViolationDetected
        );
    }

    #[test]
    fn calendar_arbitrage_for_decreasing_thetas() {
        let ssvi = Ssvi::new(
            -0.3,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.calendar_assessment(&[0.09, 0.04, 0.01]).status(),
            ArbitrageStatus::ViolationDetected
        );
    }

    #[test]
    fn calendar_free_at_zero_rho() {
        let ssvi = Ssvi::new(
            0.0,
            Phi::modified_power_law(0.5, 0.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        assert_eq!(
            ssvi.calendar_assessment(&[0.04]).status(),
            ArbitrageStatus::NoViolationDetected
        );
    }
}