Gamma for a basket option in Python - Finite Differences vs. AAD Autograd library using Heaviside Approximation

Question

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 1

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

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,).