Quantitative finance

// docs: quant#black-scholes-merton-closed-form-european-call
//! Backs the BSM example on the quant catalog page.

use stochastic_rs::quant::pricing::bsm::BSMCoc;
use stochastic_rs::quant::pricing::bsm::BSMPricer;
use stochastic_rs::quant::types::OptionType;
use stochastic_rs::traits::ModelPricer;

#[test]
fn bsm_price_and_greeks() {
  // Model state only: volatility + cost-of-carry convention.
  let model = BSMPricer::new(/* v */ 0.2, BSMCoc::Bsm1973);
  // Query: (s, k, r, q, tau).
  let (s, k, r, q, tau) = (100.0, 100.0, 0.05, 0.0, 1.0);

  let call = model.price_call(s, k, r, q, tau);
  let put = model.price_put(s, k, r, q, tau);
  let delta = model.delta(s, k, r, q, tau, OptionType::Call);
  let vega = model.vega(s, k, r, q, tau);

  assert!(call > 0.0 && put > 0.0);
  assert!((0.0..=1.0).contains(&delta));
  assert!(vega > 0.0);
}

import stochastic_rs as srs

# The Python wrapper still bundles the query into the constructor — the
# option type included — so the accessors below take no arguments.
pricer = srs.BSMPricer(s=100, v=0.2, k=100, r=0.05, tau=1.0, option_type="call", q=0.0)
call, put = pricer.call_put()
print("call =", call)
print("put  =", put)

print(f"delta={pricer.delta():.4f}, gamma={pricer.gamma():.4f}, "
      f"vega={pricer.vega():.4f}, theta={pricer.theta():.4f}, "
      f"rho={pricer.rho():.4f}")

// docs: quant#quanto-option-fixed-exchange-rate-foreign-equity
//! Backs the quanto example on the quant catalog page.

use stochastic_rs::quant::pricing::bsm::BSMCoc;
use stochastic_rs::quant::pricing::bsm::BSMPricer;
use stochastic_rs::quant::pricing::quanto::QuantoPricer;
use stochastic_rs::traits::ModelPricer;

#[test]
fn quanto_price_forward_and_merton_reduction() {
  // Model state: asset vol, FX vol, correlation, foreign rate, fixed rate.
  let model = QuantoPricer::new(0.2, 0.12, 0.3, 0.05, 1.5);
  // Query: (s, k, r = domestic rate, q, tau) — the Haug §2.13.4 inputs.
  let (s, k, r, q, tau) = (100.0, 105.0, 0.08, 0.04, 0.5);

  let (call, put) = model.call_put(s, k, r, q, tau);
  assert!((call - 5.2936847941).abs() < 1e-4);
  assert!((put - 12.2976985036).abs() < 1e-4);
  assert_eq!(model.price_call(s, k, r, q, tau), call);

  // The quanto forward carries the adjusted drift r_f − q − ρ σ_S σ_E.
  assert!((model.forward(s, q, tau) - 150.2101470686).abs() < 1e-9);

  // Without correlation and with equal rates it is E_p times Merton (1973).
  let plain = QuantoPricer::new(0.2, 0.12, 0.0, r, 1.5);
  let merton = BSMPricer::new(0.2, BSMCoc::Merton1973);
  assert!(
    (plain.price_call(s, k, r, q, tau) - 1.5 * merton.price_call(s, k, r, q, tau)).abs() < 1e-12
  );
}

import stochastic_rs as srs

# Foreign equity option paid in domestic currency at the fixed rate 1.5;
# r is the domestic rate, r_f the foreign one, rho = corr(asset, FX).
pricer = srs.QuantoPricer(s=100, v=0.2, k=105, r=0.08, tau=0.5, r_f=0.05, v_fx=0.12, rho=0.3, fixed_rate=1.5, q=0.04)
call, put = pricer.call_put()
print(call, put)          # 5.2937, 12.2977 on the Reiner / Haug §2.13.4 inputs
print(pricer.forward())   # 150.21 = E_p · S · exp((r_f − q − ρ σ_S σ_E) τ)

// docs: quant#dual-curve-bootstrap-ois-discounting-tenor-forecasting
//! Backs the dual-curve bootstrap example on the quant catalog page.

use ndarray::Array1;
use stochastic_rs::quant::curves::DiscountCurve;
use stochastic_rs::quant::curves::InterpolationMethod;
use stochastic_rs::quant::curves::MultiCurve;
use stochastic_rs::quant::curves::dual_curve::ForecastInstrument;
use stochastic_rs::quant::curves::dual_curve::bootstrap_forecast;

