regit-svi 2.0.0

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

//! Joint SSVI surface calibration (Gatheral & Jacquier 2014, Sections 4-5).
//!
//! SSVI calibration fits the whole surface at once and reports its analytic
//! sufficient butterfly and necessary-and-sufficient calendar evidence:
//!
//! 1. **ATM term structure.** Validate the caller-supplied non-decreasing
//!    `theta_t` curve.
//! 2. **Joint objective.** Use deterministic multi-start local minimization of
//!    the total weighted residual over `rho` and the `phi` parameters, with
//!    `w(k, theta)` from the SSVI form.
//! 3. **Constraints.** The Theorem 4.1 and 4.2 inequalities are included in
//!    the objective. The calibration-integrity layer performs final hard
//!    feasibility checks before a constrained success is emitted.
//!
//! The outer search is 2-D (Heston-like `phi`, parameters `rho, lambda`) or
//! 3-D (power-law `phi`, parameters `rho, eta, gamma`), solved with the
//! Nelder-Mead simplex.
//!
//! # References
//!
//! - Gatheral, J. & Jacquier, A., "Arbitrage-free SVI volatility surfaces",
//!   *Quantitative Finance* 14(1):59-71 (2014), Sections 4-5.

use core::cell::Cell;

use crate::calibration::config::{
    ConstraintMode, InitializationPolicy, SurfaceCalibrationConfig, validate_effective_quotes,
};
use crate::calibration::report::{
    ParameterizationEvidence, ResidualDiagnostics, SurfaceCalibrationReport,
    SurfaceParameterMargins, TerminationReason,
};
use crate::error::{CalibrationError, ParamError};
use crate::market::quote::Quote;
use crate::market::units::{Maturity, TotalVariance};
use crate::no_arb::evidence::ArbitrageStatus;
use crate::numerics::nelder_mead;
use crate::surface::ssvi::{Phi, Ssvi};

/// Which family of `phi` smoothing function to calibrate.
///
/// # Examples
///
/// ```
/// use regit_svi::calibration::surface::PhiFamily;
///
/// let f = PhiFamily::PowerLaw;
/// assert_eq!(f, PhiFamily::PowerLaw);
/// ```
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum PhiFamily {
    /// Heston-like `phi`, with the single parameter `lambda`.
    Heston,
    /// Power-law `phi`, with parameters `eta` and `gamma`.
    PowerLaw,
}

/// One calibration maturity: a time to expiry, its ATM total variance, and
/// the slice of quotes observed at that maturity.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::market::quote::Quote;
/// use regit_svi::calibration::surface::SsviMaturity;
///
/// let quotes = vec![Quote::new(0.0, 0.04, 1.0)?];
/// let mat = SsviMaturity::new(1.0, 0.04, quotes)?;
/// assert!((mat.theta().get() - 0.04).abs() < 1e-15);
/// # Ok(())
/// # }
/// ```
#[derive(Debug, Clone, PartialEq)]
pub struct SsviMaturity {
    maturity: Maturity,
    theta: TotalVariance,
    quotes: Vec<Quote>,
}

impl SsviMaturity {
    /// Creates one validated SSVI calibration maturity.
    ///
    /// # Errors
    ///
    /// Returns [`ParamError::EmptyCollection`] if the quote set is empty,
    /// [`ParamError::NonFinite`] or [`ParamError::NonPositiveMaturity`] for an
    /// invalid maturity, and [`ParamError::NonFinite`],
    /// [`ParamError::NegativeTotalVariance`], or
    /// [`ParamError::NonPositiveTheta`] for an invalid theta.
    pub fn new(t: f64, theta: f64, quotes: Vec<Quote>) -> Result<Self, ParamError> {
        if quotes.is_empty() {
            return Err(ParamError::EmptyCollection {
                name: "SSVI maturity quotes",
            });
        }
        let maturity = Maturity::new(t)?;
        let theta = TotalVariance::new(theta)?;
        if theta.get() <= 0.0 {
            return Err(ParamError::NonPositiveTheta { theta: theta.get() });
        }
        Ok(Self {
            maturity,
            theta,
            quotes,
        })
    }

    /// Time to expiry in years.
    #[must_use]
    pub const fn maturity(&self) -> Maturity {
        self.maturity
    }
    /// ATM total variance for this maturity.
    #[must_use]
    pub const fn theta(&self) -> TotalVariance {
        self.theta
    }
    /// Quotes used at this maturity.
    #[must_use]
    pub fn quotes(&self) -> &[Quote] {
        &self.quotes
    }
}

