Basket Option Gamma: AAD, Finite Differences, and Payoff Smoothing
Article Quant Q&A · Author: Matt
Summary
The document compares automatic differentiation (AAD) with finite differences (FD) for estimating gamma on a simulated basket option. It uses a smoothed exercise indicator around the strike so derivatives can be computed through the payoff, and reports that the sample AAD deltas and gammas are close to the FD estimates. It also compares the smoothed-payoff valuation with a maximum-based payoff calculation, which gives similar reported prices and FD Greeks in this example.
The response explains that elementwise gradient calls can sum derivative contributions across inputs, so applying them twice may not return each asset's diagonal second derivative. A full Hessian followed by extracting its diagonal is suggested for asset-specific gamma. The comparison is specific to the chosen simulation, smoothing, bump size, and random draws; closeness in this example does not establish general accuracy. The discussion also cautions about Python tuple syntax and integer division when implementing derivative calculations.
Key ideas
- Smoothing a discontinuous option payoff can make pathwise automatic differentiation feasible for Greek estimates.
- In the reported simulation, AAD and finite-difference deltas and gammas are close, but this is evidence from one setup only.
- Nested elementwise gradients may sum second-derivative contributions instead of returning the desired diagonal gamma.
- A Hessian can be used to extract each underlying's own second derivative directly.
- Greek comparisons depend on implementation details such as smoothing, bump size, and simulated paths.
Tags
Full text
# Gamma for a basket option in Python - Finite Differences vs. AAD Autograd library using Heaviside Approximation
# Gamma for a basket option in Python - Finite Differences vs. AAD Autograd library using Heaviside Approximation
I have been trying to use the Heaviside Approximation for a simple basket option so that I can solve for Gammas with AAD (Adjoint Automatic Differentiation). This routine smooths the payoff function vs. a max function, so the payoff is continuous rather than discrete - it's a smooth version of the step function. I believe I have properly called the Heaviside Approximation function, but could have made a small error. You can see the line in the code: `payoff = HeavisideApprox(payoff, np.sum(F0*weights*vols*np.sqrt(T)* 0.05))* payoff`
My questions are 1) have I applied Heaviside correctly? 2) Is the difference in Gammas between AAD and FD acceptable? 3) I had to multiply the AAD Gammas by the weights to match FD, and I'm not sure if any other Greeks need that done as well (2nd order maybe). Any suggestions are appreciated! Code follows:
```
import autograd.numpy as np
from autograd import elementwise_grad as egrad # egrad handles multiple sets of inputs
def Dirac(x):
return np.sin((1+np.sign(x))/2*np.pi)
def HeavisideApprox(x, epsilon):
p =0.5*(np.sign(x+epsilon) - np.sign(x - epsilon)) - 1/2* Dirac(x+epsilon) + 1/2* Dirac(x - epsilon)
p = p* (0.5+ 0.5* np.sin(x/2/epsilon* np.pi))
p = p+ (1+ np.sign(x - epsilon))/2 - 0.5* Dirac(x - epsilon)
return p
def MC_basket(F0, vols, K, T, IR, corr, weights, trials, steps, maximum=False):
# simple Euro option if steps=1, or (forward) Asian with more steps
np.random.seed(123)
eigvals, eigvects = np.linalg.eigh(corr)
eigvals[eigvals < 0.00001] = 0.00001
V = np.real(eigvects.dot(np.sqrt(np.diag(eigvals))))
randnums = np.random.randn(steps, F0.shape[0], trials) #3D Asian
randnumscorr = V.dot(randnums).T
paths = F0 *np.exp((-0.5*vols*vols)*T + randnumscorr*vols*np.sqrt(T))
paths = np.mean(paths, axis=1) # average aross steps #3D Asian
if maximum == False:
payoff = np.sum(paths*weights, axis=1)-K
payoff = HeavisideApprox(payoff, np.sum(F0*weights*vols*np.sqrt(T)* 0.05))* payoff
return np.mean(payoff)* np.exp(-IR*T)
else:
# normal max function, does not work with AAD Gamma
payoff = np.mean(np.maximum(np.sum(paths*weights, axis=1)-K,0)) * np.exp(-IR*T)
return payoff
def Greeks(F0, vols, K, T, IR, corr, weights, trials, steps=None):
# setup gradient functions to evaluate with AAD
gradient_func = egrad(MC_basket, (0, 1, 3, 4)) # tuple specifies inputs we want to differentiate
gradient_func2 = egrad(egrad(MC_basket, (0))) # second derivative egrad(egrad())
# solve for Greeks
delta, vega, theta, rho = gradient_func(F0, vols, K, T, IR, corr, weights, trials, steps)
gamma = weights * gradient_func2(F0, vols, K, T, IR, corr, weights, trials, steps)
# standardize Greeks for normal output
theta /= -365
vega /= 100
rho /= 100
return delta, gamma, vega, theta, rho
corr = np.array([[1,0.9,0.8],[0.9,1,0.94],[0.8,0.9,1]])
F0 = np.array([12.0,10.0,8.0])
weights = np.array([0.6, 0.3, 0.1])
vols = np.array([0.1,0.2,0.3])
K = 12
T = 0.5
IR = 0.03
steps = 1
trials = 100000
MCb = MC_basket(F0, vols, K, T, IR, corr, weights, trials, steps)
delta, gamma, vega, theta, rho = Greeks(F0, vols, K, T, IR, corr, weights, trials, steps)
print("\n[HEAVISIDE]\n\nMC basket price:", round(MCb,3), "\ndeltas AAD:", np.round(delta,3), "\ngammas AAD:", np.round(gamma,3))#, "\nvegas:", np.round(vega,3), "\ntheta:", np.round(theta,3), "\nrho:", np.round(rho,3))
bump = 0.005
# calculate FD delta and gamma using Heaviside for comparison with AAD
gammaFD = np.zeros((len(F0),))
deltaFD = np.zeros((len(F0),))
for i in range(0,F0.shape[0]):
MCb_up = np.zeros((len(F0),))
MCb_down = np.zeros((len(F0),))
MCb_up[i] = bump
MCb_down[i] = -bump
# percentage bumps
MCb_up[i] = MC_basket((MCb_up+1)*F0, vols, K, T, IR, corr, weights, trials, steps)
MCb_down[i] = MC_basket((MCb_down+1)*F0, vols, K, T, IR, corr, weights, trials, steps)
gammaFD[i] = (MCb_up[i] - 2*MCb + MCb_down[i]) / (bump**2*F0[i]**2)
deltaFD[i] = (MCb_up[i] - MCb_down[i]) / (2*bump*F0[i])
print("\ndeltas FD:",deltaFD.round(3))
print("gammas FD:", gammaFD.round(3))
# calculate FD delta and gamma using the max payoff function
MCb = MC_basket(F0, vols, K, T, IR, corr, weights, trials, steps, maximum=True)
gammaFD = np.zeros((len(F0),))
deltaFD = np.zeros((len(F0),))
for i in range(0,F0.shape[0]):
MCb_up = np.zeros((len(F0),))
MCb_down = np.zeros((len(F0),))
MCb_up[i] = bump
MCb_down[i] = -bump
# percentage bumps
MCb_up[i] = MC_basket((MCb_up+1)*F0, vols, K, T, IR, corr, weights, trials, steps, maximum=True)
MCb_down[i] = MC_basket((MCb_down+1)*F0, vols, K, T, IR, corr, weights, trials, steps, maximum=True)
gammaFD[i] = (MCb_up[i] - 2*MCb + MCb_down[i]) / (bump**2*F0[i]**2)
deltaFD[i] = (MCb_up[i] - MCb_down[i]) / (2*bump*F0[i])
print("\n[MAX]\n\nMC basket price:",MCb.round(3))
print("deltas FD:",deltaFD.round(3))
print("gammas FD:", gammaFD.round(3))
```
Outputs:
```
[HEAVISIDE]
MC basket price: 0.109
deltas AAD: [0.11 0.061 0.022]
gammas AAD: [0.09 0.024 0.003]
deltas FD: [0.11 0.061 0.022]
gammas FD: [0.087 0.025 0.003]
[MAX]
MC basket price: 0.109
deltas FD: [0.11 0.061 0.022]
gammas FD: [0.087 0.025 0.003]
```
## Answer by kwinto (score 5, accepted)
https://quant.stackexchange.com/a/75079
1) have I applied Heaviside correctly?
I'd say yes, as the result matches with the final difference calculation. Although, I'm puzzled why you use so overcomplicated smoothing for Heaviside function. In principle, any smooth approximation should work. I get the same result with
```
def HeavisideApprox(x, epsilon):
return 0.5 + 0.5 * np.tanh(100*x)
```
For more possible choices, see for example
https://github.com/norse/norse/blob/main/norse/torch/functional/threshold.py
2) Is the difference in Gammas between AAD and FD acceptable?
This very much depends on the use case. In general, AAD numbers should be more accurate than FD on top of Monte-Carlo.
3) I had to multiply the AAD Gammas by the weights to match FD, and I'm not sure if any other Greeks need that done as well (2nd order maybe).
It looks like `egrad(egrad(MC_basket, (0)))` returns a sum of second order derivatives:
```
Help on function elementwise_grad in module autograd.wrap_util:
elementwise_grad(fun, argnum=0, *nary_op_args, **nary_op_kwargs)
Returns a function that computes the sum of each column of the Jacobian of
`fun`, in one pass. If the Jacobian is diagonal, then this is the diagonal
of the Jacobian.
```
Hence, `gamma[0]` contains
$$ \frac{d^2 V}{d F_1^2} + \frac{d^2 V}{d F_1 d F_2} + \frac{d^2 V}{d F_1 d F_3} $$
After multiplying by `weights[0]` it contains what you need: $\frac{d^2 V}{d F_1^2}$, since deltas are proportional to weights:
$$ \frac{1}{w_1} \frac{d V}{d F_1} = \frac{1}{w_2} \frac{d V}{d F_2} = \frac{1}{w_3} \frac{d V}{d F_3} $$
To avoid explicitly multiplying by `weights`, you may consider to use
```
gamma = hessian(MC_basket, 0)(F0, vols, K, T, IR, corr, weights, trials, steps).diagonal()
```
Extra Comments
- Be careful with `1/2`, this results to `0` in Python and most other languages. Instead, prefer `1./2` or `0.5` to stay on the safe side.
- In Python `(0)` results in `0`. A one-element tuple is defined as `(0,)`.Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.