#[test]
fn dual_curve_bootstrap_recovers_the_tenor_basis() {
  // Exogenous OIS discount curve (flat 2 %) and a 3M tenor trading 50 bp above it.
  let ois = DiscountCurve::from_zero_rates(
    &Array1::from_vec(vec![0.25, 1.0, 5.0, 10.0]),
    &Array1::from_vec(vec![0.02; 4]),
    InterpolationMethod::LogLinearOnDiscountFactors,
  );
  let tenor_df = |t: f64| (-0.025_f64 * t).exp();

  // Market quotes implied by that tenor curve: 3M deposit, 3×6 FRA, par swaps
  // against 3M with annual fixed and quarterly floating payments.
  let mut quotes = vec![
    ForecastInstrument::Deposit {
      maturity: 0.25,
      rate: (1.0 / tenor_df(0.25) - 1.0) / 0.25,
    },
    ForecastInstrument::Fra {
      start: 0.25,
      end: 0.5,
      rate: (tenor_df(0.25) / tenor_df(0.5) - 1.0) / 0.25,
    },
  ];
  for years in [1_usize, 2, 3, 5] {
    let fixed_times: Vec<f64> = (1..=years).map(|i| i as f64).collect();
    let float_times: Vec<f64> = (1..=4 * years).map(|i| 0.25 * i as f64).collect();
    let float_pv: f64 = float_times
      .iter()
      .scan(0.0, |prev, &t| {
        let leg = ois.discount_factor(t) * (tenor_df(*prev) / tenor_df(t) - 1.0);
        *prev = t;
        Some(leg)
      })
      .sum();
    let annuity: f64 = fixed_times
      .iter()
      .scan(0.0, |prev, &t| {
        let leg = (t - *prev) * ois.discount_factor(t);
        *prev = t;
        Some(leg)
      })
      .sum();
    quotes.push(ForecastInstrument::Swap {
      rate: float_pv / annuity,
      fixed_times,
      float_times,
    });
  }

  // Forecast curve against OIS discounting: pseudo-discount factors come back exactly.
  let forecast = bootstrap_forecast(
    &quotes,
    &ois,
    InterpolationMethod::LogLinearOnDiscountFactors,
  );
  for t in [0.25, 0.5, 1.0, 2.0, 3.0, 5.0] {
    assert!((forecast.discount_factor(t) - tenor_df(t)).abs() < 1e-9);
  }

  // The multi-curve container reads the 50 bp tenor basis back off the two curves.
  let mut multi = MultiCurve::new(ois);
  multi.add_forecast("3M", forecast);
  let basis = multi.basis_spread("3M", 1.0, 1.25).unwrap();
  assert!((basis - 0.005).abs() < 2e-4, "basis {basis}");
}

import numpy as np
import stochastic_rs as srs

# OIS discount curve from par OIS quotes (one fixed-leg schedule per quote) …
ois = srs.DiscountCurve.bootstrap_ois([[1.0], [1.0, 2.0], [1.0, 2.0, 3.0]], [0.0202, 0.0204, 0.0206], interp="log_df")
# … and a 3M forecast curve bootstrapped against it from tenor quotes.
forecast = srs.DiscountCurve.bootstrap_forecast(
    ois,
    deposits=[(0.25, 0.0251)],
    fras=[(0.25, 0.5, 0.0252)],
    swaps=[(0.0254, [1.0], [0.25, 0.5, 0.75, 1.0]), (0.0256, [1.0, 2.0], [0.25 * i for i in range(1, 9)])],
    interp="log_df",
)
multi = srs.MultiCurve(ois)
multi.add_forecast("3M", forecast)
print(multi.basis_spread("3M", 1.0, 1.25))               # ≈ 0.005: the tenor trades ~50 bp over OIS
print(multi.fair_swap_rate("3M", [0.0, 0.5, 1.0, 1.5, 2.0]))

// docs: quant#callable-and-puttable-bonds-on-the-hullwhite-tree
//! Backs the callable bond example on the quant catalog page.

use stochastic_rs::quant::lattice::CallableBondSpec;
use stochastic_rs::quant::lattice::HullWhiteTree;
use stochastic_rs::quant::lattice::HullWhiteTreeModel;
use stochastic_rs::quant::lattice::price_callable_bond;

#[test]
fn callable_and_puttable_bonds_bracket_the_straight_bond() {
  // Hull–White tree: r₀ = 4 %, a = 0.3, θ = 4 %, σ = 1 %; monthly steps to the 3-year maturity.
  let tree = HullWhiteTree::new(HullWhiteTreeModel::new(0.04, 0.3, 0.04, 0.01), 3.0, 36);
  // 6 % annual coupon, face 100, callable and puttable at par after years 1 and 2.
  let bond = CallableBondSpec::new(100.0, 0.06, vec![1.0, 2.0, 3.0]);
  let callable = bond.clone().with_calls(vec![(1.0, 100.0), (2.0, 100.0)]);
  let puttable = bond.clone().with_puts(vec![(1.0, 100.0), (2.0, 100.0)]);

  let straight = price_callable_bond(&tree.tree, &tree.model, &bond);
  let called = price_callable_bond(&tree.tree, &tree.model, &callable);
  let put = price_callable_bond(&tree.tree, &tree.model, &puttable);

  // The issuer's call lowers the price, the holder's put raises it.
  assert!(called.price < straight.price && put.price > straight.price);
  assert!(called.call_value > 0.0 && put.put_value > 0.0);
  assert_eq!(straight.price, straight.straight_price);
}

