optionstratlib-pricing 0.22.0

/******************************************************************************
   Author: Joaquín Béjar García
   Email: jb@taunais.com
   Date: 15/8/24
******************************************************************************/

use crate::error::VolatilityError;
use num_traits::FromPrimitive;
#[cfg(test)]
use num_traits::ToPrimitive;
use optionstratlib_core::constants::MIN_VOLATILITY;
use optionstratlib_core::model::Positive;
use optionstratlib_core::model::decimal::{
    d_add, d_div, d_mul, d_sqrt, d_sub, d_sum, decimal_to_f64, finite_decimal, p_sqrt,
};
use optionstratlib_core::model::{ExpirationDate, OptionStyle, OptionType, Options, Side};
use optionstratlib_core::utils::time::TimeFrame;
use rand::Rng;
use rand_distr::{Distribution, StandardNormal};
use rust_decimal::{Decimal, RoundingStrategy};
use tracing::instrument;

use crate::pricing::OptionPricing;
use crate::pricing::black_scholes_model::european_price_band;
use crate::pricing::constants::IV_TOLERANCE;
use crate::pricing::iv_solver::{
    IV_BRACKET_MAX, IvOutcome, iv_initial_guess, long_residual, solve_implied_volatility,
};
#[cfg(test)]
use optionstratlib_core::pos_or_panic;

/// Calculates the constant volatility from a series of returns.
///
/// # Arguments
///
/// * `returns` - A slice of Decimal values representing the returns.
///
/// # Returns
///
/// The calculated volatility as a Decimal.
///
/// # Errors
///
/// Returns `VolatilityError::NumericalFailure` when the length of
/// `returns` cannot be represented as a `Decimal` or when the
/// final `variance.sqrt()` fails (overflow).
///
/// When `returns.len() < 2` the function returns `Ok(Positive::ZERO)`
/// rather than an error. The count and the final square root are never
/// negative, so the [`VolatilityError::PositiveError`] their `Positive`
/// conversions could report does not occur in practice.
pub fn constant_volatility(returns: &[Decimal]) -> Result<Positive, VolatilityError> {
    let n_dec =
        Decimal::from_usize(returns.len()).ok_or_else(|| VolatilityError::NumericalFailure {
            reason: format!(
                "constant_volatility: returns.len() {} not representable as Decimal",
                returns.len()
            ),
        })?;
    // A length is never negative: the `?` cannot fire.
    let n = Positive::new_decimal(n_dec)?;

    if n < Decimal::TWO {
        return Ok(Positive::ZERO);
    }

    // Sample mean and sample variance. Every step is checked:
    //   * the raw sum is folded with `d_sum`;
    //   * the mean projection uses `d_div` with banker's rounding;
    //   * each centred deviation is built with `d_sub`;
    //   * the square is built with `d_mul` (no `.powi(2)`, which is
    //     unchecked and can saturate on adversarial inputs);
    //   * the sum of squares is folded inline so we do not allocate a
    //     temporary `Vec`.
    let total = d_sum(returns, "volatility::constant::total")?;
    let mean = d_div(total, n.to_dec(), "volatility::constant::mean")?;
    let mut sq_total = Decimal::ZERO;
    for &r in returns {
        let centred = d_sub(r, mean, "volatility::constant::centred")?;
        let centred_sq = d_mul(centred, centred, "volatility::constant::centred_sq")?;
        sq_total = d_add(sq_total, centred_sq, "volatility::constant::sq_total")?;
    }
    let degrees_of_freedom = d_sub(n.to_dec(), Decimal::ONE, "volatility::constant::dof")?;
    let variance = d_div(
        sq_total,
        degrees_of_freedom,
        "volatility::constant::variance",
    )?;

    let std_dev = d_sqrt(variance, "volatility::constant::std_dev").map_err(|_| {
        VolatilityError::NumericalFailure {
            reason: "constant_volatility: sqrt(variance) failed (overflow)".to_string(),
        }
    })?;
    // A square root is never negative: the `?` cannot fire.
    Ok(Positive::new_decimal(std_dev)?)
}

/// Calculates historical volatility using a moving window approach.
///
/// # Arguments
///
/// * `returns` - A slice of Decimal values representing the returns.
/// * `window_size` - The size of the moving window.
///
/// # Returns
///
/// A vector of Decimal values representing the historical volatility for each window.
///
/// # Errors
///
/// Returns [`VolatilityError::NumericalFailure`] when `window_size`
/// is zero: `slice::windows` has no meaning at that width and
/// panics rather than yielding an empty iterator.
///
/// Propagates any `VolatilityError::NumericalFailure` raised by
/// `constant_volatility` on an individual window.
///
/// When `window_size > returns.len()` the function returns
/// `Ok(vec![])` instead of an error, so callers that need a
/// non-empty result must validate the inputs themselves.
pub fn historical_volatility(
    returns: &[Decimal],
    window_size: usize,
) -> Result<Vec<Positive>, VolatilityError> {
    if window_size == 0 {
        return Err(VolatilityError::NumericalFailure {
            reason: "historical_volatility: window_size must be non-zero".to_string(),
        });
    }
    match rolling_historical_volatility(returns, window_size) {
        Some(volatilities) => Ok(volatilities),
        // An intermediate left the `Decimal` range: the per-window two-pass
        // estimate reports the error it always reported, on the same window.
        None => returns
            .windows(window_size)
            .map(constant_volatility)
            .collect(),
    }
}

/// [`historical_volatility`] in `O(n)`: one pass of rolling sums instead of a
/// two-pass [`constant_volatility`] over every window (#859).
///
/// Each window's sample variance is `(S2 - S1² / w) / (w - 1)` with
/// `S1 = Σ (r - K)` and `S2 = Σ (r - K)²`, kept up to date by adding the
/// return that enters the window and subtracting the one that leaves. `K` is
/// the first window's mean: shifting the data by a value near its mean is the
/// standard way to keep `S2 - S1² / w` from cancelling (Chan, Golub and
/// LeVeque, 1983). The result can differ from the two-pass estimate in the
/// last digits of a `Decimal`, never more than the tests' `1e-10` relative
/// tolerance on the fixtures.
///
/// A window whose rolling variance is not positive (a rounding residue on a
/// window of equal returns) takes the two-pass estimate, so no value is
/// invented. Returns `None`, and the caller falls back to the two-pass
/// estimate for every window, when an intermediate leaves the `Decimal`
/// range or a window cannot be evaluated.
fn rolling_historical_volatility(returns: &[Decimal], window_size: usize) -> Option<Vec<Positive>> {
    let first = returns.get(..window_size)?;
    if window_size < 2 {
        // `constant_volatility` reports zero for fewer than two returns.
        return Some(vec![Positive::ZERO; returns.len()]);
    }
    let width = Decimal::from_usize(window_size)?;
    let dof = width.checked_sub(Decimal::ONE)?;
    let shift = d_div(d_sum(first, "").ok()?, width, "").ok()?;

    let mut s1 = Decimal::ZERO;
    let mut s2 = Decimal::ZERO;
    for &r in first {
        let x = r.checked_sub(shift)?;
        s1 = s1.checked_add(x)?;
        s2 = s2.checked_add(x.checked_mul(x)?)?;
    }

    let count = returns.len().checked_sub(window_size)?.checked_add(1)?;
    let mut volatilities = Vec::new();
    volatilities.try_reserve_exact(count).ok()?;

    let window_volatility = |s1: Decimal, s2: Decimal, window: &[Decimal]| -> Option<Positive> {
        let centred = s2.checked_sub(s1.checked_mul(s1)?.checked_div(width)?)?;
        let variance = centred.checked_div(dof)?;
        if variance > Decimal::ZERO {
            let std_dev = d_sqrt(variance, "").ok()?;
            Positive::new_decimal(std_dev).ok()
        } else {
            constant_volatility(window).ok()
        }
    };

    volatilities.push(window_volatility(s1, s2, first)?);
    let entering = returns.get(window_size..)?;
    for (start, (&new, &old)) in (1usize..).zip(entering.iter().zip(returns)) {
        let x_new = new.checked_sub(shift)?;
        let x_old = old.checked_sub(shift)?;
        s1 = s1.checked_add(x_new)?.checked_sub(x_old)?;
        s2 = s2
            .checked_add(x_new.checked_mul(x_new)?)?
            .checked_sub(x_old.checked_mul(x_old)?)?;
        let window = returns.get(start..start.checked_add(window_size)?)?;
        volatilities.push(window_volatility(s1, s2, window)?);
    }
    Some(volatilities)
}

/// Calculates EWMA (Exponentially Weighted Moving Average) volatility.
///
/// # Arguments
///
/// * `returns` - A slice of Decimal values representing the returns.
/// * `lambda` - The decay factor (typically 0.94 for daily data).
///
/// # Returns
///
/// A vector of Decimal values representing the EWMA volatility.
///
/// # Errors
///
/// Returns `VolatilityError::NumericalFailure` when `returns` is
/// empty or when any `variance.sqrt()` call fails (negative operand
/// reported as a generic overflow failure).
///
/// Each running std-dev is a square root, never negative, so the
/// [`VolatilityError::PositiveError`] its `Positive` conversion could
/// report does not occur in practice; a negative variance is already
/// caught by the `sqrt` check above.
pub fn ewma_volatility(
    returns: &[Decimal],
    lambda: Decimal,
) -> Result<Vec<Positive>, VolatilityError> {
    let first_return = returns
        .first()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "ewma_volatility: returns slice is empty".to_string(),
        })?;
    // `first_return^2` is built with a checked multiplication because
    // `.powi(2)` on `Decimal` silently saturates on overflow.
    let mut variance = d_mul(
        *first_return,
        *first_return,
        "volatility::ewma::initial_variance",
    )?;
    let initial_std_dev = d_sqrt(variance, "volatility::ewma::initial_std_dev").map_err(|_| {
        VolatilityError::NumericalFailure {
            reason: "ewma_volatility: sqrt(initial variance) failed (overflow)".to_string(),
        }
    })?;
    // Square roots are never negative: this `?` and the one in the loop
    // cannot fire.
    let mut volatilities = vec![Positive::new_decimal(initial_std_dev)?];

    for &return_value in returns.iter().skip(1) {
        // EWMA variance recursion: v' = λ·v + (1 - λ) · r².
        // `r²` is built via `d_mul(r, r, ..)` so that a saturating
        // squaring cannot feed a silently capped value into the
        // innovation term.
        let persistent = d_mul(lambda, variance, "volatility::ewma::persistent")?;
        let innovation_weight = d_sub(Decimal::ONE, lambda, "volatility::ewma::innovation_weight")?;
        let return_sq = d_mul(return_value, return_value, "volatility::ewma::return_sq")?;
        let innovation = d_mul(innovation_weight, return_sq, "volatility::ewma::innovation")?;
        variance = d_add(persistent, innovation, "volatility::ewma::variance")?;
        let std_dev = d_sqrt(variance, "volatility::ewma::std_dev").map_err(|_| {
            VolatilityError::NumericalFailure {
                reason: "ewma_volatility: sqrt(variance) failed (overflow)".to_string(),
            }
        })?;
        volatilities.push(Positive::new_decimal(std_dev)?);
    }

    Ok(volatilities)
}

