optionstratlib 0.15.3

/******************************************************************************
   Author: Joaquín Béjar García
   Email: jb@taunais.com
   Date: 11/8/24
******************************************************************************/
use crate::constants::{TRADING_DAYS, ZERO};
use crate::error::greeks::GreeksError;
use crate::greeks::utils::{big_n, d1, d2, n};
use crate::model::types::{OptionStyle, OptionType};
use crate::{Options, Side};
use positive::Positive;
use pretty_simple_display::{DebugPretty, DisplaySimple};
use rust_decimal::{Decimal, MathematicalOps};
use serde::{Deserialize, Serialize};
use utoipa::ToSchema;

/// Represents a complete set of option Greeks, which measure the sensitivity of an option's
/// price to various market factors.
///
/// Option Greeks are essential metrics in options trading and risk management, each quantifying
/// how the theoretical value of an option changes with respect to different parameters.
///
/// ## Fields
///
/// Each field represents a specific Greek measure:
///
/// * `delta`: Measures the rate of change in the option price relative to changes in the underlying asset price
/// * `gamma`: Measures the rate of change of delta in relation to changes in the underlying asset price
/// * `theta`: Measures the rate of change in the option price with respect to time decay (time sensitivity)
/// * `vega`: Measures the rate of change in the option price with respect to changes in implied volatility
/// * `rho`: Measures the rate of change in the option price with respect to the risk-free interest rate
/// * `rho_d`: Measures the rate of change in the option price with respect to the dividend yield
/// * `alpha`: Represents a measure of an option's excess return relative to what would be predicted by models
/// * `vanna`: Measures the rate of change of delta in relation to changes in implied volatility
/// * `vomma`: Measures the rate of change of vega in relation to changes in implied volatility
/// * `veta`: Measures the rate of change of vega in relation to changes in time
/// * `charm`: Measures the rate of change of delta in relation to changes in time
/// * `color`: Measures the rate of change of gamma in relation to changes in time
///
/// These metrics help traders understand and manage the various dimensions of risk in option positions.
#[derive(PartialEq, Clone, DebugPretty, DisplaySimple, ToSchema, Serialize)]
pub struct Greek {
    /// Measures sensitivity to changes in the underlying asset's price (first derivative)
    pub delta: Decimal,
    /// Measures the rate of change in delta (second derivative of the option price)
    pub gamma: Decimal,
    /// Measures the time decay of an option's value (sensitivity to the passage of time)
    pub theta: Decimal,
    /// Measures sensitivity to changes in implied volatility
    pub vega: Decimal,
    /// Measures sensitivity to changes in the risk-free interest rate
    pub rho: Decimal,
    /// Measures sensitivity to changes in the dividend yield
    pub rho_d: Decimal,
    /// Measures the option's theoretical value not explained by other Greeks
    pub alpha: Decimal,
    /// Measures the rate of change of delta in relation to changes in implied volatility
    pub vanna: Decimal,
    /// Measures the rate of change of vega in relation to changes in implied volatility
    pub vomma: Decimal,
    /// Measures the rate of change of vega in relation to changes in time
    pub veta: Decimal,
    /// Measures the rate of change of delta in relation to changes in time
    pub charm: Decimal,
    /// Measures the rate of change of gamma in relation to changes in time
    pub color: Decimal,
}

/// A struct representing a snapshot of the Greeks, financial measures used to assess risk and
/// sensitivity of derivative instruments such as options.
///
/// The Greeks provide insights into how various factors, such as price movement, time decay,
/// or volatility, affect the theoretical value of derivatives. This struct supports serialization
/// and deserialization for storage or communication purposes, and implements common traits like
/// `Debug`, `Clone`, and `PartialEq`.
#[derive(DebugPretty, DisplaySimple, Clone, PartialEq, Serialize, Deserialize, ToSchema)]
pub struct GreeksSnapshot {
    /// Measures sensitivity to changes in the underlying asset's price (first derivative)
    pub delta: Decimal,
    /// Measures the rate of change in delta (second derivative of the option price)
    pub gamma: Decimal,
    /// Measures the time decay of an option's value (sensitivity to the passage of time)
    pub theta: Decimal,
    /// Measures sensitivity to changes in implied volatility
    pub vega: Decimal,
    /// Measures sensitivity to changes in the risk-free interest rate
    pub rho: Option<Decimal>,
    /// Measures sensitivity to changes in the dividend yield
    pub rho_d: Option<Decimal>,
    /// Measures the option's theoretical value not explained by other Greeks
    pub alpha: Option<Decimal>,
    /// Measures the rate of change of delta in relation to changes in implied volatility
    pub vanna: Decimal,
    /// Measures the rate of change of vega in relation to changes in implied volatility
    pub vomma: Decimal,
    /// Measures the rate of change of vega in relation to changes in time
    pub veta: Decimal,
    /// Measures the rate of change of delta in relation to changes in time
    pub charm: Decimal,
    /// Measures the rate of change of gamma in relation to changes in time
    pub color: Decimal,
}

/// Trait that provides option Greeks calculation functionality for financial instruments.
///
/// The `Greeks` trait enables implementing types to calculate option sensitivity metrics
/// (Greeks) across multiple option positions. Any type that can provide access to a collection
/// of options can implement this trait to gain the ability to calculate aggregate Greek values.
///
/// This trait uses a composition approach where implementation only requires defining the
/// `get_options()` method, while default implementations for all Greek calculations are provided.
///
/// # Greek Calculations
///
/// The trait provides calculations for:
/// - Delta: Sensitivity to changes in the underlying asset's price
/// - Gamma: Rate of change of delta (acceleration of price movement)
/// - Theta: Time decay of option value
/// - Vega: Sensitivity to changes in volatility
/// - Rho: Sensitivity to changes in interest rates
/// - Rho_d: Sensitivity to changes in dividend yield
/// - Alpha: Ratio between gamma and theta
/// - Vanna: Rate of change of delta in relation to changes in implied volatility
/// - Vomma: Rate of change of vega in relation to changes in implied volatility
/// - Veta: Rate of change of vega in relation to changes in time
/// - Charm: Rate of change of delta in relation to changes in time
/// - Color: Rate of change of gamma in relation to changes in time
///
/// # Usage
///
/// Implementers only need to provide the `get_options()` method which returns a vector of
/// references to option contracts. The trait will handle aggregating the Greek values across
/// all options in the collection.
///
/// # Errors
///
/// Methods return `Result<T, GreeksError>` to handle various calculation errors that may
/// occur during Greek computations.
pub trait Greeks {
    /// Returns a vector of references to the option contracts for which Greeks will be calculated.
    ///
    /// This is the only method that must be implemented by types adopting this trait.
    /// All other methods have default implementations based on this method.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if there is an issue retrieving the options.
    fn get_options(&self) -> Result<Vec<&Options>, GreeksError>;

    /// Calculates and returns all Greeks as a single `Greek` struct.
    ///
    /// This method provides a convenient way to obtain all Greek values at once.
    /// It calls each individual Greek calculation method and compiles the results.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if any individual Greek calculation fails.
    fn greeks(&self) -> Result<Greek, GreeksError> {
        let delta = self.delta()?;
        let gamma = self.gamma()?;
        let theta = self.theta()?;
        let vega = self.vega()?;
        let rho = self.rho()?;
        let rho_d = self.rho_d()?;
        let alpha = self.alpha()?;
        let vanna = self.vanna()?;
        let vomma = self.vomma()?;
        let veta = self.veta()?;
        let charm = self.charm()?;
        let color = self.color()?;
        Ok(Greek {
            delta,
            gamma,
            theta,
            vega,
            rho,
            rho_d,
            alpha,
            vanna,
            vomma,
            veta,
            charm,
            color,
        })
    }

    /// Calculates the aggregate delta value for all options.
    ///
    /// Delta measures the rate of change in an option's price with respect to
    /// changes in the underlying asset's price.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or delta calculation fails.
    fn delta(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut delta_value = Decimal::ZERO;
        for option in options {
            delta_value += delta(option)?;
        }
        Ok(delta_value)
    }

    /// Calculates the aggregate gamma value for all options.
    ///
    /// Gamma measures the rate of change of delta with respect to
    /// changes in the underlying asset's price.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or gamma calculation fails.
    fn gamma(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut gamma_value = Decimal::ZERO;
        for option in options {
            gamma_value += gamma(option)?;
        }
        Ok(gamma_value)
    }

    /// Calculates the aggregate theta value for all options.
    ///
    /// Theta measures the rate of change of the option price with respect to time,
    /// also known as time decay.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or theta calculation fails.
    fn theta(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut theta_value = Decimal::ZERO;
        for option in options {
            theta_value += theta(option)?;
        }
        Ok(theta_value)
    }

    /// Calculates the aggregate vega value for all options.
    ///
    /// Vega measures the sensitivity of the option price to changes in
    /// the volatility of the underlying asset.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or vega calculation fails.
    fn vega(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut vega_value = Decimal::ZERO;
        for option in options {
            vega_value += vega(option)?;
        }
        Ok(vega_value)
    }

    /// Calculates the aggregate rho value for all options.
    ///
    /// Rho measures the sensitivity of the option price to changes in
    /// the risk-free interest rate.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or rho calculation fails.
    fn rho(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut rho_value = Decimal::ZERO;
        for option in options {
            rho_value += rho(option)?;
        }
        Ok(rho_value)
    }

    /// Calculates the aggregate rho_d value for all options.
    ///
    /// Rho_d measures the sensitivity of the option price to changes in
    /// the dividend yield of the underlying asset.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or rho_d calculation fails.
    fn rho_d(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut rho_d_value = Decimal::ZERO;
        for option in options {
            rho_d_value += rho_d(option)?;
        }
        Ok(rho_d_value)
    }

    /// Calculates the aggregate alpha value for all options.
    ///
    /// Alpha represents the ratio between gamma and theta, providing insight into
    /// the option's risk/reward efficiency with respect to time decay.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or alpha calculation fails.
    fn alpha(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut alpha_value = Decimal::ZERO;
        for option in options {
            alpha_value += alpha(option)?;
        }
        Ok(alpha_value)
    }

    /// Calculates the aggregate vanna value for all options.
    ///
    /// Vanna measures the sensitivity of the option delta to changes in
    /// the volatility of the underlying asset
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or vanna calculation fails.
    fn vanna(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut vanna_value = Decimal::ZERO;
        for option in options {
            vanna_value += vanna(option)?;
        }
        Ok(vanna_value)
    }

    /// Calculates the aggregate vomma value for all options.
    ///
    /// Vomma measures the sensitivity of the option vega to changes in
    /// the volatility of the underlying asset
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or vomma calculation fails.
    fn vomma(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut vomma_value = Decimal::ZERO;
        for option in options {
            vomma_value += vomma(option)?;
        }
        Ok(vomma_value)
    }

    /// Calculates the aggregate veta value for all options.
    ///
    /// Veta measures the sensitivity of the option vega in relation
    /// to changes in time
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or veta calculation fails.
    fn veta(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut veta_value = Decimal::ZERO;
        for option in options {
            veta_value += veta(option)?;
        }
        Ok(veta_value)
    }

    /// Calculates the aggregate charm value for all options.
    ///
    /// Charm measures the rate of change of the option delta with respect to time,
    /// also known as delta decay.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or charm calculation fails.
    fn charm(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut charm_value = Decimal::ZERO;
        for option in options {
            charm_value += charm(option)?;
        }
        Ok(charm_value)
    }

    /// Calculates the aggregate color value for all options.
    ///
    /// Color measures the rate of change of the option gamma with respect to time,
    /// also known as gamma decay.
    ///
    /// # Errors
    ///
    /// Returns a `GreeksError` if the options can't be retrieved or color calculation fails.
    fn color(&self) -> Result<Decimal, GreeksError> {
        let options = self.get_options()?;
        let mut color_value = Decimal::ZERO;
        for option in options {
            color_value += color(option)?;
        }
        Ok(color_value)
    }
}

