optionstratlib-simulation 0.22.1

//! The Ornstein-Uhlenbeck path generator.
//!
//! A mean-reverting path is a simulation kernel, not a volatility model: it
//! reports [`crate::error::SimulationError`] and names nothing from the
//! pricing layer. `volatility::generate_ou_process` re-exports it, so the
//! 0.21 path is unchanged.

use crate::error::SimulationError;
use crate::simulation::path_buffer;
use optionstratlib_core::model::Positive;
use optionstratlib_core::model::decimal::{d_add, d_mul, decimal_normal_sample_with, p_sqrt};
use optionstratlib_core::model::utils::sub_floor_zero;
use rand::Rng;
use rust_decimal::Decimal;

/// 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
///
/// # Randomness
///
/// The normal increments come from the thread RNG, so two calls differ. A
/// reproducible OU path is a [`crate::simulation::WalkType::MeanReverting`]
/// walk with [`crate::simulation::WalkParams::seed`] set.
///
/// # Errors
///
/// Returns [`SimulationError::PositiveError`] when `sqrt(dt)` or the mean
/// reversion term `theta * (mu - x) * dt` leaves the representable `Positive`
/// range, and [`SimulationError::Decimal`] when the diffusion term or the
/// accumulated level leaves the representable `Decimal` range.
/// Returns [`SimulationError::InvalidParameters`] when a path of `steps`
/// points cannot be allocated.
///
/// There is no limit to return in their place: the process level is the state
/// being simulated, so a step that cannot be represented ends the path rather
/// than continuing from a stand-in value.
///
/// # Examples
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use optionstratlib_core::model::Positive;
/// use optionstratlib_core::pos_or_panic;
/// use optionstratlib_simulation::simulation::ou::generate_ou_process;
///
/// # fn main() -> Result<(), optionstratlib_simulation::error::SimulationError> {
/// // 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
/// )?;
/// # Ok(())
/// # }
/// ```
pub fn generate_ou_process(
    x0: Positive,
    mu: Positive,
    theta: Positive,
    volatility: Positive,
    dt: Positive,
    steps: usize,
) -> Result<Vec<Positive>, SimulationError> {
    ou_path(x0, mu, theta, volatility, dt, steps, &mut rand::rng())
}

/// Body of [`generate_ou_process`], drawing its normal samples from `rng`.
///
/// The seeded walk kernels (`MeanReverting`, and the volatility path of
/// `Custom`) call it with the walk's generator so the whole path comes from
/// one stream.
///
/// # Errors
///
/// Same as [`generate_ou_process`].
pub(crate) fn ou_path<R: Rng + ?Sized>(
    x0: Positive,
    mu: Positive,
    theta: Positive,
    volatility: Positive,
    dt: Positive,
    steps: usize,
    rng: &mut R,
) -> Result<Vec<Positive>, SimulationError> {
    let sqrt_dt = p_sqrt(&dt, "volatility::utils::generate_ou_process")?;
    let mut x = x0.to_dec();
    let mut result = path_buffer(steps, 0)?;
    result.push(Positive::new_decimal(x).unwrap_or(Positive::ZERO));

    for _ in 1..steps {
        // Z√dt
        let dw = d_mul(
            decimal_normal_sample_with(rng),
            sqrt_dt.to_dec(),
            "volatility::ou::dw",
        )?;
        // θ(μ−x)dt
        let drift = theta
            .checked_mul(&sub_floor_zero(mu, &x))?
            .checked_mul(&dt)?
            .to_dec();
        // σ·Z√dt (Decimal)
        let diffusion = d_mul(dw, volatility.to_dec(), "volatility::ou::diffusion")?;
        // OU process step
        x = d_add(
            x,
            d_add(drift, diffusion, "volatility::ou::step")?,
            "volatility::ou::step",
        )?;
        x = x.max(Decimal::ZERO); // no negative values
        result.push(Positive::new_decimal(x).unwrap_or(Positive::ZERO));
    }

    Ok(result)
}

#[cfg(test)]
mod tests_generate_ou_process {
    use super::*;
    use optionstratlib_core::pos_or_panic;
    use optionstratlib_core::utils::deterministic_rng;
    use rust_decimal_macros::dec;

    /// Seed of the simulated paths, so every run draws the same stream.
    const SEED: u64 = 685;

    #[test]
    fn test_process_length() {
        let steps = 500;
        // Unseeded on purpose: covers the thread-RNG entry point, whose length is stream-independent.
        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,
        )
        .unwrap();
        assert_eq!(process.len(), steps);
    }

    #[test]
    fn test_all_values_positive() {
        let process = ou_path(
            Positive::ONE,
            pos_or_panic!(1.5),
            pos_or_panic!(0.2),
            pos_or_panic!(0.3),
            pos_or_panic!(0.01),
            1000,
            &mut deterministic_rng(SEED),
        )
        .unwrap();

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

    #[test]
    fn test_mean_reversion_tendency() {
        let process = ou_path(
            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,
            &mut deterministic_rng(SEED),
        )
        .unwrap();

        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}");
    }
}