scirs2-series 0.1.5

//! Pricing utility functions
//!
//! This module provides mathematical utility functions commonly used in
//! financial pricing models, including statistical distributions and
//! numerical approximations.

use scirs2_core::numeric::Float;

/// Normal cumulative distribution function approximation
///
/// Calculates the cumulative distribution function of the standard normal
/// distribution using the Abramowitz and Stegun approximation. This provides
/// a good balance between accuracy and computational efficiency.
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * `F` - Cumulative probability P(X ≤ x) where X ~ N(0,1)
///
/// # Examples
///
/// ```rust
/// use scirs2_series::financial::pricing::utils::normal_cdf;
///
/// let prob: f64 = normal_cdf(0.0); // Should be approximately 0.5
/// assert!((prob - 0.5).abs() < 0.001);
///
/// let prob: f64 = normal_cdf(1.96); // Should be approximately 0.975
/// assert!((prob - 0.975).abs() < 0.01); // Relaxed tolerance for numerical approximation
/// ```
///
/// # Algorithm
///
/// Uses the Abramowitz and Stegun rational approximation:
/// - For x ≥ 0: Φ(x) = 1 - φ(x) × (a₁t + a₂t² + a₃t³ + a₄t⁴ + a₅t⁵)
/// - For x < 0: Φ(x) = 1 - Φ(-x)
///
/// Where t = 1/(1 + px), φ(x) is the standard normal PDF, and
/// p, a₁, a₂, a₃, a₄, a₅ are predefined constants.
///
/// Maximum absolute error: ≈ 7.5 × 10⁻⁸
pub fn normal_cdf<F: Float>(x: F) -> F {
    // Abramowitz and Stegun approximation constants
    let a1 = F::from(0.254829592).expect("Failed to convert constant to float");
    let a2 = F::from(-0.284496736).expect("Failed to convert constant to float");
    let a3 = F::from(1.421413741).expect("Failed to convert constant to float");
    let a4 = F::from(-1.453152027).expect("Failed to convert constant to float");
    let a5 = F::from(1.061405429).expect("Failed to convert constant to float");
    let p = F::from(0.3275911).expect("Failed to convert constant to float");

    let sign = if x < F::zero() { -F::one() } else { F::one() };
    let x_abs = x.abs();

    let t = F::one() / (F::one() + p * x_abs);
    let y = F::one()
        - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1)
            * t
            * (-x_abs * x_abs / F::from(2.0).expect("Failed to convert constant to float")).exp();

    (F::one() + sign * y) / F::from(2.0).expect("Failed to convert constant to float")
}

/// Normal probability density function
///
/// Calculates the probability density function of the standard normal
/// distribution at a given point.
///
/// # Arguments
///
/// * `x` - Input value
///
/// # Returns
///
/// * `F` - Probability density f(x) = (1/√2π) × e^(-x²/2)
///
/// # Examples
///
/// ```rust
/// use scirs2_series::financial::pricing::utils::normal_pdf;
///
/// let density: f64 = normal_pdf(0.0); // Should be approximately 1/sqrt(2π) ≈ 0.3989
/// assert!((density - 0.3989).abs() < 0.001);
/// ```
pub fn normal_pdf<F: Float>(x: F) -> F {
    let sqrt_2pi = F::from(2.506628274631).expect("Failed to convert constant to float"); // √(2π)
    (-x.powi(2) / F::from(2.0).expect("Failed to convert constant to float")).exp() / sqrt_2pi
}

