optionstratlib 0.17.0

use crate::Options;
use crate::error::PricingError;
use crate::model::decimal::{d_div, d_mul, d_sub, finite_decimal};
use crate::pricing::utils::wiener_increment;
use num_traits::{FromPrimitive, ToPrimitive};
use positive::Positive;
use rust_decimal::{Decimal, MathematicalOps};
use std::num::NonZeroUsize;
use tracing::instrument;

/// This function performs Monte Carlo simulation to price an option.
///
/// # Arguments
///
/// * `option` - An `Options` struct containing the option's parameters, such as underlying price, strike price, risk-free rate, implied volatility, and expiration date.
/// * `steps` - An integer indicating the number of time steps in the simulation.
/// * `simulations` - An integer indicating the number of Monte Carlo simulations to run.
///
/// # Returns
///
/// * A floating-point number representing the estimated price of the option.
///
/// # Description
///
/// The function follows the below steps:
///
/// 1. Calculate the time increment `dt` based on the number of steps.
/// 2. Initialize a sum variable `payoff_sum` to accumulate the payoffs from each simulation.
/// 3. Loop through the number of simulations:
///     - For each simulation, initialize the stock price `st` to the underlying price.
///     - Loop through the number of steps:
///         - Calculate a Wiener process increment `w`.
///         - Update the stock price `st` using the discrete approximation of the geometric Brownian motion model.
///     - Calculate the payoff of the option for this simulation (for a call option, this is `max(st - strike_price, 0)`).
///     - Add the payoff to the `payoff_sum`.
/// 4. Return the average payoff discounted to its present value.
///
/// # Errors
///
/// Returns `PricingError::ExpirationDate` when the option's
/// expiration cannot be converted to a positive year fraction, and
/// `PricingError::MethodError` when the GBM discretisation
/// encounters a non-finite value (e.g. volatility overflow) or when
/// the terminal payoff averaging produces a non-representable
/// `Decimal`.
#[instrument(skip(option), fields(
    steps = steps.get(),
    simulations = simulations.get(),
    strike = %option.strike_price,
    style = ?option.option_style,
    side = ?option.side,
))]
pub fn monte_carlo_option_pricing(
    option: &Options,
    steps: NonZeroUsize,       // Number of time steps per path
    simulations: NonZeroUsize, // Number of Monte Carlo simulations
) -> Result<Decimal, PricingError> {
    let steps_raw = steps.get();
    let simulations_raw = simulations.get();
    let dt = option.expiration_date.get_years()? / steps_raw as f64;
    let mut payoff_sum = 0.0;

    for _ in 0..simulations_raw {
        let mut st = option.underlying_price.to_dec();
        for _ in 0..steps_raw {
            let w = wiener_increment(dt.to_dec())?;
            st *=
                Decimal::ONE + option.risk_free_rate * dt + option.implied_volatility.to_dec() * w;
        }
        // Calculate the payoff for a call option
        let payoff_dec = d_sub(
            st,
            option.strike_price.to_dec(),
            "pricing::monte_carlo::gbm::payoff",
        )?
        .max(Decimal::ZERO);
        let payoff: f64 = payoff_dec.to_f64().ok_or_else(|| {
            PricingError::non_finite("pricing::monte_carlo::gbm::payoff_cast", f64::NAN)
        })?;
        if !payoff.is_finite() {
            return Err(PricingError::non_finite(
                "pricing::monte_carlo::gbm::payoff",
                payoff,
            ));
        }
        payoff_sum += payoff;
    }
    // Average value of the payoffs discounted to present value.
    // Guard every `f64` boundary against NaN / ±∞ so saturation on
    // the rate, the discount exponent, or the final average surfaces
    // a tagged `PricingError::NonFinite` instead of silently collapsing
    // to `Decimal::ZERO` through the `f2d!` cast.
    let rate_f64 = option.risk_free_rate.to_f64().ok_or_else(|| {
        PricingError::non_finite("pricing::monte_carlo::rate_f64::cast", f64::NAN)
    })?;
    if !rate_f64.is_finite() {
        return Err(PricingError::non_finite(
            "pricing::monte_carlo::rate_f64",
            rate_f64,
        ));
    }
    let years = option.expiration_date.get_years()?.to_f64();
    if !years.is_finite() {
        return Err(PricingError::non_finite(
            "pricing::monte_carlo::years",
            years,
        ));
    }
    let discount = (-rate_f64 * years).exp();
    if !discount.is_finite() {
        return Err(PricingError::non_finite(
            "pricing::monte_carlo::discount",
            discount,
        ));
    }
    let average_payoff = (payoff_sum / simulations_raw as f64) * discount;
    if !average_payoff.is_finite() {
        return Err(PricingError::non_finite(
            "pricing::monte_carlo::average_payoff",
            average_payoff,
        ));
    }
    finite_decimal(average_payoff).ok_or_else(|| {
        PricingError::non_finite("pricing::monte_carlo::average_payoff::cast", average_payoff)
    })
}