/// Calculates the delta of an option.
///
/// The delta measures the sensitivity of an option's price to changes in the price of the
/// underlying asset. It is calculated differently for call and put options. For options
/// with zero implied volatility, the delta is determined based on whether the option is
/// in-the-money or out-of-the-money.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing all the relevant parameters for the calculation:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The annualized risk-free interest rate.
///   - `expiration_date`: The time to expiration of the option, in years.
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///   - `option_style`: The style of the option (Call or Put).
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated delta value.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculations fail.
///
/// # Calculation Details
///
/// - If `implied_volatility == 0`, the delta is determined based on whether the option is
///   in-the-money or out-of-the-money:
///   - Call Option:
///     - In-the-money: Delta = `sign`
///     - Out-of-the-money: Delta = 0
///   - Put Option:
///     - In-the-money: Delta = `-sign`
///     - Out-of-the-money: Delta = 0
/// - For options with non-zero implied volatility, the delta is calculated as:
///   - Call Option:
///     \[ \Delta_{\text{call}} = \text{sign} \cdot N(d1) \cdot e^{-qT} \]
///   - Put Option:
///     \[ \Delta_{\text{put}} = \text{sign} \cdot (N(d1) - 1) \cdot e^{-qT} \]
///     Where:
///     - \(N(d1)\): The cumulative distribution function (CDF) of the standard normal distribution evaluated at \(d1\).
///     - \(q\): The dividend yield.
///     - \(T\): Time to expiration.
///
/// - The delta is adjusted by multiplying it by the option quantity.
///
/// # Errors
///
/// - `GreeksError`: If the calculation of \(d1\) or the standard normal CDF (`big_n`) fails.
///
/// # Example
///
/// ```rust
/// use rust_decimal::Decimal;
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::constants::ZERO;
/// use optionstratlib::greeks::delta;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic, Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: Positive::ZERO,
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "AAPL".to_string(),
///     exotic_params: None,
/// };
///
/// match delta(&option) {
///     Ok(result) => info!("Delta: {}", result),
///     Err(e) => error!("Error calculating delta: {:?}", e),
/// }
/// ```
pub fn delta(option: &Options) -> Result<Decimal, GreeksError> {
    if !matches!(option.option_type, OptionType::European) {
        return crate::greeks::numerical::numerical_delta(option);
    }
    let expiration_date = option.expiration_date.get_years()?;

    // For an option when the time to expiration is zero (i.e., at the moment of expiration),
    // the delta takes discrete values based solely on whether the option is In-The-Money (ITM) or
    // Out-of-The-Money (OTM):
    //
    // For a Call option:
    //
    // - **Delta = 1.0** if ITM (underlying price > strike price)
    // - **Delta = 0.0** if OTM (underlying price < strike price)
    //
    // For a Put option:
    //
    // - **Delta = -1.0** if ITM (underlying price < strike price)
    // - **Delta = 0.0** if OTM (underlying price > strike price)
    //
    // In both cases, when the underlying price is exactly equal to the strike price (At-The-Money,
    // ATM), technically, the delta would be **0.5 for Calls** and **-0.5 for Puts**, although this
    // scenario is less common in practice.
    //
    // This happens because at expiration, the option effectively becomes a direct position in the
    // underlying asset (**delta = 1 or -1**) if it is ITM, or has no value (**delta = 0**) if it is OTM.
    if expiration_date == Decimal::ZERO {
        return match (
            &option.option_style,
            &option.side,
            &option.strike_price,
            &option.underlying_price,
        ) {
            // Call Options
            (OptionStyle::Call, Side::Long, strike, price) if price > strike => Ok(Decimal::ONE),
            (OptionStyle::Call, Side::Long, _, _) => Ok(Decimal::ZERO),
            (OptionStyle::Call, Side::Short, strike, price) if price > strike => Ok(-Decimal::ONE),
            (OptionStyle::Call, Side::Short, _, _) => Ok(Decimal::ZERO),

            // Put Options
            (OptionStyle::Put, Side::Long, strike, price) if price < strike => Ok(-Decimal::ONE),
            (OptionStyle::Put, Side::Long, _, _) => Ok(Decimal::ZERO),
            (OptionStyle::Put, Side::Short, strike, price) if price < strike => Ok(Decimal::ONE),
            (OptionStyle::Put, Side::Short, _, _) => Ok(Decimal::ZERO),
        };
    }

    let dividend_yield: Positive = option.dividend_yield;

    let sign = if option.is_long() {
        Decimal::ONE
    } else {
        Decimal::NEGATIVE_ONE
    };
    if option.implied_volatility == ZERO {
        return match option.option_style {
            OptionStyle::Call => {
                if option.underlying_price >= option.strike_price {
                    Ok(sign) // Delta is 1 for Call in-the-money
                } else {
                    Ok(Decimal::ZERO) // Delta is 0 for Call out-of-the-money
                }
            }
            OptionStyle::Put => {
                if option.underlying_price <= option.strike_price {
                    Ok(sign * Decimal::NEGATIVE_ONE) // Delta is -1 for Put in-the-money
                } else {
                    Ok(Decimal::ZERO) // Delta is 0 for Put out-of-the-money
                }
            }
        };
    }

    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;

    let div_date = (-expiration_date.to_dec() * dividend_yield).exp();
    let delta = match option.option_style {
        OptionStyle::Call => sign * big_n(d1)? * div_date,
        OptionStyle::Put => sign * (big_n(d1)? - Decimal::ONE) * div_date,
    };
    let delta: Decimal = delta.clamp(Decimal::NEGATIVE_ONE, Decimal::ONE);
    let quantity: Decimal = option.quantity.into();
    Ok(delta * quantity)
}

/// Computes the gamma of an option.
///
/// Gamma measures the rate of change of the option's delta with respect to changes in the underlying
/// asset's price. It is a second-order derivative of the option price and provides insight into the
/// sensitivity of delta to movements in the underlying price.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years.
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated gamma value.
/// - `Err(GreeksError)`: Returns an error if the computation of `d1` or the probability density function `n(d1)` fails.
///
/// # Calculation
///
/// Gamma is calculated using the formula:
///
/// ```math
/// \Gamma = \frac{e^{-qT} \cdot N'(d1)}{S \cdot \sigma \cdot \sqrt{T}}
/// ```
///
/// Where:
/// - \(N'(d1)\): The standard normal probability density function (PDF) evaluated at \(d1\).
/// - \(S\): The price of the underlying asset.
/// - \(\sigma\): The implied volatility of the option.
/// - \(T\): The time to expiration in years.
/// - \(q\): The dividend yield of the underlying asset.
///
/// ### Steps:
/// 1. Compute \(d1\) using the `d1` function.
/// 2. Evaluate \(N'(d1)\) using the `n` function.
/// 3. Apply the gamma formula, accounting for the effect of the dividend yield \(e^{-qT}\).
/// 4. Multiply the result by the option's quantity.
///
/// # Edge Cases
///
/// - If the implied volatility (\(\sigma\)) is zero, gamma is returned as `0`.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::gamma;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match gamma(&option) {
///     Ok(result) => info!("Gamma: {}", result),
///     Err(e) => error!("Error calculating gamma: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - This function assumes that the dividend yield \(q\) and the time to expiration \(T\) are
///   provided in consistent units.
/// - If the implied volatility or time to expiration is very small, the result may be close to 0,
///   as gamma becomes negligible in those cases.
pub fn gamma(option: &Options) -> Result<Decimal, GreeksError> {
    if !matches!(option.option_type, OptionType::European) {
        return crate::greeks::numerical::numerical_gamma(option);
    }
    if option.implied_volatility == ZERO {
        return Ok(Decimal::ZERO);
    }
    let expiration_date: Positive = option.expiration_date.get_years()?;
    if expiration_date == Decimal::ZERO {
        // At expiration, gamma is 0 for all cases
        return Ok(Decimal::ZERO);
    }

    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;

    let dividend_yield: Decimal = option.dividend_yield.into();
    let underlying_price: Decimal = option.underlying_price.into();
    let implied_volatility: Positive = option.implied_volatility;

    let gamma: Decimal = (expiration_date.to_dec() * -dividend_yield).exp() * n(d1)?
        / (underlying_price * implied_volatility * expiration_date.sqrt().to_dec());

    let quantity: Decimal = option.quantity.into();
    Ok(gamma * quantity)
}

/// Computes the Theta of an option.
///
/// Theta measures the sensitivity of the option's price to time decay, indicating the rate
/// at which the value of the option decreases as the expiration date approaches. This is
/// particularly important in options trading, as Theta reflects the "time decay" of the
/// option's extrinsic value.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `option_style`: The style of the option (Call or Put).
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated Theta value for the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d1`, `d2`, or `n`).
///
/// # Formula
///
/// The Theta is calculated using the Black-Scholes model. The formula differs for call and put options:
///
/// **Call Options:**
///
/// ```math
/// \Theta_{\text{call}} =
/// -\frac{S \cdot \sigma \cdot e^{-qT} \cdot n(d1)}{2 \sqrt{T}}
/// - r \cdot K \cdot e^{-rT} \cdot N(d2)
/// + q \cdot S \cdot e^{-qT} \cdot N(d1)
/// ```
///
/// **Put Options:**
///
/// ```math
/// \Theta_{\text{put}} =
/// -\frac{S \cdot \sigma \cdot e^{-qT} \cdot n(d1)}{2 \sqrt{T}}
/// + r \cdot K \cdot e^{-rT} \cdot N(-d2)
/// - q \cdot S \cdot e^{-qT} \cdot N(-d1)
/// ```
///
/// Where:
/// - \( S \): Underlying price
/// - \( \sigma \): Implied volatility
/// - \( T \): Time to expiration (in years)
/// - \( r \): Risk-free rate
/// - \( q \): Dividend yield
/// - \( K \): Strike price
/// - \( N(d1) \): Cumulative distribution function (CDF) of the standard normal distribution at \( d1 \).
/// - \( n(d1) \): Probability density function (PDF) of the standard normal distribution at \( d1 \).
///
/// # Calculation Steps
/// 1. Compute \( d1 \) and \( d2 \) using the `d1` and `d2` functions.
/// 2. Calculate the common term:
///    ```math
///    \text{common\_term} = -\frac{S \cdot \sigma \cdot e^{-qT} \cdot n(d1)}{2 \sqrt{T}}
///    ```
/// 3. Apply the corresponding formula for Call or Put options, accounting for the effect of
///    dividends (\( e^{-qT} \)) and risk-free rate (\( e^{-rT} \)).
/// 4. Multiply the resulting Theta by the quantity of options.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::theta;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match theta(&option) {
///     Ok(result) => info!("Theta: {}", result),
///     Err(e) => error!("Error calculating Theta: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - A positive Theta means the option gains value as time passes (rare and usually for short positions).
/// - A negative Theta is typical for long positions, as the option loses extrinsic value over time.
/// - If the implied volatility is zero, Theta may be close to zero for far-out-of-the-money options.
pub fn theta(option: &Options) -> Result<Decimal, GreeksError> {
    let t = option.expiration_date.get_years()?;
    if t == Decimal::ZERO {
        return Ok(Decimal::ZERO);
    }

    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        t,
        option.implied_volatility,
    )?;
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        t,
        option.implied_volatility,
    )?;

    let s = option.underlying_price.to_dec();
    let k = option.strike_price.to_dec();
    let r = option.risk_free_rate;
    let q = option.dividend_yield.to_dec();
    let sigma = option.implied_volatility.to_dec();

    // Common term using n
    let common_term = -(s * n(d1)? * sigma) / (Decimal::TWO * t.sqrt());

    // Pre-calculate discount factors
    let exp_minus_rt = (-r * t).exp();
    let exp_minus_qt = (-q * t).exp();

    let theta = match option.option_style {
        OptionStyle::Call => {
            common_term - r * k * exp_minus_rt * big_n(d2)? + q * s * exp_minus_qt * big_n(d1)?
        }
        OptionStyle::Put => {
            common_term + r * k * exp_minus_rt * big_n(-d2)? - q * s * exp_minus_qt * big_n(-d1)?
        }
    };

    // Adjust for quantity and convert to daily value
    Ok((theta * option.quantity.to_dec()) / Decimal::from(365))
}