/// The result of an SSVI surface calibration.
#[derive(Debug, Clone, PartialEq)]
pub struct SsviCalibration {
    ssvi: Ssvi,
    thetas: Vec<TotalVariance>,
    report: SurfaceCalibrationReport,
}

impl SsviCalibration {
    /// Fitted SSVI model.
    #[must_use]
    pub const fn ssvi(&self) -> Ssvi {
        self.ssvi
    }
    /// Validated ATM variance knots in maturity order.
    #[must_use]
    pub fn thetas(&self) -> &[TotalVariance] {
        &self.thetas
    }
    /// Weighted root-mean-square residual recomputed from the final surface.
    #[must_use]
    pub const fn rmse(&self) -> f64 {
        self.report.residuals().rmse()
    }
    /// Termination, residual, and feasibility evidence.
    #[must_use]
    pub const fn report(&self) -> SurfaceCalibrationReport {
        self.report
    }
}

/// Calibrates an SSVI surface jointly across maturities.
///
/// The `theta_t` term structure is taken from the `theta` field of each
/// supplied maturity; the global parameters `(rho, phi-params)` are fitted by
/// a multi-started Nelder-Mead search. In constrained mode infeasible
/// candidates lie outside the objective domain and a final hard check guards
/// every successful return.
///
/// # Errors
///
/// - [`CalibrationError::EmptyQuotes`] if no maturities or no quotes are
///   supplied.
/// - [`CalibrationError::AllWeightsZero`] if every fitting weight is zero.
/// - [`CalibrationError::InsufficientEffectiveQuotes`] if the positive-weight
///   surface design has too few sufficiently distinct strikes.
/// - [`CalibrationError::InsufficientThetaLevels`] if a power-law fit has
///   fewer than two distinct positive ATM-variance levels.
/// - [`CalibrationError::Param`] if maturity/theta ordering, aggregate quote
///   weight, or fitted parameters are invalid or non-finite.
/// - [`CalibrationError::DidNotConverge`] if no feasible surface was found.
/// - [`CalibrationError::Infeasible`] if constrained mode cannot establish the
///   global butterfly envelope or calendar conditions.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::surface::ssvi::{Phi, Ssvi};
/// use regit_svi::market::quote::Quote;
/// use regit_svi::calibration::surface::{calibrate, PhiFamily, SsviMaturity};
///
/// // Generate a synthetic surface from a known SSVI and recover it.
/// let truth = Ssvi::new(-0.3, Phi::modified_power_law(0.5, 0.5)?)?;
/// let ks = [-0.3, -0.15, 0.0, 0.15, 0.3];
/// let mats: Vec<SsviMaturity> = [(0.5, 0.02), (1.0, 0.04), (2.0, 0.07)]
///     .iter()
///     .map(|&(t, theta)| {
///         let quotes = ks
///             .iter()
///             .map(|&k| Quote::new(k, truth.total_variance(k, theta), 1.0))
///             .collect::<Result<Vec<_>, _>>()?;
///         SsviMaturity::new(t, theta, quotes)
///     })
///     .collect::<Result<_, _>>()?;
/// let fit = calibrate(&mats, PhiFamily::PowerLaw)?;
/// assert!(fit.rmse() < 1e-3);
/// # Ok(())
/// # }
/// ```
pub fn calibrate(
    maturities: &[SsviMaturity],
    family: PhiFamily,
) -> Result<SsviCalibration, CalibrationError> {
    calibrate_with_config(maturities, family, SurfaceCalibrationConfig::default())
}

