optionstratlib 0.17.0

/******************************************************************************
   Author: Joaquín Béjar García
   Email: jb@taunais.com
   Date: 15/8/24
******************************************************************************/
// Scoped allow: bulk migration of unchecked `[]` indexing to
// `.get().ok_or_else(..)` tracked as follow-ups to #341. The existing
// call sites are internal to this file and audited for invariant-bound
// indices (fixed-length buffers, just-pushed slices, etc.).
#![allow(clippy::indexing_slicing)]

use crate::constants::{MAX_VOLATILITY, MIN_VOLATILITY};
use crate::error::VolatilityError;
use crate::model::decimal::{
    d_add, d_div, d_mul, d_sub, d_sum, decimal_normal_sample, finite_decimal,
};
use crate::utils::time::TimeFrame;
use crate::{ExpirationDate, OptionStyle, OptionType, Options, Side};
use num_traits::{FromPrimitive, ToPrimitive};
use positive::Positive;
use rand::random;
use rayon::prelude::*;
use rust_decimal::{Decimal, MathematicalOps};
use tracing::instrument;

#[cfg(test)]
use positive::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. Conversion of the final std-dev into
/// `Positive` is clamped to `Positive::ZERO`, so `PositiveError` is
/// never surfaced.
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()
            ),
        })?;
    let n = Positive::new_decimal(n_dec).unwrap_or(Positive::ZERO);

    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 variance = d_div(
        sq_total,
        (n - Decimal::ONE).to_dec(),
        "volatility::constant::variance",
    )?;

    let std_dev = variance
        .sqrt()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "constant_volatility: sqrt(variance) failed (overflow)".to_string(),
        })?;
    Ok(Positive::new_decimal(std_dev).unwrap_or(Positive::ZERO))
}

/// 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
///
/// Propagates any `VolatilityError::NumericalFailure` raised by
/// `constant_volatility` on an individual window.
///
/// # Panics
///
/// The underlying `slice::windows(window_size)` call panics when
/// `window_size == 0`. 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> {
    returns
        .windows(window_size)
        .map(constant_volatility)
        .collect()
}

/// 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 converted via
/// `Positive::new_decimal(...).unwrap_or(Positive::ZERO)`, so
/// `PositiveError` is never surfaced — negative intermediates are
/// 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 = variance
        .sqrt()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "ewma_volatility: sqrt(initial variance) failed (overflow)".to_string(),
        })?;
    let mut volatilities = vec![Positive::new_decimal(initial_std_dev).unwrap_or(Positive::ZERO)];

    for &return_value in &returns[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 = variance
            .sqrt()
            .ok_or_else(|| VolatilityError::NumericalFailure {
                reason: "ewma_volatility: sqrt(variance) failed (overflow)".to_string(),
            })?;
        volatilities.push(Positive::new_decimal(std_dev).unwrap_or(Positive::ZERO));
    }

    Ok(volatilities)
}

/// Calculates the implied volatility of an option given its market price.
///
/// This function uses the Newton-Raphson method to iteratively approximate the implied
/// volatility that corresponds to the observed market price of the option. The implied
/// volatility is updated within the `Options` struct provided as a mutable reference.
///
/// # Parameters
/// - `market_price`: The observed market price of the option.
/// - `options`: A mutable reference to an `Options` struct, which should contain the necessary
///   methods and fields such as `implied_volatility`, `calculate_price_black_scholes()`, and `vega()`.
/// - `max_iterations`: The maximum number of iterations allowed for the Newton-Raphson method.
///
/// # Returns
/// The function returns the estimated implied volatility of the option.
///
/// # Remarks
/// - If the price difference between the calculated and market price is within the tolerated threshold (`TOLERANCE`),
///   the current implied volatility is returned.
/// - The function ensures that the implied volatility stays positive.
///
/// # Errors
///
/// The implementation is a parallel grid search over
/// `100 * max_iterations` candidate volatilities rather than a
/// Newton–Raphson iteration. It returns
/// `VolatilityError::IvNotFound` when the best candidate is the
/// lower boundary `1 / (100 * max_iterations)` (the search never
/// improved), `VolatilityError::NoValidVolatility` when every grid
/// point failed the Black–Scholes evaluation, and
/// `VolatilityError::PositiveError` when the boundary `Positive`
/// conversion fails. Black–Scholes errors on individual candidates
/// are discarded by the parallel filter rather than propagated.
#[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> {
    let base_option = options.clone();
    let iterations = 100 * max_iterations;
    let result = (1..iterations)
        .into_par_iter()
        .map(|i| {
            let mut option = base_option.clone();
            option.side = Side::Long; // Ensure the option is long
            let iv = Positive::new(i as f64 / iterations as f64).unwrap_or(Positive::ZERO);
            option.implied_volatility = iv;

            match option.calculate_price_black_scholes() {
                Ok(price) => {
                    let diff = (price - market_price.to_dec()).abs();
                    Some((iv, diff))
                }
                Err(_) => None,
            }
        })
        .filter_map(|x| x) // Remove errors
        .min_by(|a, b| a.1.partial_cmp(&b.1).unwrap_or(std::cmp::Ordering::Equal));

    match result {
        Some((best_iv, _)) => {
            let iv = best_iv.clamp(*MIN_VOLATILITY, MAX_VOLATILITY);
            if iv == Positive::new(1f64 / iterations as f64)? {
                Err(VolatilityError::IvNotFound)
            } else {
                Ok(iv)
            }
        }
        None => Err(VolatilityError::NoValidVolatility),
    }
}