/// Calculates the implied volatility of an option given its market price.
///
/// Returns the volatility at which one long unit of `options` prices
/// `market_price` under Black-Scholes. `options` is read, not modified.
///
/// # Parameters
/// - `market_price`: The observed market price of the option.
/// - `options`: The option whose volatility is implied; its side, quantity
///   and current `implied_volatility` are ignored.
/// - `max_iterations`: The maximum number of solver steps.
///
/// # Returns
/// The implied volatility, in `[MIN_VOLATILITY, 5]`.
///
/// # Method
///
/// The volatility is found by Newton's method on the Black–Scholes vega,
/// kept inside the bracket `[MIN_VOLATILITY, 5]` with a bisection fallback
/// (`pricing::iv_solver`, #859 P4); `max_iterations` caps its steps.
/// It used to be a parallel grid search over `100 * max_iterations`
/// candidates, which returned the nearest grid point instead of the root.
///
/// # Errors
///
/// Returns `VolatilityError::IvNotFound` when a European target lies outside
/// the no-arbitrage band of Black–Scholes prices, below
/// `max(S e^(-qT) - K e^(-rT), 0)` (call) / `max(K e^(-rT) - S e^(-qT), 0)`
/// (put) or above `S e^(-qT)` / `K e^(-rT)`, widened by `1e-5` (#652), when
/// the target lies outside the prices of the bracket, or when the solver
/// does not converge within `max_iterations` steps;
/// `VolatilityError::NoValidVolatility` when the option has expired (at
/// `T = 0` the price is the intrinsic value whatever the volatility, #843)
/// or `max_iterations` is not positive; `VolatilityError::Options` or
/// `VolatilityError::DecimalError` when a price, the band or the arithmetic
/// fails.
#[instrument(skip(options), fields(
    market_price = %market_price,
    strike = %options.strike_price,
    max_iterations,
))]
pub fn implied_volatility(
    market_price: Positive,
    options: &mut Options,
    max_iterations: i64,
) -> Result<Positive, VolatilityError> {
    // At expiry the price is the intrinsic value for every volatility
    // (#843), so no volatility is implied.
    if options.time_to_expiration()?.is_zero() {
        return Err(VolatilityError::NoValidVolatility);
    }
    if max_iterations <= 0 {
        return Err(VolatilityError::NoValidVolatility);
    }
    let base_option = options.clone();
    let target = market_price.to_dec();
    // A European target outside the no-arbitrage band has no implied
    // volatility (#652).
    if matches!(base_option.option_type, OptionType::European) {
        let (lower, upper) = european_price_band(&base_option)
            .map_err(optionstratlib_core::error::OptionsError::from)?;
        if target < d_sub(lower, IV_TOLERANCE, "volatility::iv::band::lower")?
            || target > d_add(upper, IV_TOLERANCE, "volatility::iv::band::upper")?
        {
            return Err(VolatilityError::IvNotFound);
        }
    }
    // A cap above `u32::MAX` iterations is no cap the solver can reach.
    let cap = u32::try_from(max_iterations).unwrap_or(u32::MAX);
    let guess = iv_initial_guess(&base_option, target).unwrap_or(IV_BRACKET_MAX);
    match solve_implied_volatility(
        *MIN_VOLATILITY,
        IV_BRACKET_MAX,
        guess,
        cap,
        |sigma, with_slope| long_residual(&base_option, target, sigma, with_slope),
    )? {
        IvOutcome::Found(iv) => Ok(iv),
        IvOutcome::BelowBracket | IvOutcome::AboveBracket | IvOutcome::NoConvergence(_) => {
            Err(VolatilityError::IvNotFound)
        }
    }
}

/// Calculates the implied volatility (IV) of an option given its parameters.
///
/// # Parameters
///
/// * `option_price` - A `Positive` value representing the current price of the option.
/// * `strike` - A `Positive` value representing the strike price of the option.
/// * `option_style` - An `OptionStyle` enum indicating the style of the option (e.g., European, American).
/// * `underlying_price` - A `Positive` value representing the current price of the underlying asset.
/// * `days` - A `Positive` value indicating the number of days to expiration for the option.
/// * `symbol` - A `String` representing the symbol of the financial instrument for the option.
///
/// # Returns
///
/// * `Ok(Positive)` - The calculated implied volatility as a `Positive` value, if the computation is successful.
/// * `Err(Box<dyn Error>)` - An error if the implied volatility cannot be calculated due to invalid inputs or other reasons.
///
/// # Errors
///
/// This function will return an error if:
/// * The inputs do not meet the required constraints.
/// * The implied volatility calculation fails to converge within the set iteration limit.
///
/// # Notes
///
/// This function internally creates an `Options` object with the given parameters,
/// and calls the `implied_volatility` function with the option data. The solver
/// is capped at 100 steps.
///
/// Ensure that all input parameters are valid and conform to the expected types
/// and ranges for meaningful results.
pub fn calculate_iv(
    option_price: Positive,
    strike: Positive,
    option_style: OptionStyle,
    underlying_price: Positive,
    days: Positive,
    symbol: String,
) -> Result<Positive, VolatilityError> {
    let mut option = Options::new(
        OptionType::European,
        Side::Long,
        symbol,
        strike,
        ExpirationDate::Days(days),
        Positive::ZERO,
        Positive::ONE,
        underlying_price,
        Decimal::ZERO,
        option_style,
        Positive::ZERO,
        None,
    );
    // 100 solver steps: Newton converges in a handful, and the cap leaves
    // room for the bisection fallback (the grid this replaced read 10 as
    // 1 000 candidates, #859 P4).
    implied_volatility(option_price, &mut option, 100)
}

/// Calculates GARCH(1,1) volatility (simplified).
///
/// # Arguments
///
/// * `returns` - A slice of Decimal values representing the returns.
/// * `omega`, `alpha`, `beta` - GARCH(1,1) parameters.
///
/// # Returns
///
/// A vector of Decimal values representing the GARCH(1,1) volatility.
///
/// # Errors
///
/// - Propagates [`optionstratlib_core::error::DecimalError::Overflow`] from the underlying
///   checked arithmetic helpers (`d_mul`, `d_add`) wrapped as
///   [`VolatilityError::DecimalError`] via the `#[from]` cascade.
///   This happens when any of the four monetary products
///   (`r^2`, `α · r²`, `β · v_prev`, `ω + α · r² + β · v_prev`)
///   exceeds the representable `Decimal` range, or when the
///   squared initial return itself overflows.
/// - Returns [`VolatilityError::NumericalFailure`] when a recursion
///   step produces a negative variance (e.g. pathological GARCH
///   parameters with `omega` or `beta * v_prev` negative). This is
///   a real diagnostic now instead of being silently clamped to
///   `Positive::ZERO` as in earlier implementations.
/// - Returns [`VolatilityError::NumericalFailure`] when `returns` is
///   empty: the recursion has no seed variance to start from.
/// - Returns [`VolatilityError::PositiveError`] if a running
///   square root (`p_sqrt`) ever rejects its operand. A variance is
///   non-negative by construction here, so this arm is a guard
///   rather than a reachable outcome.
pub fn garch_volatility(
    returns: &[Decimal],
    omega: Decimal,
    alpha: Decimal,
    beta: Decimal,
) -> Result<Vec<Positive>, VolatilityError> {
    // Helper: project a variance `Decimal` onto `Positive`, mapping
    // the rejection (negative variance) onto a readable
    // `VolatilityError::NumericalFailure` instead of silently
    // clamping to zero.
    let to_positive_variance =
        |value: Decimal, context: &'static str| -> Result<Positive, VolatilityError> {
            Positive::new_decimal(value).map_err(|_| VolatilityError::NumericalFailure {
                reason: format!("garch_volatility: negative variance at {context} ({value})"),
            })
        };

    // Seed the recursion with `r_0^2` via `d_mul` so a saturating
    // squaring cannot feed a silently capped value into the GARCH
    // update below. `first()` is what makes the seed total: indexing
    // `returns[0]` aborted on an empty slice.
    let first_return = returns
        .first()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "garch_volatility: returns slice is empty".to_string(),
        })?;
    let seed = d_mul(
        *first_return,
        *first_return,
        "volatility::garch::initial_variance",
    )?;
    let mut variance = to_positive_variance(seed, "initial_variance")?;
    let mut volatilities = vec![p_sqrt(&variance, "volatility::utils::garch_volatility")?];
    for &return_value in returns.iter().skip(1) {
        // GARCH(1, 1) variance recursion:
        // v' = ω + α · r² + β · v_prev.
        let return_sq = d_mul(return_value, return_value, "volatility::garch::return_sq")?;
        let alpha_r2 = d_mul(alpha, return_sq, "volatility::garch::alpha_r2")?;
        let beta_v = d_mul(beta, variance.to_dec(), "volatility::garch::beta_v")?;
        let persistent = d_add(omega, alpha_r2, "volatility::garch::omega_plus_alpha")?;
        let next_variance = d_add(persistent, beta_v, "volatility::garch::variance")?;
        variance = to_positive_variance(next_variance, "variance")?;
        volatilities.push(p_sqrt(&variance, "volatility::utils::garch_volatility")?);
    }
    Ok(volatilities)
}

/// Draws one Wiener increment `dW ~ N(0, dt)` for the Heston Euler step: a
/// standard normal from `rng` scaled by `dt_sqrt = sqrt(dt)` (#742).
#[inline]
fn heston_wiener_increment<R: Rng + ?Sized>(
    rng: &mut R,
    dt_sqrt: f64,
) -> Result<Decimal, VolatilityError> {
    let z: f64 = StandardNormal.sample(rng);
    let dw = z * dt_sqrt;
    finite_decimal(dw).ok_or_else(|| VolatilityError::non_finite("volatility::heston::dw", dw))
}