/// Computes the vega of an option.
///
/// Vega measures the sensitivity of the option's price to changes in the implied volatility
/// of the underlying asset. It quantifies the expected change in the option's price for a
/// 1% change in the implied volatility. Vega is particularly important for understanding
/// how an option's value is affected by market conditions that alter volatility.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the necessary parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The annualized risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///   - `option_style`: The style of the option (e.g., European).
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed vega value of the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d1` or `big_n`).
///
/// # Formula
///
/// Vega is computed using the Black-Scholes model formula:
///
/// ```math
/// \text{Vega} = S \cdot e^{-qT} \cdot n(d1) \cdot \sqrt{T}
/// ```
///
/// Where:
/// - \( S \): The price of the underlying asset.
/// - \( q \): The dividend yield of the underlying asset.
/// - \( T \): Time to expiration in years.
/// - \( n(d1) \): The probability density function (PDF) of the standard normal distribution at \( d1 \).
/// - \( d1 \): A parameter calculated using the Black-Scholes model.
///
/// # Calculation Steps
///
/// 1. Compute \( d1 \) using the `d1` function.
/// 2. Calculate the exponential factor \( e^{-qT} \), which accounts for the effect of dividends.
/// 3. Evaluate \( n(d1) \), the PDF of the standard normal distribution at \( d1 \).
/// 4. Multiply the underlying price, the exponential factor, \( n(d1) \), and the square root of time to expiration.
/// 5. Multiply the result by the quantity of options to adjust for position size.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::vega;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
///
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match vega(&option) {
///     Ok(result) => info!("Vega: {}", result),
///     Err(e) => error!("Error calculating Vega: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - Vega is usually highest for at-the-money options and decreases as the option moves deeper
///   in-the-money or out-of-the-money.
/// - For shorter time to expiration, Vega is smaller as the sensitivity to volatility diminishes.
/// - A positive Vega indicates that an increase in implied volatility will increase the option's value.
pub fn vega(option: &Options) -> Result<Decimal, GreeksError> {
    let expiration_date: Positive = option.expiration_date.get_years()?;
    if expiration_date == Decimal::ZERO {
        // At expiration, volatility has no impact on option price
        return Ok(Decimal::ZERO);
    }
    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;

    let dividend_yield: Positive = option.dividend_yield;
    let underlying_price: Decimal = option.underlying_price.to_dec();

    let vega: Decimal = underlying_price
        * (-expiration_date.to_dec() * dividend_yield).exp()
        * n(d1)?
        * expiration_date.sqrt()
        / Decimal::ONE_HUNDRED; // percentage of change in volatility

    let quantity: Decimal = option.quantity.into();
    Ok(vega * quantity)
}

/// Computes the rho of an options contract.
///
/// Rho measures the sensitivity of the option's price to changes in the risk-free interest rate.
/// It quantifies the expected change in the option's price for a 1% change in the risk-free rate.
/// This metric is useful for understanding how interest rate fluctuations affect the value of
/// options contracts.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following fields:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The annualized risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `option_style`: The style of the option (`Call` or `Put`).
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed rho value for the options contract.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d2` or `big_n`).
///
/// # Formula
///
/// The rho is calculated differently for Call and Put options, as follows:
///
/// **Call Options:**
///
/// ```math
/// \rho_{\text{call}} = K \cdot T \cdot e^{-rT} \cdot N(d2)
/// ```
///
/// **Put Options:**
///
/// ```math
/// \rho_{\text{put}} = -K \cdot T \cdot e^{-rT} \cdot N(-d2)
/// ```
///
/// Where:
/// - \( K \): The strike price of the option.
/// - \( T \): The time to expiration (in years).
/// - \( r \): The risk-free interest rate.
/// - \( N(d2) \): The cumulative distribution function (CDF) of the standard normal distribution evaluated at \( d2 \).
/// - \( e^{-rT} \): The discount factor for the risk-free rate.
///
/// # Calculation Steps
///
/// 1. Compute \( d2 \) using the `d2` function.
/// 2. Calculate the discount factor \( e^{-rT} \).
/// 3. Evaluate \( N(d2) \) or \( N(-d2) \), depending on the option style.
/// 4. Multiply the strike price, time to expiration, discount factor, and \( N(d2) \) or \( N(-d2) \).
/// 5. Multiply the result by the option's quantity.
///
/// # Edge Cases
///
/// - If the discount factor (\( e^{-rT} \)) is zero, the rho is returned as zero.
/// - If \( N(d2) \) or \( N(-d2) \) is zero, the rho is returned as zero.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::rho;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
///
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match rho(&option) {
///     Ok(result) => info!("Rho: {}", result),
///     Err(e) => error!("Error calculating rho: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - Rho is typically higher for options with longer time to expiration, as they are more
///   sensitive to changes in the risk-free rate.
/// - Call options have positive rho values, as an increase in interest rates increases their value.
/// - Put options have negative rho values, as an increase in interest rates decreases their value.
pub fn rho(option: &Options) -> Result<Decimal, GreeksError> {
    // Get time to expiration first and validate
    let t = option.expiration_date.get_years()?;
    if t == Decimal::ZERO {
        return Ok(Decimal::ZERO);
    }

    // Use existing d2 function
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        t,
        option.implied_volatility,
    )?;

    let k = option.strike_price.to_dec();
    let r = option.risk_free_rate;

    // Calculate discount factor once
    let e_rt = (-r * t).exp();

    // Calculate base rho without sign
    let base_rho = k * t * e_rt;

    // Calculate final rho based on option type
    let rho = match option.option_style {
        OptionStyle::Call => {
            let n_d2 = big_n(d2)?;
            base_rho * n_d2
        }
        OptionStyle::Put => {
            let n_minus_d2 = big_n(-d2)?;
            -base_rho * n_minus_d2
        }
    };

    // Adjust for quantity and convert to basis points
    Ok((rho * option.quantity.to_dec()) / Decimal::from(100))
}

/// Computes the sensitivity of the option price to changes in the dividend yield (Rho_d).
///
/// This function calculates how the price of an option changes with respect to variations
/// in the dividend yield of the underlying asset. This metric, often referred to as "dividend rho",
/// is essential for understanding the impact of dividends on the option's value.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant fields:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///   - `option_style`: The style of the option (`Call` or `Put`).
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated dividend sensitivity (`Rho_d`) value for the options contract.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d1` or `big_n`).
///
/// # Formula
///
/// The dividend sensitivity is calculated differently for Call and Put options:
///
/// **Call Options:**
///
/// ```math
/// \rho_d^{\text{call}} = -T \cdot S \cdot e^{-qT} \cdot N(d1)
/// ```
///
/// **Put Options:**
///
/// ```math
/// \rho_d^{\text{put}} = T \cdot S \cdot e^{-qT} \cdot N(-d1)
/// ```
///
/// Where:
/// - \( T \): Time to expiration (in years).
/// - \( S \): Price of the underlying asset.
/// - \( q \): Dividend yield.
/// - \( N(d1) \): The cumulative distribution function (CDF) of the standard normal distribution evaluated at \( d1 \).
/// - \( d1 \): A parameter calculated using the Black-Scholes model.
///
/// # Calculation Steps
///
/// 1. Compute \( d1 \) using the `d1` function.
/// 2. Evaluate the exponential factor \( e^{-qT} \), which accounts for the dividend yield.
/// 3. Calculate \( N(d1) \) or \( N(-d1) \), depending on the option style.
/// 4. Use the appropriate formula for Call or Put options.
/// 5. Multiply the result by the option's quantity to adjust for position size.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::rho_d;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic, Positive};
///
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match rho_d(&option) {
///     Ok(result) => info!("Dividend Rho (Rho_d): {}", result),
///     Err(e) => error!("Error calculating Rho_d: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - **Call Options**: A higher dividend yield decreases the price of the call option,
///   leading to a negative dividend sensitivity.
/// - **Put Options**: A higher dividend yield increases the price of the put option,
///   leading to a positive dividend sensitivity.
/// - This calculation assumes that dividends are continuously compounded at the dividend yield rate.
/// - \( Rho_d \) is generally more significant for options with longer times to expiration.
pub fn rho_d(option: &Options) -> Result<Decimal, GreeksError> {
    let expiration_date: Positive = option.expiration_date.get_years()?;
    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let dividend_yield: Positive = option.dividend_yield;
    let underlying_price: Decimal = option.underlying_price.to_dec();

    let rhod = match option.option_style {
        OptionStyle::Call => {
            -expiration_date.to_dec()
                * underlying_price
                * (-expiration_date.to_dec() * dividend_yield).exp()
                * big_n(d1)?
        }
        OptionStyle::Put => {
            expiration_date.to_dec()
                * underlying_price
                * (-expiration_date.to_dec() * dividend_yield).exp()
                * big_n(-d1)?
        }
    };

    let quantity: Decimal = option.quantity.into();
    Ok(rhod * quantity / Decimal::from(100))
}

pub fn alpha(option: &Options) -> Result<Decimal, GreeksError> {
    let gamma = gamma(option)?;
    let theta = theta(option)?;
    match (gamma, theta) {
        (val, _) if val == Decimal::ZERO => Ok(Decimal::ZERO),
        (_, val) if val == Decimal::ZERO => Ok(Decimal::MAX),
        _ => Ok(gamma / theta),
    }
}

/// Computes the vanna of an option.
///
/// Vanna measures the rate of change of delta in relation to changes in implied volatility.
/// It is a second-order derivative of the option value and can be useful to help the trader
/// to anticipate changes to the effectiveness of a delta-hedge as volatility changes.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years.
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated vanna value.
/// - `Err(GreeksError)`: Returns an error if the computation of `d2` or the probability density function `n(d1)` fails.
///
/// # Calculation
///
/// Vanna is calculated using the formula:
///
/// ```math
/// \text{Vanna} = -e^{-qT} \cdot N'(d1) \cdot \frac {d2}{\sigma}
/// ```
///
/// Where:
/// - \(N'(d1)\): The standard normal probability density function (PDF) evaluated at \(d1\).
/// - \(\sigma\): The implied volatility of the option.
/// - \(T\): The time to expiration in years.
/// - \(q\): The dividend yield of the underlying asset.
///
/// ### Steps:
/// 1. Compute \(d1\) using the `d1` function.
/// 2. Compute \(d2\) using the `d2` function.
/// 3. Evaluate \(N'(d1)\) using the `n` function.
/// 4. Compute the effect of the dividend yield \(-e^{-qT}\).
/// 5. Apply the vanna formula and divide by the implied_volatility (\(\sigma\)).
/// 6. Multiply the result by the option's quantity.
///
/// # Edge Cases
///
/// - If the implied volatility (\(\sigma\)) is zero, vanna is returned as `0`.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::vanna;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match vanna(&option) {
///     Ok(result) => info!("Vanna: {}", result),
///     Err(e) => error!("Error calculating vanna: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - This function assumes that the dividend yield \(q\) and the time to expiration \(T\) are
///   provided in consistent units.
pub fn vanna(option: &Options) -> Result<Decimal, GreeksError> {
    if option.implied_volatility == ZERO {
        return Ok(Decimal::ZERO);
    }

    let expiration_date: Positive = option.expiration_date.get_years()?;
    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;

    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let dividend_yield: Decimal = option.dividend_yield.into();
    let implied_volatility: Positive = option.implied_volatility;
    let n_d1: Decimal = n(d1)?;
    let e_rt: Decimal = -(expiration_date.to_dec() * -dividend_yield).exp();

    let vanna: Decimal = e_rt * n_d1 * (d2 / implied_volatility);

    let quantity: Decimal = option.quantity.into();
    Ok(vanna * quantity)
}

/// Computes the vomma of an option.
///
/// Vomma (aka volga, vega convexity or DvegaDvol) measures the second order
/// sensitivity to volatility. Is the second derivative of the option value
/// with respect to the volatility. Stated in another way, vomma measures
/// the rate of change to vega as volatility changes.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the necessary parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The annualized risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///   - `option_style`: The style of the option (e.g., European).
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed vomma value of the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d1` or `d2`).
///
/// # Formula
///
/// Vomma is computed using the Black-Scholes model formula:
///
/// ```math
/// \text{Vomma} = {Vega} \cdot \frac{d1 \cdot d2}{\sigma}
/// ```
///
/// Where:
/// - \( Vega \): The option's Vega value.
/// - \( d1 \): A parameter calculated using the Black-Scholes model.
/// - \( d2 \): A parameter calculated using the Black-Scholes model.
/// - \(\sigma\): The implied volatility of the option.
///
/// # Calculation Steps
///
/// 1. Compute \( Vega \) using the `vega` function.
/// 2. Compute \( d1 \) using the `d1` function.
/// 3. Compute \( d2 \) using the `d2` function.
/// 4. Compute Vomma and multiply the result by the quantity of options to adjust for position size.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::vomma;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
///
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match vomma(&option) {
///     Ok(result) => info!("Vomma: {}", result),
///     Err(e) => error!("Error calculating Vomma: {:?}", e),
/// }
/// ```
/// # Notes
///
/// Options far out-of-the money have the highest Vomma.
/// If you are long options you typically want to have as high positive Vomma
/// as possible. If short options, you typically want negative Vomma.
/// Positive Vomma tells you that you will earn more for every percentage point
/// increase in volatility, and if implied volatility is falling you will lose
/// less and less.
/// If you think the implied volatility will be volatile in the short term
/// you should typically try to find options with high Vomma.
pub fn vomma(option: &Options) -> Result<Decimal, GreeksError> {
    let expiration_date: Positive = option.expiration_date.get_years()?;
    if expiration_date == Decimal::ZERO {
        // At expiration, volatility has no impact on option price
        return Ok(Decimal::ZERO);
    }
    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let vega = vega(option)?;
    let implied_volatility: Positive = option.implied_volatility;

    let vomma: Decimal = vega * (d1 * d2 / implied_volatility);

    let quantity: Decimal = option.quantity.into();
    Ok(vomma * quantity)
}

/// Computes the veta of an option.
///
/// Veta measures the rate of change to vega as time changes.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the necessary parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The annualized risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `quantity`: The quantity of the options.
///   - `option_style`: The style of the option (e.g., European).
///
/// # Returns
///
/// - `Ok(Decimal)`: The computed veta value of the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails (e.g., in `d1` or `d2`).
///
/// # Formula
///
/// Veta is computed using the Black-Scholes model formula:
///
/// ```math
/// \text{Veta} = {-Vega} \left [q+\frac {(r-q)d1}{\sigma\sqrt T} -\frac {1 + d1 d2}{2T} \right]
/// ```
///
/// Where:
/// - \( Vega \): The option's Vega value.
/// - \( d1 \): A parameter calculated using the Black-Scholes model.
/// - \( d2 \): A parameter calculated using the Black-Scholes model.
/// - \(\sigma\): The implied volatility of the option.
///
/// # Calculation Steps
///
/// 1. Compute \( Vega \) using the `vega` function.
/// 2. Compute \( d1 \) using the `d1` function.
/// 3. Compute \( d2 \) using the `d2` function.
/// 4. Compute Veta and multiply the result by the quantity of options to adjust for position size.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::veta;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
///
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match veta(&option) {
///     Ok(result) => info!("Veta: {}", result),
///     Err(e) => error!("Error calculating Veta: {:?}", e),
/// }
/// ```
/// # Notes
///
/// - It is common practice to divide the mathematical result of veta by 100 times
///   the number of days per year to reduce the value to the percentage change in
///   vega per one day.
pub fn veta(option: &Options) -> Result<Decimal, GreeksError> {
    let expiration_date: Positive = option.expiration_date.get_years()?;
    if expiration_date == Decimal::ZERO {
        // At expiration, volatility has no impact on option price
        return Ok(Decimal::ZERO);
    }
    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        expiration_date,
        option.implied_volatility,
    )?;
    let vega = vega(option)?;
    let implied_volatility: Positive = option.implied_volatility;
    let dividend_yield: Decimal = option.dividend_yield.into();
    let risk_free_rate: Decimal = option.risk_free_rate;
    let add1 =
        (risk_free_rate - dividend_yield) * d1 / (implied_volatility * expiration_date.sqrt());
    let add2 = (Decimal::ONE + d1 * d2) / (Decimal::TWO * expiration_date);

    let veta: Decimal = -vega * (dividend_yield + add1 - add2);
    // It is common practice to divide the mathematical result of veta by
    // 100 times the number of days per year to reduce the value to the
    // percentage change in vega per one day
    let veta_adj: Decimal = veta / (TRADING_DAYS * Decimal::ONE_HUNDRED);

    let quantity: Decimal = option.quantity.into();
    Ok(veta_adj * quantity)
}

/// Computes the Charm of an option.
///
/// Charm, also known as DdeltaDtime or Delta decay, measures the sensitivity of
/// the option's delta to time decay. The mathematical result of the formula for
/// charm is expressed in delta per year. It is useful to divide this by the
/// number of days per year to arrive at the delta decay per day. This usage is
/// fairly accurate when the number of days remaining until option expiration is
/// large. When an option nears expiration, charm itsel may change quickly,
/// rendering full day estimate of delta decay inaccurate.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `option_style`: The style of the option (Call or Put).
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated Charm value for the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails
///
/// # Formula
///
/// The Charm is calculated using the Black-Scholes model. The formula differs
/// for call and put options:
///
/// **Call Options:**
///
/// ```math
/// \text{Charm}_{\text{Call}} =
/// qe^{-q\tau}N(d_{\text{1}})
/// -e^{-q\tau}n(d_{\text{1}})
/// \frac{2(r-q)\tau-d_{\text{2}}\sigma\sqrt{\tau}}
/// {2\tau\sigma\sqrt{\tau}}
/// ```
///
/// **Put Options:**
///
/// ```math
/// \text{Charm}_{\text{Put}} =
/// -qe^{-q\tau}N(-d_{\text{1}})
/// -e^{-q\tau}n(d_{\text{1}})
/// \frac{2(r-q)\tau-d_{\text{2}}\sigma\sqrt{\tau}}
/// {2\tau\sigma\sqrt{\tau}}
/// ```
///
/// Where:
/// - \( S \): Underlying price
/// - \( \sigma \): Implied volatility
/// - \( \tau \): Time to expiration (in years)
/// - \( r \): Risk-free rate
/// - \( q \): Dividend yield
/// - \( K \): Strike price
/// - \( N(d1) \): Cumulative distribution function (CDF) of the standard normal
///   distribution at \( d1 \).
/// - \( n(d1) \): Probability density function (PDF) of the standard normal
///   distribution at \( d1 \).
///
/// # Calculation Steps
/// 1. Compute \( d1 \) and \( d2 \) using the `d1` and `d2` functions.
/// 2. Calculate the common term:
///    ```math
///    \text{common\_term} =
///    \frac{2(r-q)\tau-d_{\text{2}}\sigma\sqrt{\tau}}
///    {2\tau\sigma\sqrt{\tau}}
///    ```
/// 3. Apply the corresponding formula for Call or Put options, accounting for
///    the effect of dividends (\( e^{-q\tau} \)).
/// 4. Multiply the resulting Charm by the quantity of options.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::charm;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic, Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match charm(&option) {
///     Ok(result) => info!("Charm: {}", result),
///     Err(e) => error!("Error calculating Charm: {:?}", e),
/// }
/// ```
///
/// # Notes
///
/// - With zero DTE Charm can be considered as zero.
/// - Charm effects are more pronounced near expiration.
pub fn charm(option: &Options) -> Result<Decimal, GreeksError> {
    let tau = option.expiration_date.get_years()?;
    // if DTE is zero we can assume Charm is also zero
    if tau == Decimal::ZERO {
        return Ok(Decimal::ZERO);
    }

    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        tau, // expiration date
        option.implied_volatility,
    )?;
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        tau, // expiration date
        option.implied_volatility,
    )?;
    let r = option.risk_free_rate;
    let q = option.dividend_yield.to_dec();
    let sigma = option.implied_volatility;
    let exp_minus_qt = (-q * tau).exp();
    let common_term = (Decimal::TWO * (r - q) * tau - d2 * sigma * tau.sqrt())
        / (Decimal::TWO * tau * sigma * tau.sqrt());
    let charm = match option.option_style {
        OptionStyle::Call => {
            (q * exp_minus_qt * big_n(d1)?) - (exp_minus_qt * n(d1)? * common_term)
        }
        OptionStyle::Put => {
            (-q * exp_minus_qt * big_n(-d1)?) - (exp_minus_qt * n(d1)? * common_term)
        }
    };
    // Adjust for quantity and convert to daily value
    Ok((charm * option.quantity) / Decimal::from(365))
}

/// Computes the Color of an option.
///
/// Color, also known as DgammaDtime or Gamma decay, measures the sensitivity of
/// the option's gamma to time decay. The mathematical result of the formula for
/// color is expressed in gamma per year. It is useful to divide this by the
/// number of days per year to arrive at the gamma decay per day. This usage is
/// fairly accurate when the number of days remaining until option expiration is
/// large. When an option nears expiration, color itsel may change quickly,
/// rendering full days estimate of gamma decay inaccurate.
///
/// # Parameters
///
/// - `option: &Options`
///   A reference to an `Options` struct containing the following relevant parameters:
///   - `underlying_price`: The current price of the underlying asset.
///   - `strike_price`: The strike price of the option.
///   - `risk_free_rate`: The risk-free interest rate.
///   - `expiration_date`: The time to expiration in years (provides `get_years` method).
///   - `implied_volatility`: The implied volatility of the option.
///   - `dividend_yield`: The dividend yield of the underlying asset.
///   - `option_style`: The style of the option (Call or Put).
///   - `quantity`: The quantity of the options.
///
/// # Returns
///
/// - `Ok(Decimal)`: The calculated Color value for the option.
/// - `Err(GreeksError)`: Returns an error if any intermediate calculation fails
///
/// # Formula
///
/// The Color is calculated using the Black-Scholes model. The formula is the
/// same for call and put options:
///
/// ```math
/// \text{Color} = -e^{-q\tau}
/// \frac{n(d_{\text{1}})}
/// {2S\tau\sigma\sqrt{\tau}}
/// \left[
/// 2q\tau+1+
/// \frac{2(r-q)\tau-d_{\text{2}}\sigma\sqrt{\tau}}
/// {\sigma\sqrt{\tau}}d_{\text{1}}
/// \right]
/// ```
///
/// Where:
/// - \( S \): Underlying price
/// - \( \sigma \): Implied volatility
/// - \( \tau \): Time to expiration (in years)
/// - \( r \): Risk-free rate
/// - \( q \): Dividend yield
/// - \( n(d1) \): Probability density function (PDF) of the standard normal
///   distribution at \( d1 \).
///
/// # Calculation Steps
/// 1. Compute \( d1 \) and \( d2 \) using the `d1` and `d2` functions.
/// 2. Compute \( n(d1) \) using the `n` function.
/// 3. Apply the corresponding Color formula, accounting for
///    the effect of dividends (\( e^{-q\tau} \)).
/// 4. Multiply the resulting Color value by the quantity of options.
///
/// # Example
///
/// ```rust
/// use rust_decimal_macros::dec;
/// use tracing::{error, info};
/// use optionstratlib::greeks::color;
/// use optionstratlib::{ExpirationDate, Options};
/// use optionstratlib::model::types::{ OptionStyle, OptionType, Side};
/// use positive::{pos_or_panic,Positive};
/// let option = Options {
///     option_type: OptionType::European,
///     side: Side::Long,
///     underlying_price: Positive::HUNDRED,
///     strike_price: pos_or_panic!(95.0),
///     risk_free_rate: dec!(0.05),
///     expiration_date: ExpirationDate::Days(pos_or_panic!(30.0)),
///     implied_volatility: pos_or_panic!(0.2),
///     dividend_yield: pos_or_panic!(0.01),
///     quantity: Positive::ONE,
///     option_style: OptionStyle::Call,
///     underlying_symbol: "".to_string(),
///     exotic_params: None,
/// };
///
/// match color(&option) {
///     Ok(result) => info!("Color: {}", result),
///     Err(e) => error!("Error calculating Color: {:?}", e),
/// }
/// ```
///
/// # Notes
/// - Color is generally negative for long options and positive for short options.
/// - Color will be more pronounced as expiration date approaches.
/// - When volatility increases Color sensitivity decrease.
/// - Deep ITM and OTM options have negligible Color.
pub fn color(option: &Options) -> Result<Decimal, GreeksError> {
    let tau = option.expiration_date.get_years()?;
    // if DTE is zero we can assume Color is also zero
    if tau == Decimal::ZERO {
        return Ok(Decimal::ZERO);
    }

    let d1 = d1(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        tau, // expiration date
        option.implied_volatility,
    )?;
    let d2 = d2(
        option.underlying_price,
        option.strike_price,
        option.risk_free_rate,
        tau, // expiration date
        option.implied_volatility,
    )?;
    let r = option.risk_free_rate;
    let s = option.underlying_price;
    let q = option.dividend_yield.to_dec();
    let sigma = option.implied_volatility;
    let exp_minus_qt = (-q * tau).exp();
    let factor1 = n(d1)? / (Decimal::TWO * s * tau * sigma * tau.sqrt());
    let numerator = (Decimal::TWO * (r - q) * tau) - (d2 * sigma * tau.sqrt());
    let denominator = sigma * tau.sqrt();
    let factor2 = (Decimal::TWO * q * tau) + Decimal::ONE + ((numerator / denominator) * d1);
    let color = (-exp_minus_qt * factor1 * factor2 * option.quantity) / Decimal::from(365);
    Ok(color)
}

#[cfg(test)]
pub mod tests_delta_equations {
    use super::*;
    use crate::constants::ZERO;
    use crate::model::types::{OptionStyle, Side};
    use crate::model::utils::create_sample_option;
    use crate::strategies::DELTA_THRESHOLD;
    use positive::constants::DAYS_IN_A_YEAR;

    use crate::{ExpirationDate, assert_decimal_eq};
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::{Positive, pos_or_panic};
    use rust_decimal_macros::dec;
    use tracing::info;

    #[test]
    fn test_delta_no_volatility_itm() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value.to_f64().unwrap(), 1.0, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_otm() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(110.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, ZERO, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_itm_put() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(150.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, -1.0, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_otm_put() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(160.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, ZERO, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_itm_short() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Short,
            pos_or_panic!(150.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, -1.0, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_otm_short() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Short,
            pos_or_panic!(110.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, ZERO, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_itm_put_short() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Short,
            pos_or_panic!(150.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, 1.0, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_no_volatility_otm_put_short() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Short,
            pos_or_panic!(160.0),
            Positive::ONE,
            pos_or_panic!(150.0),
            Positive::ZERO,
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility: {}", delta_value);
        assert_relative_eq!(delta_value, ZERO, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_deep_in_the_money_call() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0),
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.20),
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Deep ITM Call Delta: {}", delta_value);
        assert_relative_eq!(delta_value, 0.9991784198733309, epsilon = 1e-8);
    }

    #[test]
    fn test_delta_deep_out_of_the_money_call() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0),
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.20),
        );
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Deep OTM Call Delta: {}", delta_value);
        assert_relative_eq!(delta_value, 2.0418256951423236e-33, epsilon = 1e-4);
    }

    #[test]
    fn test_delta_at_the_money_put() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.20),
        );
        let delta_value = delta(&option).unwrap();
        info!("ATM Put Delta: {}", delta_value);
        assert_decimal_eq!(delta_value, dec!(-0.459658497), DELTA_THRESHOLD);
    }

    #[test]
    fn test_delta_short_term_high_volatility() {
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.50),
        );
        option.expiration_date = ExpirationDate::Days(pos_or_panic!(7.0));
        let delta_value = delta(&option).unwrap().to_f64().unwrap();
        info!("Short-term High Vol Call Delta: {}", delta_value);
        assert_relative_eq!(delta_value, 0.519229469584234, epsilon = 1e-4);
    }

    #[test]
    fn test_delta_long_term_low_volatility() {
        let mut option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.10),
        );
        option.expiration_date = ExpirationDate::Days(DAYS_IN_A_YEAR);
        let delta_value = delta(&option).unwrap();
        info!("Long-term Low Vol Put Delta: {}", delta_value);
        assert_decimal_eq!(delta_value, dec!(-0.2882625996), DELTA_THRESHOLD);
    }

    #[test]
    fn test_delta_long_almost_zero_time_to_maturity() {
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Short,
            pos_or_panic!(21637.0),
            Positive::ONE,
            pos_or_panic!(21825.0),
            pos_or_panic!(0.219),
        );
        option.expiration_date = ExpirationDate::Days(Positive::ONE);
        let delta_value = delta(&option).unwrap();
        info!("Long-term Low Vol Put Delta: {}", delta_value);
        assert_decimal_eq!(delta_value, dec!(-0.230544), DELTA_THRESHOLD);
    }
}

