Skip to main content

hypercall_margin/
black_scholes.rs

1//! Black-Scholes option pricing with Edgeworth expansion for skew and kurtosis.
2//!
3//! Provides vanilla BS pricing, moment-adjusted pricing, and Greeks computation.
4//! All functions are pure: no state, no IO, deterministic output.
5
6use crate::types::OptionType;
7
8fn normal_cdf(x: f64) -> f64 {
9    0.5 * (1.0 + libm::erf(x / std::f64::consts::SQRT_2))
10}
11
12fn normal_pdf(x: f64) -> f64 {
13    (-0.5 * x * x).exp() / (2.0 * std::f64::consts::PI).sqrt()
14}
15
16fn adjusted_cdf(x: f64, skew: f64, excess_kurtosis: f64) -> f64 {
17    if skew == 0.0 && excess_kurtosis == 0.0 {
18        return normal_cdf(x);
19    }
20
21    let base_cdf = normal_cdf(x);
22    let pdf = normal_pdf(x);
23    let skew_adjustment = skew * pdf * (x * x - 1.0) / 6.0;
24    let kurtosis_adjustment = excess_kurtosis * pdf * (x * x * x * x - 6.0 * x * x + 3.0) / 24.0;
25
26    (base_cdf + skew_adjustment + kurtosis_adjustment).clamp(0.0, 1.0)
27}
28
29pub fn black_scholes(
30    option_type: &OptionType,
31    spot: f64,
32    strike: f64,
33    time_to_expiry: f64,
34    risk_free_rate: f64,
35    volatility: f64,
36) -> f64 {
37    black_scholes_with_moments(
38        option_type,
39        spot,
40        strike,
41        time_to_expiry,
42        risk_free_rate,
43        volatility,
44        0.0,
45        0.0,
46    )
47}
48
49#[allow(clippy::too_many_arguments)]
50pub fn black_scholes_with_moments(
51    option_type: &OptionType,
52    spot: f64,
53    strike: f64,
54    time_to_expiry: f64,
55    risk_free_rate: f64,
56    volatility: f64,
57    skew: f64,
58    excess_kurtosis: f64,
59) -> f64 {
60    if time_to_expiry <= 0.0 {
61        return match option_type {
62            OptionType::Call => (spot - strike).max(0.0),
63            OptionType::Put => (strike - spot).max(0.0),
64        };
65    }
66
67    let d1 = ((spot / strike).ln()
68        + (risk_free_rate + 0.5 * volatility * volatility) * time_to_expiry)
69        / (volatility * time_to_expiry.sqrt());
70    let d2 = d1 - volatility * time_to_expiry.sqrt();
71
72    // Edgeworth expansion can produce negative prices for large moment adjustments.
73    // Clamp to the no-arbitrage bounds: call in [0, spot], put in [0, K*exp(-rT)].
74    let discount = (-risk_free_rate * time_to_expiry).exp();
75    match option_type {
76        OptionType::Call => {
77            let raw = spot * adjusted_cdf(d1, skew, excess_kurtosis)
78                - strike * discount * adjusted_cdf(d2, skew, excess_kurtosis);
79            raw.clamp(0.0, spot)
80        }
81        OptionType::Put => {
82            let raw = strike * discount * adjusted_cdf(-d2, skew, excess_kurtosis)
83                - spot * adjusted_cdf(-d1, skew, excess_kurtosis);
84            raw.clamp(0.0, strike * discount)
85        }
86    }
87}
88
89pub fn calculate_greeks(
90    option_type: &OptionType,
91    spot: f64,
92    strike: f64,
93    time_to_expiry: f64,
94    risk_free_rate: f64,
95    volatility: f64,
96) -> (f64, f64, f64, f64, f64) {
97    if time_to_expiry <= 0.0 {
98        let intrinsic = match option_type {
99            OptionType::Call => (spot - strike).max(0.0),
100            OptionType::Put => (strike - spot).max(0.0),
101        };
102        let delta = match option_type {
103            OptionType::Call => {
104                if spot > strike {
105                    1.0
106                } else {
107                    0.0
108                }
109            }
110            OptionType::Put => {
111                if spot < strike {
112                    -1.0
113                } else {
114                    0.0
115                }
116            }
117        };
118        return (intrinsic, delta, 0.0, 0.0, 0.0);
119    }
120
121    let sqrt_t = time_to_expiry.sqrt();
122    let d1 = ((spot / strike).ln()
123        + (risk_free_rate + 0.5 * volatility * volatility) * time_to_expiry)
124        / (volatility * sqrt_t);
125    let d2 = d1 - volatility * sqrt_t;
126    let discount_factor = (-risk_free_rate * time_to_expiry).exp();
127
128    let price = match option_type {
129        OptionType::Call => spot * normal_cdf(d1) - strike * discount_factor * normal_cdf(d2),
130        OptionType::Put => strike * discount_factor * normal_cdf(-d2) - spot * normal_cdf(-d1),
131    };
132
133    let delta = match option_type {
134        OptionType::Call => normal_cdf(d1),
135        OptionType::Put => -normal_cdf(-d1),
136    };
137
138    let gamma = normal_pdf(d1) / (spot * volatility * sqrt_t);
139    let vega = spot * normal_pdf(d1) * sqrt_t / 100.0;
140    let theta = match option_type {
141        OptionType::Call => {
142            (-spot * normal_pdf(d1) * volatility / (2.0 * sqrt_t)
143                - risk_free_rate * strike * discount_factor * normal_cdf(d2))
144                / 365.0
145        }
146        OptionType::Put => {
147            (-spot * normal_pdf(d1) * volatility / (2.0 * sqrt_t)
148                + risk_free_rate * strike * discount_factor * normal_cdf(-d2))
149                / 365.0
150        }
151    };
152
153    (price, delta, gamma, vega, theta)
154}