Solution: Computing and Interpreting Gamma Across Strikes and Expiries

Exercise: Computing and Interpreting Gamma Across Strikes and Expiries

Computation

import numpy as np
from scipy.stats import norm

def gamma_bs(S, K, T, r, sigma, q=0):
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    return np.exp(-q*T) * norm.pdf(d1) / (S*sigma*np.sqrt(T))

def delta_call(S, K, T, r, sigma, q=0):
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    return np.exp(-q*T) * norm.cdf(d1)

# Part 1: Gamma vs strike
S, T, r, sigma = 100, 0.5, 0.05, 0.2
for K in [80, 90, 100, 110, 120]:
    g = gamma_bs(S, K, T, r, sigma)
    print(f"K={K}: gamma={g:.5f}")
# K=80:  gamma=0.00394
# K=90:  gamma=0.01770
# K=100: gamma=0.02733
# K=110: gamma=0.02247
# K=120: gamma=0.01110

# Gamma maximised around K ≈ 100 (ATM). Note max is slightly ITM (K=100 here, but
# theoretically exact max is at K where d1 = 0, i.e. K = S*exp((r+sigma^2/2)*T) ≈ 103.5).

# Part 2: Gamma vs time to expiry
K = 100
for T in [1.0, 0.5, 0.1, 0.01]:
    g = gamma_bs(S, K, T, r, sigma)
    print(f"T={T}: gamma={g:.4f}")
# T=1.0:  gamma=0.0190
# T=0.5:  gamma=0.0273
# T=0.1:  gamma=0.0629
# T=0.01: gamma=0.1995   (huge!)

# Gamma grows as 1/sqrt(T-t) near ATM. At T=0.01 (about 2.5 trading days), gamma
# has blown up almost 10x its 1-year value.

# Part 3: FD check
T = 0.5
K = 100
h = 0.01
fd = (delta_call(S+h, K, T, r, sigma) - delta_call(S-h, K, T, r, sigma)) / (2*h)
cf = gamma_bs(S, K, T, r, sigma)
print(f"closed-form gamma: {cf:.6f}")
print(f"FD gamma:         {fd:.6f}")
# closed-form gamma: 0.027326
# FD gamma:         0.027326   (agrees perfectly)

# Part 4: Put gamma
# P = K exp(-rT) Phi(-d2) - S exp(-qT) Phi(-d1)
# dP/dS = -exp(-qT) Phi(-d1) + density terms that cancel
# d^2P/dS^2 = exp(-qT) phi(-d1) * (1/(S sigma sqrt T)) * sign_of_derivative...
# = same as call. Verify by FD.

def put_bs(S, K, T, r, sigma, q=0):
    d1 = (np.log(S/K) + (r - q + 0.5*sigma**2)*T) / (sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    return K*np.exp(-r*T)*norm.cdf(-d2) - S*np.exp(-q*T)*norm.cdf(-d1)

S = 100
put_gamma_fd = (put_bs(S+h, K, T, r, sigma) - 2*put_bs(S, K, T, r, sigma) + put_bs(S-h, K, T, r, sigma)) / h**2
print(f"put gamma (FD):  {put_gamma_fd:.6f}")
print(f"call gamma (CF): {cf:.6f}")
# put gamma (FD):  0.027326
# call gamma (CF): 0.027326   (equal, as expected)

Interpretation

Gamma peaks ATM. The $\phi(d_1)$ factor is maximised when $d_1 = 0$ , i.e. at $K = Se^{(r-q+\sigma^2/2)T}$ (slightly above spot for positive rates and dividends).
Gamma decays fast moving away from ATM. At $K = 80$ or $K = 120$ , gamma is an order of magnitude smaller than ATM. Deep-ITM and deep-OTM options behave nearly linearly in spot.
Gamma explodes near expiry for ATM options. As $T - t \to 0$ , ATM gamma $\to \infty$ — the payoff's kink at $K$ becomes infinitely sharp. This is the source of pin risk on the last trading day.
Put gamma = call gamma. Put-call parity forces the second derivatives to match; only the first derivatives (deltas) differ by $-e^{-q(T-t)}$ .

Takeaways

Gamma is a moneyness indicator for ATM-ness. Unlike delta (monotone, moneyness = ITM-ness), gamma peaks ATM and decays outward.
Near-expiry ATM options have explosive gamma. This is why options traders pay close attention to "the gamma book" in the final days before expiry; small stock moves can force large delta rebalances.
Call and put gamma coincide. You can use them interchangeably for vol trading — a long straddle (long call + long put) has double the gamma, not a mix of positive and negative.
Finite-difference gamma is reliable. For models without closed-form Greeks, numerical differentiation converges quickly and is the practical computation method (e.g. for local-vol or stochastic-vol models).