#[cfg(test)]
pub mod tests_gamma_equations {
    use super::*;
    use crate::model::types::{OptionStyle, Side};
    use crate::model::utils::create_sample_option;
    use positive::constants::DAYS_IN_A_YEAR;

    use crate::ExpirationDate;
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::pos_or_panic;
    use tracing::info;

    #[test]
    fn test_gamma_deep_in_the_money_call() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0),
            Positive::ONE,
            pos_or_panic!(120.0),
            pos_or_panic!(0.2),
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Deep ITM Call Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.000016049457791525, epsilon = 1e-8);
    }

    #[test]
    fn test_gamma_deep_out_of_the_money_call() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0),
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.20),
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Deep OTM Call Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.0, epsilon = 1e-34);
    }

    #[test]
    fn test_gamma_at_the_money_put() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.20),
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("ATM Put Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.06917076441486919, epsilon = 1e-8);
    }

    #[test]
    fn test_gamma_short_term_high_volatility() {
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.50),
        );
        option.expiration_date = ExpirationDate::Days(pos_or_panic!(7.0));
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Short-term High Vol Call Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.05753657912620555, epsilon = 1e-8);
    }

    #[test]
    fn test_gamma_long_term_low_volatility() {
        let mut option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.10),
        );
        option.expiration_date = ExpirationDate::Days(DAYS_IN_A_YEAR);
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Long-term Low Vol Put Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.033953150664723986, epsilon = 1e-8);
    }

    #[test]
    fn test_gamma_zero_volatility() {
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            Positive::ZERO,
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Zero Volatility Call Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_gamma_extreme_high_volatility() {
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Short,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(5.0),
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Extreme High Volatility Put Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.002146478293943308, epsilon = 1e-8);
    }
}

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

    use crate::{ExpirationDate, OptionType};
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::pos_or_panic;
    use tracing::info;

    #[test]
    fn test_50_vol_10() {
        let option = Options::new(
            OptionType::European,
            Side::Long,
            "XYZ".parse().unwrap(),
            pos_or_panic!(50.0),
            ExpirationDate::Days(pos_or_panic!(365.0)),
            pos_or_panic!(0.10),
            Positive::ONE,
            pos_or_panic!(50.0),
            Decimal::ZERO,
            OptionStyle::Call,
            Positive::ZERO,
            None,
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.0796887828189609, epsilon = 1e-8);
    }

    #[test]
    fn test_50_vol_5() {
        let option = Options::new(
            OptionType::European,
            Side::Long,
            "XYZ".parse().unwrap(),
            pos_or_panic!(50.0),
            ExpirationDate::Days(pos_or_panic!(365.0)),
            pos_or_panic!(0.05),
            Positive::ONE,
            pos_or_panic!(50.0),
            Decimal::ZERO,
            OptionStyle::Call,
            Positive::ZERO,
            None,
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.15952705216736393, epsilon = 1e-8);
    }

    #[test]
    fn test_50_vol_20() {
        let option = Options::new(
            OptionType::European,
            Side::Long,
            "XYZ".parse().unwrap(),
            pos_or_panic!(50.0),
            ExpirationDate::Days(pos_or_panic!(365.0)),
            pos_or_panic!(0.2),
            Positive::ONE,
            pos_or_panic!(50.0),
            Decimal::ZERO,
            OptionStyle::Call,
            Positive::ZERO,
            None,
        );
        let gamma_value = gamma(&option).unwrap().to_f64().unwrap();
        info!("Gamma: {}", gamma_value);
        assert_relative_eq!(gamma_value, 0.03969525474873078, epsilon = 1e-8);
    }
}