/// Simulates stochastic volatility using the Heston model (simplified).
///
/// # Arguments
///
/// * `kappa` - Mean reversion speed.
/// * `theta` - Long-term variance.
/// * `xi` - Volatility of volatility.
/// * `v0` - Initial variance.
/// * `dt` - Time step.
/// * `steps` - Number of simulation steps.
/// * `rng` - The generator every variance shock is drawn from, one
///   standard-normal draw per step after the first, scaled by `sqrt(dt)` to
///   the Wiener increment `dW ~ N(0, dt)`. A seeded generator such as
///   [`optionstratlib_core::utils::deterministic_rng`] makes the path
///   reproducible; pass `&mut rand::rng()` to draw from the thread-local RNG.
///
/// # Returns
///
/// A vector of Decimal values representing the simulated volatility.
///
/// # Errors
///
/// Returns `VolatilityError::NumericalFailure` when the initial
/// `sqrt(dt)` call overflows, when that square root is not
/// representable as `f64`, or when the per-step `dw` sample cannot
/// be converted back into `Decimal`.
///
/// Returns [`VolatilityError::DecimalError`] when a step of the
/// Euler discretisation `v' = v + κ(θ − v)dt + ξ√v dW` leaves the
/// representable `Decimal` range. The raw operators this used to be
/// written with aborted instead; there is no limit to substitute,
/// because the variance is the state being simulated and a step that
/// cannot be represented ends the path.
///
/// Variance is floored at zero on each step (`v = v.max(Decimal::ZERO)`),
/// the documented full-truncation floor of the Euler scheme, so its
/// `Positive` conversion cannot fail.
pub fn simulate_heston_volatility<R: Rng + ?Sized>(
    kappa: Decimal,
    theta: Decimal,
    xi: Decimal,
    v0: Decimal,
    dt: Decimal,
    steps: usize,
    rng: &mut R,
) -> Result<Vec<Positive>, VolatilityError> {
    let mut v = v0.max(Decimal::ZERO);
    // `v` is floored at zero here and on every step: the `?` cannot fire.
    let mut v_pos = Positive::new_decimal(v)?;
    let mut volatilities = vec![p_sqrt(
        &v_pos,
        "volatility::utils::simulate_heston_volatility",
    )?];
    let dt_sqrt = d_sqrt(dt, "volatility::heston::sqrt_dt").map_err(|_| {
        VolatilityError::NumericalFailure {
            reason: "simulate_heston_volatility: sqrt(dt) failed (overflow)".to_string(),
        }
    })?;
    let dt_sqrt_f64 = decimal_to_f64(dt_sqrt).map_err(|_| VolatilityError::NumericalFailure {
        reason: "simulate_heston_volatility: sqrt(dt) not representable as f64".to_string(),
    })?;
    for _ in 1..steps {
        let dw = heston_wiener_increment(rng, dt_sqrt_f64)?;
        let sqrt_v = p_sqrt(&v_pos, "volatility::utils::simulate_heston_volatility")?.to_dec();
        // Euler step of dv = κ(θ − v)dt + ξ√v dW, every factor checked:
        // the raw operator form aborted with `Multiplication overflowed`
        // on an extreme `dt` or an extreme initial variance.
        let reversion = d_sub(theta, v, "volatility::heston::reversion")?;
        let drift = d_mul(
            d_mul(kappa, reversion, "volatility::heston::kappa_reversion")?,
            dt,
            "volatility::heston::drift",
        )?;
        let diffusion = d_mul(
            d_mul(xi, sqrt_v, "volatility::heston::xi_sqrt_v")?,
            dw,
            "volatility::heston::diffusion",
        )?;
        v = d_add(
            v,
            d_add(drift, diffusion, "volatility::heston::dv")?,
            "volatility::heston::variance",
        )?;
        v = v.max(Decimal::ZERO); // Ensure variance doesn't become negative
        v_pos = Positive::new_decimal(v)?;
        volatilities.push(p_sqrt(
            &v_pos,
            "volatility::utils::simulate_heston_volatility",
        )?);
    }
    Ok(volatilities)
}

/// Calculates bounds for uncertain volatility.
///
/// # Arguments
///
/// * `option` - The option for which to calculate bounds.
/// * `min_volatility` - The minimum possible volatility.
/// * `max_volatility` - The maximum possible volatility.
///
/// # Returns
///
/// A tuple `(lower_bound, upper_bound)` of Black–Scholes values of the
/// position, signed by `option.side` like every other price (#646, #648):
/// a short position's bounds are the long bounds negated and swapped,
/// `(-long_upper, -long_lower)`, so both are non-positive and still
/// ordered. The long bounds are the long prices at `min_volatility` and
/// `max_volatility`, in that order; a vanilla price rises with volatility,
/// so with `min_volatility <= max_volatility` they are ordered too.
///
/// # Errors
///
/// Propagates any [`VolatilityError::Options`] returned by the
/// underlying Black–Scholes evaluation at the minimum or maximum
/// volatility bound — typically
/// `OptionsError::PricingError` with an
/// `PricingError::ExpirationDate` inner cause.
pub fn uncertain_volatility_bounds(
    option: &Options,
    min_volatility: Positive,
    max_volatility: Positive,
) -> Result<(Decimal, Decimal), VolatilityError> {
    // Both bounds are priced long and the side is applied once (#715). A
    // short Black–Scholes price is negative, and the `Positive` bounds this
    // returned clamped it to zero.
    let long_price = |volatility: Positive| -> Result<Decimal, VolatilityError> {
        let mut long = option.clone();
        long.side = Side::Long;
        long.implied_volatility = volatility;
        Ok(long.calculate_price_black_scholes()?)
    };
    let lower = long_price(min_volatility)?;
    let upper = long_price(max_volatility)?;

    Ok(match option.side {
        Side::Long => (lower, upper),
        Side::Short => (-upper, -lower),
    })
}

/// Annualizes a volatility value from a specific timeframe.
///
/// # Arguments
///
/// * `volatility` - The volatility value to annualize
/// * `timeframe` - The timeframe of the input volatility
///
/// # Returns
///
/// The annualized volatility as Decimal
///
/// # Formula
///
/// The annualization is performed using the square root of time rule:
/// annualized_vol = vol * sqrt(periods_per_year)
///
/// # Examples
///
/// ```
/// use rust_decimal_macros::dec;
/// use optionstratlib_core::pos_or_panic;
/// use optionstratlib_core::utils::time::TimeFrame;
/// use optionstratlib_pricing::volatility::{annualized_volatility};
/// let daily_vol = pos_or_panic!(0.01); // 1% daily volatility
/// let annual_vol = annualized_volatility(daily_vol, TimeFrame::Day);
/// // annual_vol ≈ 0.1587 or about 15.87%
/// ```
///
/// # Errors
///
/// Returns [`VolatilityError::PositiveError`] when the scaling factor
/// multiplied by the base volatility cannot be represented as a
/// `Positive` — for example `Positive::MAX` annualised from
/// `TimeFrame::Microsecond`, which used to abort with
/// `Positive arithmetic overflow in mul`.
#[inline]
pub fn annualized_volatility(
    volatility: Positive,
    timeframe: TimeFrame,
) -> Result<Positive, VolatilityError> {
    let scale_factor = p_sqrt(
        &timeframe.periods_per_year(),
        "volatility::utils::annualized_volatility",
    )?;
    Ok(volatility.checked_mul(&scale_factor)?)
}

/// De-annualizes a volatility value to a specific timeframe.
///
/// # Arguments
///
/// * `annual_volatility` - The annualized volatility value
/// * `timeframe` - The target timeframe
///
/// # Returns
///
/// The de-annualized volatility as Decimal
///
/// # Formula
///
/// The de-annualization is performed using:
/// timeframe_vol = annual_vol / sqrt(periods_per_year)
///
/// # Examples
///
/// ```
/// use rust_decimal_macros::dec;
/// use optionstratlib_core::pos_or_panic;
/// use optionstratlib_core::utils::time::TimeFrame;
/// use optionstratlib_pricing::volatility::{de_annualized_volatility};
/// let annual_vol = pos_or_panic!(0.20); // 20% annual volatility
/// let daily_vol = de_annualized_volatility(annual_vol, TimeFrame::Day);
/// // daily_vol ≈ 0.0126 or about 1.26%
/// ```
///
/// # Errors
///
/// Returns [`VolatilityError::InvalidTime`] when the timeframe reports
/// zero periods per year — `TimeFrame::Custom(Positive::ZERO)` — since
/// the square-root-of-time rule has no de-annualised value there. The
/// guard and its message mirror the one [`adjust_volatility`] applies
/// to its target timeframe, so the two entry points agree on where the
/// domain ends.
///
/// Returns [`VolatilityError::PositiveError`] when the rescaling
/// produces a value that violates the `Positive` invariant or
/// overflows, typically due to division by an extremely small
/// number of periods.
#[inline]
pub fn de_annualized_volatility(
    annual_volatility: Positive,
    timeframe: TimeFrame,
) -> Result<Positive, VolatilityError> {
    let periods = timeframe.periods_per_year();
    if periods.is_zero() {
        return Err(VolatilityError::InvalidTime {
            time: periods,
            reason: format!(
                "Cannot de-annualize volatility to timeframe with zero periods per year: {timeframe:?}"
            ),
        });
    }
    let scale_factor = p_sqrt(&periods, "volatility::utils::de_annualized_volatility")?;
    // Explicit rounding rather than the dependency's default, matching the
    // strategy `d_div` applies crate-wide, so a de-annualised volatility
    // rounds the same way whichever path computes it.
    Ok(annual_volatility
        .checked_div_with_strategy(&scale_factor, RoundingStrategy::MidpointNearestEven)?)
}

/// Adjusts volatility between different timeframes using the square root of time rule
///
/// # Arguments
/// * `volatility` - The volatility to adjust
/// * `from_frame` - The original timeframe of the volatility
/// * `to_frame` - The target timeframe for the volatility
///
/// # Returns
/// The adjusted volatility for the target timeframe
///
/// # Example
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use optionstratlib_core::pos_or_panic;
/// use optionstratlib_core::utils::TimeFrame;
/// use optionstratlib_pricing::volatility::adjust_volatility;
/// let daily_vol = pos_or_panic!(0.2); // 20% daily volatility
/// let minute_vol = adjust_volatility(daily_vol, TimeFrame::Day, TimeFrame::Minute)?;
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns [`VolatilityError::InvalidTime`] when `to_frame` reports zero
/// periods per year, and [`VolatilityError::PositiveError`] when the
/// period ratio or the rescaled volatility leaves the representable
/// `Positive` range — for example `TimeFrame::Microsecond` into
/// `TimeFrame::Custom(1e-27)`, which used to abort with
/// `Positive arithmetic overflow in div`.
pub fn adjust_volatility(
    volatility: Positive,
    from_frame: TimeFrame,
    to_frame: TimeFrame,
) -> Result<Positive, VolatilityError> {
    if from_frame == to_frame {
        return Ok(volatility);
    }

    let from_periods = from_frame.periods_per_year();
    let to_periods = to_frame.periods_per_year();

    if from_periods == to_periods {
        return Ok(volatility);
    }

    // Check for division by zero
    if to_periods.is_zero() {
        return Err(VolatilityError::InvalidTime {
            time: to_periods,
            reason: format!(
                "Cannot adjust volatility to timeframe with zero periods per year: {to_frame:?}"
            ),
        });
    }

    // Scale factor is the square root of `from_periods / to_periods`, taken
    // as `sqrt(from) / sqrt(to)` rather than `sqrt(from / to)`. The two are
    // equal in exact arithmetic and not in this one: the ratio is rounded to
    // 28 decimal places before the root can rescue it, so a `Year` source
    // against a `Custom(Positive::MAX)` target rounds `1 / MAX` — about
    // `1.26e-29` — to zero and the volatility silently collapses, even
    // though the root of that ratio is around `1.1e-15` and perfectly
    // representable. Rooting first keeps both operands far from the scale
    // limit. The division that follows chooses `MidpointNearestEven`, the
    // strategy `d_div` applies to every `Decimal` division in the crate.
    let scale_factor = p_sqrt(&from_periods, "volatility::utils::adjust_volatility")?
        .checked_div_with_strategy(
            &p_sqrt(&to_periods, "volatility::utils::adjust_volatility")?,
            RoundingStrategy::MidpointNearestEven,
        )?;

    Ok(volatility.checked_mul(&scale_factor)?)
}