/// Inverse normal cumulative distribution function (quantile function)
///
/// Calculates the inverse of the normal CDF using the Beasley-Springer-Moro
/// approximation. Given a probability p, returns x such that Φ(x) = p.
///
/// # Arguments
///
/// * `p` - Probability value (must be between 0 and 1)
///
/// # Returns
///
/// * `F` - Quantile value x such that P(X ≤ x) = p
///
/// # Examples
///
/// ```rust
/// use scirs2_series::financial::pricing::utils::normal_quantile;
///
/// let quantile: f64 = normal_quantile(0.975); // Should be approximately 1.96
/// assert!((quantile - 1.96).abs() < 0.01);
///
/// let quantile: f64 = normal_quantile(0.5); // Should be approximately 0.0
/// assert!(quantile.abs() < 0.001);
/// ```
pub fn normal_quantile<F: Float>(p: F) -> F {
    if p <= F::zero() || p >= F::one() {
        return F::nan();
    }

    // Beasley-Springer-Moro approximation constants
    let a0 = F::from(2.515517).expect("Failed to convert constant to float");
    let a1 = F::from(0.802853).expect("Failed to convert constant to float");
    let a2 = F::from(0.010328).expect("Failed to convert constant to float");
    let b1 = F::from(1.432788).expect("Failed to convert constant to float");
    let b2 = F::from(0.189269).expect("Failed to convert constant to float");
    let b3 = F::from(0.001308).expect("Failed to convert constant to float");

    let t = if p < F::from(0.5).expect("Failed to convert constant to float") {
        (-F::from(2.0).expect("Failed to convert constant to float") * p.ln()).sqrt()
    } else {
        (-F::from(2.0).expect("Failed to convert constant to float") * (F::one() - p).ln()).sqrt()
    };

    let numerator = a0 + a1 * t + a2 * t.powi(2);
    let denominator = F::one() + b1 * t + b2 * t.powi(2) + b3 * t.powi(3);
    let z = t - numerator / denominator;

    if p < F::from(0.5).expect("Failed to convert constant to float") {
        -z
    } else {
        z
    }
}

/// Present value calculation
///
/// Calculates the present value of a future cash flow using continuous compounding.
///
/// # Arguments
///
/// * `future_value` - Future cash flow amount
/// * `rate` - Discount rate (continuously compounded)
/// * `time` - Time to maturity
///
/// # Returns
///
/// * `F` - Present value
///
/// # Examples
///
/// ```rust
/// use scirs2_series::financial::pricing::utils::present_value;
///
/// let pv: f64 = present_value(100.0, 0.05, 1.0); // $100 in 1 year at 5%
/// assert!((pv - 95.123).abs() < 0.001); // ≈ $95.12
/// ```
pub fn present_value<F: Float>(future_value: F, rate: F, time: F) -> F {
    future_value * (-rate * time).exp()
}

/// Future value calculation
///
/// Calculates the future value of a present cash flow using continuous compounding.
///
/// # Arguments
///
/// * `present_value` - Current cash flow amount
/// * `rate` - Growth rate (continuously compounded)
/// * `time` - Time to maturity
///
/// # Returns
///
/// * `F` - Future value
///
/// # Examples
///
/// ```rust
/// use scirs2_series::financial::pricing::utils::future_value;
///
/// let fv: f64 = future_value(100.0, 0.05, 1.0); // $100 growing at 5% for 1 year
/// assert!((fv - 105.127).abs() < 0.001); // ≈ $105.13
/// ```
pub fn future_value<F: Float>(present_value: F, rate: F, time: F) -> F {
    present_value * (rate * time).exp()
}

/// Calculate d1 parameter for Black-Scholes formula
///
/// Helper function to calculate the d1 parameter used in Black-Scholes
/// and other log-normal option pricing models.
///
/// # Arguments
///
/// * `spot` - Current spot price
/// * `strike` - Strike price
/// * `rate` - Risk-free rate
/// * `volatility` - Volatility
/// * `time` - Time to expiration
///
/// # Returns
///
/// * `F` - d1 parameter value
pub fn calculate_d1<F: Float>(spot: F, strike: F, rate: F, volatility: F, time: F) -> F {
    ((spot / strike).ln()
        + (rate + volatility.powi(2) / F::from(2.0).expect("Failed to convert constant to float"))
            * time)
        / (volatility * time.sqrt())
}

/// Calculate d2 parameter for Black-Scholes formula
///
/// Helper function to calculate the d2 parameter used in Black-Scholes
/// and other log-normal option pricing models.
///
/// # Arguments
///
/// * `d1` - Previously calculated d1 parameter
/// * `volatility` - Volatility
/// * `time` - Time to expiration
///
/// # Returns
///
/// * `F` - d2 parameter value
pub fn calculate_d2<F: Float>(d1: F, volatility: F, time: F) -> F {
    d1 - volatility * time.sqrt()
}