/// Estimates the price of a financial option using the Monte Carlo simulation method.
///
/// # Parameters
/// - `option`: A reference to an `Options` object that represents the option being evaluated.
///   This object contains necessary details such as risk-free rate, dividend yield,
///   expiration date, and the payoff calculation logic.
/// - `final_prices`: A slice of `Positive` values representing the simulated final prices of
///   the underlying asset at the option's expiration. The length of this slice
///   corresponds to the number of simulations.
///
/// # Returns
/// - `Result<Positive, PricingError>`: Returns a `Positive` value encapsulating the estimated
///   option price calculated using the Monte Carlo method. If an error occurs during intermediate
///   calculations (e.g., getting the expiration year), it returns a `PricingError`.
///
/// # How it Works
/// 1. The number of simulations is determined by the length of the `final_prices` slice. If the slice
///    is empty, the function immediately returns a price of `Positive::ZERO`.
/// 2. Calculates the effective discount factor based on the risk-free rate, dividend yield, and
///    time to expiration. This factor is used to discount future payoffs to their present value.
/// 3. For each simulated final price in the `final_prices` slice:
///    - Compute the payoff using the `option.payoff_at_price` method.
///    - Accumulate the total payoff across all simulations.
/// 4. Compute the average payoff by dividing the total payoff by the number of simulations.
///    The average payoff is then discounted using the calculated discount factor.
/// 5. Return the discounted average payoff as the estimated option price.
///
/// # Errors
/// - Propagates any error from `option.expiration_date.get_years()?` (e.g.,
///   invalid expiration date).
/// - Returns `PricingError::method_error` when `num_simulations` cannot be
///   represented as a `Decimal` (effectively unreachable for valid `usize`
///   inputs but surfaced explicitly for completeness).
/// - Per-simulation `option.payoff_at_price(...)` failures fall back silently
///   to `Decimal::ZERO` for that simulation; the function does not panic.
///
/// This function assumes that the `Options` struct and `Positive` type
/// are implemented elsewhere in the codebase and provide necessary functionality (e.g., payoff calculation).
pub fn price_option_monte_carlo(
    option: &Options,
    final_prices: &[Positive],
) -> Result<Positive, PricingError> {
    // The number of simulations is the length of the final prices vector
    let num_simulations = final_prices.len();

    if num_simulations == 0 {
        return Ok(Positive::ZERO);
    }

    // Calculate total discount factor (risk-free rate adjusted for dividends)
    let effective_rate = option.risk_free_rate - option.dividend_yield;
    let discount_factor = (-effective_rate * option.expiration_date.get_years()?).exp();

    // Calculate payoff for each final price and sum them
    let total_payoff: Decimal = final_prices
        .iter()
        .map(|&final_price| {
            option
                .payoff_at_price(&final_price)
                .unwrap_or(Decimal::ZERO)
        })
        .sum();

    // Average payoff discounted to present value. Both the mean and the
    // discounting are fused monetary flows, so they go through the checked
    // helpers; `d_div` applies banker's rounding at `DIV_DEFAULT_SCALE`.
    let n_dec = Decimal::from_usize(num_simulations).ok_or_else(|| {
        PricingError::method_error(
            "price_option_monte_carlo",
            &format!("num_simulations not representable as Decimal: {num_simulations}"),
        )
    })?;
    let mean_payoff = d_div(total_payoff, n_dec, "pricing::monte_carlo::mean")?;
    let avg_payoff = d_mul(discount_factor, mean_payoff, "pricing::monte_carlo::price")?;
    Ok(Positive::new_decimal(avg_payoff.abs()).unwrap_or(Positive::ZERO))
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::constants::ZERO;
    use crate::model::types::{OptionStyle, OptionType, Side};
    use crate::{ExpirationDate, assert_decimal_eq, f2du};
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::{Positive, pos_or_panic};
    use rust_decimal::MathematicalOps;
    use rust_decimal_macros::dec;

    fn create_test_option() -> Options {
        Options {
            option_type: OptionType::European,
            side: Side::Long,
            underlying_symbol: "TEST".to_string(),
            strike_price: Positive::HUNDRED,
            expiration_date: ExpirationDate::Days(DAYS_IN_A_YEAR), // 1 year
            implied_volatility: pos_or_panic!(0.2),
            quantity: Positive::ONE,
            underlying_price: Positive::HUNDRED,
            risk_free_rate: dec!(0.05),
            option_style: OptionStyle::Call,
            dividend_yield: Positive::ZERO,
            exotic_params: None,
        }
    }

    #[test]
    fn test_monte_carlo_option_pricing_at_the_money() {
        let option = create_test_option();
        let price = monte_carlo_option_pricing(&option, crate::nz!(252), crate::nz!(1000)).unwrap();
        // The price should be close to the Black-Scholes price for these parameters
        let expected_price = dec!(9.100); // Calculated using Black-Scholes
        assert_decimal_eq!(price, expected_price, dec!(5));
    }

    #[test]
    fn test_monte_carlo_option_pricing_zero_volatility() {
        let mut option = create_test_option();
        option.implied_volatility = Positive::ZERO;
        let price = monte_carlo_option_pricing(&option, crate::nz!(25), crate::nz!(100)).unwrap();
        let expected_price = f64::max(
            (option.underlying_price - option.strike_price * (-option.risk_free_rate).exp()).into(),
            ZERO,
        );
        assert_decimal_eq!(price, f2du!(expected_price).unwrap(), dec!(0.1));
    }

    #[test]
    fn test_monte_carlo_option_pricing_high_volatility() {
        let mut option = create_test_option();
        option.implied_volatility = pos_or_panic!(0.5);
        let price = monte_carlo_option_pricing(&option, crate::nz!(252), crate::nz!(100)).unwrap();
        // The price should be higher with higher volatility
        assert!(price > dec!(10.0));
    }

    #[test]
    fn test_monte_carlo_option_pricing_short_expiration() {
        let mut option = create_test_option();
        option.expiration_date = ExpirationDate::Days(pos_or_panic!(30.0)); // 30 days
        let price = monte_carlo_option_pricing(&option, crate::nz!(30), crate::nz!(100)).unwrap();
        // The price should be lower for a shorter expiration
        assert!(price < dec!(5.0));
    }

    #[test]
    fn test_monte_carlo_option_pricing_consistency() {
        let option = create_test_option();
        let _price1 =
            monte_carlo_option_pricing(&option, crate::nz!(100), crate::nz!(100)).unwrap();
        let _price2 =
            monte_carlo_option_pricing(&option, crate::nz!(100), crate::nz!(100)).unwrap();
        // Two runs should produce similar results
        // assert_relative_eq!(price1, price2,  0.05);
    }
}