/// Adjusts annualized volatility for use in random walk simulations with a specific dt.
///
/// This function converts an annualized volatility to the appropriate volatility for the
/// base timeframe used in the simulation. The random walk implementation will then scale
/// this by sqrt(dt) internally, so this function only performs the de-annualization step.
///
/// **Important**: The random walk implementations (Brownian, GeometricBrownian, etc.) already
/// multiply the volatility by `sqrt(dt)` internally, so this function should NOT do that scaling.
/// It only converts from annual volatility to the base timeframe volatility.
///
/// # Arguments
///
/// * `annual_volatility` - The annualized volatility (e.g., 0.20 for 20% annual volatility)
/// * `_dt` - The time step size (unused, kept for API compatibility)
/// * `_dt_timeframe` - The timeframe that dt represents (unused, kept for API compatibility)
/// * `dt_base_timeframe` - The base timeframe for the simulation (e.g., `TimeFrame::Day`)
///
/// # Returns
///
/// The de-annualized volatility for the base timeframe as `Positive`
///
/// # Formula
///
/// The adjustment follows:
/// ```text
/// volatility_base = annual_volatility / sqrt(periods_per_year)
/// ```
///
/// The random walk will then apply: `volatility_base * sqrt(dt)` internally.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use optionstratlib_core::{pos_or_panic, model::Positive};
/// use optionstratlib_core::utils::time::{TimeFrame, convert_time_frame};
/// use optionstratlib_pricing::volatility::volatility_for_dt;
///
/// let annual_vol = pos_or_panic!(0.20); // 20% annual volatility
/// let days = pos_or_panic!(7.0);
/// let dt = convert_time_frame(Positive::ONE / days, &TimeFrame::Minute, &TimeFrame::Day)?;
///
/// // Get volatility for daily timeframe (random walk will scale by sqrt(dt))
/// let vol_for_walk = volatility_for_dt(
///     annual_vol,
///     dt,
///     TimeFrame::Minute,
///     TimeFrame::Day
/// )?;
/// # Ok(())
/// # }
/// ```
///
/// # Errors
///
/// Returns an error if the de-annualization fails
#[inline]
pub fn volatility_for_dt(
    annual_volatility: Positive,
    _dt: Positive,
    _dt_timeframe: TimeFrame,
    dt_base_timeframe: TimeFrame,
) -> Result<Positive, VolatilityError> {
    // Only de-annualize to the base timeframe
    // The random walk implementation will multiply by sqrt(dt) internally
    de_annualized_volatility(annual_volatility, dt_base_timeframe)
}

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

    /// A representable scale factor must not collapse to zero on its way to
    /// being computed. `Year` against `Custom(Positive::MAX)` divides one by
    /// `Decimal::MAX`, which is about `1.26e-29` and rounds to zero at the
    /// 28-place scale — but the square root of that ratio is around
    /// `1.1e-15`, which is representable, so the answer exists and the
    /// intermediate was the only thing that did not.
    #[test]
    fn test_adjust_volatility_survives_an_extreme_target_timeframe() {
        let adjusted = adjust_volatility(
            pos_or_panic!(0.2),
            TimeFrame::Year,
            TimeFrame::Custom(Positive::MAX),
        )
        .expect("the conversion is representable");

        assert!(
            adjusted > Positive::ZERO,
            "a representable scale factor collapsed to zero: {adjusted}"
        );
    }
    use optionstratlib_core::assert_pos_relative_eq;

    #[test]
    fn test_annualize_daily_volatility() {
        let daily_vol = pos_or_panic!(0.01); // 1% daily volatility
        let annual_vol = annualized_volatility(daily_vol, TimeFrame::Day).unwrap();
        let expected = daily_vol * pos_or_panic!(252.0f64.sqrt());
        assert_pos_relative_eq!(annual_vol, expected, pos_or_panic!(1e-10));
    }

    #[test]
    fn test_deannualize_annual_volatility() {
        let annual_vol = pos_or_panic!(0.20); // 20% annual volatility
        let daily_vol = de_annualized_volatility(annual_vol, TimeFrame::Day).unwrap();
        let expected = pos_or_panic!(0.01259881576697424);
        assert_pos_relative_eq!(daily_vol, expected, pos_or_panic!(1e-10));
    }

    #[test]
    fn test_custom_timeframe() {
        let custom_periods = Positive::HUNDRED;
        let vol = pos_or_panic!(0.05);
        let annual_vol = annualized_volatility(vol, TimeFrame::Custom(custom_periods)).unwrap();
        let expected = vol * custom_periods.sqrt();
        assert_pos_relative_eq!(annual_vol, expected, pos_or_panic!(1e-10));
    }

    #[test]
    fn test_conversion_roundtrip() {
        let original_vol = pos_or_panic!(0.15);
        let annualized = annualized_volatility(original_vol, TimeFrame::Day).unwrap();
        let roundtrip = de_annualized_volatility(annualized, TimeFrame::Day).unwrap();
        assert_pos_relative_eq!(original_vol, roundtrip, pos_or_panic!(1e-10));
    }

    #[test]
    fn test_different_timeframes() {
        let daily_vol = pos_or_panic!(0.01);
        let weekly_vol = annualized_volatility(daily_vol, TimeFrame::Day).unwrap();
        let monthly_vol = de_annualized_volatility(weekly_vol, TimeFrame::Month).unwrap();
        assert!(monthly_vol > daily_vol); // Monthly vol should be higher than daily
    }
}

#[cfg(test)]
mod tests_constant_volatility {
    use super::*;
    use optionstratlib_core::assert_pos_relative_eq;
    use optionstratlib_core::constants::ZERO;
    use rust_decimal_macros::dec;

    #[test]
    fn test_constant_volatility_single_value() {
        let returns = [dec!(0.05)];
        let result = constant_volatility(&returns).unwrap();
        assert_eq!(result, ZERO);
    }

    #[test]
    fn test_constant_volatility_identical_values() {
        let returns = [dec!(0.02), dec!(0.02), dec!(0.02), dec!(0.02)];
        let result = constant_volatility(&returns).unwrap();
        assert_eq!(result, ZERO);
    }

    #[test]
    fn test_constant_volatility_varying_values() {
        let returns = [dec!(0.01), dec!(0.03), dec!(0.02), dec!(0.04)];
        let result = constant_volatility(&returns).unwrap();
        assert_pos_relative_eq!(
            result,
            pos_or_panic!(0.012909944487358056),
            pos_or_panic!(1e-10)
        );
    }
}

#[cfg(test)]
mod tests_historical_volatility {
    use super::*;
    use optionstratlib_core::assert_pos_relative_eq;
    use rust_decimal_macros::dec;

    #[test]
    fn test_historical_volatility_empty_returns() {
        let returns: [Decimal; 0] = [];
        let result = historical_volatility(&returns, 3).unwrap();
        assert!(result.is_empty());
    }

    #[test]
    fn test_historical_volatility_single_value() {
        let returns = [dec!(0.02)];
        let result = historical_volatility(&returns, 3).unwrap();
        assert!(result.is_empty());
    }

    #[test]
    fn test_historical_volatility_insufficient_data() {
        let returns = [dec!(0.01), dec!(0.02)];
        let result = historical_volatility(&returns, 3).unwrap();
        assert!(result.is_empty());
    }

    #[test]
    fn test_historical_volatility_exact_window() {
        let returns = [dec!(0.01), dec!(0.02), dec!(0.03)];
        let result = historical_volatility(&returns, 3).unwrap();
        assert_eq!(result.len(), 1);
        assert_pos_relative_eq!(result[0], pos_or_panic!(0.01), pos_or_panic!(1e-10));
    }

    #[test]
    fn test_historical_volatility_larger_window() {
        let returns = [dec!(0.01), dec!(0.02), dec!(0.03), dec!(0.04)];
        let result = historical_volatility(&returns, 3).unwrap();
        assert_eq!(result.len(), 2);
        assert_pos_relative_eq!(result[0], pos_or_panic!(0.01), pos_or_panic!(1e-10));
        assert_pos_relative_eq!(result[1], pos_or_panic!(0.01), pos_or_panic!(1e-10));
    }
}

#[cfg(test)]
mod tests_ewma_volatility {
    use super::*;
    use optionstratlib_core::assert_pos_relative_eq;
    use rust_decimal::MathematicalOps;
    use rust_decimal_macros::dec;

    #[test]
    fn test_ewma_volatility_single_return() {
        let returns = [dec!(0.02)];
        let lambda = dec!(0.94);
        let result = ewma_volatility(&returns, lambda).unwrap();
        assert_eq!(result.len(), 1);
        assert_eq!(result[0], pos_or_panic!(0.02)); // The volatility is simply the return itself
    }

    #[test]
    fn test_ewma_volatility_constant_returns() {
        let returns = [dec!(0.02), dec!(0.02), dec!(0.02), dec!(0.02)];
        let lambda = dec!(0.94);
        let single_value = dec!(0.02);
        let result = ewma_volatility(&returns, lambda).unwrap();
        assert_eq!(result.len(), 4);

        // Test the EWMA calculation
        let expected = [
            pos_or_panic!(0.02),
            Positive::try_from(
                (lambda * single_value.powi(2) + (Decimal::ONE - lambda) * single_value.powi(2))
                    .sqrt()
                    .unwrap(),
            )
            .unwrap_or_default(),
            Positive::try_from(
                (lambda * single_value.powi(2) + (Decimal::ONE - lambda) * single_value.powi(2))
                    .sqrt()
                    .unwrap()
                    .powi(2),
            )
            .unwrap_or_default()
            .sqrt(),
            Positive::try_from(
                (lambda * single_value.powi(2) + (Decimal::ONE - lambda) * single_value.powi(2))
                    .sqrt()
                    .unwrap()
                    .powi(2)
                    .sqrt()
                    .unwrap()
                    .powi(2),
            )
            .unwrap_or_default()
            .sqrt(),
        ];

        for (res, exp) in result.iter().zip(expected.iter()) {
            assert_pos_relative_eq!(*res, *exp, pos_or_panic!(1e-10));
        }
    }

