blacksholes.rs

![feature(const_fn_floating_point_arithmetic)]

fn rand() -> f64 {
    const RAND_SEED: u64 = 69;
    // https://en.wikipedia.org/wiki/Permuted_congruential_generator
    thread_local! {
        static INIT: std::cell::Cell<u64> = std::cell::Cell::new(RAND_SEED);
    }
    INIT.with(|cell| {
        let s = cell.get();
        cell.set(
            s.wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407),
        );
        f64::from((((s ^ (s >> 18)) >> 27) as u32).rotate_right((s >> 59) as u32))
    })
}

#[inline]
const fn float_max(a: f64, b: f64) -> f64 {
    if a >= b {
        a
    } else {
        b
    }
}

fn gaussian_box_muller() -> f64 {
    loop {
        let x = 2.0f64 * (rand() / std::u32::MAX as f64) - 1.0f64;
        let y = 2.0f64 * (rand() / std::u32::MAX as f64) - 1.0f64;
        let euclid_sq = x * x + y * y;
        if euclid_sq <= 1.0 {
            break x * (-2.0f64 * euclid_sq.ln() / euclid_sq).sqrt();
        }
    }
}

/// Black-Scholes European Option
#[derive(Debug, Clone, Copy, PartialEq)]
struct EuropeanOption {
    /// underlying price
    asset_price: f64,
    /// strike price
    strike: f64,
    /// risk free rate
    rfr: f64,
    /// asset volatility
    volatility: f64,
}

impl EuropeanOption {
    fn monte_carlo_call_price(&self, num_sims: u32, expiry: f64) -> f64 {
        let Self {
            asset_price,
            strike,
            rfr,
            volatility: sigma,
        } = self;
        let uprice_adjust = asset_price * (expiry * (rfr - 0.5 * sigma * sigma)).exp();
        let mut payoff_sum = 0.0;

        for _ in 0..num_sims {
            let uprice_cur =
                uprice_adjust * ((sigma * sigma * expiry).sqrt() * gaussian_box_muller()).exp();
            payoff_sum += float_max(uprice_cur - strike, 0.0);
        }
        payoff_sum / f64::from(num_sims) * (-rfr * expiry).exp()
    }

    fn monte_carlo_put_price(&self, num_sims: u32, expiry: f64) -> f64 {
        let Self {
            asset_price,
            strike,
            rfr,
            volatility: sigma,
        } = self;
        let uprice_adjust = asset_price * (expiry * (rfr - 0.5 * sigma * sigma)).exp();
        let mut payoff_sum = 0.0;

        for _ in 0..num_sims {
            let uprice_cur =
                uprice_adjust * ((sigma * sigma * expiry).sqrt() * gaussian_box_muller()).exp();
            payoff_sum += float_max(strike - uprice_cur, 0.0);
        }
        payoff_sum / f64::from(num_sims) * (-rfr * expiry).exp()
    }
}

fn main() {
    const OPTION: EuropeanOption = EuropeanOption {
        asset_price: 100.0, // Option price
        strike: 100.0,      // Strike price
        rfr: 0.05,          // Risk-free rate
        volatility: 0.2,    // Volatility of the underlying (20%)
    };

    const SIMULATIONS: u32 = 10_000; // Number of simulated asset paths
    const T: f64 = 1.0; // One year until expiry
    let call = OPTION.monte_carlo_call_price(SIMULATIONS, T);
    let put = OPTION.monte_carlo_put_price(SIMULATIONS, T);

    println!(
        r#"Parameters:
- Paths: {paths},
- Maturity (years): {maturity},
- Asset: {option:#?}"#,
        paths = SIMULATIONS,
        maturity = T,
        option = OPTION,
    );
    println!("Call Price: {}", call);
    println!("Put Price: {}", put);
}