# 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

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):
# 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)
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("deltas FD:",deltaFD.round(3))
print("gammas FD:", gammaFD.round(3))


Outputs:

[HEAVISIDE]

deltas FD: [0.11  0.061 0.022]
gammas FD: [0.087 0.025 0.003]

[MAX]

deltas FD: [0.11  0.061 0.022]
gammas FD: [0.087 0.025 0.003]


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:

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()


1. 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.
2. In Python (0) results in 0. A one-element tuple is defined as (0,).