#[cfg(test)]
pub mod tests_vega_equation {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::{OptionType, Side};
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    fn create_test_option(
        underlying_price: Positive,
        strike_price: Positive,
        implied_volatility: Positive,
        dividend_yield: Positive,
        expiration_in_days: Positive,
    ) -> Options {
        Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            strike_price,
            ExpirationDate::Days(expiration_in_days),
            implied_volatility,
            Positive::ONE, // Quantity
            underlying_price,
            dec!(0.05), // Risk-free rate
            OptionStyle::Call,
            dividend_yield,
            None, // No exotic params for this test
        )
    }

    #[test]
    fn test_vega_atm() {
        let option = create_test_option(
            Positive::HUNDRED,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO,
            DAYS_IN_A_YEAR,
        );
        let vega = vega(&option).unwrap().to_f64().unwrap();
        let expected_vega = 0.3752403469;
        assert!(
            (vega - expected_vega).abs() < 1e-5,
            "Vega ATM test failed: expected {expected_vega}, got {vega}"
        );
    }

    #[test]
    fn test_vega_otm() {
        let option = create_test_option(
            pos_or_panic!(90.0),
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO,
            DAYS_IN_A_YEAR,
        );
        let vega = vega(&option).unwrap().to_f64().unwrap();
        let expected_vega = 0.35347991;
        assert!(
            (vega - expected_vega).abs() < 1e-5,
            "Vega OTM test failed: expected {expected_vega}, got {vega}"
        );
    }

    #[test]
    fn test_vega_short_expiration() {
        let option = create_test_option(
            Positive::HUNDRED,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO,
            Positive::ONE,
        );
        let vega = vega(&option).unwrap().to_f64().unwrap();
        let expected_vega = 0.020878089;
        assert!(
            (vega - expected_vega).abs() < 1e-5,
            "Vega short expiration test failed: expected {expected_vega}, got {vega}"
        );
    }

    #[test]
    fn test_vega_with_dividends() {
        let option = create_test_option(
            Positive::HUNDRED,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(0.03),
            Positive::ONE,
        );
        let vega = vega(&option).unwrap().to_f64().unwrap();
        let expected_vega = 0.0208763735;
        assert!(
            (vega - expected_vega).abs() < 1e-5,
            "Vega with dividends test failed: expected {expected_vega}, got {vega}"
        );
    }

    #[test]
    fn test_vega_itm() {
        let option = create_test_option(
            pos_or_panic!(110.0),
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO,
            Positive::ONE,
        );
        let vega = vega(&option).unwrap().to_f64().unwrap();
        let expected_vega = 0.0;
        assert!(
            (vega - expected_vega).abs() < 1e-5,
            "Vega ITM test failed: expected {expected_vega}, got {vega}"
        );
    }
}

#[cfg(test)]
pub mod tests_rho_equations {
    use super::*;
    use crate::model::types::{OptionStyle, OptionType, Side};
    use crate::{ExpirationDate, assert_decimal_eq};
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    fn create_test_option(style: OptionStyle) -> 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),
            implied_volatility: pos_or_panic!(0.2),
            quantity: Positive::ONE,
            underlying_price: Positive::HUNDRED,
            risk_free_rate: dec!(0.05),
            option_style: style,
            dividend_yield: Positive::ZERO,
            exotic_params: None,
        }
    }

    #[test]
    fn test_rho_call_option() {
        let option = create_test_option(OptionStyle::Call);
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, 0.532324815464, epsilon = 1e-8);
    }

    #[test]
    fn test_rho_put_option() {
        let option = create_test_option(OptionStyle::Put);
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, -0.41890460905, epsilon = 1e-8);
    }

    #[test]
    fn test_rho_zero_time_to_expiry() {
        let mut option = create_test_option(OptionStyle::Call);
        option.expiration_date = ExpirationDate::Days(Positive::ZERO);
        let result = rho(&option).is_ok();
        assert!(result);
        assert_decimal_eq!(rho(&option).unwrap(), Decimal::ZERO, dec!(1e-8));
    }

    #[test]
    fn test_rho_zero_risk_free_rate() {
        let mut option = create_test_option(OptionStyle::Call);
        option.risk_free_rate = dec!(0.0);
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, 0.460172162, epsilon = 1e-8);
    }

    #[test]
    fn test_rho_deep_out_of_money_call() {
        let mut option = create_test_option(OptionStyle::Call);
        option.strike_price = pos_or_panic!(1000.0);
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_rho_deep_out_of_money_put() {
        let mut option = create_test_option(OptionStyle::Put);
        option.strike_price = Positive::ONE;
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, 0.0, epsilon = 1e-8);
    }

    #[test]
    fn test_rho_high_volatility() {
        let mut option = create_test_option(OptionStyle::Call);
        option.implied_volatility = Positive::ONE;
        let result = rho(&option).unwrap().to_f64().unwrap();
        assert_relative_eq!(result, 0.3104386883, epsilon = 0.0001);
    }
}

#[cfg(test)]
pub mod tests_theta_long_equations {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::Side;
    use crate::model::utils::create_sample_option;
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;

    #[test]
    fn test_theta_call_option() {
        // Create a sample call option
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(155.0), // strike price
            pos_or_panic!(0.20),  // implied volatility
        );

        // Expected theta value for a call option (precomputed or from known source)
        let expected_theta = -0.0561725050;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_put_option() {
        // Create a sample put option
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(145.0), // strike price
            pos_or_panic!(0.25),  // implied volatility
        );

        // Expected theta value for a put option (precomputed or from known source)
        let expected_theta = -0.055928204732;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_call_option_near_expiry() {
        // Create a sample call option near expiry
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(150.0), // strike price
            pos_or_panic!(0.15),  // implied volatility
        );
        option.expiration_date = ExpirationDate::Days(Positive::ONE); // Option close to expiry

        // Expected theta value for a near-expiry call option (precomputed)
        let expected_theta = -0.24315788969;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_put_option_far_from_expiry() {
        // Create a sample put option far from expiry
        let mut option = create_sample_option(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(140.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(130.0), // strike price
            pos_or_panic!(0.30),  // implied volatility
        );
        option.expiration_date = ExpirationDate::Days(DAYS_IN_A_YEAR); // Option far from expiry

        // Expected theta value for a far-expiry put option (precomputed)
        let expected_theta = -0.0139607780805;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }
}

#[cfg(test)]
pub mod tests_theta_short_equations {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::Side;
    use crate::model::utils::create_sample_option;
    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;