/// Calibrates an SSVI surface with explicit limits and feasibility policy.
///
/// # Errors
///
/// Returns the same variants as [`calibrate`]:
/// [`CalibrationError::EmptyQuotes`], [`CalibrationError::AllWeightsZero`],
/// [`CalibrationError::InsufficientEffectiveQuotes`],
/// [`CalibrationError::InsufficientThetaLevels`], [`CalibrationError::Param`],
/// [`CalibrationError::DidNotConverge`], and
/// [`CalibrationError::Infeasible`].
#[allow(clippy::too_many_lines)] // Keeps validation, hard-feasible search, and final audit in source order.
pub fn calibrate_with_config(
    maturities: &[SsviMaturity],
    family: PhiFamily,
    config: SurfaceCalibrationConfig,
) -> Result<SsviCalibration, CalibrationError> {
    if maturities.is_empty() || maturities.iter().all(|m| m.quotes.is_empty()) {
        return Err(CalibrationError::EmptyQuotes);
    }
    for (index, pair) in maturities.windows(2).enumerate() {
        if pair[1].maturity <= pair[0].maturity {
            return Err(CalibrationError::Param(ParamError::NotStrictlyIncreasing {
                name: "maturity",
                index: index + 1,
                previous: pair[0].maturity.get(),
                value: pair[1].maturity.get(),
            }));
        }
        if pair[1].theta < pair[0].theta {
            return Err(CalibrationError::Param(ParamError::DecreasingAtmVariance {
                index: index + 1,
                previous: pair[0].theta.get(),
                value: pair[1].theta.get(),
            }));
        }
    }
    let all_quotes: Vec<Quote> = maturities
        .iter()
        .flat_map(|maturity| maturity.quotes.iter().copied())
        .collect();
    let minimum = match family {
        PhiFamily::Heston => 2,
        PhiFamily::PowerLaw => 3,
    };
    let (usable, distinct) = validate_effective_quotes(&all_quotes, minimum, 1e-12)?;
    let total_weight: f64 = maturities
        .iter()
        .flat_map(|m| m.quotes.iter())
        .map(|q| q.weight)
        .sum();
    if total_weight <= 0.0 {
        return Err(CalibrationError::AllWeightsZero);
    }

    let thetas: Vec<f64> = maturities.iter().map(|m| m.theta.get()).collect();
    if family == PhiFamily::PowerLaw {
        let mut distinct_theta = 1_usize;
        let mut previous = thetas[0];
        for theta in thetas.iter().copied().skip(1) {
            if theta - previous > 64.0 * f64::EPSILON * (1.0 + theta.max(previous)) {
                distinct_theta = distinct_theta.saturating_add(1);
                previous = theta;
            }
        }
        if distinct_theta < 2 {
            return Err(CalibrationError::InsufficientThetaLevels {
                got: distinct_theta,
                need: 2,
            });
        }
    }

    let evaluations = Cell::new(0_usize);
    // In constrained mode infeasible points are outside the objective domain.
    let objective = |p: &[f64]| -> f64 {
        evaluations.set(evaluations.get().saturating_add(1));
        let Some(ssvi) = surface_from_params(p, family) else {
            return f64::INFINITY;
        };
        let mut cost = 0.0;
        for mat in maturities {
            for q in &mat.quotes {
                if q.weight <= 0.0 {
                    continue;
                }
                let model = ssvi.total_variance(q.k, mat.theta.get());
                let r = model - q.w;
                cost += q.weight * r * r;
            }
        }
        if config.constraint_mode() == ConstraintMode::Constrained
            && (ssvi.global_butterfly_assessment().status() != ArbitrageStatus::NoViolationDetected
                || ssvi.calendar_assessment(&thetas).status()
                    != ArbitrageStatus::NoViolationDetected)
        {
            return f64::INFINITY;
        }
        cost
    };

    // Multi-start seeds; rho coordinates are unconstrained via tanh.
    let rho_seeds = [-0.6_f64, -0.2, 0.0, 0.2, 0.6];
    let mut best_obj = f64::INFINITY;
    let mut best_params: Vec<f64> = Vec::new();
    let mut selected_start = None;
    let mut selected_termination = crate::numerics::OptimizerTermination::IterationLimit;
    let mut total_iterations = 0_usize;
    let mut start_index = 0_usize;
    let rho_count = if config.initialization() == InitializationPolicy::DeterministicSingleStart {
        1
    } else {
        rho_seeds.len()
    };

    match family {
        PhiFamily::Heston => {
            for &rho0 in &rho_seeds[..rho_count] {
                for &lambda0 in &[0.5_f64, 1.0, 2.0, 5.0] {
                    let correlation = rho0.abs();
                    let floor = (1.0 + correlation) / 4.0;
                    let start = [atanh_clamped(rho0), (lambda0 - floor).max(0.05).ln()];
                    let res = nelder_mead(
                        objective,
                        &start,
                        config.tolerance(),
                        config.outer_iterations(),
                    );
                    total_iterations = total_iterations.saturating_add(res.iterations);
                    let eligible =
                        res.converged || config.constraint_mode() == ConstraintMode::BestEffort;
                    if eligible && res.fx < best_obj {
                        best_obj = res.fx;
                        best_params = res.x;
                        selected_start = Some(start_index);
                        selected_termination = res.termination;
                    }
                    start_index = start_index.saturating_add(1);
                }
            }
        }
        PhiFamily::PowerLaw => {
            for &rho0 in &rho_seeds[..rho_count] {
                for &eta0 in &[0.3_f64, 0.6, 1.0] {
                    for &gamma0 in &[0.15_f64, 0.3, 0.45] {
                        let eta_max = power_eta_max(rho0, gamma0);
                        let eta_fraction = (eta0 / eta_max).clamp(1e-6, 1.0 - 1e-6);
                        let start = [
                            atanh_clamped(rho0),
                            logit(eta_fraction),
                            logit(2.0 * gamma0),
                        ];
                        let res = nelder_mead(
                            objective,
                            &start,
                            config.tolerance(),
                            config.outer_iterations(),
                        );
                        total_iterations = total_iterations.saturating_add(res.iterations);
                        let eligible =
                            res.converged || config.constraint_mode() == ConstraintMode::BestEffort;
                        if eligible && res.fx < best_obj {
                            best_obj = res.fx;
                            best_params = res.x;
                            selected_start = Some(start_index);
                            selected_termination = res.termination;
                        }
                        start_index = start_index.saturating_add(1);
                    }
                }
            }
        }
    }

    let ssvi =
        surface_from_params(&best_params, family).ok_or(CalibrationError::DidNotConverge {
            iterations: total_iterations,
            residual: best_obj,
        })?;

    let butterfly = ssvi.global_butterfly_assessment();
    let calendar = ssvi.calendar_assessment(&thetas);
    if config.constraint_mode() == ConstraintMode::Constrained {
        if butterfly.status() != ArbitrageStatus::NoViolationDetected {
            return Err(CalibrationError::Infeasible {
                condition: "global SSVI butterfly envelope",
                margin: butterfly.margin(),
            });
        }
        if calendar.status() != ArbitrageStatus::NoViolationDetected {
            return Err(CalibrationError::Infeasible {
                condition: "SSVI calendar conditions",
                margin: calendar.margin(),
            });
        }
    }

    // Residual diagnostics are recomputed after every parameter transform.
    let mut residual = 0.0;
    let mut max_absolute = 0.0_f64;
    for mat in maturities {
        for q in &mat.quotes {
            if q.weight <= 0.0 {
                continue;
            }
            let r = ssvi.total_variance(q.k, mat.theta.get()) - q.w;
            residual += q.weight * r * r;
            max_absolute = max_absolute.max(r.abs());
        }
    }
    let rmse = (residual / total_weight).sqrt();
    let residuals = ResidualDiagnostics::new(
        residual,
        rmse,
        max_absolute,
        usable,
        distinct,
        all_quotes.len().saturating_sub(usable),
        total_weight,
    );
    let termination = TerminationReason::from(selected_termination);
    let (parameterization, phi_scale, phi_shape) = match family {
        PhiFamily::Heston => (
            ParameterizationEvidence::new(
                "rho=(1-margin)*tanh; lambda=(1+|rho|)/4+exp",
                "no projection; infeasible candidates excluded from the objective domain",
            ),
            ssvi.phi().heston_lambda().unwrap_or(f64::NAN) - (1.0 + ssvi.rho().abs()) / 4.0,
            f64::INFINITY,
        ),
        PhiFamily::PowerLaw => {
            let (eta, gamma) = ssvi
                .phi()
                .modified_power_law_parameters()
                .unwrap_or((f64::NAN, f64::NAN));
            (
                ParameterizationEvidence::new(
                    "rho=(1-margin)*tanh; gamma=0.5*logistic; eta=envelope*logistic",
                    "no projection; infeasible candidates excluded from the objective domain",
                ),
                eta.min(power_eta_max(ssvi.rho(), gamma) - eta),
                gamma.min(0.5 - gamma),
            )
        }
    };
    let report = SurfaceCalibrationReport::new(
        total_iterations,
        evaluations.get(),
        termination,
        residuals,
        butterfly,
        calendar,
        "hard-feasible SSVI / Nelder–Mead",
        match family {
            PhiFamily::Heston => rho_count * 4,
            PhiFamily::PowerLaw => rho_count * 9,
        },
        selected_start,
        config.tolerance(),
        parameterization,
        SurfaceParameterMargins::new(
            1.0 - ssvi.rho().abs(),
            phi_scale,
            phi_shape,
            butterfly.margin(),
            calendar.margin(),
        ),
    );

    Ok(SsviCalibration {
        ssvi,
        thetas: maturities.iter().map(|maturity| maturity.theta).collect(),
        report,
    })
}

