Option Gamma via MC using JAX

[G10VA]

Hi all, I'm new about Jax and I'm trying to use it for a finance project: compute option greeks in the MC setting.
I'm able to recover the first order Greeks but I'm getting wrong result from the higher order ones. Here my code:

import jax
import jax.scipy
from jax.scipy import special as Sfn
from jax import numpy

S = 100.0
K = 100.0
r = 0.01
sigma = 0.05
T = 1.0

from jax import random
key = random.PRNGKey(1701)

def BlackMC(S, K, r, sigma, T):
    Z = jax.random.normal(key, shape=(1000000,))
    W_T = numpy.sqrt(T)*Z
    prices = S * numpy.exp((r - 0.5 * numpy.square(sigma)) * T + sigma * W_T)
    payoffs = numpy.maximum(prices - K,0) #torch.max(prices - K, torch.zeros(1000000))
    return numpy.mean(payoffs) * numpy.exp(-r * T)

value = BlackMC(S, K, r, sigma, T)

DeltaMC, RhoMC, VegaMC, mThetaMC = [jax.grad(BlackMC, argnums=x) for x in [0, 2, 3, 4] ]


"""
Expected values

Delta: 0.5884349942207336
Rho: 56.32249450683594
Vega: 38.89812469482422
Theta: 1.5356781482696533
"""


print("Value: ", value)
print("Delta: ",DeltaMC(S, K, r, sigma, T))
print("Rho: ",RhoMC(S, K, r, sigma, T))
print("Vega: ",VegaMC(S, K, r, sigma, T))
print("mTheta: ",mThetaMC(S, K, r, sigma, T))

"""
Expected value

BS Gamma: 0.0777941569685936
"""

GammaMC = jax.grad(DeltaMC, argnums=0) 

print("Gamma: ", GammaMC(S, K, r, sigma, T))

I get 0 as answer for Gamma. Can someone provide help?
As said I'm quite inexperienced about this topic so I'm sure it sia just a trivial error from my side.

[jakevdp]

It looks like this is working as expected, given how you've defined your functions – GammaMC is the gradient of DeltaMC with respect to the first argument, and DeltaMC is the gradient of BlackMC with respect to the first argument.

So put together, you are taking the second derivative of your function with respect to the argument S, and because your function is linear in S, the second derivative is zero.

Does that make sense?

[frskplis]

Your Black Bcholes MC function is in fact linear, therefore second derivative wrt to price is 0. This is in fact unfortunate, as analytical gamma is definitely not zero for vanilla call. Please read this topic for a workaround (i.e. change the max function to heavisde) and the linked paper:
link

[G10VA]

Hi @jakevdp and @frskplis, thank you both for checking inside the post and suggesting the answer to the problem.
I assumed that the problem was related to a wrong implementation of the AAD while the origin was elsewhere.

Thank you again for your help.
It has been heavily appreciated.
Best regards.