hypercall/shared/

black_scholes_utils.rs

use crate::types::OptionType;
use hypercall_types::Greeks;

fn normal_cdf(x: f64) -> f64 {
    0.5 * (1.0 + libm::erf(x / std::f64::consts::SQRT_2))
}

fn normal_pdf(x: f64) -> f64 {
    (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt()
}

// Adjusted CDF using a simplified approach for skew and kurtosis
fn adjusted_cdf(x: f64, skew: f64, excess_kurtosis: f64) -> f64 {
    if skew == 0.0 && excess_kurtosis == 0.0 {
        return normal_cdf(x);
    }

    // Get base normal CDF value
    let base_cdf = normal_cdf(x);

    // Apply skew adjustment using Edgeworth expansion approximation
    // Positive skew increases right tail probability
    let pdf = normal_pdf(x);
    let skew_adjustment = skew * pdf * (x * x - 1.0) / 6.0;

    // Apply kurtosis adjustment
    // Positive excess kurtosis increases tail probabilities
    let kurtosis_adjustment = excess_kurtosis * pdf * (x * x * x * x - 6.0 * x * x + 3.0) / 24.0;

    // Apply adjustments and ensure result is in [0, 1]
    #[allow(clippy::manual_clamp)]
    (base_cdf + skew_adjustment + kurtosis_adjustment)
        .max(0.0)
        .min(1.0)
}

// Adjusted PDF for skew and kurtosis
#[allow(dead_code)]
fn adjusted_pdf(x: f64, skew: f64, excess_kurtosis: f64) -> f64 {
    if skew == 0.0 && excess_kurtosis == 0.0 {
        return normal_pdf(x);
    }

    let n = normal_pdf(x);
    let adjustment = 1.0
        + (skew / 6.0) * (x * x * x - 3.0 * x)
        + (excess_kurtosis / 24.0) * (x * x * x * x - 6.0 * x * x + 3.0)
        + (skew * skew / 72.0)
            * (x * x * x * x * x * x - 15.0 * x * x * x * x + 45.0 * x * x - 15.0);

    n * adjustment.max(0.0) // Ensure non-negative PDF
}

pub fn black_scholes(
    option_type: &OptionType,
    spot: f64,
    strike: f64,
    time_to_expiry: f64,
    risk_free_rate: f64,
    volatility: f64,
) -> f64 {
    black_scholes_with_moments(
        option_type,
        spot,
        strike,
        time_to_expiry,
        risk_free_rate,
        volatility,
        0.0,
        0.0,
    )
}

pub fn black_scholes_with_moments(
    option_type: &OptionType,
    spot: f64,
    strike: f64,
    time_to_expiry: f64,
    risk_free_rate: f64,
    volatility: f64,
    skew: f64,
    excess_kurtosis: f64,
) -> f64 {
    if time_to_expiry <= 0.0 {
        return match option_type {
            OptionType::Call => (spot - strike).max(0.0),
            OptionType::Put => (strike - spot).max(0.0),
        };
    }

    let d1 = ((spot / strike).ln()
        + (risk_free_rate + 0.5 * volatility * volatility) * time_to_expiry)
        / (volatility * time_to_expiry.sqrt());
    let d2 = d1 - volatility * time_to_expiry.sqrt();

    match option_type {
        OptionType::Call => {
            spot * adjusted_cdf(d1, skew, excess_kurtosis)
                - strike
                    * (-risk_free_rate * time_to_expiry).exp()
                    * adjusted_cdf(d2, skew, excess_kurtosis)
        }
        OptionType::Put => {
            strike
                * (-risk_free_rate * time_to_expiry).exp()
                * adjusted_cdf(-d2, skew, excess_kurtosis)
                - spot * adjusted_cdf(-d1, skew, excess_kurtosis)
        }
    }
}

pub fn calculate_greeks(
    option_type: &OptionType,
    spot: f64,
    strike: f64,
    time_to_expiry: f64,
    risk_free_rate: f64,
    volatility: f64,
) -> Greeks {
    if time_to_expiry <= 0.0 {
        let intrinsic = match option_type {
            OptionType::Call => (spot - strike).max(0.0),
            OptionType::Put => (strike - spot).max(0.0),
        };
        let delta = match option_type {
            OptionType::Call => {
                if spot > strike {
                    1.0
                } else {
                    0.0
                }
            }
            OptionType::Put => {
                if spot < strike {
                    -1.0
                } else {
                    0.0
                }
            }
        };
        return Greeks {
            delta,
            gamma: 0.0,
            theta: 0.0,
            vega: 0.0,
            rho: 0.0,
            implied_vol: volatility,
            theoretical_price: intrinsic,
            market_mid_price: None,
        };
    }

    let sqrt_t = time_to_expiry.sqrt();
    let d1 = ((spot / strike).ln()
        + (risk_free_rate + 0.5 * volatility * volatility) * time_to_expiry)
        / (volatility * sqrt_t);
    let d2 = d1 - volatility * sqrt_t;
    let discount_factor = (-risk_free_rate * time_to_expiry).exp();

    let price = match option_type {
        OptionType::Call => spot * normal_cdf(d1) - strike * discount_factor * normal_cdf(d2),
        OptionType::Put => strike * discount_factor * normal_cdf(-d2) - spot * normal_cdf(-d1),
    };

    let delta = match option_type {
        OptionType::Call => normal_cdf(d1),
        OptionType::Put => -normal_cdf(-d1),
    };

    let gamma = normal_pdf(d1) / (spot * volatility * sqrt_t);
    let vega = spot * normal_pdf(d1) * sqrt_t / 100.0; // Divided by 100 for 1% vol move
    let theta = match option_type {
        OptionType::Call => {
            (-spot * normal_pdf(d1) * volatility / (2.0 * sqrt_t)
                - risk_free_rate * strike * discount_factor * normal_cdf(d2))
                / 365.0
        }
        OptionType::Put => {
            (-spot * normal_pdf(d1) * volatility / (2.0 * sqrt_t)
                + risk_free_rate * strike * discount_factor * normal_cdf(-d2))
                / 365.0
        }
    };

    let rho = match option_type {
        OptionType::Call => time_to_expiry * strike * discount_factor * normal_cdf(d2) / 100.0,
        OptionType::Put => -time_to_expiry * strike * discount_factor * normal_cdf(-d2) / 100.0,
    };

    Greeks {
        delta,
        gamma,
        theta,
        vega,
        rho,
        implied_vol: volatility,
        theoretical_price: price,
        market_mid_price: None,
    }
}

// Calculate implied volatility using Newton-Raphson method
pub fn calculate_implied_volatility(
    option_type: &OptionType,
    spot: f64,
    strike: f64,
    time_to_expiry: f64,
    risk_free_rate: f64,
    market_price: f64,
    initial_vol: Option<f64>,
) -> Option<f64> {
    if time_to_expiry <= 0.0 {
        return None;
    }

    // Check if price is within valid bounds
    let intrinsic = match option_type {
        OptionType::Call => (spot - strike).max(0.0),
        OptionType::Put => (strike - spot).max(0.0),
    };

    if market_price < intrinsic {
        return None;
    }

    let mut vol = initial_vol.unwrap_or(0.3); // Start with 30% if no initial guess
    let max_iterations = 50;
    let tolerance = 1e-6;

    for _ in 0..max_iterations {
        let greeks = calculate_greeks(
            option_type,
            spot,
            strike,
            time_to_expiry,
            risk_free_rate,
            vol,
        );
        let price_diff = greeks.theoretical_price - market_price;

        if price_diff.abs() < tolerance {
            return Some(vol);
        }

        // Newton-Raphson update
        let vega_full = greeks.vega * 100.0; // Convert back to full vega
        if vega_full.abs() < 1e-10 {
            break;
        }

        vol -= price_diff / vega_full;

        // Keep volatility in reasonable bounds
        #[allow(clippy::manual_clamp)]
        {
            vol = vol.max(0.001).min(5.0);
        }
    }

    None
}