Debugging Monte Carlo Greeks for a European Call
Article Quant Q&A · Author: Hasek
Summary
The document investigates why a Monte Carlo European call valuation can appear close to its Black–Scholes price while its finite-difference delta and vega disagree sharply. The responses identify implementation issues, including a volatility argument that is not actually used in the path simulation, accidental dependence on a global variable, and seed handling that does not reliably produce comparable simulations. They also emphasize discounting simulated payoffs and accounting for dividends when calculating analytical Greeks.
A revised example simulates the terminal asset value directly under geometric Brownian motion and uses matched random seeds for bumped valuations, reducing noise in finite differences. Reported comparisons show the revised Monte Carlo price and Greeks close to analytic values in that example. These are illustrative results, not a general convergence guarantee; finite-difference estimates remain sensitive to bump size and simulation noise. The document’s sample Greek snippets also warrant careful review before reuse, since bump conventions must be applied consistently.
Key ideas
- The simulation must use the passed volatility parameter consistently rather than a global value.
- Discount simulated option payoffs before comparing them with Black–Scholes prices.
- Dividend yield affects the analytical call delta and should be included in comparisons.
- Using the same random draws for bumped prices can reduce Monte Carlo noise in finite-difference Greeks.
- Finite-difference estimates depend on correct bump conventions and remain subject to simulation error.
Tags
Full text
# Monte Carlo Greeks aren't matching Black-Scholes ones
# Monte Carlo Greeks aren't matching Black-Scholes ones
I'm implementing a Monte Carlo simulation of Geometric Brownian Motion for having a hands on model for pricing and hedging exotic barrier and autocallable payoffs. The idea is to test my implementation against Black-Scholes valuation of vanilla European call to make sure the price and greeks are more or less in line with analytical ones.
```
Monte-Carlo price 0.0034602543023745684 %
Monte-Carlo probability of expiration ITM 74.77669999999999 %
Monte-Carlo cashflow 34.60254302374568 USD
Black-Scholes price 0.0032218232666278717 %
Black-Scholes probability of expiration ITM 74.75673753241554 %
Black-Scholes cashflow 32.218232666278716 USD
Monte-Carlo delta 0.8232844460563724
Black-Scholes delta 0.7687961324513397
Monte-Carlo vega -4.402817163679625
Black-Scholes vega 63.063541724176005
```
Despite the price being close enough to analytical the greeks seems to be completely off on one million simulations with negative vega making no sense for a long European call. How did numerical vega ended up being negative here? Is there anything wrong with the delta as well not converging closer to the analytical one?
Please see my code below. Any help is appreciated.
```
import numpy as np
from scipy.stats import norm
import math
import datetime
notional = 1000000
coupon = 0.10 * notional
current_spot = 294.
initial_spot = 250
barrier = 1.10 * initial_spot
rate = 0.12
div = 0.05
sigma = 0.19
valuationDate = datetime.date(2025, 2, 17)
maturityDate = datetime.date(2025, 9, 22)
maturity = (maturityDate - valuationDate).days / 365
n_of_steps = (maturityDate - valuationDate).days
def MonteCarlo(S, r, d, vol, T, n_of_steps, n_of_sims, seed=np.random.seed(2025)):
seed
dt = float(T) / n_of_steps
paths = np.zeros((n_of_steps + 1, n_of_sims), np.float64)
paths[0] = S
for t in range(1, n_of_steps + 1):
rand = np.random.standard_normal(n_of_sims)
paths[t] = paths[t - 1] * np.exp((r - d - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * rand)
return paths
def callMonteCarlo(paths, strike):
price = np.mean(np.clip(paths[-1] - strike, a_min=0, a_max=None))
return price
def callBlackScholes(spot, strike, time, vol, rate=0., div=0.):
d1 = (math.log(spot / strike) + (rate - div + vol**2 / 2) * time) / (vol * math.sqrt(time))
d2 = d1 - vol * math.sqrt(time)
price = math.exp(- div * time) * norm.cdf(d1) * spot - norm.cdf(d2) * strike * math.exp(- rate * time)
return price
def callProbabilityInTheMoney(paths, strike):
n = len(paths[-1])
m = sum(s >= strike for s in paths[-1])
return float(m / n)
def N_d2(spot, strike, time, vol, rate=0., div=0.):
d1 = (math.log(spot / strike) + (rate - div + vol**2 / 2) * time) / (vol * math.sqrt(time))
d2 = d1 - vol * math.sqrt(time)
return norm.cdf(d2)
def analyticalCallDelta(spot, strike, time, vol, rate=0., div=0.):
d1 = (math.log(spot / strike) + (rate - div + vol**2 / 2) * time) / (vol * math.sqrt(time))
delta = math.exp(- div * time) * norm.cdf(d1)
return delta
def analyticalVega(spot, strike, time, vol, rate=0., div=0.):
d1 = (math.log(spot / strike) + (rate - div + vol**2 / 2) * time) / (vol * math.sqrt(time))
vega = spot * math.exp(- div * time) * norm.pdf(d1) * math.sqrt(time)
return vega
def numericalCallDelta(K, S, r, d, vol, T, n_of_steps, n_of_sims, seed=np.random.seed(2025)):
S_up = S + S / 100.
S_down = S - S / 100.
seed
paths_up = MonteCarlo(S_up, r, d, vol, T, n_of_steps, n_of_sims)
paths_down = MonteCarlo(S_down, r, d, vol, T, n_of_steps, n_of_sims)
price_up = callMonteCarlo(paths_up, K)
price_down = callMonteCarlo(paths_down, K)
delta = (price_up - price_down) / (S_up - S_down)
return delta
def numericalVega(K, S, r, d, vol, T, n_of_steps, n_of_sims, seed=np.random.seed(2025)):
vol_up = vol + 0.01
vol_down = vol - 0.01
seed
paths_up = MonteCarlo(S, r, d, vol_up, T, n_of_steps, n_of_sims)
paths_down = MonteCarlo(S, r, d, vol_down, T, n_of_steps, n_of_sims)
price_up = callMonteCarlo(paths_up, K)
price_down = callMonteCarlo(paths_down, K)
vega = (price_up - price_down) / (vol_up - vol_down)
return vega
seed = np.random.seed(2025)
paths = MonteCarlo(S=current_spot, r=rate, d=div, vol=sigma, T=maturity, n_of_steps=n_of_steps, n_of_sims=1000000, seed=seed)
priceMC = math.exp(- rate * maturity) * callMonteCarlo(paths=paths, strike=barrier) / notional
priceBS = math.exp(- rate * maturity) * callBlackScholes(spot=current_spot, strike=barrier, time=maturity, vol=sigma, rate=rate, div=div) / notional
print("Monte-Carlo price {0} %".format(100 * priceMC))
print("Monte-Carlo probability of expiration ITM {0} %".format(100 * monteCarloProbability))
print("Monte-Carlo cashflow {0} USD".format(priceMC * notional))
print()
print("Black-Scholes price {0} %".format(100 * priceBS))
print("Black-Scholes probability of expiration ITM {0} %".format(100 * blackScholesProbability))
print("Black-Scholes cashflow {0} USD".format(priceBS * notional))
print()
deltaMC = numericalCallDelta(K=barrier, S=current_spot, r=rate, d=div, vol=sigma, T=maturity, n_of_steps=n_of_steps, n_of_sims=1000000)
deltaBS = analyticalCallDelta(spot=current_spot, strike=barrier, time=maturity, vol=sigma, rate=rate, div=div)
print("Monte-Carlo delta {0}".format(deltaMC))
print("Black-Scholes delta {0}".format(deltaBS))
print()
vegaMC = numericalVega(K=barrier, S=current_spot, r=rate, d=div, vol=sigma, T=maturity, n_of_steps=n_of_steps, n_of_sims=1000000)
vegaBS = analyticalVega(spot=current_spot, strike=barrier, time=maturity, vol=sigma, rate=rate, div=div)
print("Monte-Carlo vega {0}".format(vegaMC))
print("Black-Scholes vega {0}".format(vegaBS))
```
## Answer by Andrea (score 1, accepted)
https://quant.stackexchange.com/a/81953
- Don't use global variables, not even in small scripts, they hide issues. `vol` is not used in `MonteCarlo()`.
- `seed=np.random.seed(2025)` as a default argument is only evaluated once. Make sure your `delta price` is 0 for a 0 bump, or you get extra MC noise.
## Answer by Pedro (score 5)
https://quant.stackexchange.com/a/81955
I think you had a few issues in your code, including discounting. This code works (notice you can sample the final spot in one step, you do not have to create unnecessary numerical noise)
```
import numpy as np
from scipy.stats import norm
import datetime
def black_scholes_mc_call(
spot0: float,
rate: float,
dividend: float,
sigma: float,
strike: float,
maturity: float,
num_paths: int,
seeds: tuple[int, ...],
) -> float:
partials = np.zeros(len(seeds))
for i, seed in enumerate(seeds):
normals = np.random.RandomState(seed).standard_normal(num_paths)
s_maturity = (
spot0
* np.exp((rate - dividend) * maturity)
* np.exp(-maturity * sigma**2 / 2)
* np.exp(sigma * np.sqrt(maturity) * normals)
)
val = np.exp(-maturity * rate) * np.mean(np.maximum(s_maturity - strike, 0))
partials[i] = val
return np.mean(partials)
def black_scholes_call(
spot: float,
strike: float,
time: float,
vol: float,
rate: float,
div: float,
) -> float:
fwd = spot * np.exp((rate - div) * time)
z = math.log(fwd / strike)
term_vol = vol * math.sqrt(time)
d1 = z / term_vol + term_vol / 2
d2 = z / term_vol - term_vol / 2
df = np.exp(-time * rate)
price = df * (fwd * norm.cdf(d1) - norm.cdf(d2) * strike)
return price
current_spot = 5000.0
initial_spot = 5000.0
rate = 0.05
barrier = (1 + rate) * initial_spot
dividend = 0.08
sigma = 0.15
valuationDate = datetime.date(2025, 1, 1)
maturityDate = datetime.date(2026, 1, 1)
maturity = (maturityDate - valuationDate).days / 365
priceMC = black_scholes_mc_call(
spot0=current_spot,
rate=rate,
dividend=dividend,
sigma=sigma,
maturity=maturity,
num_paths=int(1e7),
strike=barrier,
seeds=(1,2,3,4,5),
)
priceBS = black_scholes_call(
spot=current_spot,
strike=barrier,
time=maturity,
vol=sigma,
rate=rate,
div=dividend,
)
print(f"Monte-Carlo price {priceMC}")
print(f"Analytic price {priceBS}")
```
I get
```
Monte-Carlo price 136.68000878498842
Analytic price 136.71161616372243
```
For the Greeks, you have to be careful about the dividends, among other things:
```
def delta(
spot: float,
strike: float,
maturity: float,
volatility: float,
rate: float,
dividend: float,
):
fwd = spot * np.exp((rate - dividend) * maturity)
z = math.log(fwd / strike)
d1 = z / volatility * math.sqrt(maturity) + volatility * math.sqrt(maturity) / 2
return np.exp(-maturity * dividend) * norm.cdf(d1)
def delta_mc(
spot0: float,
strike: float,
bump: float,
rate: float,
div: float,
vol: float,
maturity: float,
num_paths: int,
seeds: tuple[int, ...],
):
price_up = black_scholes_mc_call(
spot0=spot0 + bump,
rate=rate,
dividend=div,
sigma=vol,
strike=strike,
maturity=maturity,
num_paths=num_paths,
seeds=seeds,
)
price_down = black_scholes_mc_call(
spot0=spot0,
rate=rate,
dividend=div,
sigma=vol,
strike=strike,
maturity=maturity,
num_paths=num_paths,
seeds=seeds,
)
delta = (price_up - price_down) / bump
return delta
```
Running
```
delta = delta(
spot=initial_spot,
strike=barrier,
maturity=maturity,
volatility=sigma,
rate=rate,
dividend=dividend,
)
deltamc = delta_mc(
spot0=initial_spot,
strike=barrier,
bump=1,
maturity=maturity,
vol=sigma,
rate=rate,
div=dividend,
num_paths=int(1e7),
seeds=(1, 2, 3, 4,5),
)
delta, deltamc
```
gives me
```
(0.3011747308109762, 0.30133159022952327)
```
which is acceptable.
Finally for vegas
```
def vega_mc(
spot0: float,
strike: float,
rate: float,
dividend: float,
vol: float,
vol_bump: float,
maturity: float,
num_paths: int,
seeds: tuple[int, ...],
):
price_up = black_scholes_mc_call(
spot0,
rate,
dividend,
vol + vol_bump,
strike,
maturity,
num_paths,
seeds,
)
price_down = black_scholes_mc_call(
spot0,
rate,
dividend,
vol,
strike,
maturity,
num_paths,
seeds,
)
return (price_up - price_down) / vol_bump
```
and running
```
vegaa = vega(initial_spot, barrier, maturity, sigma, rate, dividend)
vegamc = vega_mc(
spot0=initial_spot,
strike=barrier,
rate=rate,
dividend=dividend,
vol=sigma,
vol_bump=sigma*0.005,
maturity=maturity,
num_paths=1000000,
seeds=(1, 2, 3, 4),
)
vegaa, vegamc
```
I get
```
(1663.8411096873679, 1661.890544055268)
```
which again is in the ballpark. Here are some plots, with 200000 paths per strike: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.