    #[test]
    fn test_ewma_volatility_varying_returns() {
        let returns = [dec!(0.01), dec!(0.02), dec!(0.03), dec!(0.04)];
        let lambda = dec!(0.94);
        let result = ewma_volatility(&returns, lambda).unwrap();
        assert_eq!(result.len(), 4);

        // Verify general properties rather than exact values due to complexity of calculation
        assert!(result.iter().all(|&x| x > Positive::ZERO));
        assert!(
            result
                .iter()
                .zip(result.iter().skip(1))
                .all(|(a, b)| *b >= *a)
        );
    }

    #[test]
    fn test_ewma_volatility_low_lambda() {
        let returns = [dec!(0.01), dec!(0.02), dec!(0.03), dec!(0.04)];
        let lambda = dec!(0.5); // Low lambda means faster decay
        let result = ewma_volatility(&returns, lambda).unwrap();
        assert_eq!(result.len(), 4);

        // Verify more weight on recent values
        let last = result.last().unwrap();
        assert!(*last > *result.first().unwrap());
    }

    #[test]
    fn test_ewma_volatility_high_lambda() {
        let returns = [dec!(0.01), dec!(0.02), dec!(0.03), dec!(0.04)];
        let lambda = dec!(0.99); // High lambda means slower decay
        let result = ewma_volatility(&returns, lambda).unwrap();
        assert_eq!(result.len(), 4);

        // Verify less weight on recent values
        let differences: Vec<_> = result.windows(2).map(|w| w[1] - w[0]).collect();
        assert!(differences.iter().all(|&d| d.to_dec() < dec!(0.01)));
    }
}

#[cfg(test)]
mod tests_implied_volatility {
    use super::*;
    use crate::greeks::Greeks;
    use optionstratlib_core::assert_decimal_eq;
    use optionstratlib_core::assert_pos_relative_eq;
    use optionstratlib_core::constants::{MAX_VOLATILITY, MIN_VOLATILITY};
    use optionstratlib_core::model::ExpirationDate;
    use optionstratlib_core::model::types::{OptionStyle, OptionType, Side};
    use rust_decimal_macros::dec;
    use tracing::error;

    fn create_test_option() -> Options {
        Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            Positive::HUNDRED,                         // strike price
            ExpirationDate::Days(pos_or_panic!(30.0)), // 30 days to expiration
            pos_or_panic!(0.2),                        // initial implied volatility
            Positive::ONE,                             // quantity
            Positive::HUNDRED,                         // underlying price (ATM)
            dec!(0.05),                                // risk-free rate
            OptionStyle::Call,                         // call option
            Positive::ZERO,                            // no dividend yield
            None,                                      // no exotic params
        )
    }

    #[test]
    fn test_implied_volatility_max_iterations() {
        let mut option = create_test_option();
        let market_price = pos_or_panic!(5.0);

        // One solver step does not converge. The grid search replaced in
        // #859 P4 read `1` as 100 candidates and always returned one; the
        // solver reports that it did not find the root instead.
        let result = implied_volatility(market_price, &mut option, 1);
        assert!(
            matches!(result, Err(VolatilityError::IvNotFound)),
            "{result:?}"
        );

        let result = implied_volatility(market_price, &mut option, 100);
        let iv = result.unwrap();
        assert!(iv >= *MIN_VOLATILITY && iv <= MAX_VOLATILITY);

        let result = calculate_iv(
            market_price,
            Positive::HUNDRED,
            OptionStyle::Call,
            Positive::HUNDRED,
            pos_or_panic!(30.0),
            "TEST".to_string(),
        );
        assert!(result.is_ok());

        let iv = result.unwrap();
        assert!(iv >= *MIN_VOLATILITY && iv <= MAX_VOLATILITY);
        // The grid search replaced in #859 P4 returned the grid point 0.437;
        // the solver returns the volatility that reprices 5.0.
        assert_pos_relative_eq!(iv, pos_or_panic!(0.4374518792594474), pos_or_panic!(1e-12));
    }

    #[test]
    fn test_implied_volatility_zero_dte() {
        let iv = pos_or_panic!(0.25);
        let mut option = Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            Positive::HUNDRED,
            ExpirationDate::Days(pos_or_panic!(0.5)),
            iv,
            Positive::ONE,
            Positive::HUNDRED,
            dec!(0.0),
            OptionStyle::Call,
            Positive::ZERO,
            None,
        );

        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(0.502), dec!(0.002));
        assert_decimal_eq!(gamma, dec!(0.431), dec!(0.002));
        assert_decimal_eq!(vega, dec!(0.015), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-0.369), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.001), dec!(0.002));
        assert_decimal_eq!(vanna, dec!(0.0073826791123684657342819348), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(-0.0000012641), dec!(0.0000000001));
        assert_decimal_eq!(veta, dec!(0.0002138), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(-0.0018456), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(-0.4311576), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(0.369), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001)); // More lenient tolerance
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }

        option.implied_volatility = iv;
        option.option_style = OptionStyle::Put;
        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(-0.498), dec!(0.002));
        assert_decimal_eq!(gamma, dec!(0.431), dec!(0.002));
        assert_decimal_eq!(vega, dec!(0.015), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-0.369), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.001), dec!(0.002));
        assert_decimal_eq!(vanna, dec!(0.0073826791123684657342819348), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(-0.0000012641), dec!(0.0000000001));
        assert_decimal_eq!(veta, dec!(0.0002138), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(-0.0018456), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(-0.4311576), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(0.369), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001)); // More lenient tolerance
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }
    }

    #[test]
    fn test_implied_volatility_zero_dte_real() {
        let iv = pos_or_panic!(0.356831);
        let mut option = Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            pos_or_panic!(20600.0),
            ExpirationDate::Days(pos_or_panic!(0.52)),
            iv,
            Positive::ONE,
            pos_or_panic!(21049.88),
            dec!(0.0),
            OptionStyle::Call,
            pos_or_panic!(0.05),
            None,
        );

        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(0.946), dec!(0.002));
        assert_decimal_eq!(gamma, dec!(0.00038), dec!(0.002));
        assert_decimal_eq!(vega, dec!(0.8735285922867751407644169906), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-27.244276204109831381546317238), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.277), dec!(0.002));
        assert_decimal_eq!(
            vanna,
            dec!(-0.4905170891818293972168486173),
            dec!(0.0000001)
        );
        assert_decimal_eq!(
            vomma,
            dec!(6.2569425433351718484053752895),
            dec!(0.00000001)
        );
        assert_decimal_eq!(veta, dec!(0.0434650505282011247018063025), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1695474991264431915978503078), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005864547013714062753140932), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(454.917), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }

        option.implied_volatility = iv;
        option.option_style = OptionStyle::Put;
        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(-0.0541909521971111839091560979), dec!(0.001));
        assert_decimal_eq!(gamma, dec!(0.0), dec!(0.001));
        assert_decimal_eq!(vega, dec!(0.8735285922867751407644169906), dec!(0.001));
        assert_decimal_eq!(theta, dec!(-30.127616013684725380529216795), dec!(0.001));
        assert_decimal_eq!(rho, dec!(-0.016), dec!(0.001));
        assert_decimal_eq!(
            vanna,
            dec!(-0.4905170891818293972168486173),
            dec!(0.0000001)
        );
        assert_decimal_eq!(
            vomma,
            dec!(6.2569425433351718484053752895),
            dec!(0.00000001)
        );
        assert_decimal_eq!(veta, dec!(0.0434650505282011247018063025), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1694105225826541109000514178), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005864547013714062753140932), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(6.537), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }
    }

    #[test]
    fn test_implied_volatility_zero_dte_real_short() {
        let iv = pos_or_panic!(0.356831);
        let mut option = Options::new(
            OptionType::European,
            Side::Short,
            "TEST".to_string(),
            pos_or_panic!(20600.0),
            ExpirationDate::Days(pos_or_panic!(0.52)),
            iv,
            Positive::ONE,
            pos_or_panic!(21049.88),
            dec!(0.0),
            OptionStyle::Call,
            pos_or_panic!(0.05),
            None,
        );

        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        // Short position, so every greek carries the negative of the equivalent
        // long one; see #428. Theta is positive because the short collects decay.
        assert_decimal_eq!(delta, dec!(-0.946), dec!(0.002));
        assert_decimal_eq!(gamma, dec!(-0.00038), dec!(0.002));
        assert_decimal_eq!(vega, dec!(-0.8735285922867751407644169906), dec!(0.002));
        assert_decimal_eq!(theta, dec!(27.244276204109831381546317238), dec!(0.002));
        assert_decimal_eq!(rho, dec!(-0.277), dec!(0.002));
        assert_decimal_eq!(vanna, dec!(0.4905170891818293972168486173), dec!(0.0000001));
        assert_decimal_eq!(
            vomma,
            dec!(-6.2569425433351718484053752895),
            dec!(0.00000001)
        );
        assert_decimal_eq!(veta, dec!(-0.0434650505282011247018063025), dec!(0.0000001));
        assert_decimal_eq!(
            charm,
            dec!(-0.1695474991264431915978503078),
            dec!(0.0000001)
        );
        assert_decimal_eq!(
            color,
            dec!(-0.0005864547013714062753140932),
            dec!(0.0000001)
        );

        let market_price = option.calculate_price_black_scholes().unwrap().abs();
        assert_decimal_eq!(market_price, dec!(454.917), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }

        option.implied_volatility = iv;
        option.option_style = OptionStyle::Put;
        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        // Still the same short position, now as a put; see #428.
        assert_decimal_eq!(delta, dec!(0.0541909521971111839091560979), dec!(0.001));
        assert_decimal_eq!(gamma, dec!(0.0), dec!(0.001));
        assert_decimal_eq!(vega, dec!(-0.8735285922867751407644169906), dec!(0.001));
        assert_decimal_eq!(theta, dec!(30.127616013684725380529216795), dec!(0.001));
        assert_decimal_eq!(rho, dec!(0.016), dec!(0.001));
        assert_decimal_eq!(vanna, dec!(0.4905170891818293972168486173), dec!(0.0000001));
        assert_decimal_eq!(
            vomma,
            dec!(-6.2569425433351718484053752895),
            dec!(0.00000001)
        );
        assert_decimal_eq!(veta, dec!(-0.0434650505282011247018063025), dec!(0.0000001));
        assert_decimal_eq!(
            charm,
            dec!(-0.1694105225826541109000514178),
            dec!(0.0000001)
        );
        assert_decimal_eq!(
            color,
            dec!(-0.0005864547013714062753140932),
            dec!(0.0000001)
        );

        let market_price = option.calculate_price_black_scholes().unwrap().abs();
        assert_decimal_eq!(market_price, dec!(6.537), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }
    }

    #[test]
    fn test_implied_volatility_zero_dte_real_put() {
        let iv = pos_or_panic!(0.356831);
        let mut option = Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            pos_or_panic!(23325.0),
            ExpirationDate::Days(pos_or_panic!(0.29)),
            iv,
            Positive::ONE,
            pos_or_panic!(24118.5),
            dec!(0.0),
            OptionStyle::Put,
            Positive::ZERO,
            None,
        );

        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(-0.00043), dec!(0.002));
        assert_decimal_eq!(gamma, dec!(0.0000064), dec!(0.00001));
        assert_decimal_eq!(vega, dec!(0.0105), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-0.64997), dec!(0.002));
        assert_decimal_eq!(rho, dec!(-0.000083), dec!(0.002));
        assert_decimal_eq!(
            vanna,
            dec!(-0.0144632755644302251546856722),
            dec!(0.0000001)
        );
        assert_decimal_eq!(vomma, dec!(0.3275301270), dec!(0.0000000001));
        assert_decimal_eq!(veta, dec!(0.0031824), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.0088981), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0001111), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(0.027), dec!(0.002)); // too small to be accurate

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(pos_or_panic!(1.5), &mut option, 10);

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, pos_or_panic!(0.528), pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }
    }

    #[test]
    fn test_implied_volatility_zero_dte_real_call() {
        let iv = pos_or_panic!(0.356831);
        let mut option = Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            pos_or_panic!(23325.0),
            ExpirationDate::Days(pos_or_panic!(0.32)),
            iv,
            Positive::ONE,
            pos_or_panic!(24103.00),
            dec!(0.0),
            OptionStyle::Call,
            Positive::ZERO,
            None,
        );

        let delta = option.delta().unwrap();
        let gamma = option.gamma().unwrap();
        let vega = option.vega().unwrap();
        let theta = option.theta().unwrap();
        let rho = option.rho().unwrap();
        let vanna = option.vanna().unwrap();
        let vomma = option.vomma().unwrap();
        let veta = option.veta().unwrap();
        let charm = option.charm().unwrap();
        let color = option.color().unwrap();
        assert_decimal_eq!(delta, dec!(0.999), dec!(0.001));
        assert_decimal_eq!(gamma, dec!(0.00001), dec!(0.00001));
        assert_decimal_eq!(vega, dec!(0.02255), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-1.25733), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.204), dec!(0.002));
        assert_decimal_eq!(
            vanna,
            dec!(-0.0274529936439301097044284653),
            dec!(0.0000001)
        );
        assert_decimal_eq!(vomma, dec!(0.6094669204), dec!(0.0000000001));
        assert_decimal_eq!(veta, dec!(0.0054321), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.0153063), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0001675), dec!(0.0000001));

        let market_price = option.calculate_price_black_scholes().unwrap();
        assert_decimal_eq!(market_price, dec!(778.065), dec!(0.002));

        // For very short-term options, use a more lenient tolerance
        option.implied_volatility = pos_or_panic!(0.4); // Start with different IV
        let iv_result = implied_volatility(
            Positive::new_decimal(market_price).unwrap(),
            &mut option,
            10,
        );

        // For zero DTE options, expect either convergence or a reasonable approximation
        match iv_result {
            Ok(iv_aprox) => {
                // If it converges, it should be close to the original IV
                assert_pos_relative_eq!(iv_aprox, iv, pos_or_panic!(0.001));
            }
            Err(_) => {
                // If it doesn't converge, that's also acceptable for zero DTE options
                // as they have extremely low vega and are numerically challenging
                error!("Non-convergence is acceptable for zero DTE options");
            }
        }
    }
}

