Why Automatic Differentiation Can Miss Monte Carlo Gamma
Article Quant Q&A · Author: Matt
Summary
This document investigates why second-order automatic differentiation gives nearly zero Gamma for a Monte Carlo Black call pricer, even though the analytical Black formula produces a plausible result. The example differentiates a simulated payoff containing a maximum function. Away from the exercise boundary, that payoff is piecewise linear in the underlying, so pathwise second derivatives are zero; the curvature of the option value comes from paths crossing the kink, which ordinary pathwise differentiation does not capture directly.
The accepted answer proposes replacing the hard payoff indicator with a smoothed Heaviside approximation and uses antithetic normal draws. Its reported Monte Carlo AAD Greeks are compared with finite-difference estimates and analytical Black Greeks, illustrating the effect in that example. Smoothing introduces an approximation whose bandwidth can affect bias and stability, and finite-difference estimates also depend on bump size and sampling noise. The document offers one workaround, not a general guarantee for all payoffs or AAD implementations.
Key ideas
- A call payoff’s kink can make pathwise second derivatives miss the option’s Gamma.
- The example’s hard maximum function is the source identified for the near-zero Monte Carlo Gamma.
- The accepted answer smooths the payoff boundary with a Heaviside approximation.
- Antithetic sampling is used alongside smoothing in the revised example.
- Results are compared with finite differences and analytical Black Greeks, but smoothing and sampling choices can affect estimates.
Tags
Full text
# use either jax or autograd functions, same format
# 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 import grad # from jax import grad
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):
# setup gradient functions to evaluate with AAD
gradient_func = grad(MC_call_price, (0)) # tuple specifies inputs we want to differentiate
gradient_func2 = grad(gradient_func) # second derivative
# 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):
# setup gradient functions to evaluate with AAD
gradient_func = grad(BS_call_price, (0)) # tuple specifies inputs we want to differentiate
gradient_func2 = grad(gradient_func) # second derivative
# 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
AAD: Delta 0.7497560731718061 Gamma 5.204170427930421e-18
FD: Delta 0.7497588343497341 Gamma 0.013656582353860358
BS Option price 14.880707205664976
AAD: Delta 0.7496697714173796 Gamma 0.01499062307583833
```
Thanks for any observations or suggestions on how to fix Gamma! I've tried the same programming model for more advanced options (baskets) and the Gamma is always close to 0, so I imagine I am missing something.
## Answer by Yong-guang Gong (score 3, accepted)
https://quant.stackexchange.com/a/74940
I think the issues is because of the payoff function. You should replace the maximum() with HeavisideApprox(). Read this paper: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=1626547
Here's the codes:
```
# use either jax or autograd functions, same format
import autograd.numpy as np # import jax.numpy as np
from autograd import grad # from jax import grad
from autograd.scipy.stats import norm # from jax.scipy.stats import norm
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_call_price(F, vol, K, T, IR, steps, trials):
np.random.seed(123)
dt = T/steps
rands = np.random.normal(size=(int(steps),int(trials)));
rands = np.hstack((rands,-rands))
paths = np.sum(((IR - 0.5*vol**2)*dt + vol*np.sqrt(dt) *rands ), axis=0)
payoff = F*np.exp(paths)-K
payoff = HeavisideApprox(payoff, F* vol *np.sqrt(T)* 0.05)* payoff
return np.mean(payoff)* np.exp(-IR*T)
def AAD_Greeks(F, vol, K, T, IR, steps, trials):
# setup gradient functions to evaluate with AAD
gradient_func = grad(MC_call_price, (0)) # tuple specifies inputs we want to differentiate
gradient_func2 = grad(gradient_func,(0)) # second derivative
gradient_func3 = grad(gradient_func,(1)) # second derivative
gradient_func4 = grad(MC_call_price,(1)) #
gradient_func5 = grad(gradient_func4,(1)) # second derivative
# solve for Greeks
delta = gradient_func(F, vol, K, T, IR, steps, trials)
gamma = gradient_func2(F, vol, K, T, IR, steps, trials)
vanna = gradient_func3(F, vol, K, T, IR, steps, trials)
vega = gradient_func4(F, vol, K, T, IR, steps, trials)
volga = gradient_func5(F, vol, K, T, IR, steps, trials)
return delta, gamma, vanna, vega, volga
def BS_call_price(F, vol, K, T, IR):
b = np.exp(-IR*T)
x1 = np.log(F/(b*K))/vol/np.sqrt(T) + 0.5*vol*np.sqrt(T)
x2 = x1 - vol* np.sqrt(T)
z1 = F * norm.cdf(x1)
z2 = b* K* norm.cdf(x2)
return z1 - z2
def AAD_Greeks2(F, vol, K, T, IR):
# setup gradient functions to evaluate with AAD
gradient_func = grad(BS_call_price, (0)) # tuple specifies inputs we want to differentiate
gradient_func2 = grad(gradient_func,(0)) # second derivative
gradient_func3 = grad(gradient_func,(1)) # second derivative
gradient_func4 = grad(BS_call_price,(1)) #
gradient_func5 = grad(gradient_func4,(1)) # second derivative
# solve for Greeks
delta = gradient_func(F, vol, K, T, IR)
gamma = gradient_func2(F, vol, K, T, IR)
vanna = gradient_func3(F, vol, K, T, IR)
vega = gradient_func4(F, vol, K, T, IR)
volga = gradient_func5(F, vol, K, T, IR)
return delta, gamma, vanna, vega, volga
if __name__ == '__main__':
F = 100.0
vol = 0.2
K = 100.
T = 0.5
IR = 0.03
steps = 100
trials = 10000
bump = 0.001
optionprice = MC_call_price(F, vol, K, T, IR, steps, trials)
delta, gamma, vanna, vega, volga = AAD_Greeks(F, vol, K, T, IR, steps, trials)
print("\nMC Option price", optionprice,"\n", "\nAAD: Delta", delta, "Gamma", gamma, "vanna", vanna, "vega", vega, "volga", volga)
up1 = MC_call_price(F*(1+bump), vol, K, T, IR, steps, trials)
down1 = MC_call_price(F*(1-bump), vol, K, T, IR, steps, trials)
FDdelta = (up1-down1)/(2*bump*F)
FDgamma = (up1 - 2*optionprice + down1)/(bump**2*F**2)
up2 = MC_call_price(F, vol+bump, K, T, IR, steps, trials)
down2 = MC_call_price(F, vol-bump, K, T, IR, steps, trials)
FDvega = (up2-down2)/(2*bump)
FDvolga = (up2 - 2*optionprice + down2)/(bump**2)
up3 = MC_call_price(F*(1+bump), vol+bump, K, T, IR, steps, trials)
up4 = MC_call_price(F*(1+bump), vol-bump, K, T, IR, steps, trials)
down3 = MC_call_price(F*(1-bump), vol+bump, K, T, IR, steps, trials)
down4 = MC_call_price(F*(1-bump), vol-bump, K, T, IR, steps, trials)
FDvanna = (up3 +down4 - up4- down3)/F/bump/bump/4
print("FD: Delta", FDdelta, "Gamma", FDgamma, "Vanna", FDvanna, "Vega", FDvega, "Volga", FDvolga)
BSoptionprice = BS_call_price(F, vol, K, T, IR)
BSdelta, BSgamma, BSVanna, BSVega, BSVolga = AAD_Greeks2(F, vol, K, T, IR)
print("\n\nBS Option price", BSoptionprice,"\n", "\nAAD: Delta", BSdelta, "Gamma", BSgamma, "Vanna", BSVanna, "Vega", BSVega, "Volga", BSVolga)
```
All works with AAD, the outputs:
```
MC Option price 6.294541301567644
AAD: Delta 0.5703929155658839 Gamma 0.02985300460851462 vanna -0.09991144416894931 vega 27.335608708711906 volga 0.7481238019846563
FD: Delta 0.5703875381820334 Gamma 0.029732742335131235 Vanna -0.09715704481827457 Vega 27.33560790448486 Volga 0.7617713926322267
BS Option price 6.371027942167473
AAD: Delta 0.5701581024006674 Gamma 0.027772131739916557 Vanna -0.0694303293497937 Vega 27.772131739916563 Volga 0.8678791168724274
```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.