Python Implementation Guide¶
This section provides a practical guide to implementing the SABR model in Python, connecting the mathematical theory developed in earlier sections to working code. The implementation covers the Hagan implied volatility formula (with numerical stability safeguards), the calibration procedure, Monte Carlo simulation, and a complete usage workflow. The companion file python_implementation.py contains the full runnable code; this section explains the design decisions, highlights the numerical pitfalls, and demonstrates typical usage patterns.
Learning Objectives
By the end of this section, you will be able to:
- Implement the Hagan implied volatility formula with proper edge case handling
- Structure a SABR calibration module with ATM-first parameter determination
- Implement Monte Carlo simulation with exact volatility sampling
- Validate implementations against known analytical results
- Use the module for a complete pricing and calibration workflow
Module Design¶
Architecture¶
The implementation is organized around a single class SABRModel that encapsulates the model parameters and provides methods for pricing, calibration, and simulation:
```python class SABRModel: """SABR model for option pricing and calibration."""
def __init__(self, F: float, beta: float = 0.5):
self.F = F # Forward price
self.beta = beta # CEV exponent (fixed by convention)
# Core methods:
# _hagan_implied_vol(K, T, alpha, rho, nu) -> float
# call_price_hagan(K, T, r, alpha, rho, nu) -> float
# simulate_paths(T, alpha, rho, nu, ...) -> arrays
# price_monte_carlo(K, T, r, alpha, rho, nu, ...) -> (float, float)
```
A separate function calibrate_sabr(F, T, market_data, ...) handles calibration, since calibration is a workflow that uses the model rather than being an intrinsic property of it.
Design Decisions¶
Forward-based: The model takes the forward \(F\) as input (not the spot). This is natural for the SABR model, which is formulated in forward space.
Beta fixed at construction: The CEV exponent \(\beta\) is set once and not changed, reflecting the convention that \(\beta\) is fixed by the market.
Parameters passed to methods: The SABR parameters \((\alpha, \rho, \nu)\) are passed to each method rather than stored, because the same forward may be priced with different parameter sets during calibration.
Hagan Formula Implementation¶
The Core Function¶
The Hagan implied volatility formula is the most critical function to implement correctly. The key challenge is handling the edge cases that arise when \(K \approx F\) (ATM), \(\nu \approx 0\) (CEV limit), or \(|\rho| \approx 1\).
```python def _hagan_implied_vol(self, K, T, alpha, rho, nu): F, beta = self.F, self.beta
# ATM case: use simplified formula
if abs(F - K) < 1e-8 * F:
fk = F * F # F*K = F^2 at ATM
fk_beta = fk ** ((1 - beta) / 2)
correction = (
(1 - beta) ** 2 * alpha ** 2 / (24 * fk ** (1 - beta))
+ rho * beta * nu * alpha / (4 * fk_beta)
+ (2 - 3 * rho ** 2) * nu ** 2 / 24
)
return alpha / F ** (1 - beta) * (1 + correction * T)
# General case
fk = F * K
fk_beta = fk ** ((1 - beta) / 2)
log_fk = np.log(F / K)
# z and x(z) computation
z = (nu / alpha) * fk_beta * log_fk
# ...
```
Numerical Stability: The z/x(z) Ratio¶
The function \(x(z) = \ln((\sqrt{1 - 2\rho z + z^2} + z - \rho)/(1-\rho))\) and the ratio \(z/x(z)\) require careful handling:
When \(|z| < 10^{-7}\): Use the Taylor expansion:
python
if abs(z) < 1e-7:
zxz = 1.0 # z/x(z) -> 1 as z -> 0
else:
sqrt_term = np.sqrt(1 - 2 * rho * z + z ** 2)
x = np.log((sqrt_term + z - rho) / (1 - rho))
zxz = z / x if abs(x) > 1e-10 else 1.0
When \(|\rho|\) is near 1: The denominator \(1 - \rho\) in the logarithm approaches zero. Check for this:
python
if abs(1 - rho) < 1e-10:
# rho ≈ 1: use limiting formula
zxz = z / np.log(1 / (1 - z)) if z < 0.999 else 1.0
When the argument of the logarithm is non-positive: This can happen for extreme parameter combinations. Guard against it:
python
arg = (sqrt_term + z - rho) / (1 - rho)
if arg <= 0:
return alpha / F ** (1 - beta) # Fall back to backbone
Calibration Implementation¶
ATM-First Procedure¶
The calibration follows the two-step procedure described in the calibration section:
```python def calibrate_sabr(F, T, market_vols, beta=0.5): """ Calibrate SABR to market implied volatilities.
Args:
F: Forward price
T: Time to expiry
market_vols: dict {strike: implied_vol}
beta: CEV exponent (fixed)
Returns:
dict with calibrated (alpha, rho, nu)
"""
model = SABRModel(F, beta)
atm_vol = market_vols.get(F, None)
def alpha_from_atm(rho, nu):
"""Solve for alpha given ATM vol, rho, nu."""
# Newton's method on the ATM Hagan formula
alpha = atm_vol * F ** (1 - beta) # Initial guess
for _ in range(10):
vol = model._hagan_implied_vol(F, T, alpha, rho, nu)
dvol = ... # derivative w.r.t. alpha
alpha -= (vol - atm_vol) / dvol
return alpha
def objective(params):
rho, nu = params
alpha = alpha_from_atm(rho, nu)
error = 0
for K, vol_mkt in market_vols.items():
vol_model = model._hagan_implied_vol(K, T, alpha, rho, nu)
error += (vol_model - vol_mkt) ** 2
return error
# Optimize over (rho, nu)
from scipy.optimize import minimize
result = minimize(objective, [-0.3, 0.4],
bounds=[(-0.999, 0.999), (0.01, 2.0)])
rho_opt, nu_opt = result.x
alpha_opt = alpha_from_atm(rho_opt, nu_opt)
return {"alpha": alpha_opt, "rho": rho_opt, "nu": nu_opt}
```
Bounds and Constraints¶
The optimization enforces \(\rho \in (-0.999, 0.999)\) and \(\nu \in (0.01, 2.0)\) using bounds. These prevent the optimizer from exploring unphysical parameter regions where the Hagan formula is unreliable.
Monte Carlo Implementation¶
Exact Volatility, Euler Forward¶
The recommended simulation scheme uses exact sampling for the lognormal volatility and Euler discretization for the forward:
```python def simulate_paths(self, T, alpha, rho, nu, num_paths=10000, num_steps=100, seed=None): if seed is not None: rng = np.random.default_rng(seed) else: rng = np.random.default_rng()
dt = T / num_steps
F = np.full(num_paths, self.F)
sigma = np.full(num_paths, alpha)
for _ in range(num_steps):
# Independent normals
Z1 = rng.standard_normal(num_paths)
Z2 = rng.standard_normal(num_paths)
# Correlated increments
W1 = Z1
W2 = rho * Z1 + np.sqrt(1 - rho ** 2) * Z2
# Exact volatility step (lognormal)
sigma *= np.exp(-0.5 * nu ** 2 * dt + nu * np.sqrt(dt) * W2)
# Euler forward step with absorbing boundary
dF = sigma * np.maximum(F, 0) ** self.beta * np.sqrt(dt) * W1
F = np.maximum(F + dF, 0)
return F, sigma
```
Key Implementation Details¶
Vectorization: The entire simulation is vectorized over paths using NumPy arrays. No Python loops over paths.
Absorbing boundary: np.maximum(F, 0) in the diffusion coefficient and np.maximum(F + dF, 0) after the step ensure the forward stays non-negative.
Exact volatility: sigma *= np.exp(...) is the exact lognormal evolution, avoiding the negative-volatility issue.
Random number generation: Uses np.random.default_rng() (the modern NumPy generator) instead of the deprecated np.random.seed().
Validation¶
Against Known Results¶
The implementation should be validated against:
-
CEV limit (\(\nu = 0\)): The Hagan formula should reduce to \(\sigma_B = \alpha / F^{1-\beta}\) at the money.
-
Black--Scholes limit (\(\beta = 1\), \(\nu = 0\)): The Hagan formula should return \(\sigma_B = \alpha\) for all strikes.
-
Symmetry (\(\rho = 0\), \(\beta = 1\)): The smile should be symmetric around ATM.
-
Monte Carlo vs Hagan: For European options, the MC price should converge to the Hagan price as the number of paths and time steps increase.
```python def validate(): model = SABRModel(F=0.03, beta=0.5) alpha, rho, nu = 0.035, -0.3, 0.4
# Test 1: ATM Hagan vs MC
vol_hagan = model._hagan_implied_vol(0.03, 1.0, alpha, rho, nu)
price_hagan = model.call_price_hagan(0.03, 1.0, 0.0, alpha, rho, nu)
F_terminal, _ = model.simulate_paths(1.0, alpha, rho, nu,
num_paths=500000, num_steps=200)
payoffs = np.maximum(F_terminal - 0.03, 0)
price_mc = np.mean(payoffs)
print(f"Hagan price: {price_hagan:.6f}")
print(f"MC price: {price_mc:.6f} +/- {np.std(payoffs)/np.sqrt(len(payoffs)):.6f}")
```
Typical Validation Results¶
| Test | Expected | Actual | Pass |
|---|---|---|---|
| ATM vol (\(\nu=0\)) = \(\alpha/F^{0.5}\) | 20.21% | 20.21% | Yes |
| MC vs Hagan (ATM, 500K paths) | Agreement within 1 bps | 0.3 bps difference | Yes |
| Symmetric smile (\(\rho=0\), \(\beta=1\)) | \(\sigma(K) = \sigma(F^2/K)\) | Verified | Yes |
Usage Workflow¶
Complete Example¶
A typical workflow using the SABR module:
```python
1. Set up model¶
model = SABRModel(F=0.035, beta=0.5)
2. Market data¶
market_vols = { 0.015: 0.248, # -200 bps 0.025: 0.221, # -100 bps 0.035: 0.202, # ATM 0.045: 0.192, # +100 bps 0.055: 0.187, # +200 bps }
3. Calibrate¶
params = calibrate_sabr(F=0.035, T=5.0, market_vols=market_vols, beta=0.5) alpha, rho, nu = params["alpha"], params["rho"], params["nu"] print(f"Calibrated: alpha={alpha:.4f}, rho={rho:.3f}, nu={nu:.3f}")
4. Compute smile at fine strikes¶
strikes = np.linspace(0.01, 0.06, 50) vols = [model._hagan_implied_vol(K, 5.0, alpha, rho, nu) for K in strikes]
5. Price an exotic via MC¶
F_paths, sigma_paths = model.simulate_paths(5.0, alpha, rho, nu, num_paths=100000, num_steps=500)
... compute path-dependent payoff¶
```
Summary¶
The Python implementation of the SABR model centers on three components: the Hagan implied volatility formula (with edge case handling for ATM, small \(z\), and extreme \(\rho\)), the calibration procedure (ATM-first with Levenberg--Marquardt or Nelder--Mead for the smile), and Monte Carlo simulation (exact lognormal volatility with Euler forward stepping). The implementation is validated against CEV limits, Black--Scholes limits, symmetry checks, and Monte Carlo convergence. The companion file python_implementation.py provides the complete runnable code.
Further Reading¶
- Le Floc'h, F. (2014). Fast and accurate analytic basis point SABR. SSRN preprint.
- Hagan, P. et al. (2002). Managing smile risk. Wilmott Magazine, 1, 84--108.
- VanderPlas, J. (2016). Python Data Science Handbook. O'Reilly, Chapters 2--4 (NumPy and SciPy).
Exercises¶
Exercise 1. The Hagan formula involves the ratio \(z/x(z)\) where \(x(z) = \ln((\sqrt{1 - 2\rho z + z^2} + z - \rho)/(1 - \rho))\). Compute \(z/x(z)\) numerically for \(\rho = -0.3\) and \(z = 0.5\). Then compute it for \(z = 10^{-8}\) and verify that the Taylor expansion \(z/x(z) \approx 1\) is appropriate for small \(|z|\).
Solution to Exercise 1
For \(\rho = -0.3\) and \(z = 0.5\), we compute step by step:
This ratio exceeds 1, which means the implied volatility is boosted above the backbone level --- the smile effect.
For \(z = 10^{-8}\), we can verify the Taylor expansion. The function \(x(z)\) near \(z = 0\) satisfies \(x(0) = 0\) and \(x'(0) = 1\) (by L'Hopital's rule or direct expansion). The Taylor expansion gives
so that
For \(z = 10^{-8}\), the correction \(\rho z/2 \approx 1.5 \times 10^{-9}\) is well below machine epsilon (\(\approx 10^{-16}\) for double precision), so \(z/x(z) \approx 1\) to full floating-point accuracy. Direct computation of \(x(z)\) for such small \(z\) would involve \(\ln(1 + \epsilon)\) where \(\epsilon \approx 10^{-8}\), which suffers from catastrophic cancellation. The Taylor approximation \(z/x(z) = 1\) is therefore both accurate and numerically necessary for \(|z| < 10^{-7}\).
Exercise 2. The ATM-first calibration determines \(\alpha\) by solving \(\sigma_B^{\text{Hagan}}(\alpha, \rho, \nu) = \sigma_B^{\text{ATM}}\) using Newton's method. Write the Newton iteration formula explicitly:
For \(F = 0.03\), \(\beta = 0.5\), \(\sigma_B^{\text{ATM}} = 0.22\), \(\rho = -0.25\), and \(\nu = 0.35\), compute the initial guess \(\alpha_0 = \sigma_B^{\text{ATM}} \cdot F^{1-\beta}\).
Solution to Exercise 2
The ATM Hagan formula (simplified for \(K = F\)) is
Define \(f(\alpha) = \sigma_B^{\mathrm{ATM}}(\alpha) - \sigma_B^{\mathrm{market}}\). Newton's iteration is
The derivative is
Note that \(\partial/\partial\alpha\) of \(\alpha^2/(24 F^{2(1-\beta)})\) gives \(2\alpha/(24 F^{2(1-\beta)})\), so the \(\alpha^2\) term in the bracket contributes with a factor of 3 instead of 1 when differentiating the product \(\alpha \times [1 + \cdots \alpha^2 \cdots]\).
For the initial guess with \(F = 0.03\), \(\beta = 0.5\), and \(\sigma_B^{\mathrm{ATM}} = 0.22\):
This initial guess corresponds to setting the bracket to 1 (ignoring the \(O(T)\) correction). It is typically within 5--10% of the true \(\alpha\), so Newton's method converges in 3--5 iterations.
Exercise 3. In the Monte Carlo implementation, the volatility is simulated exactly: \(\sigma_{n+1} = \sigma_n \exp(-\frac{1}{2}\nu^2\Delta t + \nu\sqrt{\Delta t}\,W_2)\). Show that this is the exact solution to \(d\sigma_t = \nu\sigma_t\,dW_t^{(2)}\) over the interval \([t_n, t_{n+1}]\). Why does this eliminate the need for a positivity floor on \(\sigma\)?
Solution to Exercise 3
The SDE for the volatility is \(d\sigma_t = \nu\sigma_t\,dW_t^{(2)}\), which is a geometric Brownian motion with zero drift. By Ito's lemma, applying \(f(\sigma) = \ln\sigma\):
Integrating from \(t_n\) to \(t_{n+1} = t_n + \Delta t\):
Since \(W_{t_{n+1}}^{(2)} - W_{t_n}^{(2)} \sim \mathcal{N}(0, \Delta t) = \sqrt{\Delta t}\,Z\) where \(Z \sim \mathcal{N}(0,1)\):
This is exactly the update used in the implementation: sigma *= np.exp(-0.5 * nu**2 * dt + nu * np.sqrt(dt) * W2).
This formula is the exact solution of the SDE over the interval \([t_n, t_{n+1}]\), not an Euler approximation. As a consequence:
- There is no discretization error in the volatility path.
- Since \(\sigma_{t_{n+1}} = \sigma_{t_n} \exp(\cdot)\) and the exponential function is always positive, \(\sigma_{t_{n+1}} > 0\) automatically whenever \(\sigma_{t_n} > 0\). No positivity floor is needed.
By contrast, an Euler scheme \(\sigma_{t_{n+1}} = \sigma_{t_n} + \nu\sigma_{t_n}\sqrt{\Delta t}\,Z\) could in principle produce negative values for large negative \(Z\) (though the probability is exponentially small for reasonable \(\Delta t\)). The exact scheme eliminates this possibility entirely.
Exercise 4. The implementation uses np.maximum(F + dF, 0) to enforce the absorbing boundary. Discuss the bias this introduces: does it overestimate or underestimate the absorption probability? How would you implement the more accurate chi-squared-based absorption test described in the Monte Carlo section?
Solution to Exercise 4
The np.maximum(F + dF, 0) absorbing boundary introduces a downward bias in the forward and therefore overestimates the absorption probability. Here is why:
The true SABR dynamics have an absorbing boundary at \(F = 0\) (for \(\beta < 1\)). The exact absorption occurs when a continuous path first hits zero. In the Euler scheme, the forward is updated discretely:
The np.maximum(..., 0) clamp sets \(F_{n+1} = 0\) whenever the discrete step overshoots to a negative value. However, the continuous path might have:
- Crossed zero and returned positive between \(t_n\) and \(t_{n+1}\): The Euler scheme absorbs the path, but the continuous path would survive. This overestimates absorption.
- Hit zero between \(t_n\) and \(t_{n+1}\) but the Euler step remains positive: The scheme misses the absorption entirely.
The net effect for typical parameters (small \(\Delta t\), moderate \(\beta\)) is that the simple clamp overestimates absorption, leading to option prices that are slightly too low (since absorbed paths contribute zero payoff).
A more accurate approach uses the chi-squared bridge technique. For the CEV process with \(d\tilde{F} = \sigma \tilde{F}^\beta\,dW\), the transition density is related to the noncentral chi-squared distribution. Conditional on \(F_n\) and \(F_{n+1}\) (both positive), the probability that the continuous path hit zero in between can be computed exactly using the formula
which has a known closed-form expression involving noncentral chi-squared CDFs (see Chen and Glasserman, 2008). Implementation steps:
- Simulate \(F_{n+1}\) from the noncentral chi-squared transition density (not Euler).
- Draw a uniform \(U \sim \text{Uniform}(0,1)\).
- If \(U <\) the bridge absorption probability, set \(F_{n+1} = 0\) (absorbed).
- Otherwise, keep \(F_{n+1}\) as simulated.
This approach is exact for the conditional absorption probability and eliminates the discretization bias.
Exercise 5. Design a validation test for the SABR implementation that checks the Black-Scholes limit (\(\beta = 1\), \(\nu = 0\)). In this limit, the SABR model reduces to geometric Brownian motion with constant volatility \(\alpha\). Write a test that compares the Hagan formula output, the Monte Carlo price, and the exact Black-Scholes price for an ATM call with \(F = 0.05\), \(T = 1\), and \(\alpha = 0.20\). What tolerances should you use for each comparison?
Solution to Exercise 5
In the Black--Scholes limit (\(\beta = 1\), \(\nu = 0\)), the SABR dynamics reduce to
which is geometric Brownian motion with constant volatility \(\alpha\). The Hagan formula should give \(\sigma_B(K) = \alpha\) for all strikes, and option prices should match the exact Black formula.
Test design:
Setup: \(F = 0.05\), \(T = 1\), \(\alpha = 0.20\), \(\beta = 1\), \(\rho = 0\), \(\nu = 0\) (or \(\nu = 10^{-6}\) to avoid division-by-zero issues), \(K = F = 0.05\) (ATM call).
```python import numpy as np from scipy.stats import norm
def black_call(F, K, T, sigma): """Exact Black formula for a call on a forward.""" d1 = (np.log(F / K) + 0.5 * sigma**2 * T) / (sigma * np.sqrt(T)) d2 = d1 - sigma * np.sqrt(T) return F * norm.cdf(d1) - K * norm.cdf(d2)
model = SABRModel(F=0.05, beta=1.0) alpha, rho, nu = 0.20, 0.0, 1e-6 K, T = 0.05, 1.0
Test 1: Hagan implied vol == alpha¶
vol_hagan = model._hagan_implied_vol(K, T, alpha, rho, nu) assert abs(vol_hagan - alpha) < 1e-6, f"Hagan vol {vol_hagan} != alpha {alpha}"
Test 2: Hagan price == Black price¶
price_hagan = model.call_price_hagan(K, T, 0.0, alpha, rho, nu) price_black = black_call(0.05, 0.05, 1.0, 0.20) assert abs(price_hagan - price_black) < 1e-8, "Hagan price != Black price"
Test 3: MC price ≈ Black price¶
F_terminal, _ = model.simulate_paths(T, alpha, rho, nu, num_paths=500000, num_steps=200) payoffs = np.maximum(F_terminal - K, 0) price_mc = np.mean(payoffs) se = np.std(payoffs) / np.sqrt(len(payoffs)) assert abs(price_mc - price_black) < 3 * se, "MC price outside 3-sigma" ```
Tolerances:
- Hagan vs exact \(\alpha\): Tolerance \(10^{-6}\) (the error arises only from the \(\nu = 10^{-6}\) approximation and floating-point arithmetic).
- Hagan price vs Black price: Tolerance \(10^{-8}\) (both are deterministic calculations using the same formula in this limit).
- MC price vs Black price: Use a statistical tolerance of 3 standard errors. For 500,000 paths with \(F = 0.05\), \(\sigma = 0.20\), \(T = 1\), the ATM call price is approximately \(F \cdot \sigma\sqrt{T}/(2\sqrt{2\pi}) \approx 0.004\), and the standard error is roughly \(0.004/\sqrt{500000} \approx 5.7 \times 10^{-6}\), so the tolerance is about \(1.7 \times 10^{-5}\).
Exercise 6. The calibration function uses bounds \(\rho \in (-0.999, 0.999)\) and \(\nu \in (0.01, 2.0)\). Explain why each bound is necessary. What numerical problems arise if \(\rho\) is allowed to reach exactly \(\pm 1\)? What happens to the Hagan formula if \(\nu = 0\), and why is a small positive lower bound used instead?
Solution to Exercise 6
Bounds on \(\rho \in (-0.999, 0.999)\):
The correlation parameter \(\rho\) appears in the Hagan formula in several places that become singular as \(|\rho| \to 1\):
- The function \(x(z) = \ln((\sqrt{1-2\rho z + z^2} + z - \rho)/(1-\rho))\) has the denominator \((1-\rho)\). When \(\rho = 1\), this denominator is zero and the logarithm diverges for any \(z > 0\). The ratio \(z/x(z)\) requires a separate limiting formula.
- The geodesic distance on the SABR manifold involves the factor \(1/(1-\rho^2)\), which diverges as \(|\rho| \to 1\). Geometrically, at \(|\rho| = 1\) the two Brownian motions are perfectly (anti-)correlated, and the two-dimensional diffusion degenerates to a one-dimensional process on a curve within the \((F, \sigma)\) plane.
- The Cholesky decomposition used in Monte Carlo (\(W_2 = \rho W_1 + \sqrt{1-\rho^2}\,Z_2\)) has \(\sqrt{1-\rho^2} = 0\) at \(|\rho| = 1\), so the second Brownian motion loses all independent randomness.
Using bounds \(\pm 0.999\) keeps the optimizer away from these singularities while still allowing near-perfect correlation when the data demands it.
Bounds on \(\nu \in (0.01, 2.0)\):
Lower bound: When \(\nu = 0\), the SABR model reduces to the CEV model with constant \(\sigma = \alpha\). The Hagan formula is valid in this limit, but the \(z/x(z)\) factor involves the ratio \(0/0\) (since \(z = (\nu/\alpha)(FK)^{(1-\beta)/2}\ln(F/K)\) is proportional to \(\nu\), and \(x(z) \to 0\) as \(z \to 0\)). While the limit \(z/x(z) \to 1\) is well-defined analytically, the numerical computation of \(x(z) = \ln(\cdots)\) for \(z \approx 0\) suffers from catastrophic cancellation. Using \(\nu \geq 0.01\) avoids this numerical issue. If the true best-fit \(\nu\) is smaller than 0.01, the model is effectively CEV and should be calibrated as such.
Upper bound: The Hagan formula is a short-time asymptotic expansion with error \(O(\nu^2 T)\). For \(\nu > 2\), the expansion parameter \(\nu^2 T\) exceeds 4 even for \(T = 1\), making the approximation unreliable. The upper bound prevents the optimizer from exploring parameter regions where the Hagan formula produces meaningless results. In practice, \(\nu > 1\) is already unusual for most markets.