Why Autodiff Misses Monte Carlo Greeks for Barrier Options

# Barrier option Greeks using AD


# Barrier option Greeks using AD












I am trying to price a Down-and-Out Barrier option using Monte Carlo and get the Greeks using autodiff as provided by PyTorch. However, comparing the output to bumping, I get vastly different values for the Greeks, and have verified this by bumping using QuantLib also. The price is correct and agrees with QL. This is my code:

```
import torch
import time

torch.set_default_dtype(torch.double)
torch.set_default_device('cpu')
torch.manual_seed(1234)

def simulate_barrier(spot, strike, barrier, vol, rate, expiry, numPaths: int, numSteps: int, normals):
    dt = expiry / numSteps
    times = torch.linspace(start=0., end=expiry, steps=numSteps+1)

    W = torch.concatenate([torch.zeros(numPaths, 1), torch.cumsum(normals, dim=1)], dim=1)
    gbm = spot * torch.exp((rate - vol**2/2)*times + vol*torch.sqrt(dt)*W)

    # (Paul Glasserman.1997. A Continuity correction for Discrete Barrier Options)
    alive = torch.all(gbm > (barrier + (barrier*0.5826 * vol * torch.sqrt(dt))), dim=1)
    payoffs = torch.where(alive, torch.max(gbm[:, -1] - strike, torch.tensor(0.)), torch.tensor(0.))
    return payoffs.mean() * torch.exp(-rate*expiry)

if __name__ == "__main__":
    spot = torch.tensor(100., requires_grad=True)
    strike = torch.tensor(100., requires_grad=True)
    rfr = torch.tensor(0.01, requires_grad=True)
    vol = torch.tensor(0.2, requires_grad=True)
    barrier = torch.tensor(90.)
    expiry = torch.tensor(1., requires_grad=True)

    numPaths = 100000
    numSteps = 252
    print(f"Paths:      {numPaths}")
    print(f"Time steps: {numSteps}")

    normals = torch.randn((numPaths, numSteps))

    delta1 = simulate_barrier(spot+1, strike, barrier, vol, rfr, expiry, numPaths, numSteps, normals)
    delta2 = simulate_barrier(spot-1, strike, barrier, vol, rfr, expiry, numPaths, numSteps, normals)

    theta1 = simulate_barrier(spot, strike, barrier, vol, rfr, expiry+0.01, numPaths, numSteps, normals)
    theta2 = simulate_barrier(spot, strike, barrier, vol, rfr, expiry-0.01, numPaths, numSteps, normals)
    
    rho1 = simulate_barrier(spot, strike, barrier, vol, rfr+0.001, expiry, numPaths, numSteps, normals)
    rho2 = simulate_barrier(spot, strike, barrier, vol, rfr-0.001, expiry, numPaths, numSteps, normals)
    
    vega1 = simulate_barrier(spot, strike, barrier, vol+0.01, rfr, expiry, numPaths, numSteps, normals)
    vega2 = simulate_barrier(spot, strike, barrier, vol-0.01, rfr, expiry, numPaths, numSteps, normals)
    
    start = time.time()
    price = simulate_barrier(spot, strike, barrier, vol, rfr, expiry, numPaths, numSteps, normals)
    price.backward()
    end = time.time()

    print(f"\nPrice: {price.item():.2f}")
    print(f"Delta: {spot.grad.item():.4f}")
    print(f"Bump Delta: {(delta1-delta2)/2:.4f}")
    print(f"Theta: {expiry.grad.item():.4f}")
    print(f"Bump Theta: {(theta1-theta2)/0.02:.4f}")
    print(f"Rho:   {rfr.grad.item():.4f}")
    print(f"Bump Rho: {(rho1-rho2)/0.002:.4f}")
    print(f"Vega:  {vol.grad.item():.4f}")
    print(f"Bump Vega: {(vega1-vega2)/0.02:.4f}")
    print(f"\nExecution time: {(end-start):.4f} seconds")
```

With the output:

```
Paths:      100000
Time steps: 252

Price: 6.89
Delta: 0.4027
Bump Delta: 0.7187
Theta: 4.0707
Bump Theta: 1.9855
Rho:   33.3845
Bump Rho: 42.1327
Vega:  33.3415
Bump Vega: 15.3068

Execution time: 0.8077 seconds
```

I cannot seem to find the issue here, and cannot tell whether it is to do with the nature of the Barrier option and difficulties in simulating it or a mistake in my reasoning.

Any thoughts/advice/comparison would be appreciated. Thanks!

## Answer by Andrea (score 2)

https://quant.stackexchange.com/a/81099

It is the nature of the barrier.

You cannot apply AD to discontinuous payoffs, so even a digital or an "asset or nothing" option will show problems.

But the question is a lot more fundamental: what do you actually mean by "derivative"?

If you plot the MonteCarlo price of a digital with a fixed number of paths and very very high numeric precision, you will get a piecewise flat function, and so its derivative is 0 almost everywhere (and AD gives you this).

But what you really mean by derivative is different: you want an approximation of the real one, i.e. infinite number of paths.

And with discontinuous payoff you cannot apply the Lebesgue dominated convergence theorem (or equivalent) and flip it with the sum (for infinite paths).