/// Reconstructs an [`Ssvi`] from the unconstrained outer-search vector.
///
/// `rho` is mapped through `tanh`; the `phi` parameters through `exp`
/// (positivity) or the logistic (open interval `(0, 1)` for `gamma`).
fn surface_from_params(p: &[f64], family: PhiFamily) -> Option<Ssvi> {
    const RHO_MARGIN: f64 = 32.0 * f64::EPSILON;
    let rho = (1.0 - RHO_MARGIN) * p.first()?.tanh();
    let phi = match family {
        PhiFamily::Heston => {
            let lambda = (1.0 + rho.abs()) / 4.0 + p.get(1)?.exp();
            Phi::heston(lambda).ok()?
        }
        PhiFamily::PowerLaw => {
            let gamma = 0.5 * logistic(*p.get(2)?);
            let eta = power_eta_max(rho, gamma) * logistic(*p.get(1)?);
            Phi::modified_power_law(eta, gamma).ok()?
        }
    };
    Ssvi::new(rho, phi).ok()
}

fn power_eta_max(rho: f64, gamma: f64) -> f64 {
    let correlation_factor = 1.0 + rho.abs();
    let x = 1.0 - 2.0 * gamma;
    let h_gamma = if x.abs() <= 64.0 * f64::EPSILON {
        1.0
    } else {
        x.powf(x) / (2.0 - 2.0 * gamma).powf(2.0 - 2.0 * gamma)
    };
    (4.0 / correlation_factor).min(2.0 / (correlation_factor * h_gamma).sqrt())
}