    #[test]
    fn test_theta_short_call_option() {
        // Create a sample short call option
        let option = create_sample_option(
            OptionStyle::Call,
            Side::Short,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(155.0), // strike price
            pos_or_panic!(0.20),  // implied volatility
        );

        // Expected theta value for a short call option (precomputed or from known source)
        let expected_theta = -0.05617250509;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_short_put_option() {
        // Create a sample short put option
        let option = create_sample_option(
            OptionStyle::Put,
            Side::Short,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(145.0), // strike price
            pos_or_panic!(0.25),  // implied volatility
        );

        // Expected theta value for a short put option (precomputed or from known source)
        let expected_theta = -0.05592820473;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_short_call_option_near_expiry() {
        // Create a sample short call option near expiry
        let mut option = create_sample_option(
            OptionStyle::Call,
            Side::Short,
            pos_or_panic!(150.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(150.0), // strike price
            pos_or_panic!(0.15),  // implied volatility
        );
        option.expiration_date = ExpirationDate::Days(Positive::ONE); // Option close to expiry

        // Expected theta value for a short near-expiry call option (precomputed)
        let expected_theta = -0.2431578896;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }

    #[test]
    fn test_theta_short_put_option_far_from_expiry() {
        // Create a sample short put option far from expiry
        let mut option = create_sample_option(
            OptionStyle::Put,
            Side::Short,
            pos_or_panic!(140.0), // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(130.0), // strike price
            pos_or_panic!(0.30),  // implied volatility
        );
        option.expiration_date = ExpirationDate::Days(DAYS_IN_A_YEAR); // Option far from expiry

        // Expected theta value for a far-expiry short put option (precomputed)
        let expected_theta = -0.01396077;

        // Compute the theta value using the function
        let calculated_theta = theta(&option).unwrap().to_f64().unwrap();

        // Assert the calculated theta is close to the expected value
        assert_relative_eq!(calculated_theta, expected_theta, epsilon = 1e-8);
    }
}

#[cfg(test)]
mod tests_greeks_trait {
    use super::*;
    use crate::model::types::{OptionStyle, OptionType, Side};
    use crate::{ExpirationDate, assert_decimal_eq};
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    // A simple struct for testing the Greeks trait
    struct TestOptionCollection {
        options: Vec<Options>,
    }

    impl Greeks for TestOptionCollection {
        fn get_options(&self) -> Result<Vec<&Options>, GreeksError> {
            Ok(self.options.iter().collect())
        }
    }

    // Helper function to create a test option
    fn create_test_option(side: Side, style: OptionStyle, quantity: Positive) -> Options {
        Options::new(
            OptionType::European,
            side,
            "TEST".to_string(),
            Positive::HUNDRED, // strike_price
            ExpirationDate::Days(pos_or_panic!(30.0)),
            pos_or_panic!(0.2), // implied_volatility
            quantity,
            Positive::HUNDRED, // underlying_price
            dec!(0.05),        // risk_free_rate
            style,
            pos_or_panic!(0.01), // dividend_yield
            None,                // exotic_params
        )
    }

    #[test]
    fn test_greeks_single_option() {
        let option = create_test_option(Side::Long, OptionStyle::Call, Positive::ONE);
        let collection = TestOptionCollection {
            options: vec![option],
        };

        let greeks = collection.greeks().unwrap();

        // Test each greek value
        assert_decimal_eq!(greeks.delta, dec!(0.539519922), dec!(0.000001));
        assert_decimal_eq!(greeks.gamma, dec!(0.069170764), dec!(0.000001));
        assert_decimal_eq!(greeks.theta, dec!(-0.04351001), dec!(0.000001));
        assert_decimal_eq!(greeks.vega, dec!(0.1137053), dec!(0.000001));
        assert_decimal_eq!(greeks.rho, dec!(0.04233121458), dec!(0.000001));
        assert_decimal_eq!(greeks.rho_d, dec!(-0.04434410), dec!(0.000001));
        assert_decimal_eq!(greeks.vanna, dec!(-0.08527902), dec!(0.000001));
        assert_decimal_eq!(greeks.vomma, dec!(0.00245323), dec!(0.000001));
        assert_decimal_eq!(greeks.veta, dec!(0.00002720), dec!(0.000001));
    }

    #[test]
    fn test_greeks_multiple_options() {
        let option1 = create_test_option(Side::Long, OptionStyle::Call, Positive::ONE);
        let option2 = create_test_option(Side::Short, OptionStyle::Put, Positive::ONE);
        let collection = TestOptionCollection {
            options: vec![option1, option2],
        };

        let greeks = collection.greeks().unwrap();

        // Test aggregated greek values
        assert!(
            greeks.delta.abs() > dec!(0.0),
            "Delta should be non-zero for multiple options"
        );
        assert!(
            greeks.gamma.abs() > dec!(0.0),
            "Gamma should be non-zero for multiple options"
        );
        assert!(
            greeks.theta.abs() > dec!(0.0),
            "Theta should be non-zero for multiple options"
        );
        assert!(
            greeks.vega.abs() > dec!(0.0),
            "Vega should be non-zero for multiple options"
        );
        assert!(
            greeks.rho.abs() > dec!(0.0),
            "Rho should be non-zero for multiple options"
        );
        assert!(
            greeks.rho_d.abs() > dec!(0.0),
            "Rho_d should be non-zero for multiple options"
        );
        assert!(
            greeks.vanna.abs() > dec!(0.0),
            "Vanna should be non-zero for multiple options"
        );
        assert!(
            greeks.vomma.abs() > dec!(0.0),
            "Vomma should be non-zero for multiple options"
        );
        assert!(
            greeks.veta.abs() > dec!(0.0),
            "Veta should be non-zero for multiple options"
        );
    }

    #[test]
    fn test_greeks_simple_validation() {
        let option = Options::new(
            OptionType::European,
            Side::Long,
            "AAPL".to_string(),
            pos_or_panic!(155.0),
            ExpirationDate::Days(pos_or_panic!(30.0)),
            pos_or_panic!(0.20),
            Positive::ONE,
            pos_or_panic!(150.0),
            dec!(0.05),
            OptionStyle::Call,
            pos_or_panic!(0.00),
            None,
        );

        let greeks = option.greeks().unwrap();

        assert_decimal_eq!(greeks.delta, dec!(0.3186329), dec!(0.000001));
        assert_decimal_eq!(greeks.gamma, dec!(0.0415044), dec!(0.000001));
        assert_decimal_eq!(greeks.theta, dec!(-0.0574808), dec!(0.000001));
        assert_decimal_eq!(greeks.vega, dec!(0.15350973), dec!(0.000001));
        assert_decimal_eq!(greeks.rho, dec!(0.03786580), dec!(0.000001));
        assert_decimal_eq!(greeks.rho_d, dec!(-0.03928351), dec!(0.000001));
        assert_decimal_eq!(greeks.vanna, dec!(0.94393865), dec!(0.000001));
        assert_decimal_eq!(greeks.vomma, dec!(0.19140525), dec!(0.000001));
        assert_decimal_eq!(greeks.veta, dec!(0.00004880), dec!(0.000001));
    }

    #[test]
    fn test_greeks_zero_quantity() {
        let option = create_test_option(Side::Long, OptionStyle::Call, Positive::ZERO);
        let collection = TestOptionCollection {
            options: vec![option],
        };

        let greeks = collection.greeks().unwrap();

        // All greeks should be zero for zero quantity
        assert_eq!(greeks.delta, dec!(0.0));
        assert_eq!(greeks.gamma, dec!(0.0));
        assert_eq!(greeks.theta, dec!(0.0));
        assert_eq!(greeks.vega, dec!(0.0));
        assert_eq!(greeks.rho, dec!(0.0));
        assert_eq!(greeks.rho_d, dec!(0.0));
        assert_eq!(greeks.vanna, dec!(0.0));
        assert_eq!(greeks.vomma, dec!(0.0));
        assert_eq!(greeks.veta, dec!(0.0));
    }

    #[test]
    fn test_greeks_opposing_positions() {
        let option1 = Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            pos_or_panic!(50.0), // strike_price
            ExpirationDate::Days(pos_or_panic!(365.0)),
            pos_or_panic!(0.2), // implied_volatility
            Positive::ONE,
            pos_or_panic!(50.0), // underlying_price
            dec!(0.05),          // risk_free_rate
            OptionStyle::Call,
            pos_or_panic!(0.01), // dividend_yield
            None,                // exotic_params
        );
        let option2 = Options::new(
            OptionType::European,
            Side::Short,
            "TEST".to_string(),
            pos_or_panic!(50.0), // strike_price
            ExpirationDate::Days(pos_or_panic!(365.0)),
            pos_or_panic!(0.2), // implied_volatility
            Positive::ONE,
            pos_or_panic!(50.0), // underlying_price
            dec!(0.05),          // risk_free_rate
            OptionStyle::Call,
            pos_or_panic!(0.01), // dividend_yield
            None,                // exotic_params
        );
        let collection = TestOptionCollection {
            options: vec![option1, option2],
        };

        let greeks = collection.greeks().unwrap();

        // Opposing positions should mostly cancel out
        assert_decimal_eq!(greeks.delta, Decimal::ZERO, dec!(0.000001));
        assert_decimal_eq!(greeks.gamma, dec!(0.0743013), dec!(0.000001));
        assert_decimal_eq!(greeks.vega, dec!(0.37150664), dec!(0.000001));
        assert_decimal_eq!(greeks.rho, dec!(0.532324815), dec!(0.000001));
        assert_decimal_eq!(greeks.vanna, dec!(-0.55725996), dec!(0.000001));
        assert_decimal_eq!(greeks.vomma, dec!(0.09752049), dec!(0.000001));
        assert_decimal_eq!(greeks.veta, dec!(0.00000657), dec!(0.000001));
    }

    #[test]
    fn test_individual_greek_methods() {
        let option1 = create_test_option(Side::Long, OptionStyle::Call, Positive::ONE);
        let option2 = create_test_option(Side::Short, OptionStyle::Put, Positive::ONE);
        let collection = TestOptionCollection {
            options: vec![option1, option2],
        };

        // Test each individual greek method
        let delta = collection.delta().unwrap();
        let gamma = collection.gamma().unwrap();
        let theta = collection.theta().unwrap();
        let vega = collection.vega().unwrap();
        let rho = collection.rho().unwrap();
        let rho_d = collection.rho_d().unwrap();
        let vanna = collection.vanna().unwrap();
        let vomma = collection.vomma().unwrap();
        let veta = collection.veta().unwrap();

        // Verify each value is non-zero (actual values depend on input parameters)
        assert!(delta.abs() > dec!(0.0), "Delta calculation failed");
        assert!(gamma.abs() > dec!(0.0), "Gamma calculation failed");
        assert!(theta.abs() > dec!(0.0), "Theta calculation failed");
        assert!(vega.abs() > dec!(0.0), "Vega calculation failed");
        assert!(rho.abs() > dec!(0.0), "Rho calculation failed");
        assert!(rho_d.abs() > dec!(0.0), "Rho_d calculation failed");
        assert!(vanna.abs() > dec!(0.0), "Vanna calculation failed");
        assert!(vomma.abs() > dec!(0.0), "Vomma calculation failed");
        assert!(veta.abs() > dec!(0.0), "Veta calculation failed");
    }

    #[test]
    fn test_empty_option_collection() {
        let collection = TestOptionCollection { options: vec![] };

        // All greeks should be zero for empty collection
        let greeks = collection.greeks().unwrap();
        assert_eq!(greeks.delta, dec!(0.0));
        assert_eq!(greeks.gamma, dec!(0.0));
        assert_eq!(greeks.theta, dec!(0.0));
        assert_eq!(greeks.vega, dec!(0.0));
        assert_eq!(greeks.rho, dec!(0.0));
        assert_eq!(greeks.rho_d, dec!(0.0));
        assert_eq!(greeks.vanna, dec!(0.0));
        assert_eq!(greeks.vomma, dec!(0.0));
        assert_eq!(greeks.veta, dec!(0.0));
    }

    #[test]
    fn test_greeks_with_different_expirations() {
        let mut option1 = create_test_option(Side::Long, OptionStyle::Call, Positive::ONE);
        let mut option2 = create_test_option(Side::Long, OptionStyle::Call, Positive::ONE);

        // Set different expiration dates
        option1.expiration_date = ExpirationDate::Days(pos_or_panic!(30.0));
        option2.expiration_date = ExpirationDate::Days(pos_or_panic!(60.0));

        let collection = TestOptionCollection {
            options: vec![option1, option2],
        };

        let greeks = collection.greeks().unwrap();

        // Verify values are calculated correctly for different expirations
        assert!(greeks.delta.abs() > dec!(0.0));
        assert!(greeks.gamma.abs() > dec!(0.0));
        assert!(greeks.theta.abs() > dec!(0.0));
        assert!(greeks.vega.abs() > dec!(0.0));
        assert!(greeks.rho.abs() > dec!(0.0));
        assert!(greeks.rho_d.abs() > dec!(0.0));
        assert!(greeks.vanna.abs() > dec!(0.0));
        assert!(greeks.vomma.abs() > dec!(0.0));
        assert!(greeks.veta.abs() > dec!(0.0));
    }
}

#[cfg(test)]
pub mod tests_vanna_equation {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::{OptionType, Side};
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    fn create_test_option(
        underlying_price: Positive,
        strike_price: Positive,
        implied_volatility: Positive,
        dividend_yield: Positive,
        expiration_in_days: Positive,
    ) -> Options {
        Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            strike_price,
            ExpirationDate::Days(expiration_in_days),
            implied_volatility,
            Positive::ONE, // Quantity
            underlying_price,
            dec!(0.05), // Risk-free rate
            OptionStyle::Call,
            dividend_yield,
            None, // No exotic params for this test
        )
    }

    #[test]
    fn test_vanna_atm() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            DAYS_IN_A_YEAR,     // expiration_in_days
        );
        let vanna = vanna(&option).unwrap().to_f64().unwrap();
        let expected_vanna = -0.28143026;
        assert!(
            (vanna - expected_vanna).abs() < 1e-5,
            "Vega ATM test failed: expected {expected_vanna}, got {vanna}"
        );
    }

