there is a fair number of questions on gradients out there, but I haven't been able to fix my problem. In a nutshell: Trying to run a Monte Carlo simulation and get pathwise differentials. There are a few tutorials out there, but they do run into the same problem as my own code. I have boiled it down into the toy example below. The second order derivative called gamma is wrong however (I'll highlight it below). So something is wrong with my nested gradient tape. Having the tape watch the variable has no effect actually. There must be something I am not aware of here. Any hint much appreciated.
edit: I have cross checked this with a deterministic function and the code works just fine. Second derivatives are calculated correctly using gradient tape. So no idea why it doesn't work on a Monte Carlo.
import numpy as np
import pandas as pd
import tensorflow as tf
from pprint import pprint
DTYPE = tf.float32
SEED = 3232
S0 = tf.Variable(100, dtype=DTYPE)
strike = tf.Variable(110, dtype=DTYPE)
time_to_expiry = tf.Variable(1, dtype=DTYPE)
implied_vol = tf.Variable(0.3, dtype=DTYPE)
v = dict(S0=S0,strike=strike,time_to_expiry=time_to_expiry,implied_vol=implied_vol)
@tf.function
def brownian(S0, dt, sigma, mu, dw):
dt_sqrt = tf.math.sqrt(dt)
shock = sigma * dt_sqrt * dw
drift = (mu - (sigma ** 2) / 2)
bm = tf.math.exp(drift * dt + shock)
out = S0 * tf.math.cumprod(bm, axis=1)
return out
@tf.function
def pricer_montecarlo(S0, strike, time_to_expiry, implied_vol, dw):
sigma = implied_vol
T = time_to_expiry
r = tf.constant(0.0,dtype=DTYPE)
K = strike
dt = T / dw.shape[1]
st = brownian(S0, dt, sigma, r, dw)
payout = tf.math.maximum(st[:, -1] - K, 0)
npv = tf.exp(-r * T) * tf.reduce_mean(payout)
return npv
def calculate_montecarlo(greeks=True):
nsims = 10**7
nobs = 2
dw = tf.random.normal((nsims, nobs), seed=SEED)
out = dict()
if greeks:
with tf.GradientTape() as g2:
g2.watch(v['S0'])
with tf.GradientTape() as g1:
g1.watch(v['S0'])
npv = pricer_montecarlo(**v, dw=dw)
dv = g1.gradient(npv, v)
g2.watch(dv)
d2v = g2.gradient(dv['S0'], v)
out["dv"] = {k: v.numpy() for k, v in dv.items()}
out["d2v"] = {k: v.numpy() for k, v in d2v.items()}
else:
npv = pricer_montecarlo(**v, dw=dw).numpy()
out["npv"] = npv.numpy()
return out
out = calculate_montecarlo()
pprint(out)
from py_vollib import black_scholes
from py_vollib.black_scholes.greeks import analytical
print('npv='+str(black_scholes.black_scholes('c', 100, 110, 1, 0, 0.3)))
print('dv S0='+str(analytical.delta('c', 100, 110, 1, 0, 0.3)))
print('d2v S0='+str(analytical.gamma('c', 100, 110, 1, 0, 0.3)))
print('dv implied_vol='+str(analytical.vega('c', 100, 110, 1, 0, 0.3)))
print('dv time_to_expiry='+str(analytical.theta('c', 100, 110, 1, 0, 0.3)))
Output: 'd2v': {'S0': 0.0..., should be close to 0.013112390443974165 (plus some stochastic noise).
{'d2v': {'S0': 0.0,
'implied_vol': 0.3933603,
'strike': 0.0,
'time_to_expiry': 0.059004053},
'dv': {'S0': 0.43342653,
'implied_vol': 39.336025,
'strike': -0.32001525,
'time_to_expiry': 5.9004045},
'npv': 8.140971}
npv=8.141012048964207
dv S0=0.4334094123285094
d2v S0=0.013112390443974165
dv implied_vol=0.39337171331922494
dv time_to_expiry=-0.01616596082133801