Coverage for src / qsmile / core / black76.py: 98%
89 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-05-01 22:47 +0000
1"""Black76 forward option pricing and implied volatility inversion."""
3from __future__ import annotations
5import numpy as np
6from numpy.typing import ArrayLike, NDArray
7from scipy.stats import invgauss, norm
10def _validate_common(
11 forward: ArrayLike,
12 strike: ArrayLike,
13 discount_factor: ArrayLike,
14 expiry: float,
15) -> None:
16 """Validate common inputs for Black76 functions."""
17 if expiry <= 0:
18 msg = f"expiry must be positive, got {expiry}"
19 raise ValueError(msg)
20 forward_arr = np.asarray(forward)
21 strike_arr = np.asarray(strike)
22 df_arr = np.asarray(discount_factor)
23 if np.any(forward_arr <= 0):
24 msg = "forward must be positive"
25 raise ValueError(msg)
26 if np.any(strike_arr <= 0):
27 msg = "strike must be positive"
28 raise ValueError(msg)
29 if np.any(df_arr <= 0):
30 msg = "discount_factor must be positive"
31 raise ValueError(msg)
34def black76_call(
35 forward: ArrayLike,
36 strike: ArrayLike,
37 discount_factor: ArrayLike,
38 vol: ArrayLike,
39 expiry: float,
40) -> NDArray[np.float64] | np.floating:
41 """Compute Black76 call option price.
43 C = D * [F * Phi(d1) - K * Phi(d2)]
45 Parameters
46 ----------
47 forward : ArrayLike
48 Forward price. Must be positive.
49 strike : ArrayLike
50 Strike price. Must be positive.
51 discount_factor : ArrayLike
52 Discount factor. Must be positive.
53 vol : ArrayLike
54 Volatility. Must be non-negative.
55 expiry : float
56 Time to expiry in years. Must be positive.
57 """
58 _validate_common(forward, strike, discount_factor, expiry)
59 vol_arr = np.asarray(vol, dtype=np.float64)
60 if np.any(vol_arr < 0):
61 msg = "vol must be non-negative"
62 raise ValueError(msg)
64 F = np.asarray(forward, dtype=np.float64)
65 K = np.asarray(strike, dtype=np.float64)
66 D = np.asarray(discount_factor, dtype=np.float64)
67 sqrt_t = np.sqrt(expiry)
69 # Handle zero vol case
70 zero_vol = vol_arr == 0.0
71 safe_vol = np.where(zero_vol, 1.0, vol_arr) # avoid division by zero
73 d1 = (np.log(F / K) + 0.5 * safe_vol**2 * expiry) / (safe_vol * sqrt_t)
74 d2 = d1 - safe_vol * sqrt_t
76 price = D * (F * norm.cdf(d1) - K * norm.cdf(d2))
78 # Zero vol: intrinsic value
79 intrinsic = D * np.maximum(F - K, 0.0)
80 return np.where(zero_vol, intrinsic, price)
83def black76_put(
84 forward: ArrayLike,
85 strike: ArrayLike,
86 discount_factor: ArrayLike,
87 vol: ArrayLike,
88 expiry: float,
89) -> NDArray[np.float64] | np.floating:
90 """Compute Black76 put option price.
92 P = D * [K * Phi(-d2) - F * Phi(-d1)]
94 Parameters
95 ----------
96 forward : ArrayLike
97 Forward price. Must be positive.
98 strike : ArrayLike
99 Strike price. Must be positive.
100 discount_factor : ArrayLike
101 Discount factor. Must be positive.
102 vol : ArrayLike
103 Volatility. Must be non-negative.
104 expiry : float
105 Time to expiry in years. Must be positive.
106 """
107 _validate_common(forward, strike, discount_factor, expiry)
108 vol_arr = np.asarray(vol, dtype=np.float64)
109 if np.any(vol_arr < 0):
110 msg = "vol must be non-negative"
111 raise ValueError(msg)
113 F = np.asarray(forward, dtype=np.float64)
114 K = np.asarray(strike, dtype=np.float64)
115 D = np.asarray(discount_factor, dtype=np.float64)
116 sqrt_t = np.sqrt(expiry)
118 zero_vol = vol_arr == 0.0
119 safe_vol = np.where(zero_vol, 1.0, vol_arr)
121 d1 = (np.log(F / K) + 0.5 * safe_vol**2 * expiry) / (safe_vol * sqrt_t)
122 d2 = d1 - safe_vol * sqrt_t
124 price = D * (K * norm.cdf(-d2) - F * norm.cdf(-d1))
126 intrinsic = D * np.maximum(K - F, 0.0)
127 return np.where(zero_vol, intrinsic, price)
130def black76_implied_vol(
131 price: float,
132 forward: float,
133 strike: float,
134 discount_factor: float,
135 expiry: float,
136 *,
137 is_call: bool,
138 tol: float = 1e-12,
139 max_vol: float = 10.0,
140) -> float:
141 """Invert Black76 to recover implied volatility.
143 Uses the explicit closed-form solution of Schadner (2026), "An Explicit
144 Solution to Black-Scholes Implied Volatility" (arXiv:2604.24480), which
145 expresses implied volatility as a direct transform of the option price
146 via the inverse Gaussian quantile function -- no root finding required.
148 For a call with normalized price ``c = C / (D F)`` and forward
149 log-moneyness ``k = log(K/F)``::
151 sigma = 2 / sqrt(T * F_IG^{-1}((1-c)/m; 2/|k|, 1))
153 where ``m = 1`` if ``K > F`` and ``m = K/F`` if ``K < F``. At ``k = 0``
154 the formula collapses to ``sigma = (2/sqrt(T)) * Phi^{-1}((c+1)/2)``.
155 The put case follows from put-call parity.
157 Parameters
158 ----------
159 price : float
160 Observed option price.
161 forward : float
162 Forward price. Must be positive.
163 strike : float
164 Strike price. Must be positive.
165 discount_factor : float
166 Discount factor. Must be positive.
167 expiry : float
168 Time to expiry in years. Must be positive.
169 is_call : bool
170 True for call, False for put.
171 tol : float
172 Tolerance for arbitrage-bound checks and the intrinsic-value
173 short-circuit (returns ``0.0``).
174 max_vol : float
175 Retained for backward compatibility. The closed-form solution does
176 not perform a search, so this argument is unused.
177 """
178 del max_vol # unused; closed-form solution does not search
179 _validate_common(forward, strike, discount_factor, expiry)
181 # No-arbitrage bounds
182 if is_call:
183 intrinsic = discount_factor * max(forward - strike, 0.0)
184 upper_bound = discount_factor * forward
185 else:
186 intrinsic = discount_factor * max(strike - forward, 0.0)
187 upper_bound = discount_factor * strike
189 if price < intrinsic - tol:
190 msg = f"price {price} is below intrinsic value {intrinsic}"
191 raise ValueError(msg)
192 if price > upper_bound + tol:
193 msg = f"price {price} exceeds upper bound {upper_bound}"
194 raise ValueError(msg)
196 # Edge case: price equals intrinsic -> vol is 0
197 if price <= intrinsic + tol:
198 return 0.0
200 # Normalized price c = C / (D * F); k = log(K / F)
201 df_f = discount_factor * forward
202 k = float(np.log(strike / forward))
204 # ATM case: explicit Gaussian-quantile form (Schadner Eq. 2)
205 if abs(k) < 1e-15:
206 c = price / df_f
207 v = 2.0 * float(norm.ppf((c + 1.0) / 2.0))
208 return v / float(np.sqrt(expiry))
210 # Compute the IG quantile argument q in a numerically stable form that
211 # avoids catastrophic cancellation when the option is deep ITM (i.e. when
212 # the normalized price c is close to 1). Schadner's m factor folds the
213 # ITM case into the OTM one via parity; combined with the price/(D*F)
214 # normalization this yields:
215 # call: q = (D*F - price) / (D * min(F, K))
216 # put : q = (D*K - price) / (D * min(F, K))
217 # We additionally compute the complementary probability qc = 1 - q in a
218 # form that does not cancel, so we can use the survival-function inverse
219 # ``invgauss.isf`` whenever q > 0.5. Computing the smaller of (q, qc)
220 # accurately preserves machine precision in deep ITM/OTM regimes where
221 # the quantile lies far in the tail of the inverse Gaussian.
222 df = discount_factor
223 denom = df * min(forward, strike)
224 if is_call:
225 numer = df * forward - price
226 # qc = 1 - q; time value above intrinsic, normalized by D*K when ITM
227 qc = price / (df * forward) if strike >= forward else (price - df * (forward - strike)) / (df * strike)
228 else:
229 numer = df * strike - price
230 qc = (price - df * (strike - forward)) / (df * forward) if strike >= forward else price / (df * strike)
231 q = numer / denom
233 # Numerical guard: q must lie in (0, 1).
234 q = min(max(q, 1e-300), 1.0)
235 qc = min(max(qc, 1e-300), 1.0)
237 # scipy.stats.invgauss is parameterised by mu with shape lambda = 1,
238 # matching Schadner's F_IG(.; 2/|k|, 1). Use isf when q > 0.5 because the
239 # survival-function inversion is better-conditioned in the right tail.
240 mu = 2.0 / abs(k)
241 x = float(invgauss.isf(qc, mu)) if q > 0.5 else float(invgauss.ppf(q, mu))
242 return 2.0 / float(np.sqrt(expiry * x))