import stochastic_rs as srs

tree = srs.HullWhiteCallableBond(initial_rate=0.04, mean_reversion=0.3, theta=0.04, sigma=0.01, horizon=3.0, steps=36)
price, straight, call_value, put_value = tree.price(
    face=100.0, coupon_rate=0.06, coupon_times=[1.0, 2.0, 3.0],
    calls=[(1.0, 100.0), (2.0, 100.0)],
)
print(price, straight, call_value)   # callable < straight; call_value = straight − callable

// docs: quant#bachelier-normal-model-and-implied-normal-volatility
//! Backs the Bachelier example on the quant catalog page.

use stochastic_rs::quant::pricing::bachelier::BachelierPricer;
use stochastic_rs::quant::pricing::bachelier::normal_implied_volatility;
use stochastic_rs::quant::types::OptionType;
use stochastic_rs::traits::ModelPricer;

#[test]
fn bachelier_price_and_normal_implied_volatility() {
  // Model state: normal volatility of 20 price units per √year.
  let model = BachelierPricer::new(20.0);
  // Query: (s, k, r, q, tau) — a strike below the forward.
  let (s, k, r, q, tau) = (100.0, 95.0, 0.05, 0.02, 0.75);

  let call = model.price_call(s, k, r, q, tau);
  let put = model.price_put(s, k, r, q, tau);
  let forward = model.forward(s, r, q, tau);
  assert!((call - put - (-r * tau).exp() * (forward - k)).abs() < 1e-9);

  // The inversion recovers the normal volatility from the price.
  let implied = model.implied_volatility(call, s, k, r, q, tau, OptionType::Call);
  assert!((implied - 20.0).abs() < 1e-9);

  // Forward-space entry point for rate options quoted in normal vol (bp).
  let undiscounted = call * (r * tau).exp();
  assert!(
    (normal_implied_volatility(undiscounted, forward, k, tau, OptionType::Call) - 20.0).abs()
      < 1e-9
  );
}

import stochastic_rs as srs

pricer = srs.BachelierPricer(s=100, v=20.0, k=95, r=0.05, tau=0.75, q=0.02)
call, put = pricer.call_put()
print(call, put, pricer.forward())
print(pricer.implied_volatility(call, "call"))   # 20.0: the normal vol round-trips

// docs: quant#swaption-calibration-of-the-tree-short-rate-models
//! Backs the tree swaption calibration example on the quant catalog page.

use ndarray::Array1;
use stochastic_rs::quant::calibration::hw_swaption::SwaptionQuote;
use stochastic_rs::quant::calibration::tree_swaption::BlackKarasinskiSwaptionCalibrator;
use stochastic_rs::quant::curves::DiscountCurve;
use stochastic_rs::quant::curves::InterpolationMethod;
use stochastic_rs::quant::instruments::option::types::SwaptionDirection;
use stochastic_rs::traits::Calibrator;
use stochastic_rs::traits::ToShortRateModel;

#[test]
fn black_karasinski_fits_a_small_swaption_grid() {
  // Flat 3 % curve and two ATM payer swaptions quoted in Black-76 vol.
  let curve = DiscountCurve::from_zero_rates(
    &Array1::from_vec(vec![0.5, 1.0, 5.0, 10.0]),
    &Array1::from_vec(vec![0.03; 4]),
    InterpolationMethod::LogLinearOnDiscountFactors,
  );
  let quote = |expiry: f64, tenor: f64, black_vol: f64| SwaptionQuote {
    expiry,
    tenor,
    black_vol,
    fixed_accrual: 0.5,
    direction: SwaptionDirection::Payer,
    weight: None,
  };
  let quotes = [quote(1.0, 2.0, 0.22), quote(2.0, 2.0, 0.20)];

  // (a, σ) of the log-rate, repriced on the Black–Karasinski tree at 8 levels per year.
  let calibrator =
    BlackKarasinskiSwaptionCalibrator::new(&quotes, &curve, 1.0, 0.03, 0.03, 8).with_max_iters(200);
  let result = calibrator.calibrate(Some((0.1, 0.2))).unwrap();
  let scale = result.market_prices.iter().sum::<f64>() / 2.0;
  assert!(
    result.rmse / scale < 0.05,
    "relative rmse {}",
    result.rmse / scale
  );

  // The result plugs straight into the lattice pipeline through `ToShortRateModel`.
  let model = ToShortRateModel::to_short_rate_model(&result, 0.03, 0.03);
  assert_eq!(model.sigma, result.sigma);
}