#[cfg(test)]
mod tests_garch_volatility {
    use super::*;
    use optionstratlib_core::assert_decimal_eq;
    use optionstratlib_core::assert_pos_relative_eq;
    use rust_decimal_macros::dec;

    #[test]
    fn test_garch_single_return() {
        let returns = vec![dec!(0.02)];
        let omega = dec!(0.1);
        let alpha = dec!(0.2);
        let beta = dec!(0.7);

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();
        assert_eq!(result.len(), 1);
        assert_pos_relative_eq!(result[0], pos_or_panic!(0.02), pos_or_panic!(1e-10));
    }

    #[test]
    fn test_garch_constant_returns() {
        let returns = vec![dec!(0.02), dec!(0.02), dec!(0.02), dec!(0.02)];
        let omega = dec!(0.1);
        let alpha = dec!(0.2);
        let beta = dec!(0.7);

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();
        assert_eq!(result.len(), 4);

        // First volatility should be the absolute value of first return
        assert_pos_relative_eq!(result[0], pos_or_panic!(0.02), pos_or_panic!(1e-10));

        // For constant returns, volatility should converge
        let last_two_diff = (result[3].to_dec() - result[2].to_dec()).abs();
        assert_decimal_eq!(last_two_diff, dec!(0.0555), dec!(0.001));
    }

    #[test]
    fn test_garch_varying_returns() {
        let returns = vec![dec!(0.01), dec!(-0.02), dec!(0.03), dec!(-0.01)];
        let omega = dec!(0.1);
        let alpha = dec!(0.2);
        let beta = dec!(0.7);

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();
        assert_eq!(result.len(), 4);

        // All volatilities should be positive
        assert!(result.iter().all(|&v| v > Positive::ZERO));

        // Volatility should increase after large returns
        assert!(result[2] > result[1]); // After 0.03 return
    }

    #[test]
    fn test_garch_zero_initial_return() {
        let returns = vec![dec!(0.0), dec!(0.02), dec!(0.03)];
        let omega = dec!(0.1);
        let alpha = dec!(0.2);
        let beta = dec!(0.7);

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();
        assert_eq!(result.len(), 3);
        assert_eq!(result[0], Positive::ZERO);
    }

    #[test]
    fn test_garch_parameter_sum_one() {
        let returns = vec![dec!(0.01), dec!(0.02), dec!(0.03)];
        let omega = dec!(0.05);
        let alpha = dec!(0.15);
        let beta = dec!(0.8); // alpha + beta = 0.95 < 1

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();

        // Volatility should remain stable when alpha + beta < 1
        assert!(result.iter().all(|&v| v < Positive::ONE));
    }

    #[test]
    fn test_garch_extreme_returns() {
        let returns = vec![dec!(0.01), dec!(0.2), dec!(-0.2), dec!(0.01)]; // Large returns
        let omega = dec!(0.1);
        let alpha = dec!(0.2);
        let beta = dec!(0.7);

        let result = garch_volatility(&returns, omega, alpha, beta).unwrap();

        // Volatility should spike after large returns
        assert!(result[2] > result[1]);
        assert!(result[2] > result[0] * Positive::TWO); // At least double the initial volatility
    }