/// Inverse hyperbolic tangent with the argument clamped just inside `(-1, 1)`.
#[inline]
fn atanh_clamped(x: f64) -> f64 {
    let x = x.clamp(-0.999_999, 0.999_999);
    0.5 * ((1.0 + x) / (1.0 - x)).ln()
}

/// Logistic function mapping `R -> (0, 1)`, used for the `gamma` parameter.
#[inline]
fn logistic(x: f64) -> f64 {
    1.0 / (1.0 + (-x).exp())
}

/// Inverse logistic (logit) mapping `(0, 1) -> R`.
#[inline]
fn logit(p: f64) -> f64 {
    let p = p.clamp(1e-6, 1.0 - 1e-6);
    (p / (1.0 - p)).ln()
}

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

    /// Builds synthetic maturities from a known SSVI surface.
    fn synthetic(truth: &Ssvi, ts_thetas: &[(f64, f64)], ks: &[f64]) -> Vec<SsviMaturity> {
        ts_thetas
            .iter()
            .map(|&(t, theta)| {
                SsviMaturity::new(
                    t,
                    theta,
                    ks.iter()
                        .map(|&k| {
                            Quote::new(k, truth.total_variance(k, theta), 1.0)
                                .expect("valid test or documentation fixture")
                        })
                        .collect(),
                )
                .expect("valid test or documentation fixture")
            })
            .collect()
    }

    #[test]
    fn logistic_logit_invert() {
        for &x in &[0.05, 0.3, 0.5, 0.8, 0.95] {
            assert!((logistic(logit(x)) - x).abs() < 1e-10);
        }
    }

    #[test]
    fn rejects_empty() {
        assert!(matches!(
            calibrate(&[], PhiFamily::PowerLaw),
            Err(CalibrationError::EmptyQuotes)
        ));
    }

    #[test]
    fn rejects_unordered_or_decreasing_maturities() {
        let quotes = vec![
            Quote::new(-0.1, 0.04, 1.0).expect("valid test or documentation fixture"),
            Quote::new(0.0, 0.04, 1.0).expect("valid test or documentation fixture"),
            Quote::new(0.1, 0.04, 1.0).expect("valid test or documentation fixture"),
        ];
        let late = SsviMaturity::new(2.0, 0.06, quotes.clone())
            .expect("valid test or documentation fixture");
        let early = SsviMaturity::new(1.0, 0.04, quotes.clone())
            .expect("valid test or documentation fixture");
        assert!(matches!(
            calibrate(&[late, early], PhiFamily::PowerLaw),
            Err(CalibrationError::Param(
                ParamError::NotStrictlyIncreasing { .. }
            ))
        ));
        let high = SsviMaturity::new(1.0, 0.06, quotes.clone())
            .expect("valid test or documentation fixture");
        let low =
            SsviMaturity::new(2.0, 0.04, quotes).expect("valid test or documentation fixture");
        assert!(matches!(
            calibrate(&[high, low], PhiFamily::PowerLaw),
            Err(CalibrationError::Param(
                ParamError::DecreasingAtmVariance { .. }
            ))
        ));
    }

    #[test]
    fn power_law_requires_two_distinct_theta_levels() {
        let quotes = vec![
            Quote::new(-0.1, 0.04, 1.0).expect("valid fixture"),
            Quote::new(0.0, 0.04, 1.0).expect("valid fixture"),
            Quote::new(0.1, 0.04, 1.0).expect("valid fixture"),
        ];
        let one = SsviMaturity::new(1.0, 0.04, quotes.clone()).expect("valid maturity");
        assert!(matches!(
            calibrate(core::slice::from_ref(&one), PhiFamily::PowerLaw),
            Err(CalibrationError::InsufficientThetaLevels { got: 1, need: 2 })
        ));
        let same_theta_later = SsviMaturity::new(2.0, 0.04, quotes).expect("valid maturity");
        assert!(matches!(
            calibrate(&[one, same_theta_later], PhiFamily::PowerLaw),
            Err(CalibrationError::InsufficientThetaLevels { got: 1, need: 2 })
        ));
    }

    #[test]
    fn recovers_power_law_surface() {
        let truth = 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 ks = [-0.3, -0.15, 0.0, 0.15, 0.3];
        let mats = synthetic(&truth, &[(0.5, 0.02), (1.0, 0.04), (2.0, 0.07)], &ks);
        let fit =
            calibrate(&mats, PhiFamily::PowerLaw).expect("valid test or documentation fixture");
        assert!(fit.rmse() < 1e-3, "rmse = {}", fit.rmse());
        assert!(fit.report().selected_start().is_some());
        let margins = fit.report().margins();
        let (_, gamma) = fit
            .ssvi()
            .phi()
            .modified_power_law_parameters()
            .expect("power-law fit");
        assert!((margins.phi_shape() - gamma.min(0.5 - gamma)).abs() < 1e-15);
        assert!(margins.phi_scale() >= 0.0);
        assert!(
            fit.report()
                .parameterization()
                .repair()
                .contains("no projection")
        );
        assert!(
            (fit.ssvi().rho() - (-0.3)).abs() < 0.1,
            "rho = {}",
            fit.ssvi().rho()
        );
    }

    #[test]
    fn recovers_heston_surface() {
        let truth = Ssvi::new(
            -0.2,
            Phi::heston(1.5).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        let ks = [-0.3, -0.15, 0.0, 0.15, 0.3];
        let mats = synthetic(&truth, &[(0.5, 0.02), (1.0, 0.04), (2.0, 0.07)], &ks);
        let fit = calibrate(&mats, PhiFamily::Heston).expect("valid test or documentation fixture");
        assert!(fit.rmse() < 5e-3, "rmse = {}", fit.rmse());
    }

    #[test]
    fn fitted_surface_is_arbitrage_free() {
        let truth = Ssvi::new(
            -0.4,
            Phi::modified_power_law(0.6, 0.4).expect("valid test or documentation fixture"),
        )
        .expect("valid test or documentation fixture");
        let ks = [-0.4, -0.2, 0.0, 0.2, 0.4];
        let mats = synthetic(&truth, &[(0.25, 0.015), (1.0, 0.05), (3.0, 0.11)], &ks);
        let fit =
            calibrate(&mats, PhiFamily::PowerLaw).expect("valid test or documentation fixture");
        assert_eq!(
            fit.report().butterfly_assessment().status(),
            ArbitrageStatus::NoViolationDetected
        );
        assert_eq!(
            fit.report().calendar_assessment().status(),
            ArbitrageStatus::NoViolationDetected
        );
    }
}