    #[test]
    fn test_vanna_otm() {
        let option = create_test_option(
            pos_or_panic!(90.0), // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            Positive::ZERO,      // dividend_yield
            DAYS_IN_A_YEAR,      // expiration_in_days
        );
        let vanna = vanna(&option).unwrap().to_f64().unwrap();
        let expected_vanna = 0.73995634;
        assert!(
            (vanna - expected_vanna).abs() < 1e-5,
            "Vanna OTM test failed: expected {expected_vanna}, got {vanna}"
        );
    }

    #[test]
    fn test_vanna_short_expiration() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            Positive::ONE,      // expiration_in_days
        );
        let vanna = vanna(&option).unwrap().to_f64().unwrap();
        let expected_vanna = -0.01565856;
        assert!(
            (vanna - expected_vanna).abs() < 1e-5,
            "Vanna short expiration test failed: expected {expected_vanna}, got {vanna}"
        );
    }

    #[test]
    fn test_vanna_with_dividends() {
        let option = create_test_option(
            Positive::HUNDRED,   // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            pos_or_panic!(0.03), // dividend_yield
            Positive::ONE,       // expiration_in_days
        );
        let vanna = vanna(&option).unwrap().to_f64().unwrap();
        let expected_vanna = -0.01565728;
        assert!(
            (vanna - expected_vanna).abs() < 1e-5,
            "Vanna with dividends test failed: expected {expected_vanna}, got {vanna}"
        );
    }

    #[test]
    fn test_vanna_itm() {
        let option = create_test_option(
            pos_or_panic!(110.0), // underlying_price
            Positive::HUNDRED,    // strike_price
            pos_or_panic!(0.2),   // implied_volatility
            Positive::ZERO,       // dividend_yield
            Positive::ONE,        // expiration_in_days
        );
        let vanna = vanna(&option).unwrap().to_f64().unwrap();
        let expected_vanna = 0.0;
        assert!(
            (vanna - expected_vanna).abs() < 1e-5,
            "Vanna ITM test failed: expected {expected_vanna}, got {vanna}"
        );
    }
}

#[cfg(test)]
pub mod tests_vomma_equation {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::{OptionType, Side};
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    fn create_test_option(
        underlying_price: Positive,
        strike_price: Positive,
        implied_volatility: Positive,
        dividend_yield: Positive,
        expiration_in_days: Positive,
    ) -> Options {
        Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            strike_price,
            ExpirationDate::Days(expiration_in_days),
            implied_volatility,
            Positive::ONE, // Quantity
            underlying_price,
            dec!(0.05), // Risk-free rate
            OptionStyle::Call,
            dividend_yield,
            None, // No exotic params for this test
        )
    }

    #[test]
    fn test_vomma_atm() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            DAYS_IN_A_YEAR,     // expiration_in_days
        );
        let vomma = vomma(&option).unwrap().to_f64().unwrap();
        let expected_vomma = 0.09850059;
        assert!(
            (vomma - expected_vomma).abs() < 1e-5,
            "Vomma ATM test failed: expected {expected_vomma}, got {vomma}"
        );
    }

    #[test]
    fn test_vomma_otm() {
        let option = create_test_option(
            pos_or_panic!(90.0), // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            Positive::ZERO,      // dividend_yield
            DAYS_IN_A_YEAR,      // expiration_in_days
        );
        let vomma = vomma(&option).unwrap().to_f64().unwrap();
        let expected_vomma = 0.11774357;
        assert!(
            (vomma - expected_vomma).abs() < 1e-5,
            "Vomma OTM test failed: expected {expected_vomma}, got {vomma}"
        );
    }

    #[test]
    fn test_vomma_short_expiration() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            Positive::ONE,      // expiration_in_days
        );
        let vomma = vomma(&option).unwrap().to_f64().unwrap();
        let expected_vomma = 0.0000150150;
        assert!(
            (vomma - expected_vomma).abs() < 1e-5,
            "Vomma short expiration test failed: expected {expected_vomma}, got {vomma}"
        );
    }

    #[test]
    fn test_vomma_with_dividends() {
        let option = create_test_option(
            Positive::HUNDRED,   // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            pos_or_panic!(0.03), // dividend_yield
            Positive::ONE,       // expiration_in_days
        );
        let vomma = vomma(&option).unwrap().to_f64().unwrap();
        let expected_vomma = 0.0000150138;
        assert!(
            (vomma - expected_vomma).abs() < 1e-5,
            "Vomma with dividends test failed: expected {expected_vomma}, got {vomma}"
        );
    }

    #[test]
    fn test_vomma_itm() {
        let option = create_test_option(
            pos_or_panic!(110.0), // underlying_price
            Positive::HUNDRED,    // strike_price
            pos_or_panic!(0.2),   // implied_volatility
            Positive::ZERO,       // dividend_yield
            Positive::ONE,        // expiration_in_days
        );
        let vomma = vomma(&option).unwrap().to_f64().unwrap();
        let expected_vomma = 0.0;
        assert!(
            (vomma - expected_vomma).abs() < 1e-5,
            "Vomma ITM test failed: expected {expected_vomma}, got {vomma}"
        );
    }
}

#[cfg(test)]
pub mod tests_veta_equation {
    use super::*;
    use crate::ExpirationDate;
    use crate::model::types::{OptionType, Side};
    use num_traits::ToPrimitive;
    use positive::constants::DAYS_IN_A_YEAR;
    use positive::pos_or_panic;
    use rust_decimal_macros::dec;

    fn create_test_option(
        underlying_price: Positive,
        strike_price: Positive,
        implied_volatility: Positive,
        dividend_yield: Positive,
        expiration_in_days: Positive,
    ) -> Options {
        Options::new(
            OptionType::European,
            Side::Long,
            "TEST".to_string(),
            strike_price,
            ExpirationDate::Days(expiration_in_days),
            implied_volatility,
            Positive::ONE, // Quantity
            underlying_price,
            dec!(0.05), // Risk-free rate
            OptionStyle::Call,
            dividend_yield,
            None, // No exotic params for this test
        )
    }

    #[test]
    fn test_veta_atm() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            DAYS_IN_A_YEAR,     // expiration_in_days
        );
        let veta = veta(&option).unwrap().to_f64().unwrap();
        let expected_veta = 0.0000065332;
        assert!(
            (veta - expected_veta).abs() < 1e-5,
            "Veta ATM test failed: expected {expected_veta}, got {veta}"
        );
    }

    #[test]
    fn test_veta_otm() {
        let option = create_test_option(
            pos_or_panic!(90.0), // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            Positive::ZERO,      // dividend_yield
            DAYS_IN_A_YEAR,      // expiration_in_days
        );
        let veta = veta(&option).unwrap().to_f64().unwrap();
        let expected_veta = 0.0000081007;
        assert!(
            (veta - expected_veta).abs() < 1e-5,
            "Veta OTM test failed: expected {expected_veta}, got {veta}"
        );
    }

    #[test]
    fn test_veta_short_expiration() {
        let option = create_test_option(
            Positive::HUNDRED,  // underlying_price
            Positive::HUNDRED,  // strike_price
            pos_or_panic!(0.2), // implied_volatility
            Positive::ZERO,     // dividend_yield
            Positive::ONE,      // expiration_in_days
        );
        let veta = veta(&option).unwrap().to_f64().unwrap();
        let expected_veta = 0.0001511497;
        assert!(
            (veta - expected_veta).abs() < 1e-5,
            "Veta short expiration test failed: expected {expected_veta}, got {veta}"
        );
    }

    #[test]
    fn test_veta_with_dividends() {
        let option = create_test_option(
            Positive::HUNDRED,   // underlying_price
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.2),  // implied_volatility
            pos_or_panic!(0.03), // dividend_yield
            Positive::ONE,       // expiration_in_days
        );
        let veta = veta(&option).unwrap().to_f64().unwrap();
        let expected_veta = 0.0001511559;
        assert!(
            (veta - expected_veta).abs() < 1e-5,
            "Veta with dividends test failed: expected {expected_veta}, got {veta}"
        );
    }

    #[test]
    fn test_veta_itm() {
        let option = create_test_option(
            pos_or_panic!(110.0), // underlying_price
            Positive::HUNDRED,    // strike_price
            pos_or_panic!(0.2),   // implied_volatility
            Positive::ZERO,       // dividend_yield
            Positive::ONE,        // expiration_in_days
        );
        let veta = veta(&option).unwrap().to_f64().unwrap();
        let expected_veta = 0.0;
        assert!(
            (veta - expected_veta).abs() < 1e-5,
            "Veta ITM test failed: expected {expected_veta}, got {veta}"
        );
    }
}

#[cfg(test)]
pub mod tests_charm_equations {
    use super::*;
    use crate::model::types::{OptionStyle, Side};
    use crate::model::utils::create_sample_option_with_days;

    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::pos_or_panic;
    use tracing::info;

    #[test]
    fn test_charm_call_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,   // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Call ITM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), 0.00277350, epsilon = 1e-8);
    }

    #[test]
    fn test_charm_put_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            Positive::HUNDRED,   // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Put ITM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), -0.00392474, epsilon = 1e-8);
    }

    #[test]
    fn test_charm_call_atm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Call ATM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), -0.00045952, epsilon = 1e-8);
    }

    #[test]
    fn test_charm_put_atm() {
        let option = create_sample_option_with_days(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Put ATM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), -0.00048690, epsilon = 1e-8);
    }

    #[test]
    fn test_charm_call_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(90.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Call OTM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), -0.00401791, epsilon = 1e-8);
    }

    #[test]
    fn test_charm_put_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Put,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(90.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Put OTM Value: {}", charm_value);
        assert_relative_eq!(charm_value.to_f64().unwrap(), 0.00285007, epsilon = 1e-8);
    }
}

/// Tests for second-order volatility Greeks (Vanna, Vomma, Veta) edge cases.
///
/// These tests cover:
/// - High and low volatility environments
/// - Near expiration scenarios
/// - Extreme changes in underlying price (deep ITM/OTM)
#[cfg(test)]
pub mod tests_volatility_greeks_edge_cases {
    use super::*;
    use crate::model::types::{OptionStyle, Side};
    use crate::model::utils::create_sample_option_with_days;
    use positive::pos_or_panic;

    use tracing::info;

    // ==================== VANNA EDGE CASES ====================