    #[test]
    fn test_garch_parameters_validation() {
        let returns = vec![dec!(0.01)];

        // Test with invalid parameters (negative values)
        let result_negative = garch_volatility(&returns, dec!(-0.1), dec!(0.2), dec!(0.7));
        assert!(result_negative.is_ok()); // Should handle negative parameters

        // Test with parameters summing to more than 1
        let result_sum = garch_volatility(&returns, dec!(0.1), dec!(0.5), dec!(0.6));
        assert!(result_sum.is_ok()); // Should handle parameters summing > 1
    }
}

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

    /// Last volatility of `seeded_path(DETERMINISTIC_RNG_DEFAULT_SEED)`.
    /// Re-baselined by #742, when `dW` became a normal draw; the uniform
    /// draw it replaced pinned 0.4720671070320042200460741098, a variance
    /// of 0.22 against a long-run `theta` of 0.04.
    const PINNED_HESTON_LAST: Decimal = dec!(0.2326994513761414494617655821);
    use optionstratlib_core::utils::{DETERMINISTIC_RNG_DEFAULT_SEED, deterministic_rng};
    use rust_decimal_macros::dec;

    #[test]
    fn test_heston_basic_properties() {
        let kappa = dec!(2.0); // Mean reversion speed
        let theta = dec!(0.04); // Long-term variance
        let xi = dec!(0.3); // Volatility of variance
        let v0 = dec!(0.04); // Initial variance
        let dt = dec!(0.01); // Time step
        let steps = 100; // Number of steps

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Check basic properties
        assert_eq!(result.len(), steps);
        assert_eq!(result[0], Positive::new_decimal(v0).unwrap().sqrt());
        assert!(result.iter().all(|&x| x >= Positive::ZERO));
    }

    #[test]
    fn test_heston_zero_volatility_of_volatility() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(0.0); // No randomness in volatility
        let v0 = dec!(0.04);
        let dt = dec!(0.01);
        let steps = 50;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // With no volatility of volatility, values should move deterministically towards theta
        for i in 1..steps {
            // Verify that volatility moves monotonically towards the long-term mean
            match v0 < theta {
                true => assert!(result[i] >= result[i - 1]),
                false => assert!(result[i] <= result[i - 1]),
            }
        }
    }

    #[test]
    fn test_heston_high_mean_reversion() {
        let kappa = dec!(10.0); // High mean reversion
        let theta = dec!(0.04);
        let xi = dec!(0.1);
        let v0 = dec!(0.08); // Start away from theta
        let dt = dec!(0.01);
        let steps = 200;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // With high mean reversion, later values should be closer to sqrt(theta)
        let theta_vol = Positive::new_decimal(theta).unwrap().sqrt();
        let initial_dist = (result[0].to_dec() - theta_vol.to_dec()).abs();
        let final_dist = (result[steps - 1].to_dec() - theta_vol.to_dec()).abs();
        assert!(final_dist < initial_dist);
    }

    #[test]
    fn test_heston_high_volatility() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(1.0); // High volatility of volatility
        let v0 = dec!(0.04);
        let dt = dec!(0.01);
        let steps = 100;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // With high volatility, we should see more variation
        let variations: Vec<_> = result
            .windows(2)
            .map(|w| (w[1].to_dec() - w[0].to_dec()).abs())
            .collect();
        assert!(variations.iter().any(|&x| x > dec!(0.001)));
    }

    #[test]
    fn test_heston_zero_initial_variance() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(0.3);
        let v0 = dec!(0.0); // Start at zero
        let dt = dec!(0.01);
        let steps = 100;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Should start at zero and then increase
        assert_eq!(result[0], Positive::ZERO);
        assert!(result.iter().skip(1).any(|&x| x > Positive::ZERO));
    }

    #[test]
    fn test_heston_numerical_stability() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(0.3);
        let v0 = dec!(0.04);
        let dt = dec!(0.001); // Very small time step
        let steps = 1000;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Check that values remain finite and positive
        assert!(result.iter().all(|&x| x < Positive::HUNDRED));
        assert!(result.iter().all(|&x| x >= Positive::ZERO));
    }

    #[test]
    fn test_heston_large_time_steps() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(0.3);
        let v0 = dec!(0.04);
        let dt = dec!(0.1); // Large time step
        let steps = 10;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Should still produce valid results with large steps
        assert!(result.iter().all(|&x| x >= Positive::ZERO));
    }

    #[test]
    fn test_heston_extreme_parameters() {
        let kappa = dec!(20.0); // Very high mean reversion
        let theta = dec!(1.0); // High target variance
        let xi = dec!(2.0); // Very high vol of vol
        let v0 = dec!(0.5);
        let dt = dec!(0.01);
        let steps = 100;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Even with extreme parameters, results should be finite and non-negative
        assert!(result.iter().all(|&x| x >= Positive::ZERO));
    }

    #[test]
    fn test_heston_min_steps() {
        let kappa = dec!(2.0);
        let theta = dec!(0.04);
        let xi = dec!(0.3);
        let v0 = dec!(0.04);
        let dt = dec!(0.01);
        let steps = 1;

        let result = simulate_heston_volatility(
            kappa,
            theta,
            xi,
            v0,
            dt,
            steps,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();

        // Should work with minimum number of steps
        assert_eq!(result.len(), steps);
        assert_eq!(result[0], Positive::new_decimal(v0).unwrap().sqrt());
    }

    fn seeded_path(seed: u64) -> Vec<Positive> {
        simulate_heston_volatility(
            dec!(2.0),
            dec!(0.04),
            dec!(0.3),
            dec!(0.04),
            dec!(0.01),
            50,
            &mut deterministic_rng(seed),
        )
        .unwrap()
    }

    #[test]
    fn test_simulate_heston_volatility_same_seed_identical_path() {
        assert_eq!(seeded_path(11), seeded_path(11));
        assert_ne!(seeded_path(11), seeded_path(12));
    }

    #[test]
    fn test_simulate_heston_volatility_seeded_regression_pinned() {
        let path = seeded_path(DETERMINISTIC_RNG_DEFAULT_SEED);
        assert_eq!(path.len(), 50);
        assert_eq!(path.last().map(Positive::to_dec), Some(PINNED_HESTON_LAST));
    }

    #[test]
    fn test_heston_wiener_increment_seeded_moments_match_n0_dt() {
        // #742: dW must be N(0, dt). The uniform draw it replaced had mean
        // sqrt(dt) / 2 = 0.05 and variance dt / 12.
        const DRAWS: u32 = 100_000;
        let dt = 0.01_f64;
        let mut rng = deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED);
        let (mut sum, mut sum_sq) = (0.0_f64, 0.0_f64);
        for _ in 0..DRAWS {
            let dw = heston_wiener_increment(&mut rng, dt.sqrt())
                .unwrap()
                .to_f64()
                .unwrap();
            sum += dw;
            sum_sq += dw * dw;
        }
        let n = f64::from(DRAWS);
        let mean = sum / n;
        let variance = sum_sq / n - mean * mean;
        // Standard error of the mean is sqrt(dt / n) ~ 3.2e-4; of the
        // variance, dt * sqrt(2 / n) ~ 4.5e-5. Both bounds are ~5 sigma.
        assert!(mean.abs() < 1.6e-3, "sample mean {mean}");
        assert!((variance - dt).abs() < 2.5e-4, "sample variance {variance}");
    }

    #[test]
    fn test_simulate_heston_volatility_long_run_mean_reverts_to_theta() {
        // #742: with a symmetric dW the variance mean-reverts to theta. The
        // uniform draw biased every shock upward and settled near 0.64.
        // Feller: 2 * kappa * theta = 0.16 > xi^2 = 0.09. The normal draw
        // gives 0.0406 over these seeds.
        const PATHS: u64 = 200;
        const STEPS: usize = 1_000;
        const BURN_IN: usize = 500;
        let theta = dec!(0.04);
        let mut total = Decimal::ZERO;
        let mut count = Decimal::ZERO;
        for seed in 0..PATHS {
            let path = simulate_heston_volatility(
                dec!(2.0),
                theta,
                dec!(0.3),
                dec!(0.04),
                dec!(0.01),
                STEPS,
                &mut deterministic_rng(seed),
            )
            .unwrap();
            for vol in path.iter().skip(BURN_IN) {
                total += vol.to_dec() * vol.to_dec();
                count += Decimal::ONE;
            }
        }
        let long_run_mean = total / count;
        assert!(
            (long_run_mean - theta).abs() < dec!(0.004),
            "long-run mean variance {long_run_mean}, theta {theta}"
        );
    }

    #[test]
    fn test_simulate_heston_volatility_zero_dt_ignores_rng_draws() {
        // With dt = 0 the shock sqrt(dt) * z and the drift are zero, so the
        // path stays at sqrt(v0) whatever the generator yields.
        let path = simulate_heston_volatility(
            dec!(2.0),
            dec!(0.04),
            dec!(0.3),
            dec!(0.09),
            Decimal::ZERO,
            20,
            &mut deterministic_rng(DETERMINISTIC_RNG_DEFAULT_SEED),
        )
        .unwrap();
        assert!(path.iter().all(|v| *v == pos_or_panic!(0.3)));
    }
}

#[cfg(test)]
mod tests_uncertain_volatility_bounds {
    use super::*;
    use optionstratlib_core::assert_decimal_eq;
    use optionstratlib_core::model::types::{OptionStyle, OptionType, Side};

    use optionstratlib_core::model::ExpirationDate;
    use rust_decimal_macros::dec;

    // Helper function to create a test option
    fn create_test_option(style: OptionStyle, side: Side, strike: Positive) -> Options {
        Options::new(
            OptionType::European,
            side,
            "TEST".to_string(),
            strike,
            ExpirationDate::Days(pos_or_panic!(30.0)),
            pos_or_panic!(0.2), // Initial implied volatility
            Positive::ONE,      // Quantity
            Positive::HUNDRED,  // Underlying price
            dec!(0.05),         // Risk-free rate
            style,
            Positive::ZERO, // No dividend yield
            None,           // No exotic params
        )
    }

    #[test]
    fn test_bounds_basic_call() {
        let option = create_test_option(OptionStyle::Call, Side::Long, Positive::HUNDRED);
        let (lower, upper) =
            uncertain_volatility_bounds(&option, pos_or_panic!(0.1), pos_or_panic!(0.3)).unwrap();

        // Lower bound should be less than upper bound for a call
        assert!(lower < upper);
        // Both bounds should be positive
        assert!(lower > Decimal::ZERO);
        assert!(upper > Decimal::ZERO);
    }

    #[test]
    fn test_bounds_basic_put() {
        let option = create_test_option(OptionStyle::Put, Side::Long, Positive::HUNDRED);
        let (lower, upper) =
            uncertain_volatility_bounds(&option, pos_or_panic!(0.1), pos_or_panic!(0.3)).unwrap();

        // Lower bound should be less than upper bound for a put
        assert!(lower < upper);
        // Both bounds should be positive
        assert!(lower > Decimal::ZERO);
        assert!(upper > Decimal::ZERO);
    }

    #[test]
    fn test_bounds_same_volatility() {
        let option = create_test_option(OptionStyle::Call, Side::Long, Positive::HUNDRED);
        let vol = pos_or_panic!(0.2);
        let (lower, upper) = uncertain_volatility_bounds(&option, vol, vol).unwrap();

        // Bounds should be equal when min and max volatilities are the same
        assert_decimal_eq!(lower, upper, dec!(1e-10));
    }

    #[test]
    fn test_bounds_itm_call() {
        let itm_option = create_test_option(OptionStyle::Call, Side::Long, pos_or_panic!(90.0));
        let (lower, upper) =
            uncertain_volatility_bounds(&itm_option, pos_or_panic!(0.1), pos_or_panic!(0.3))
                .unwrap();

        // For ITM call, bounds should be above intrinsic value
        let intrinsic = dec!(10.0); // 100 - 90
        assert!(lower > intrinsic);
        assert!(upper > intrinsic);
    }

    #[test]
    fn test_bounds_otm_call() {
        let otm_option = create_test_option(OptionStyle::Call, Side::Long, pos_or_panic!(110.0));
        let (lower, upper) =
            uncertain_volatility_bounds(&otm_option, pos_or_panic!(0.1), pos_or_panic!(0.3))
                .unwrap();

        // For OTM call, both bounds should be positive but lower than strike price
        assert!(lower > Decimal::ZERO);
        assert!(upper < dec!(110.0));
    }

    #[test]
    fn test_bounds_otm_put() {
        let otm_option = create_test_option(OptionStyle::Put, Side::Long, pos_or_panic!(90.0));
        let (lower, upper) =
            uncertain_volatility_bounds(&otm_option, pos_or_panic!(0.1), pos_or_panic!(0.3))
                .unwrap();

        // For OTM put, both bounds should be positive but lower than strike price
        assert!(lower > Decimal::ZERO);
        assert!(upper < dec!(90.0));
    }

    #[test]
    fn test_bounds_extreme_volatilities() {
        let option = create_test_option(OptionStyle::Call, Side::Long, Positive::HUNDRED);
        let (lower, upper) =
            uncertain_volatility_bounds(&option, pos_or_panic!(0.01), Positive::ONE).unwrap();

        assert!(lower < upper);
        assert!(lower > Decimal::ZERO);
        assert!(upper < option.underlying_price.to_dec());
    }

    /// A short position's bounds are the long bounds negated and swapped
    /// (#715), the `short == -long` convention of #646 and #648: both are
    /// negative and still ordered. They used to clamp to zero.
    #[test]
    fn test_bounds_short_are_negated_long_bounds() {
        for style in [OptionStyle::Call, OptionStyle::Put] {
            for strike in [90.0, 100.0, 110.0] {
                let long = create_test_option(style, Side::Long, pos_or_panic!(strike));
                let short = create_test_option(style, Side::Short, pos_or_panic!(strike));
                let (long_lower, long_upper) =
                    uncertain_volatility_bounds(&long, pos_or_panic!(0.1), pos_or_panic!(0.3))
                        .unwrap();
                let (short_lower, short_upper) =
                    uncertain_volatility_bounds(&short, pos_or_panic!(0.1), pos_or_panic!(0.3))
                        .unwrap();
                assert_eq!(short_lower, -long_upper, "{style:?} K={strike}");
                assert_eq!(short_upper, -long_lower, "{style:?} K={strike}");
                assert!(short_lower < short_upper, "{style:?} K={strike}");
                assert!(short_upper < Decimal::ZERO, "{style:?} K={strike}");
            }
        }
    }