import numpy as np
import stochastic_rs as srs

curve = srs.DiscountCurve.from_zero_rates(np.array([0.5, 1.0, 5.0, 10.0]), np.full(4, 0.03), interp="log_df")
quotes = [(1.0, 2.0, 0.22, 0.5, "payer"), (2.0, 2.0, 0.20, 0.5, "payer")]   # (expiry, tenor, black_vol, accrual, direction)

bk = srs.BlackKarasinskiSwaptionCalibrator(quotes, curve, initial_rate=0.03, long_run_rate=0.03, steps_per_year=8)
print(bk.calibrate(initial_guess=(0.1, 0.2)))          # (a, sigma, rmse, converged)
hw = srs.HullWhiteSwaptionCalibrator(quotes, curve)
print(hw.calibrate())                                    # (a, sigma, rmse, converged) via Jamshidian
g2 = srs.G2ppSwaptionCalibrator(quotes, curve, initial_rate=0.03, steps_per_year=4, max_iters=200)
print(g2.calibrate())                                    # (a, b, sigma, eta, rho, rmse, converged)

// docs: quant#essvi-slices-with-arbitrage-free-interpolation
//! Backs the eSSVI example on the quant catalog page.

use stochastic_rs::quant::vol_surface::SsviParams;
use stochastic_rs::quant::vol_surface::essvi::calibrate_essvi;
use stochastic_rs::quant::vol_surface::ssvi::SsviSlice;

#[test]
fn essvi_slices_recover_a_global_ssvi_surface_and_interpolate_without_arbitrage() {
  // Market slices generated by a global SSVI surface (ρ = −0.4, η = 0.6, γ = 0.4) at four maturities.
  let params = SsviParams::new(-0.4, 0.6, 0.4);
  let maturities = [0.25, 0.5, 1.0, 2.0];
  let ks: Vec<f64> = (0..21).map(|i| -0.5 + 0.05 * i as f64).collect();
  let slices: Vec<SsviSlice<f64>> = maturities
    .iter()
    .map(|&t| SsviSlice {
      log_moneyness: ks.clone(),
      total_variance: ks
        .iter()
        .map(|&k| params.total_variance(k, 0.04 * t))
        .collect(),
      theta: 0.04 * t,
    })
    .collect();

  // Anchored (ρ, ψ) per slice, going forward in maturity within the butterfly and calendar bounds.
  let surface = calibrate_essvi(&slices, &maturities);
  assert!(surface.is_butterfly_free() && surface.is_calendar_spread_free());
  for (slice, &t) in surface.slices.iter().zip(&maturities) {
    assert!((slice.rho + 0.4).abs() < 2e-3);
    assert!((slice.total_variance(0.2) - params.total_variance(0.2, 0.04 * t)).abs() < 2e-6);
  }

  // Between the slices the parameters interpolate linearly, keeping w non-decreasing in t.
  let w_early = surface.total_variance(0.1, 0.75);
  let w_late = surface.total_variance(0.1, 1.5);
  assert!(w_early > surface.total_variance(0.1, 0.5) && w_late > w_early);
}

import numpy as np
import stochastic_rs as srs

ks = np.linspace(-0.5, 0.5, 21)
maturities = [0.25, 0.5, 1.0, 2.0]
ssvi = srs.SsviParams(rho=-0.4, eta=0.6, gamma=0.4)
slices = [(ks, np.array([ssvi.total_variance(k, 0.04 * t) for k in ks]), 0.04 * t) for t in maturities]

surface = srs.EssviSurface.calibrate(maturities, slices)     # (ks, total_variances, theta) per slice
print(surface.slices())                                      # [(t, theta, rho, psi), ...]
print(surface.is_butterfly_free(), surface.is_calendar_spread_free())
print(surface.implied_vol(0.1, 0.75))                        # interpolated between the 0.5y and 1y slices

// docs: quant#regularised-calibration
//! Backs the regularised calibration example on the quant catalog page.

use ndarray::Array1;
use stochastic_rs::quant::calibration::Regularization;
use stochastic_rs::quant::calibration::SabrCalibrator;
use stochastic_rs::quant::pricing::sabr::SabrPricer;
use stochastic_rs::quant::types::OptionType;
use stochastic_rs::traits::Calibrator;

