Pricing models, implied volatility and the surface

Notation: S spot, K strike, T years to expiry (calendar days / 365), r continuous risk-free rate, q continuous dividend yield, sigma volatility, N the standard normal CDF, phi its density.

Black-Scholes-Merton (models/bsm.py)

d1 = [ln(S/K) + (r - q + sigma^2/2) T] / (sigma sqrt T),   d2 = d1 - sigma sqrt T
call = S e^{-qT} N(d1) - K e^{-rT} N(d2)
put  = K e^{-rT} N(-d2) - S e^{-qT} N(-d1)

Greeks (call; put signs in brackets), reported in trader units:

Greek Formula Unit reported
delta e^{-qT} N(d1) [-e^{-qT} N(-d1)] per 1 unit of S
gamma e^{-qT} phi(d1) / (S sigma sqrt T) per 1 unit of S
vega S e^{-qT} phi(d1) sqrt T per 1 vol point (/100)
theta -S e^{-qT} phi(d1) sigma / (2 sqrt T) -+ r K e^{-rT} N(+-d2) +- q S e^{-qT} N(+-d1) per calendar day (/365)
rho K T e^{-rT} N(d2) [-K T e^{-rT} N(-d2)] per 1% rate (/100)
vanna -e^{-qT} phi(d1) d2 / sigma = d(delta)/d(sigma) per vol point
charm -d(delta)/dT, closed form with the (2(r-q)T - d2 sigma sqrt T) / (2T sigma sqrt T) term per calendar day
vomma vega d1 d2 / sigma per vol point squared

Rho holds q fixed. Every Greek is checked against central finite differences of the price in the test suite, and the headline numbers against Hull's examples (S=42, K=40, r=10%, sigma=20%, T=0.5 gives 4.76 / 0.81; the index option S=930, K=900, r=8%, q=3%, sigma=20%, T=2/12 gives 51.83; S=49, K=50 Greeks example gives delta 0.522, gamma 0.066, vega 12.1, theta -4.31/yr, rho 8.91).

Edge cases: T <= 0 returns intrinsic value and a step delta; sigma <= 0 returns the discounted forward intrinsic; S = 0 is clipped to avoid log(0) so spot grids may include zero.

Black-76 (models/black76.py)

For an option on a forward or future F:

call = e^{-rT} [F N(d1) - K N(d2)],   d1 = [ln(F/K) + sigma^2 T / 2] / (sigma sqrt T)

This is BSM with S = F and q = r, so the module delegates to it. The only change is rho: the forward does not move when r moves, so rho = -T x price. Hull's futures example (F=20, K=20, r=9%, sigma=25%, T=4/12) gives 1.1166 for both call and put.

CRR binomial tree (models/binomial.py)

dt = T/n,  u = exp(sigma sqrt dt),  d = 1/u,  p = (e^{(r-q) dt} - d) / (u - d)
V_node = max( e^{-r dt} [p V_up + (1-p) V_down],  intrinsic )   # American only

The tree is vectorised over a whole vector of spot values (a 2-D array of node values), so payoff curves cost one tree walk, not one per point.

Convergence. A CRR price oscillates with n depending on where the strike falls between terminal nodes. The default is the average of the n and n+1 step trees, which cancels most of that oscillation. Worked example, European put S=K=100, r=5%, sigma=20%, T=1 (BSM 5.5735):

steps single tree averaged
100 oscillates around BSM closer
800 converged within 2e-3 across strikes (tested)

Known American value: S=50, K=50, r=10%, sigma=40%, T=5/12 gives 4.284 (Hull: 4.28) against 4.076 European. An American call without dividends equals the European call, and with a high yield it exceeds it (both asserted in tests).

Greeks are bump-and-reprice: delta and gamma with a 2% spot bump (smaller bumps fall inside the node spacing and make gamma noisy), vega with +-0.5 vol points, theta with a one-day step, rho with +-1bp. Vanna, charm and vomma are not computed for trees.

The tree requires 0 < p < 1, which fails for tiny sigma relative to r sqrt(dt); the solver raises a clear error and the American IV solver starts its bracket above that limit.

Monte Carlo (models/montecarlo.py)

Terminal price S_T = S exp((r - q - sigma^2/2) T + sigma sqrt T Z). With antithetic variates each sample is the pair average 0.5 (f(Z) + f(-Z)); the standard error is computed over pairs, so it is honest about the variance reduction. Result: price, standard error, 95% CI. It exists to validate the closed forms and the tree: tests require the BSM price inside the CI.

Implied volatility (models/iv.py)

  1. Reject prices outside the no-arbitrage band: below max(discounted forward intrinsic, 0) or above S e^{-qT} (call) / K e^{-rT} (put).
  2. Deep ITM transform. Put-call parity turns an ITM quote into the OTM option of the other type (call price minus discounted forward intrinsic is the put price). This keeps the time value from drowning in float error; a price that equals intrinsic to machine precision is reported as "no time value".
  3. Safeguarded Newton-Raphson on f(sigma) = BSM(sigma) - price with the analytic vega. Start: Brenner-Subrahmanyam sigma0 = sqrt(2 pi / T) price / S, clamped to [0.01, 3]. A bracket [lo, hi] is maintained from the sign of f; a Newton step that leaves the bracket (or vega near zero, typical for deep OTM) is replaced by bisection. The loop ends on a relative price tolerance, with brentq as the last resort.
  4. American options: Brent root-finding on the binomial price (no closed-form vega).

Example: price 4.7594 for the Hull call returns 19.9997% (the input 20% less the rounding of the price). Tests recover the input vol across moneyness 0.6-1.6, expiries from 7 days to 3 years and vols from 8% to 90%, skipping only cases where the time value is below float resolution.

Volatility surface (surface.py)

Smile points. For each expiry and strike the IV comes from the out-of-the-money option (call at or above the forward, put below), because it has the most time value and the tightest relative spread; the other side is used only if that one cannot be inverted. Mid price when bid and ask are both live, else last.

Raw SVI. Total variance w = sigma_imp^2 T as a function of forward log-moneyness k = ln(K/F), F = S e^{(r-q)T}:

w(k) = a + b ( rho (k - m) + sqrt((k - m)^2 + sigma^2) )

Fitted per expiry by bounded least squares in total variance (b > 0, |rho| < 1, sigma > 0), with 12 starting points (3 rho x 2 m x 2 sigma) and the best cost kept, because SVI has flat directions that trap a single start. Noise-free data recovers all five parameters to 1e-4 relative.

Interpolation. At a fixed k, total variance is interpolated linearly in T between neighbouring fitted expiries. Linear interpolation of total variance between two calendar-arbitrage-free slices stays free of calendar arbitrage. Before the first expiry the vol is held flat in k (w scales with T); beyond the last it is likewise extrapolated flat. Strategy pricing queries the surface by strike with k computed from the as-of forward, i.e. sticky strike when spot moves in scenarios.

Arbitrage checks.

Worked example: an SVI slice with a=0, b=2, rho=0.95, m=0, sigma=0.01 (a steep right wing) gives g(k) < 0 at 119 of 241 grid points and is flagged; the example slices give none.

References