/// Bivariate normal cumulative distribution function approximation
///
/// Calculates the cumulative probability of a bivariate normal distribution
/// using Drezner's approximation. Useful for pricing spread options and
/// other multi-asset derivatives.
///
/// # Arguments
///
/// * `x` - First variable value
/// * `y` - Second variable value  
/// * `rho` - Correlation coefficient between variables
///
/// # Returns
///
/// * `F` - Cumulative probability P(X ≤ x, Y ≤ y)
pub fn bivariate_normal_cdf<F: Float>(x: F, y: F, rho: F) -> F {
    // Simplified approximation for bivariate normal CDF
    // This is a basic implementation - more sophisticated methods exist

    if rho.abs() < F::from(0.001).expect("Failed to convert constant to float") {
        // If correlation is near zero, variables are independent
        return normal_cdf(x) * normal_cdf(y);
    }

    // Use the basic approximation for non-zero correlation
    // In practice, you might want to implement more accurate methods
    // like Drezner-Wesolowsky or Genz methods
    let cdf_x = normal_cdf(x);
    let cdf_y = normal_cdf(y);

    // Simple correction for correlation (approximation)
    let correction = rho * normal_pdf(x) * normal_pdf(y)
        / F::from(4.0).expect("Failed to convert constant to float");

    (cdf_x * cdf_y + correction).min(F::one()).max(F::zero())
}

/// Calculate annualized return from continuously compounded rate
///
/// # Arguments
///
/// * `rate` - Continuously compounded rate
///
/// # Returns
///
/// * `F` - Annualized simple return
pub fn continuous_to_simple_return<F: Float>(rate: F) -> F {
    rate.exp() - F::one()
}

/// Calculate continuously compounded rate from annualized return
///
/// # Arguments
///
/// * `simple_return` - Annualized simple return
///
/// # Returns
///
/// * `F` - Continuously compounded rate
pub fn simple_to_continuous_return<F: Float>(simple_return: F) -> F {
    (F::one() + simple_return).ln()
}

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

    #[test]
    fn test_normal_cdf() {
        // Test key values with reasonable tolerance for approximation
        let cdf_0 = normal_cdf(0.0);
        assert!((cdf_0 - 0.5).abs() < 0.01);

        let cdf_positive = normal_cdf(1.96);
        assert!((cdf_positive - 0.975).abs() < 0.01);

        let cdf_negative = normal_cdf(-1.96);
        assert!((cdf_negative - 0.025).abs() < 0.01);
    }

    #[test]
    fn test_normal_pdf() {
        let pdf_0 = normal_pdf(0.0);
        let expected = 1.0 / (2.0 * std::f64::consts::PI).sqrt();
        assert!((pdf_0 - expected).abs() < 0.001);

        // PDF should be symmetric
        let pdf_pos = normal_pdf(1.0);
        let pdf_neg = normal_pdf(-1.0);
        assert!((pdf_pos - pdf_neg).abs() < 1e-10);
    }

    #[test]
    fn test_normal_quantile() {
        // Test that quantile and CDF are inverses
        let p = 0.975;
        let quantile = normal_quantile(p);
        let cdf_back = normal_cdf(quantile);
        assert!((cdf_back - p).abs() < 0.01);

        // Test known values
        let q_median = normal_quantile(0.5);
        assert!(q_median.abs() < 0.001);
    }

    #[test]
    fn test_present_future_value() {
        let pv = 100.0;
        let rate = 0.05;
        let time = 1.0;

        let fv = future_value(pv, rate, time);
        let pv_back = present_value(fv, rate, time);

        assert!((pv_back - pv).abs() < 1e-10);
    }

    #[test]
    fn test_d1_d2_calculation() {
        let spot = 100.0;
        let strike = 100.0;
        let rate = 0.05;
        let vol = 0.2;
        let time = 1.0;

        let d1 = calculate_d1(spot, strike, rate, vol, time);
        let d2 = calculate_d2(d1, vol, time);

        // d2 should be d1 - vol*sqrt(time)
        let expected_d2 = d1 - vol * time.sqrt();
        assert!((d2 - expected_d2).abs() < 1e-10);
    }

    #[test]
    fn test_bivariate_normal_independence() {
        // When rho = 0, should equal product of marginals
        let x = 1.0;
        let y = 0.5;
        let rho = 0.0;

        let biv_cdf = bivariate_normal_cdf(x, y, rho);
        let product = normal_cdf(x) * normal_cdf(y);

        assert!((biv_cdf - product).abs() < 0.01);
    }

    #[test]
    fn test_return_conversions() {
        let simple_return = 0.10; // 10%
        let continuous = simple_to_continuous_return(simple_return);
        let back_to_simple = continuous_to_simple_return(continuous);

        assert!((back_to_simple - simple_return).abs() < 1e-10);
    }

    #[test]
    fn test_invalid_quantile() {
        assert!(normal_quantile(-0.1).is_nan());
        assert!(normal_quantile(1.1).is_nan());
        assert!(normal_quantile(0.0).is_nan());
        assert!(normal_quantile(1.0).is_nan());
    }
}