#[test]
fn a_tikhonov_anchor_pulls_the_sabr_fit() {
  // Synthetic smile from SABR (α = 0.2, β = 1, ν = 0.6, ρ = −0.3), one year, seven strikes.
  let (s, r, tau) = (100.0, 0.01, 1.0);
  let strikes = vec![80.0, 90.0, 95.0, 100.0, 105.0, 110.0, 120.0];
  let pricer = SabrPricer::new(0.2, 1.0, 0.6, -0.3);
  let prices: Vec<f64> = strikes
    .iter()
    .map(|&k| pricer.call_put(s, k, r, 0.0, tau).0)
    .collect();
  let calibrator = |regularization: Option<Regularization>| {
    let mut c = SabrCalibrator::new(
      None,
      Array1::from_vec(prices.clone()),
      Array1::from_elem(strikes.len(), s),
      Array1::from_vec(strikes.clone()),
      r,
      None,
      tau,
      OptionType::Call,
      false,
    );
    c.regularization = regularization;
    c
  };

  // Plain least squares recovers ν; a heavy anchor at ν⁰ = 0.9 (weights in price² units,
  // natural order (α, ν, ρ)) pulls the fit there.
  let plain = calibrator(None).calibrate(None).unwrap();
  let pulled = calibrator(Some(Regularization::new(
    vec![0.2, 0.9, -0.3],
    vec![0.0, 1e4, 0.0],
  )))
  .calibrate(None)
  .unwrap();
  assert!((plain.nu - 0.6).abs() < 0.05 && (pulled.nu - 0.9).abs() < 0.05);
}

import stochastic_rs as srs

strikes = [80.0, 90.0, 95.0, 100.0, 105.0, 110.0, 120.0]
prices = [srs.SabrPricer(s=100.0, alpha=0.2, beta=1.0, nu=0.6, rho=-0.3, k=k, r=0.01, tau=1.0).price() for k in strikes]

plain = srs.SabrCalibrator(strikes, prices, s=100.0, r=0.01, tau=1.0).calibrate()
pulled = srs.SabrCalibrator(strikes, prices, s=100.0, r=0.01, tau=1.0,
                            regularization=([0.2, 0.9, -0.3], [0.0, 1e4, 0.0])).calibrate()
print(plain[2], pulled[2])     # nu ≈ 0.6 without the anchor, ≈ 0.9 with it

// docs: quant#heston-pde-by-adi-finite-differences
//! Backs the Heston ADI example on the quant catalog page.

use stochastic_rs::quant::pricing::heston::HestonPricer;
use stochastic_rs::quant::pricing::heston_adi::AdiScheme;
use stochastic_rs::quant::pricing::heston_adi::HestonAdiPricer;
use stochastic_rs::traits::ModelPricer;

#[test]
fn adi_solver_reprices_the_semi_analytic_heston_call() {
  // Case 1 of in 't Hout & Foulon: κ = 1.5, η = 0.04, σ = 0.3, ρ = −0.9, r_d = 2.5 %, one year.
  let (v0, kappa, eta, sigma, rho) = (0.04, 1.5, 0.04, 0.3, -0.9);
  let (s, k, r_d, r_f, tau) = (100.0, 100.0, 0.025, 0.0, 1.0);

  // Model + numerics state; the query travels as arguments (r = r_d, q = r_f).
  let pde = HestonAdiPricer::new(v0, kappa, eta, sigma, rho)
    .with_grid(100, 50, 50)
    .with_scheme(AdiScheme::ModifiedCraigSneyd);
  let adi = pde.price_call(s, k, r_d, r_f, tau);

  // Semi-analytic reference through the crate's Heston pricer.
  let analytic = HestonPricer::new(v0, rho, kappa, eta, sigma, Some(0.0))
    .call_put(s, k, r_d, r_f, tau)
    .0;
  assert!(
    (adi - analytic).abs() / analytic < 1e-2,
    "adi {adi} vs analytic {analytic}"
  );

  // The same grid prices a down-and-out call by moving the lower boundary to the barrier.
  let knocked = pde.with_barrier(90.0).price_call(s, k, r_d, r_f, tau);
  assert!(knocked > 0.0 && knocked < adi);
}

import stochastic_rs as srs

pde = srs.HestonAdiPricer(s=100.0, k=100.0, r=0.025, tau=1.0, v0=0.04, kappa=1.5, theta=0.04, sigma=0.3, rho=-0.9,
                          q=0.0, m1=100, m2=50, steps=50, scheme="mcs")
call, put = pde.call_put()
analytic = srs.HestonPricer(s=100.0, k=100.0, r=0.025, tau=1.0, v0=0.04, kappa=1.5, theta=0.04, sigma=0.3, rho=-0.9, q=0.0).price()
print(call, analytic)                     # agree to ~1 %
print(srs.HestonAdiPricer(s=100.0, k=100.0, r=0.025, tau=1.0, v0=0.04, kappa=1.5, theta=0.04, sigma=0.3, rho=-0.9,
                          barrier=90.0).price())   # down-and-out call