#[cfg(test)]
mod tests_price_option_monte_carlo {
    use super::*;
    use crate::chains::generator_positive;
    use crate::model::utils::create_sample_option;
    use crate::simulation::simulator::Simulator;
    use crate::simulation::steps::{Step, Xstep, Ystep};
    use crate::simulation::{WalkParams, WalkType, WalkTypeAble};
    use crate::utils::TimeFrame;
    use crate::utils::time::convert_time_frame;
    #[cfg(feature = "static_export")]
    use crate::visualization::Graph;
    use crate::{ExpirationDate, OptionStyle, Side};
    use positive::{Positive, assert_pos_relative_eq, pos_or_panic};
    use rust_decimal_macros::dec;

    #[test]
    fn test_empty_prices_returns_zero() {
        // Arrange
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
        );
        let empty_prices = &[];

        // Act
        let result = price_option_monte_carlo(&option, empty_prices);

        // Assert
        assert!(result.is_ok());
        assert_eq!(result.unwrap(), Positive::ZERO);
    }

    #[test]
    fn test_call_option_pricing() {
        // Arrange
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
        );

        option.risk_free_rate = dec!(0.05);
        option.dividend_yield = pos_or_panic!(0.02);
        option.expiration_date = ExpirationDate::Days(pos_or_panic!(365.0));

        // Setup simulated prices and expected payoffs
        let prices = vec![
            pos_or_panic!(110.0),
            pos_or_panic!(90.0),
            pos_or_panic!(105.0),
        ];

        // Act
        let result = price_option_monte_carlo(&option, &prices);
        assert!(result.is_ok());
        assert_pos_relative_eq!(
            result.unwrap(),
            pos_or_panic!(4.85222766),
            pos_or_panic!(0.001)
        );
    }

    #[test]
    fn test_simulation() {
        #[derive(Clone)]
        struct TestWalker;
        impl WalkTypeAble<Positive, Positive> for TestWalker {}
        let walker = Box::new(TestWalker);
        let initial_price = pos_or_panic!(1000.0);
        let days = pos_or_panic!(365.0);
        let volatility = pos_or_panic!(0.2);
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            initial_price,
            Positive::ONE,
            initial_price,
            volatility,
        );
        option.risk_free_rate = dec!(0.05);
        option.dividend_yield = pos_or_panic!(0.02);
        option.expiration_date = ExpirationDate::Days(days);

        let init_step = Step {
            x: Xstep::new(Positive::ONE, TimeFrame::Day, ExpirationDate::Days(days)),
            y: Ystep::new(0, initial_price),
        };

        let dt = convert_time_frame(Positive::ONE, &TimeFrame::Day, &TimeFrame::Year);
        let walk_params = WalkParams {
            size: 365,
            init_step,
            walk_type: WalkType::Custom {
                dt,
                drift: dec!(0.02),
                volatility,
                vov: pos_or_panic!(0.01),
                vol_speed: Default::default(),
                vol_mean: pos_or_panic!(0.2),
            },
            walker,
        };

        let Ok(simulator) = Simulator::new(
            "Test Simulator".to_string(),
            100,
            &walk_params,
            generator_positive,
        ) else {
            panic!("simulator setup failed");
        };

        #[cfg(feature = "static_export")]
        simulator
            .write_html("Draws/Simulation/simulator_test_montecarlo.html".as_ref())
            .unwrap();
        let get_last_positive_values = simulator.get_last_positive_values();

        let result = price_option_monte_carlo(&option, &get_last_positive_values);
        assert!(result.is_ok());

        let bs = option.calculate_price_black_scholes().unwrap();
        assert_pos_relative_eq!(
            result.unwrap(),
            Positive::new_decimal(bs).unwrap(),
            pos_or_panic!(10.0)
        );
    }

    // #[test]
    // fn test_put_option_pricing() {
    //     // Arrange
    //     let mut mock_option = MockOptions::new();
    //
    //     // Setup expiration date (1 year from now)
    //     let expiration = Utc::now().date_naive() + chrono::Duration::days(365);
    //     mock_option.risk_free_rate = 0.05;
    //     mock_option.dividend_yield = 0.02;
    //     mock_option.expiration_date = expiration;
    //
    //     // Setup simulated prices and expected payoffs
    //     let prices = vec![
    //         Positive::from_f64(90.0),
    //         Positive::from_f64(110.0),
    //         Positive::from_f64(95.0),
    //     ];
    //
    //     // For a put option with strike 100:
    //     // Payoffs would be: 10.0, 0.0, 5.0
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(90.0)))
    //         .returning(|_| Ok(Decimal::from_str("10.0").unwrap()));
    //
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(110.0)))
    //         .returning(|_| Ok(Decimal::from_str("0.0").unwrap()));
    //
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(95.0)))
    //         .returning(|_| Ok(Decimal::from_str("5.0").unwrap()));
    //
    //     // Act
    //     let result = price_option_monte_carlo(&mock_option, &prices);
    //
    //     // Assert
    //     assert!(result.is_ok());
    //
    //     // Expected: avg payoff = (10 + 0 + 5)/3 = 5.0
    //     // Discount factor = exp(-(0.05-0.02)*1) = exp(-0.03) ≈ 0.9704
    //     // Expected price = 5.0 * 0.9704 ≈ 4.852
    //     let expected = Positive::from_f64(4.852);
    //
    //     // Using approximate comparison due to floating-point calculations
    //     let diff = (result.unwrap().0 - expected.0).abs();
    //     assert!(diff < Decimal::from_str("0.001").unwrap(),
    //             "Expected close to {}, got {}", expected.0, result.unwrap().0);
    // }
    //
    // #[test]
    // fn test_expiration_date_error_handling() {
    //     // Arrange
    //     let mut mock_option = MockOptions::new();
    //
    //     // Setup with a problematic expiration date that will cause an error
    //     let prices = vec![Positive::from_f64(100.0)];
    //
    //     mock_option.expect_payoff_at_price()
    //         .returning(|_| Ok(Decimal::from_str("10.0").unwrap()));
    //
    //     // Make the expiration date calculation fail
    //     mock_option.expiration_date = NaiveDate::from_ymd_opt(9999, 12, 31).unwrap(); // Far future date that might cause issues
    //
    //     // Create a custom implementation for get_years that returns an error
    //     impl ExpirationDate for MockOptions {
    //         fn get_years(&self) -> Result<f64, ExpirationDateError> {
    //             Err("Invalid expiration date".unwrap_or(Positive::ZERO))
    //         }
    //     }
    //
    //     // Act
    //     let result = price_option_monte_carlo(&mock_option, &prices);
    //
    //     // Assert
    //     assert!(result.is_err());
    //     assert_eq!(result.unwrap_err().to_string(), "Invalid expiration date");
    // }
    //
    // #[test]
    // fn test_payoff_calculation_error_handling() {
    //     // Arrange
    //     let mut mock_option = MockOptions::new();
    //
    //     // Setup expiration date (1 year from now)
    //     let expiration = Utc::now().date_naive() + chrono::Duration::days(365);
    //     mock_option.risk_free_rate = 0.05;
    //     mock_option.dividend_yield = 0.02;
    //     mock_option.expiration_date = expiration;
    //
    //     let prices = vec![Positive::from_f64(100.0)];
    //
    //     // Make the payoff calculation fail
    //     mock_option.expect_payoff_at_price()
    //         .returning(|_| Err("Payoff calculation failed".into()));
    //
    //     // Act & Assert
    //     // We expect a panic since the function uses expect() on the payoff calculation
    //     std::panic::catch_unwind(|| {
    //         price_option_monte_carlo(&mock_option, &prices)
    //     }).expect_err("Expected a panic but none occurred");
    // }
    //
    // #[test]
    // fn test_large_number_of_simulations() {
    //     // Arrange
    //     let mut mock_option = MockOptions::new();
    //
    //     // Setup expiration date (1 year from now)
    //     let expiration = Utc::now().date_naive() + chrono::Duration::days(365);
    //     mock_option.risk_free_rate = 0.05;
    //     mock_option.dividend_yield = 0.02;
    //     mock_option.expiration_date = expiration;
    //
    //     // Create a large number of simulations
    //     let num_simulations = 1000;
    //     let mut prices = Vec::with_capacity(num_simulations);
    //     for _ in 0..num_simulations {
    //         prices.push(Positive::from_f64(100.0));
    //     }
    //
    //     // All payoffs are 5.0 in this test
    //     mock_option.expect_payoff_at_price()
    //         .returning(|_| Ok(Decimal::from_str("5.0").unwrap()));
    //
    //     // Act
    //     let result = price_option_monte_carlo(&mock_option, &prices);
    //
    //     // Assert
    //     assert!(result.is_ok());
    //
    //     // Expected: avg payoff = 5.0
    //     // Discount factor = exp(-(0.05-0.02)*1) = exp(-0.03) ≈ 0.9704
    //     // Expected price = 5.0 * 0.9704 ≈ 4.852
    //     let expected = Positive::from_f64(4.852);
    //
    //     // Using approximate comparison due to floating-point calculations
    //     let diff = (result.unwrap().0 - expected.0).abs();
    //     assert!(diff < Decimal::from_str("0.001").unwrap(),
    //             "Expected close to {}, got {}", expected.0, result.unwrap().0);
    // }
    //
    // #[test]
    // fn test_negative_payoffs_handled_correctly() {
    //     // Arrange
    //     let mut mock_option = MockOptions::new();
    //
    //     // Setup expiration date (1 year from now)
    //     let expiration = Utc::now().date_naive() + chrono::Duration::days(365);
    //     mock_option.risk_free_rate = 0.05;
    //     mock_option.dividend_yield = 0.02;
    //     mock_option.expiration_date = expiration;
    //
    //     let prices = vec![
    //         Positive::from_f64(100.0),
    //         Positive::from_f64(110.0),
    //         Positive::from_f64(90.0),
    //     ];
    //
    //     // Some payoffs are negative (unusual but possible in some exotic options)
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(100.0)))
    //         .returning(|_| Ok(Decimal::from_str("-5.0").unwrap()));
    //
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(110.0)))
    //         .returning(|_| Ok(Decimal::from_str("10.0").unwrap()));
    //
    //     mock_option.expect_payoff_at_price()
    //         .with(eq(Positive::from_f64(90.0)))
    //         .returning(|_| Ok(Decimal::from_str("5.0").unwrap()));
    //
    //     // Act
    //     let result = price_option_monte_carlo(&mock_option, &prices);
    //
    //     // Assert
    //     assert!(result.is_ok());
    //
    //     // Expected: avg payoff = (-5 + 10 + 5)/3 ≈ 3.333
    //     // Discount factor = exp(-(0.05-0.02)*1) = exp(-0.03) ≈ 0.9704
    //     // Expected price = 3.333 * 0.9704 ≈ 3.235
    //     // Since we take the absolute value, it should be positive
    //     let expected = Positive::from_f64(3.235);
    //
    //     // Using approximate comparison due to floating-point calculations
    //     let diff = (result.unwrap().0 - expected.0).abs();
    //     assert!(diff < Decimal::from_str("0.001").unwrap(),
    //             "Expected close to {}, got {}", expected.0, result.unwrap().0);
    // }
}