# Option pricing Greeks in Python - incorrect Gamma with MC option pricing (Black) using AAD autograd / JAX libraries - but works with closed form?

I am attempting to use AAD (Adjoint Algorithmic Differentiation) with a simple Black MC pricer, and found that the Gamma is incorrect. The output was compared to Black analytical Greeks, as well as Finite Difference Greeks, and there is a discrepancy. What confuses me is the Black analytical formula produces the correct Gamma with identical AutoGrad (or JAX) code. I don't know if anyone can spot my error or has more experience with 2nd order derivatives and AAD that can provide advice on how to fix the issue.

Here's the simple MC Black option pricer for a call option, along with the AutoGrad code to replicate the issue. MC_call_price is a simple MC, BS_call_price is the analytical one:

# use either jax or autograd functions, same format
import autograd.numpy as np # import jax.numpy as np
from autograd.scipy.stats import norm # from jax.scipy.stats import norm

def MC_call_price(F, vol, K, T, IR, steps, trials):
np.random.seed(123)
dt = T/steps
paths = np.log(F) +  np.cumsum(((IR - vol**2/2)*dt + vol*np.sqrt(dt) * \
np.random.normal(size=(int(steps),int(trials)))), axis=0)
payoff = np.mean(np.maximum(np.exp(paths)[-1]-K, 0)) * np.exp(-IR*T)
return payoff

def AAD_Greeks(F, vol, K, T, IR, steps, trials):
gradient_func = grad(MC_call_price, (0)) # tuple specifies inputs we want to differentiate
# solve for Greeks
delta = gradient_func(F, vol, K, T, IR, steps, trials)
gamma = gradient_func2(F, vol, K, T, IR, steps, trials)
return delta, gamma

def BS_call_price(F, vol, K, T, IR):
b = np.exp(-IR*T).astype('float')
x1 = np.log(F/(b*K))
x1 += (.5*(vol**2)*T)
x1 = x1/(vol*(T**.5))
z1 = norm.cdf(x1)
z1 = z1*F
x2 = np.log(F/(b*K)) - .5*(vol**2)*T
x2 = x2/(vol*(T**.5))
z2 = norm.cdf(x2)
z2 = b*K*z2
return z1 - z2

def AAD_Greeks2(F, vol, K, T, IR):
gradient_func = grad(BS_call_price, (0)) # tuple specifies inputs we want to differentiate
# solve for Greeks
delta = gradient_func(F, vol, K, T, IR)
gamma = gradient_func2(F, vol, K, T, IR)
return delta, gamma

F = 100.0
vol = 0.30
K = 90
T = 0.5
IR = 0.03
steps = 1
trials = 1000000
bump = 0.005

optionprice = MC_call_price(F, vol, K, T, IR, steps, trials)
delta, gamma = AAD_Greeks(F, vol, K, T, IR, steps, trials)
print("\nOption price", optionprice,"\n", "\nAAD: Delta", delta, "Gamma", gamma)
up = MC_call_price(F+bump, vol, K, T, IR, steps, trials)
down = MC_call_price(F-bump, vol, K, T, IR, steps, trials)
FDdelta = (up-down)/(2*bump)
FDgamma = (up - 2*optionprice + down)/(bump**2)
print("FD: Delta", FDdelta, "Gamma", FDgamma)
BSoptionprice = BS_call_price(F, vol, K, T, IR, )
BSdelta, BSgamma = AAD_Greeks2(F, vol, K, T, IR)
print("\n\nBS Option price", BSoptionprice,"\n", "\nAAD: Delta", BSdelta, "Gamma", BSgamma)


And the output showing MC doesn't produce the correct Gamma, but the BS analytical formula does with AAD, which is close to the FD method:

Option price 14.892979946909723