// docs: quant#xva-exposure-profiles-and-cva--dva--fva
//! Backs the XVA example on the quant catalog page.

use ndarray::Array1;
use stochastic_rs::quant::credit::survival_curve::HazardInterpolation;
use stochastic_rs::quant::credit::survival_curve::SurvivalCurve;
use stochastic_rs::quant::curves::DiscountCurve;
use stochastic_rs::quant::curves::InterpolationMethod;
use stochastic_rs::quant::risk::xva::cva;
use stochastic_rs::quant::risk::xva::fva;
use stochastic_rs::quant::risk::xva::irs::HullWhiteSwapExposure;
use stochastic_rs::simd_rng::Deterministic;

#[test]
fn cva_of_a_par_swap_under_hull_white() {
  // Flat 3 % curve, five-year annual payer swap of 1m notional at par.
  let curve = DiscountCurve::from_zero_rates(
    &Array1::from_vec(vec![0.5, 1.0, 5.0, 10.0]),
    &Array1::from_vec(vec![0.03; 4]),
    InterpolationMethod::LogLinearOnDiscountFactors,
  );
  let mut swap = HullWhiteSwapExposure::new(
    0.1,
    0.01,
    1_000_000.0,
    0.0,
    vec![1.0, 2.0, 3.0, 4.0, 5.0],
    1.0,
  );
  swap.fixed_rate = swap.par_rate(&curve);

  // Exposure profile on the payment dates from 4 000 Hull–White short-rate paths.
  let profile = swap.profile(&curve, 4_000, 0.95, Deterministic::new(7));
  assert!(profile.peak_epe() > 0.0 && profile.epe[4] == 0.0);

  // CVA against a 2 % flat hazard at 60 % LGD, and the symmetric FVA at a 50 bp funding spread.
  let counterparty = SurvivalCurve::from_hazard_rates(
    &Array1::from_vec(vec![1.0, 5.0, 10.0]),
    &Array1::from_vec(vec![0.02; 3]),
    HazardInterpolation::PiecewiseConstantHazard,
  );
  let cva_value = cva(&profile, &counterparty, &curve, 0.6);
  assert!(cva_value > 0.0 && cva_value < 0.01 * swap.notional);
  let _funding = fva(&profile, &curve, 0.005);
}

import numpy as np
import stochastic_rs as srs

curve = srs.DiscountCurve.from_zero_rates(np.array([0.5, 1.0, 5.0, 10.0]), np.full(4, 0.03), interp="log_df")
swap = srs.HullWhiteSwapExposure(mean_reversion=0.1, sigma=0.01, notional=1e6, payment_times=[1.0, 2.0, 3.0, 4.0, 5.0])
profile = swap.profile(curve, paths=4000, quantile=0.95, seed=7)     # at par on the curve
print(profile.epe(), profile.pfe())
print(profile.cva(hazard_rate=0.02, discount=curve, lgd=0.6), profile.fva(curve, funding_spread=0.005))

# Any MtM matrix works: here a Brownian one, whose EPE is σ√t/√(2π).
mtm = 10.0 * np.cumsum(np.random.default_rng(1).standard_normal((20000, 4)) * np.sqrt(0.5), axis=1)
print(srs.ExposureProfile.from_mtm(mtm, [0.5, 1.0, 1.5, 2.0]).epe())

// docs: quant#cds-index-and-cdo-tranches
//! Backs the CDS index and tranche example on the quant catalog page.

use chrono::NaiveDate;
use ndarray::Array1;
use stochastic_rs::quant::credit::index::CdsIndex;
use stochastic_rs::quant::credit::index::flat_survival;
use stochastic_rs::quant::credit::tranche::CdoTranche;
use stochastic_rs::quant::credit::tranche::PoolName;
use stochastic_rs::quant::curves::DiscountCurve;
use stochastic_rs::quant::curves::InterpolationMethod;

#[test]
fn index_upfront_and_tranche_spreads() {
  let discount = DiscountCurve::from_zero_rates(
    &Array1::from_vec(vec![0.5, 1.0, 5.0, 10.0]),
    &Array1::from_vec(vec![0.03; 4]),
    InterpolationMethod::LogLinearOnDiscountFactors,
  );
  let date = |y: i32, m: u32, d: u32| NaiveDate::from_ymd_opt(y, m, d).unwrap();

  // A 125-name index at a 100 bp coupon: fair spread and the ISDA upfront for a 120 bp quote.
  let index = CdsIndex::homogeneous(
    vec![flat_survival(0.02); 125],
    0.4,
    0.01,
    10_000_000.0,
    date(2026, 3, 20),
    date(2031, 6, 20),
  );
  let fair = index.fair_spread(date(2026, 3, 20), &discount);
  let upfront = index.isda_upfront(date(2026, 3, 20), &discount, 0.012, 0.4);
  assert!(fair > 0.0 && upfront > 0.0);

  // Equity and mezzanine tranches on the same pool under a 30 % Gaussian copula correlation.
  let pool: Vec<PoolName> = (0..125)
    .map(|_| PoolName {
      weight: 1.0 / 125.0,
      recovery: 0.4,
      survival: flat_survival(0.02),
    })
    .collect();
  let times: Vec<f64> = (1..=5).map(|i| i as f64).collect();
  let equity =
    CdoTranche::new(0.0, 0.03, 0.05, times.clone(), 1.0, 0.3).valuation(&pool, &discount);
  let mezz = CdoTranche::new(0.03, 0.07, 0.01, times, 1.0, 0.3).valuation(&pool, &discount);
  assert!(equity.fair_spread > mezz.fair_spread && mezz.fair_spread > 0.0);
}