    #[test]
    fn test_vanna_high_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.8), // High volatility (80%)
            pos_or_panic!(30.0),
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna High Volatility: {}", vanna_value);
        // Vanna should be smaller in absolute terms with high volatility
        assert!(vanna_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_vanna_low_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.05), // Low volatility (5%)
            pos_or_panic!(30.0),
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna Low Volatility: {}", vanna_value);
        // Vanna calculation should still work with low volatility
        assert!(vanna_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vanna_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ONE, // 1 day to expiration
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna Near Expiration: {}", vanna_value);
        // Near expiration, vanna should still be calculable
        assert!(vanna_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vanna_deep_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // Deep ITM (underlying >> strike)
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna Deep ITM: {}", vanna_value);
        // Deep ITM options have small vanna
        assert!(vanna_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_vanna_deep_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0), // Deep OTM (underlying << strike)
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna Deep OTM: {}", vanna_value);
        // Deep OTM options have small vanna
        assert!(vanna_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_vanna_zero_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            Positive::ZERO, // Zero volatility
            pos_or_panic!(30.0),
        );
        let vanna_value = vanna(&option).unwrap();
        info!("Vanna Zero Volatility: {}", vanna_value);
        // With zero volatility, vanna should be zero
        assert_eq!(vanna_value, Decimal::ZERO);
    }

    // ==================== VOMMA EDGE CASES ====================

    #[test]
    fn test_vomma_high_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.8), // High volatility (80%)
            pos_or_panic!(30.0),
        );
        let vomma_value = vomma(&option).unwrap();
        info!("Vomma High Volatility: {}", vomma_value);
        // Vomma should be calculable with high volatility
        assert!(vomma_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vomma_low_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.05), // Low volatility (5%)
            pos_or_panic!(30.0),
        );
        let vomma_value = vomma(&option).unwrap();
        info!("Vomma Low Volatility: {}", vomma_value);
        // Vomma should be calculable with low volatility
        assert!(vomma_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vomma_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ONE, // 1 day to expiration
        );
        let vomma_value = vomma(&option).unwrap();
        info!("Vomma Near Expiration: {}", vomma_value);
        // Near expiration, vomma should still be calculable
        assert!(vomma_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vomma_at_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO, // At expiration
        );
        let vomma_value = vomma(&option).unwrap();
        info!("Vomma At Expiration: {}", vomma_value);
        // At expiration, vomma should be zero
        assert_eq!(vomma_value, Decimal::ZERO);
    }

    #[test]
    fn test_vomma_deep_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0), // Deep OTM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );
        let vomma_value = vomma(&option).unwrap();
        info!("Vomma Deep OTM: {}", vomma_value);
        // Deep OTM options have highest vomma
        assert!(vomma_value.abs() < Decimal::MAX);
    }

    // ==================== VETA EDGE CASES ====================

    #[test]
    fn test_veta_high_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.8), // High volatility (80%)
            pos_or_panic!(30.0),
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta High Volatility: {}", veta_value);
        // Veta should be calculable with high volatility
        assert!(veta_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_veta_low_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.05), // Low volatility (5%)
            pos_or_panic!(30.0),
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta Low Volatility: {}", veta_value);
        // Veta should be calculable with low volatility
        assert!(veta_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_veta_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ONE, // 1 day to expiration
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta Near Expiration: {}", veta_value);
        // Near expiration, veta should still be calculable
        assert!(veta_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_veta_at_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            Positive::ZERO, // At expiration
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta At Expiration: {}", veta_value);
        // At expiration, veta should be zero
        assert_eq!(veta_value, Decimal::ZERO);
    }

    #[test]
    fn test_veta_deep_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // Deep ITM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta Deep ITM: {}", veta_value);
        // Deep ITM options have small veta
        assert!(veta_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_veta_deep_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0), // Deep OTM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta Deep OTM: {}", veta_value);
        // Deep OTM options have small veta
        assert!(veta_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_veta_long_dated_option() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(365.0), // 1 year to expiration
        );
        let veta_value = veta(&option).unwrap();
        info!("Veta Long Dated: {}", veta_value);
        // Long dated options should have calculable veta
        assert!(veta_value.abs() < Decimal::MAX);
    }

    // ==================== CHARM EDGE CASES ====================

    #[test]
    fn test_charm_high_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.8),  // High volatility (80%)
            pos_or_panic!(30.0), // expiration_days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm High Volatility: {}", charm_value);
        // Charm when is far from expiration is negligible
        assert!(charm_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_charm_low_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.05), // Low volatility (5%)
            pos_or_panic!(30.0), // expiration_days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Low Volatility: {}", charm_value);
        // Charm when is far from expiration is negligible
        // low volatility increases the value
        assert!(charm_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_charm_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED, // strike_price
            pos_or_panic!(0.2),
            Positive::ONE, // 1 day to expiration
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Near Expiration: {}", charm_value);
        // Near expiration Charm is increasing
        assert!(charm_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_charm_at_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED, // strike_price
            pos_or_panic!(0.2),
            Positive::ZERO, // At expiration
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm At Expiration: {}", charm_value);
        // At expiration, Charm should be zero
        assert_eq!(charm_value, Decimal::ZERO);
    }

    #[test]
    fn test_charm_deep_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // Deep ITM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0), // expiration_days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Deep ITM: {}", charm_value);
        // Deep ITM options far from expiration have small Charm
        assert!(charm_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_charm_deep_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0), // Deep OTM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0), // expiration_days
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Deep OTM: {}", charm_value);
        // For deep OTM options Charm should be zero
        assert_eq!(charm_value.abs(), Decimal::ZERO);
    }

    #[test]
    fn test_charm_long_dated_option() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(365.0), // 1 year to expiration
        );
        let charm_value = charm(&option).unwrap();
        info!("Charm Long Dated: {}", charm_value);
        // Long dated options should have calculable Charm
        assert!(charm_value.abs() < Decimal::ONE);
    }

    // ==================== COLOR EDGE CASES ====================

    #[test]
    fn test_color_high_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.8),  // High volatility (80%)
            pos_or_panic!(30.0), // expiration_days
        );
        let color_value = color(&option).unwrap();
        info!("Color High Volatility: {}", color_value);
        // Color when is far from expiration is negligible
        assert!(color_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_color_low_volatility() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED,   // strike_price
            pos_or_panic!(0.05), // Low volatility (5%)
            pos_or_panic!(30.0), // expiration_days
        );
        let color_value = color(&option).unwrap();
        info!("Color Low Volatility: {}", color_value);
        // Color when is far from expiration is negligible
        // low volatility increases the value
        assert!(color_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_color_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED, // strike_price
            pos_or_panic!(0.2),
            Positive::ONE, // 1 day to expiration
        );
        let color_value = color(&option).unwrap();
        info!("Color Near Expiration: {}", color_value);
        // Near expiration Color is increasing
        assert!(color_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_color_at_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED, // underlying_price
            Positive::ONE,
            Positive::HUNDRED, // strike_price
            pos_or_panic!(0.2),
            Positive::ZERO, // At expiration
        );
        let color_value = color(&option).unwrap();
        info!("Color At Expiration: {}", color_value);
        // At expiration, Color should be zero
        assert_eq!(color_value, Decimal::ZERO);
    }

    #[test]
    fn test_color_deep_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(150.0), // Deep ITM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0), // expiration_days
        );
        let color_value = color(&option).unwrap();
        info!("Color Deep ITM: {}", color_value);
        // Deep ITM options far from expiration have small Color
        assert!(color_value.abs() < Decimal::ONE);
    }

    #[test]
    fn test_color_deep_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(50.0), // Deep OTM
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0), // expiration_days
        );
        let color_value = color(&option).unwrap();
        info!("Color Deep OTM: {}", color_value);
        // For deep OTM options Color should be zero
        assert_eq!(color_value.abs(), Decimal::ZERO);
    }

    #[test]
    fn test_color_long_dated_option() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(365.0), // 1 year to expiration
        );
        let color_value = color(&option).unwrap();
        info!("Color Long Dated: {}", color_value);
        // Long dated options should have calculable Color
        assert!(color_value.abs() < Decimal::ONE);
    }

    // ==================== COMBINED SCENARIOS ====================

    #[test]
    fn test_volatility_greeks_extreme_scenario() {
        // High volatility + near expiration + ATM
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            Positive::ONE, // 100% volatility
            Positive::TWO, // 2 days to expiration
        );

        let vanna_value = vanna(&option).unwrap();
        let vomma_value = vomma(&option).unwrap();
        let veta_value = veta(&option).unwrap();
        let charm_value = charm(&option).unwrap();
        let color_value = color(&option).unwrap();

        info!("Extreme Scenario - Vanna: {}", vanna_value);
        info!("Extreme Scenario - Vomma: {}", vomma_value);
        info!("Extreme Scenario - Veta: {}", veta_value);
        info!("Extreme Scenario - Charm: {}", charm_value);
        info!("Extreme Scenario - Color: {}", color_value);

        // All should be finite values
        assert!(vanna_value.abs() < Decimal::MAX);
        assert!(vomma_value.abs() < Decimal::MAX);
        assert!(veta_value.abs() < Decimal::MAX);
        assert!(charm_value.abs() < Decimal::MAX);
        assert!(color_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_volatility_greeks_put_option() {
        let option = create_sample_option_with_days(
            OptionStyle::Put,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED,
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );

        let vanna_value = vanna(&option).unwrap();
        let vomma_value = vomma(&option).unwrap();
        let veta_value = veta(&option).unwrap();
        let charm_value = charm(&option).unwrap();
        let color_value = color(&option).unwrap();

        info!("Put Option - Vanna: {}", vanna_value);
        info!("Put Option - Vomma: {}", vomma_value);
        info!("Put Option - Veta: {}", veta_value);
        info!("Put Option - Charm: {}", charm_value);
        info!("Put Option - Color: {}", color_value);

        // All should be finite values
        assert!(vanna_value.abs() < Decimal::MAX);
        assert!(vomma_value.abs() < Decimal::MAX);
        assert!(veta_value.abs() < Decimal::MAX);
        assert!(charm_value.abs() < Decimal::MAX);
        assert!(color_value.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vanna_atm_vs_otm_comparison() {
        let atm_option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            Positive::HUNDRED, // ATM
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );

        let otm_option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,
            Positive::ONE,
            pos_or_panic!(110.0), // OTM
            pos_or_panic!(0.2),
            pos_or_panic!(30.0),
        );

        let vanna_atm = vanna(&atm_option).unwrap();
        let vanna_otm = vanna(&otm_option).unwrap();

        info!("Vanna ATM: {}", vanna_atm);
        info!("Vanna OTM: {}", vanna_otm);

        // Both should be calculable
        assert!(vanna_atm.abs() < Decimal::MAX);
        assert!(vanna_otm.abs() < Decimal::MAX);
    }

    #[test]
    fn test_vomma_smile_effect() {
        // Test vomma at different strikes to verify smile effect
        let strikes = vec![
            pos_or_panic!(90.0),
            pos_or_panic!(95.0),
            Positive::HUNDRED,
            pos_or_panic!(105.0),
            pos_or_panic!(110.0),
        ];

        for strike in strikes {
            let option = create_sample_option_with_days(
                OptionStyle::Call,
                Side::Long,
                Positive::HUNDRED,
                Positive::ONE,
                strike,
                pos_or_panic!(0.2),
                pos_or_panic!(30.0),
            );
            let vomma_value = vomma(&option).unwrap();
            info!("Vomma at strike {}: {}", strike, vomma_value);
            assert!(vomma_value.abs() < Decimal::MAX);
        }
    }
}

#[cfg(test)]
pub mod tests_color_equations {
    use super::*;
    use crate::model::types::{OptionStyle, Side};
    use crate::model::utils::create_sample_option_with_days;

    use approx::assert_relative_eq;
    use num_traits::ToPrimitive;
    use positive::pos_or_panic;
    use tracing::info;

    #[test]
    fn test_color_itm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            Positive::HUNDRED,   // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let color_value = color(&option).unwrap();
        info!("Color ITM Value: {}", color_value);
        assert_relative_eq!(color_value.to_f64().unwrap(), -0.00039105, epsilon = 1e-8);
    }

    #[test]
    fn test_color_atm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let color_value = color(&option).unwrap();
        info!("Color ATM Value: {}", color_value);
        assert_relative_eq!(color_value.to_f64().unwrap(), -0.00081635, epsilon = 1e-8);
    }

    #[test]
    fn test_color_atm_near_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(95.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(0.5),  // expiration days
        );
        let color_value = color(&option).unwrap();
        info!("Color ATM Near Expiration Value: {}", color_value);
        assert_relative_eq!(color_value.to_f64().unwrap(), -0.37822466, epsilon = 1e-8);
    }

    #[test]
    fn test_color_atm_right_before_expiration() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(95.0),  // underlying price
            Positive::ONE,        // quantity
            pos_or_panic!(95.0),  // strike price
            pos_or_panic!(0.3),   // volatility
            pos_or_panic!(0.001), // expiration days
        );
        let color_value = color(&option).unwrap();
        info!("Color ATM Right Before Expiration Value: {}", color_value);
        assert_relative_eq!(
            color_value.to_f64().unwrap(),
            -4228.45476344,
            epsilon = 1e-8
        );
    }

    #[test]
    fn test_color_otm() {
        let option = create_sample_option_with_days(
            OptionStyle::Call,
            Side::Long,
            pos_or_panic!(90.0), // underlying price
            Positive::ONE,       // quantity
            pos_or_panic!(95.0), // strike price
            pos_or_panic!(0.3),  // volatility
            pos_or_panic!(30.0), // expiration days
        );
        let color_value = color(&option).unwrap();
        info!("Color OTM Value: {}", color_value);
        assert_relative_eq!(color_value.to_f64().unwrap(), -0.00046416, epsilon = 1e-8);
    }
}