/// 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 iteration
/// limit for the IV calculation is set to 10.
///
/// 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,
    );
    implied_volatility(option_price, &mut option, 10)
}

/// 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 [`crate::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.
///
/// # Panics
///
/// Panics on `returns[0]` when `returns` is empty. Callers must
/// validate the input slice length before invoking.
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.
    let seed = d_mul(
        returns[0],
        returns[0],
        "volatility::garch::initial_variance",
    )?;
    let mut variance = to_positive_variance(seed, "initial_variance")?;
    let mut volatilities = vec![variance.sqrt()];
    for &return_value in &returns[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(variance.sqrt());
    }
    Ok(volatilities)
}

/// 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.
///
/// # 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`.
///
/// Variance is clamped to zero on each step
/// (`v = v.max(Decimal::ZERO)`) and converted via
/// `Positive::new_decimal(...).unwrap_or(Positive::ZERO)`, so
/// negative variance is silently recovered and `PositiveError` is
/// never surfaced.
pub fn simulate_heston_volatility(
    kappa: Decimal,
    theta: Decimal,
    xi: Decimal,
    v0: Decimal,
    dt: Decimal,
    steps: usize,
) -> Result<Vec<Positive>, VolatilityError> {
    let mut v = v0.max(Decimal::ZERO);
    let mut v_pos = Positive::new_decimal(v).unwrap_or(Positive::ZERO);
    let mut volatilities = vec![v_pos.sqrt()];
    let dt_sqrt_f64 = dt
        .sqrt()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "simulate_heston_volatility: sqrt(dt) failed (overflow)".to_string(),
        })?
        .to_f64()
        .ok_or_else(|| VolatilityError::NumericalFailure {
            reason: "simulate_heston_volatility: sqrt(dt) not representable as f64".to_string(),
        })?;
    for _ in 1..steps {
        let dw_f64 = random::<f64>() * dt_sqrt_f64;
        let dw = finite_decimal(dw_f64)
            .ok_or_else(|| VolatilityError::non_finite("volatility::heston::dw", dw_f64))?;
        let sqrt_v = v_pos.sqrt().to_dec();
        v += kappa * (theta - v) * dt + xi * sqrt_v * dw;
        v = v.max(Decimal::ZERO); // Ensure variance doesn't become negative
        v_pos = Positive::new_decimal(v).unwrap_or(Positive::ZERO);
        volatilities.push(v_pos.sqrt());
    }
    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 of (lower_bound, upper_bound) for the option price.
///
/// # 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<(Positive, Positive), VolatilityError> {
    // Create a clone of the option for lower bound calculation
    let mut lower_bound_option = option.clone();
    lower_bound_option.implied_volatility = min_volatility;

    // Create a clone of the option for upper bound calculation
    let mut upper_bound_option = option.clone();
    upper_bound_option.implied_volatility = max_volatility;

    // Calculate the option price with minimum volatility
    let lower_bound = Positive::new_decimal(lower_bound_option.calculate_price_black_scholes()?)
        .unwrap_or(Positive::ZERO);

    // Calculate the option price with maximum volatility
    let upper_bound = Positive::new_decimal(upper_bound_option.calculate_price_black_scholes()?)
        .unwrap_or(Positive::ZERO);

    Ok((lower_bound, upper_bound))
}

/// 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 positive::pos_or_panic;
/// use optionstratlib::utils::time::TimeFrame;
/// use optionstratlib::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` (e.g. overflow on an extreme timeframe annualisation).
#[inline]
pub fn annualized_volatility(
    volatility: Positive,
    timeframe: TimeFrame,
) -> Result<Positive, VolatilityError> {
    Ok(volatility * timeframe.periods_per_year().sqrt())
}