import numpy as np
import stochastic_rs as srs

curve = srs.DiscountCurve.from_zero_rates(np.array([0.5, 1.0, 5.0, 10.0]), np.full(4, 0.03), interp="log_df")
names = [(1 / 125, 0.4, 0.02)] * 125                     # (weight, recovery, hazard)
index = srs.CdsIndex(names, coupon=0.01, notional=1e7, effective_date="2026-03-20", maturity_date="2031-06-20")
print(index.fair_spread("2026-03-20", curve), index.isda_upfront("2026-03-20", curve, quoted_spread=0.012))

equity = srs.CdoTranche(names, attachment=0.0, detachment=0.03, spread=0.05, payment_times=[1, 2, 3, 4, 5], correlation=0.3)
print(equity.valuation(curve))                            # (protection, annuity, premium, fair_spread, upfront)
print(equity.expected_tranche_loss(5.0), equity.large_pool_expected_tranche_loss(p=1 - np.exp(-0.1), lgd=0.6))

// docs: quant#heston-fourier-pricer-with-cui-analytic-jacobian
//! Backs the Heston pricer example on the quant catalog page. The model
//! holds only its six Heston parameters; spot, strike, rate, dividend
//! yield and maturity travel to the call, so one instance prices a whole
//! strike/maturity grid.

use stochastic_rs::quant::OptionType;
use stochastic_rs::quant::pricing::heston::HestonPricer;
use stochastic_rs::traits::ModelPricer;

#[test]
fn heston_pricer_price_and_greeks() {
  let model = HestonPricer::new(
    /* v0 */ 0.04, /* rho */ -0.5, /* kappa */ 2.0, /* theta */ 0.04,
    /* sigma */ 0.3, /* lambda */ None,
  );
  let (s, k, r, q, tau) = (100.0, 100.0, 0.03, 0.0, 1.0);

  let price = model.price_call(s, k, r, q, tau);
  let g = model.greeks(s, k, r, q, tau, OptionType::Call);

  assert!(price > 0.0);
  assert!((0.0..=1.0).contains(&g.delta));
  assert!(g.vega > 0.0);
  assert!(g.vanna.is_finite());
}

import stochastic_rs as srs

pricer = srs.HestonPricer(
    s=100, v0=0.04, k=100, r=0.03, kappa=2.0, theta=0.04, sigma=0.3,
    rho=-0.5, tau=1.0, q=0.0,
)
call, put = pricer.call_put()
print(f"call={call:.4f}, put={put:.4f}")

// docs: quant#heston-calibration-to-a-market-vol-surface
//! Backs the Heston calibration example on the quant catalog page.
//! `HestonCalibrator` fits a single maturity slice directly against market
//! call prices (spot/strike/price vectors) rather than an
//! `ImpliedVolSurface` object.

use stochastic_rs::quant::calibration::heston::HestonCalibrator;
use stochastic_rs::quant::calibration::heston::HestonParams;
use stochastic_rs::quant::types::OptionType;
use stochastic_rs::traits::CalibrationResult;
use stochastic_rs::traits::Calibrator;

#[test]
fn heston_calibrator_recovers_plausible_params() {
  let s = vec![100.0; 9];
  let k = vec![80.0, 85.0, 90.0, 95.0, 100.0, 105.0, 110.0, 115.0, 120.0];
  let c_market = vec![21.5, 17.9, 14.2, 11.0, 8.2, 6.0, 4.3, 3.1, 2.2];

  let calibrator = HestonCalibrator::new(
    Some(HestonParams {
      v0: 0.04,
      kappa: 1.5,
      theta: 0.04,
      sigma: 0.5,
      rho: -0.7,
    }),
    c_market.into(),
    s.into(),
    k.into(),
    /* r */ 0.01,
    /* q */ Some(0.0),
    /* tau */ 0.5,
    OptionType::Call,
    None,
    None,
    None,
    /* record_history */ true,
  );

  let result = calibrator.calibrate(None).unwrap();
  let p = result.params();

  assert!(p.kappa > 0.0);
  assert!(p.theta > 0.0);
  assert!(p.sigma > 0.0);
  assert!((-1.0..=1.0).contains(&p.rho));
  assert!(result.rmse() >= 0.0);
}