    /// Each short bound is the Black–Scholes price of the short position
    /// at the volatility that bounds it: the most negative at the highest
    /// volatility.
    #[test]
    fn test_bounds_short_match_short_black_scholes_prices() {
        let short = create_test_option(OptionStyle::Call, Side::Short, Positive::HUNDRED);
        let at = |volatility: f64| {
            let mut option = short.clone();
            option.implied_volatility = pos_or_panic!(volatility);
            option.calculate_price_black_scholes().unwrap()
        };
        let (lower, upper) =
            uncertain_volatility_bounds(&short, pos_or_panic!(0.1), pos_or_panic!(0.3)).unwrap();
        assert_eq!(lower, at(0.3));
        assert_eq!(upper, at(0.1));
    }
}

#[cfg(test)]
mod tests_adjust_volatility {
    use super::*;
    use optionstratlib_core::assert_pos_relative_eq;

    #[test]
    fn test_same_timeframe() {
        let vol = pos_or_panic!(0.2);
        let result = adjust_volatility(vol, TimeFrame::Day, TimeFrame::Day).unwrap();
        assert_eq!(result, vol);
    }

    #[test]
    fn test_same_periods_different_timeframe() {
        // Create two custom timeframes with same periods
        let vol = pos_or_panic!(0.2);
        let result =
            adjust_volatility(vol, TimeFrame::Custom(pos_or_panic!(252.0)), TimeFrame::Day)
                .unwrap();
        assert_eq!(result, vol);
    }

    #[test]
    fn test_zero_periods() {
        let vol = pos_or_panic!(0.2);
        let result = adjust_volatility(vol, TimeFrame::Day, TimeFrame::Custom(Positive::ZERO));
        assert!(result.is_err());

        if let Err(e) = result {
            assert!(e.to_string().contains("zero periods per year"));
        } else {
            panic!("Expected error for zero periods");
        }
    }

    #[test]
    fn test_daily_to_minute() {
        let daily_vol = pos_or_panic!(0.2);
        let result = adjust_volatility(daily_vol, TimeFrame::Day, TimeFrame::Minute).unwrap();

        // For testing, we can calculate the expected value:
        // daily periods = 252
        // minute periods = 252 * 6.5 * 60 = 98280
        // scale_factor = sqrt(252/98280) ≈ 0.0506
        // expected = 0.2 * 0.0506 ≈ 0.01012
        assert_pos_relative_eq!(result, pos_or_panic!(0.0101273936), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_minute_to_daily() {
        let minute_vol = pos_or_panic!(0.01012);
        let result = adjust_volatility(minute_vol, TimeFrame::Minute, TimeFrame::Day).unwrap();

        // This should be approximately the inverse of the previous test
        assert_pos_relative_eq!(result, pos_or_panic!(0.199853), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_daily_to_hourly() {
        let daily_vol = pos_or_panic!(0.2);
        let result = adjust_volatility(daily_vol, TimeFrame::Day, TimeFrame::Hour).unwrap();

        // daily periods = 252
        // hourly periods = 252 * 6.5 = 1638
        // scale_factor = sqrt(252/1638) ≈ 0.3922
        // expected = 0.2 * 0.3922 ≈ 0.07844
        assert_pos_relative_eq!(result, pos_or_panic!(0.07844), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_monthly_to_daily() {
        let monthly_vol = pos_or_panic!(0.3);
        let result = adjust_volatility(monthly_vol, TimeFrame::Month, TimeFrame::Day).unwrap();

        // monthly periods = 12
        // daily periods = 252
        // scale_factor = sqrt(12/252) ≈ 0.218
        // expected = 0.3 * 0.218 ≈ 0.0654
        assert_pos_relative_eq!(result, pos_or_panic!(0.0654653), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_custom_timeframe() {
        let vol = pos_or_panic!(0.25);
        let custom_periods = Positive::HUNDRED;
        let result =
            adjust_volatility(vol, TimeFrame::Custom(custom_periods), TimeFrame::Day).unwrap();

        // custom periods = 100
        // daily periods = 252
        // scale_factor = sqrt(100/252) ≈ 0.629
        // expected = 0.25 * 0.629 ≈ 0.15725
        assert_pos_relative_eq!(result, pos_or_panic!(0.157485197), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_yearly_to_daily() {
        let yearly_vol = pos_or_panic!(0.4);
        let result = adjust_volatility(yearly_vol, TimeFrame::Year, TimeFrame::Day).unwrap();

        // yearly periods = 1
        // daily periods = 252
        // scale_factor = sqrt(1/252) ≈ 0.0629
        // expected = 0.4 * 0.0629 ≈ 0.02516
        assert_pos_relative_eq!(result, pos_or_panic!(0.025197631), pos_or_panic!(0.0001));
    }

    #[test]
    fn test_zero_volatility() {
        let result = adjust_volatility(Positive::ZERO, TimeFrame::Day, TimeFrame::Minute).unwrap();
        assert_eq!(result, Positive::ZERO);
    }

    #[test]
    fn test_very_small_volatility() {
        let small_vol = pos_or_panic!(0.0001);
        let result = adjust_volatility(small_vol, TimeFrame::Day, TimeFrame::Hour).unwrap();
        assert!(result > Positive::ZERO);
        assert!(result < small_vol);
    }

    #[test]
    fn test_very_large_volatility() {
        let large_vol = pos_or_panic!(10.0);
        let result = adjust_volatility(large_vol, TimeFrame::Day, TimeFrame::Minute).unwrap();
        assert!(result > Positive::ZERO);
        assert!(result < large_vol);
    }
}

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

    #[test]
    fn constant_volatility_nan_return_surfaces_decimal_error() {
        // A NaN return in the input slice feeds the `d_sub(r, mean, ..)`
        // inside constant_volatility; the current implementation treats
        // NaN as 0 at the Decimal boundary via `Decimal` defaults, but
        // this test pins the behaviour for future guard tightening.
        let returns = [Decimal::ZERO, Decimal::ZERO];
        // Happy path: no NaN, function returns zero vol.
        let v = constant_volatility(&returns).expect("finite inputs");
        assert_eq!(v, Positive::ZERO);
    }

    #[test]
    fn heston_happy_path_does_not_trip_non_finite() {
        // With finite, well-posed inputs the Heston simulator must
        // succeed without hitting the `finite_decimal(dw_f64)` guard.
        use rust_decimal_macros::dec;
        let res = simulate_heston_volatility(
            dec!(1.0),  // kappa
            dec!(0.04), // theta
            dec!(0.1),  // xi
            dec!(0.04), // v0
            dec!(0.01), // dt
            10,         // steps
            &mut optionstratlib_core::utils::deterministic_rng(
                optionstratlib_core::utils::DETERMINISTIC_RNG_DEFAULT_SEED,
            ),
        );
        assert!(res.is_ok(), "finite inputs unexpectedly failed: {res:?}");
    }
}

/// The rolling `historical_volatility` of #859 against the per-window
/// two-pass estimate it replaced.
#[cfg(test)]
mod tests_rolling_historical_volatility {
    use super::*;
    use rust_decimal_macros::dec;

    /// The previous implementation: a two-pass estimate per window.
    fn two_pass(returns: &[Decimal], window: usize) -> Result<Vec<Positive>, VolatilityError> {
        returns.windows(window).map(constant_volatility).collect()
    }

    /// Deterministic returns: `scale · ((i·7919 mod 1009) - 504) / 1009`,
    /// plus `drift · i`, so the mean moves across the series.
    fn series(len: usize, scale: Decimal, drift: Decimal) -> Vec<Decimal> {
        (0..len)
            .map(|i| {
                let k = Decimal::from((i * 7919) % 1009) - dec!(504);
                scale * k / dec!(1009) + drift * Decimal::from(i)
            })
            .collect()
    }

    fn max_relative_difference(returns: &[Decimal], window: usize) -> Decimal {
        let (rolling, reference) = match (
            historical_volatility(returns, window),
            two_pass(returns, window),
        ) {
            (Ok(a), Ok(b)) => (a, b),
            other => panic!("both must succeed: {other:?}"),
        };
        assert_eq!(rolling.len(), reference.len());
        let mut worst = Decimal::ZERO;
        for (a, b) in rolling.iter().zip(&reference) {
            let diff = (a.to_dec() - b.to_dec()).abs();
            let rel = if b.to_dec().is_zero() {
                diff
            } else {
                diff / b.to_dec()
            };
            worst = worst.max(rel);
        }
        worst
    }

    /// The largest relative difference measured on this grid is `7.5e-18`
    /// (returns of scale `1e-4` with a drift of `1e-5` per step, where the
    /// mean moves furthest from the shift). `1e-15` leaves a hundredfold
    /// margin over rounding and is five orders of magnitude below the
    /// `1e-10` an algorithmic error would breach.
    const ROLLING_TOLERANCE: Decimal = dec!(0.000000000000001);

    #[test]
    fn test_rolling_matches_two_pass_on_a_grid() {
        for scale in [dec!(0.0001), dec!(0.01), dec!(0.2), dec!(10)] {
            for drift in [dec!(0), dec!(0.00001), dec!(0.01)] {
                for window in [2usize, 3, 5, 21, 63, 252] {
                    let returns = series(1008, scale, drift);
                    let rel = max_relative_difference(&returns, window);
                    assert!(
                        rel <= ROLLING_TOLERANCE,
                        "scale {scale}, drift {drift}, window {window}: {rel}"
                    );
                }
            }
        }
    }

    #[test]
    fn test_rolling_matches_two_pass_on_the_bench_series() {
        let bench: Vec<Decimal> = (0..1008)
            .map(|i| (Decimal::from(i % 23) - dec!(11)) / dec!(997))
            .collect();
        for window in [2usize, 21, 252] {
            assert!(max_relative_difference(&bench, window) <= ROLLING_TOLERANCE);
        }
    }

    #[test]
    fn test_rolling_equal_returns_are_exactly_zero() {
        let returns = vec![dec!(0.0123); 50];
        for window in [2usize, 7, 50] {
            match historical_volatility(&returns, window) {
                Ok(volatilities) => assert!(volatilities.iter().all(|v| *v == Positive::ZERO)),
                Err(e) => panic!("equal returns must succeed: {e:?}"),
            }
        }
    }

    #[test]
    fn test_rolling_edge_windows_match_two_pass() {
        let returns = series(10, dec!(0.01), dec!(0));
        for window in [1usize, 10, 11] {
            assert_eq!(
                format!("{:?}", historical_volatility(&returns, window)),
                format!("{:?}", two_pass(&returns, window))
            );
        }
        assert!(historical_volatility(&returns, 0).is_err());
    }

    #[test]
    fn test_rolling_overflow_falls_back_to_the_two_pass_error() {
        let huge = Decimal::MAX / dec!(2);
        let returns = vec![huge, huge, huge, dec!(1)];
        assert_eq!(
            format!("{:?}", historical_volatility(&returns, 3)),
            format!("{:?}", two_pass(&returns, 3))
        );
    }
}