/// 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 positive::pos_or_panic;
/// use optionstratlib::utils::time::TimeFrame;
/// use optionstratlib::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::PositiveError`] when the rescaling
/// produces a value that violates the `Positive` invariant, typically
/// due to division rounding on an extremely small timeframe.
#[inline]
pub fn de_annualized_volatility(
    annual_volatility: Positive,
    timeframe: TimeFrame,
) -> Result<Positive, VolatilityError> {
    Ok(annual_volatility / timeframe.periods_per_year().sqrt())
}

/// 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<(), optionstratlib::error::Error> {
/// use positive::pos_or_panic;
/// use optionstratlib::utils::TimeFrame;
/// use optionstratlib::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
///
/// Propagates any [`VolatilityError::PositiveError`] surfaced by the
/// intermediate [`annualized_volatility`] or [`de_annualized_volatility`]
/// calls when either rescaling breaches the `Positive` invariant.
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 square root of (from_periods / to_periods). `Positive`
    // already implements `sqrt()`, so skip the f64 round-trip and keep the
    // computation in `Decimal` precision end-to-end.
    let scale_factor = (from_periods / to_periods).sqrt();

    Ok(volatility * 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<(), optionstratlib::error::Error> {
/// use positive::{pos_or_panic, Positive};
/// use optionstratlib::utils::time::{TimeFrame, convert_time_frame};
/// use optionstratlib::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)
}

/// Generates a mean-reverting Ornstein-Uhlenbeck process time series
///
/// This function simulates a discrete-time Ornstein-Uhlenbeck stochastic process, which is
/// commonly used in financial mathematics to model mean-reverting processes such as interest rates,
/// volatility, or commodity prices. The process follows the stochastic differential equation:
///
/// dX_t = θ(μ - X_t)dt + σdW_t
///
/// Where:
/// - θ (theta) is the speed of reversion to the mean
/// - μ (mu) is the long-term mean level
/// - σ (sigma) is the volatility or intensity of random fluctuations
/// - W_t is a Wiener process (standard Brownian motion)
///
/// # Arguments
/// * `x0` - Initial value of the process
/// * `mu` - Long-term mean the process reverts to
/// * `theta` - Speed of mean reversion (higher values cause faster reversion)
/// * `sigma` - Volatility parameter controlling the intensity of random fluctuations
/// * `dt` - Time step size for the simulation
/// * `steps` - Number of time steps to simulate
///
/// # Returns
/// A vector containing the simulated values of the Ornstein-Uhlenbeck process at each time step
///
/// # Examples
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use positive::{pos_or_panic, Positive};
/// use optionstratlib::volatility::generate_ou_process;
///
/// // Simulate an OU process with initial value 1.0, mean 1.5,
/// // reversion speed 0.1, volatility 0.2, time step 0.01, for 1000 steps
/// let process = generate_ou_process(
///     Positive::ONE,       // initial value
///     pos_or_panic!(1.5),       // long-term mean
///     pos_or_panic!(0.1),       // speed of reversion
///     pos_or_panic!(0.2),       // volatility
///     pos_or_panic!(0.01),      // time step
///     1000             // number of steps
/// );
/// ```
#[must_use]
pub fn generate_ou_process(
    x0: Positive,
    mu: Positive,
    theta: Positive,
    volatility: Positive,
    dt: Positive,
    steps: usize,
) -> Vec<Positive> {
    let sqrt_dt = dt.sqrt();
    let mut x = x0.to_dec();
    let mut result = Vec::with_capacity(steps);
    result.push(Positive::new_decimal(x).unwrap_or(Positive::ZERO));

    for _ in 1..steps {
        let dw = decimal_normal_sample() * sqrt_dt; // Z√dt
        let drift = (theta * mu.sub_or_zero(&x) * dt).to_dec(); // θ(μ−x)dt
        let diffusion = dw * volatility; // σ·Z√dt (Decimal)
        x += drift + diffusion; // OU process step
        x = x.max(Decimal::ZERO); // no negative values
        result.push(Positive::new_decimal(x).unwrap_or(Positive::ZERO));
    }

    result
}

#[cfg(test)]
mod tests_annualize_volatility {
    use super::*;
    use positive::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 crate::constants::ZERO;
    use positive::assert_pos_relative_eq;
    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 positive::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 positive::assert_pos_relative_eq;
    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::ExpirationDate;
    use crate::assert_decimal_eq;
    use crate::constants::{MAX_VOLATILITY, MIN_VOLATILITY};
    use crate::greeks::Greeks;
    use crate::model::types::{OptionStyle, OptionType, Side};
    use positive::assert_pos_relative_eq;
    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);

        // Test with very low number of iterations
        let result = implied_volatility(market_price, &mut option, 1);
        assert!(result.is_ok()); // Should still return a result

        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);
        assert_pos_relative_eq!(iv, pos_or_panic!(0.437), pos_or_panic!(1e-3));
    }

    #[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.0073826), 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.0073826), 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.866), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-26.990), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.277), dec!(0.002));
        assert_decimal_eq!(vanna, dec!(-0.4879786), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(6.2450679681), dec!(0.00000001));
        assert_decimal_eq!(veta, dec!(0.0433019), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1686671), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005877), 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.053), dec!(0.001));
        assert_decimal_eq!(gamma, dec!(0.0), dec!(0.001));
        assert_decimal_eq!(vega, dec!(0.866), dec!(0.001));
        assert_decimal_eq!(theta, dec!(-29.874), dec!(0.001));
        assert_decimal_eq!(rho, dec!(-0.016), dec!(0.001));
        assert_decimal_eq!(vanna, dec!(-0.4879786), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(6.2450679681), dec!(0.00000001));
        assert_decimal_eq!(veta, dec!(0.0433019), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1685301), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005877), 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();
        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.866), dec!(0.002));
        assert_decimal_eq!(theta, dec!(-26.990), dec!(0.002));
        assert_decimal_eq!(rho, dec!(0.277), dec!(0.002));
        assert_decimal_eq!(vanna, dec!(-0.4879786), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(6.2450679681), dec!(0.00000001));
        assert_decimal_eq!(veta, dec!(0.0433019), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1686671), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005877), 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();
        assert_decimal_eq!(delta, dec!(0.053), dec!(0.001));
        assert_decimal_eq!(gamma, dec!(0.0), dec!(0.001));
        assert_decimal_eq!(vega, dec!(0.866), dec!(0.001));
        assert_decimal_eq!(theta, dec!(-29.874), dec!(0.001));
        assert_decimal_eq!(rho, dec!(-0.016), dec!(0.001));
        assert_decimal_eq!(vanna, dec!(-0.4879786), dec!(0.0000001));
        assert_decimal_eq!(vomma, dec!(6.2450679681), dec!(0.00000001));
        assert_decimal_eq!(veta, dec!(0.0433019), dec!(0.0000001));
        assert_decimal_eq!(charm, dec!(0.1685301), dec!(0.0000001));
        assert_decimal_eq!(color, dec!(0.0005877), 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.0144632), 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.0274529), 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 crate::assert_decimal_eq;
    use positive::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::*;
    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).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).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).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).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).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).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).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).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).unwrap();

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

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

    use crate::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 > Positive::ZERO);
        assert!(upper > Positive::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 > Positive::ZERO);
        assert!(upper > Positive::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_pos_relative_eq!(lower, upper, pos_or_panic!(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 = pos_or_panic!(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 > Positive::ZERO);
        assert!(upper < pos_or_panic!(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 > Positive::ZERO);
        assert!(upper < pos_or_panic!(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 > Positive::ZERO);
        assert!(upper < option.underlying_price);
    }
}

#[cfg(test)]
mod tests_adjust_volatility {
    use super::*;
    use positive::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_generate_ou_process {
    use super::*;
    use rust_decimal_macros::dec;

    #[test]
    fn test_process_length() {
        let steps = 500;
        let process = generate_ou_process(
            Positive::ONE,
            pos_or_panic!(1.5),
            pos_or_panic!(0.1),
            pos_or_panic!(0.2),
            pos_or_panic!(0.01),
            steps,
        );
        assert_eq!(process.len(), steps);
    }

    #[test]
    fn test_all_values_positive() {
        let process = generate_ou_process(
            Positive::ONE,
            pos_or_panic!(1.5),
            pos_or_panic!(0.2),
            pos_or_panic!(0.3),
            pos_or_panic!(0.01),
            1000,
        );

        for value in process {
            assert!(
                value >= Positive::ZERO,
                "Found non-positive value: {value:?}"
            );
        }
    }

    #[test]
    fn test_mean_reversion_tendency() {
        let process = generate_ou_process(
            pos_or_panic!(0.1),
            Positive::ONE,
            Positive::ONE,       // high theta for fast reversion
            pos_or_panic!(0.01), // low volatility
            pos_or_panic!(0.01),
            1000,
        );

        let last = process.last().unwrap().to_dec();
        let diff = (last - dec!(1.0)).abs();
        assert!(diff < dec!(0.1), "Final value too far from mean: {last}");
    }
}

#[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
        );
        assert!(res.is_ok(), "finite inputs unexpectedly failed: {res:?}");
    }
}