import stochastic_rs as srs

# The calibrator takes one MarketSlice per expiry — strikes, observed
# prices, per-strike call/put flags, and tau — not a surface object.
slices = [
    srs.MarketSlice(strikes, prices, [True] * len(strikes), tau)
    for strikes, prices, tau in quotes
]
v0, kappa, theta, sigma, rho, converged, rmse = srs.HestonCalibrator(
    slices, s=100.0, r=0.03
).calibrate()
print(f"kappa={kappa:.3f}, theta={theta:.4f}, sigma={sigma:.3f}, "
      f"rho={rho:.3f}, v0={v0:.4f}")
print(f"RMSE = {rmse:.6f} (converged={converged})")

// docs: quant#heston-slv-leverage-calibration-by-the-particle-method
//! Backs the Heston SLV calibration example on the quant catalog page: a
//! vanilla call surface in, a leverage-calibrated model and its Monte Carlo
//! pricer out.

use ndarray::Array1;
use ndarray::Array2;
use stochastic_rs::quant::calibration::heston::HestonParams;
use stochastic_rs::quant::calibration::heston_slv::HestonSlvCalibrator;
use stochastic_rs::quant::pricing::fourier::HestonFourier;
use stochastic_rs::quant::pricing::slv::ParticleMethod;
use stochastic_rs::traits::CalibrationResult;
use stochastic_rs::traits::Calibrator;
use stochastic_rs::traits::ModelPricer;

#[test]
fn heston_slv_calibrates_a_vanilla_surface_and_prices_under_half_mixing() {
  let (s, r, q) = (100.0, 0.02, 0.0);
  let heston = HestonParams {
    v0: 0.04,
    kappa: 2.0,
    theta: 0.05,
    sigma: 0.4,
    rho: -0.6,
  };
  // The vanilla surface — generated from a Heston model here, market quotes
  // in practice: `calls[[j, i]]` is the call at `strikes[i]`, `maturities[j]`.
  let strikes = Array1::linspace(70.0, 140.0, 36).to_vec();
  let maturities = vec![0.25, 0.5, 0.75, 1.0];
  let model = HestonFourier {
    v0: heston.v0,
    kappa: heston.kappa,
    theta: heston.theta,
    sigma: heston.sigma,
    rho: heston.rho,
    r,
    q,
  };
  let calls = Array2::from_shape_fn((maturities.len(), strikes.len()), |(j, i)| {
    model.price_call(s, strikes[i], r, q, maturities[j])
  });

  // Half the vol-of-vol; the leverage makes up the difference so the surface
  // is reproduced. The Heston fit is pinned here — leave it out to have
  // `HestonCalibrator` fit the same quotes.
  let result = HestonSlvCalibrator::new(s, r, q, strikes, maturities, calls)
    .with_mixing(0.5)
    .with_heston_params(heston)
    .with_particle_method(
      ParticleMethod::default()
        .with_particles(5_000)
        .with_steps_per_year(50)
        .with_seed(7),
    )
    .calibrate(None)
    .unwrap();
  assert!(result.converged());
  assert!(
    result.rmse() < 1.0,
    "in-sample repricing rmse {}",
    result.rmse()
  );
  assert!(result.leverage().covers(100.0, 0.5));

  // The pricer is anchored to the calibration rates and reprices the market.
  let pricer = result
    .to_model(r, q)
    .with_paths(4_000)
    .with_steps_per_year(50);
  let call = pricer.price_call(s, 100.0, r, q, 0.5);
  let market = model.price_call(s, 100.0, r, q, 0.5);
  assert!((call - market).abs() < 1.0, "slv {call} vs market {market}");
}

import stochastic_rs as srs

# calls[j, i] = C(strikes[i], maturities[j]); eta scales the vol-of-vol and
# the leverage makes up the difference, so the surface is reproduced either way.
calibrator = srs.HestonSlvCalibrator(
    s=100.0, r=0.02, q=0.0, strikes=list(strikes), maturities=list(maturities),
    calls=calls, eta=0.5, n_particles=50_000, steps_per_year=100,
)   # method="fokker_planck" swaps the cloud for the finite-volume PDE
result = calibrator.calibrate()
kappa, theta, sigma, rho, v0, eta = result.params()
print(f"rmse = {result.rmse:.4f} (converged={result.converged})")
surface = result.leverage()             # .spots, .times, .values, .interpolate(s, t)
pricer = result.to_model(n_paths=50_000)  # anchored to the calibration (r, q)
print(pricer.price_call(100.0, 100.0, 0.02, 0.0, 0.5))