Notebook

The purpose of this project is to implement formulas for the Greeks and then use them to test various methods of doing Monte Carlo Greeks. The functions written in the first section will be used in later projects.

In [1]:

# Mount google drive and setup paths
import sys
 
# mount google drive 
from google.colab import drive
drive.mount('/content/drive/', force_remount=True)
 
# specify path name to import py file
folder_name = 'Jupyter Notebooks/Quantitative Finance'
sys.path.append("drive/My Drive/" + folder_name)

Mounted at /content/drive/

In [0]:

# Import necessary libraries
from random import random
from random import randint
from random import seed
import pandas as pd

from GBM import *
from Option import *
from Greeks import *
from finite_difference_method import *

Implementing the formulas¶

Implement the following formulas for both European call and put options

(i) the Delta (spot derivative),

(ii) the Gamma (2nd spot derivative),

(iii) the Vega (volatility derivative),

(iv) the Rho ( $r$ derivative),

(v) the Theta of an option, this is equal to minus the $T$ derivative,

The formulas are deducible by simply differentiating the Black-Scholes prices.

All implementations of greeks can be found at Option.py, class Greeks. Let $\tau = T-t.$ Since

$\begin{align*} c(t,S) & = S N(d_1) - Ke^{-r\tau} N(d_2) \quad \text{where} \\ d_1 & = \frac{\ln \left( \frac{S}{K} \right) + (r + 0.5 \sigma^2) \tau}{\sigma\sqrt{\tau}}, \\ d_2 & = \frac{\ln \left( \frac{S}{K} \right) + (r - 0.5 \sigma^2) \tau}{\sigma\sqrt{\tau}}, \\ d_1 - d_2 & = \sigma\sqrt{\tau}, \end{align*}$ by using the following facts

$\begin{align*} S N'(d_1) & = Ke^{-r\tau} N'(d_2), \\ \frac{\partial d_1}{\partial S} & = \frac{\partial d_2}{\partial S} = \frac{1}{S\sigma\sqrt{\tau}}, \\ \frac{\partial d_1}{\partial \sigma} & = - \frac{d_2}{\sigma} \quad, \frac{\partial d_2}{\partial \sigma} = - \frac{d_1}{\sigma}, \\ \frac{\partial d_1}{\partial r} & = \frac{\partial d_2}{\partial r}, \\ \end{align*}$

we have

$\begin{align*} \text{Delta} : \frac{\partial c}{\partial S} & = N(d_1), \\ \text{Gamma} : \frac{\partial^2 c}{\partial S^2} & = \frac{N'(d_1)}{S\sigma\sqrt{\tau}}, \\ \text{Vega} : \frac{\partial c}{\partial \sigma} & = S\sqrt{\tau} N'(d_1) , \\ \text{Rho} : \frac{\partial c}{\partial r} & = K\tau e^{-r\tau} N(d_2), \\ \text{Theta} : \frac{\partial c}{\partial t} & = -rKe^{-r\tau}N(d_2) - \frac{SN'(d_1)\sigma}{2\sqrt{\tau}} \\ \end{align*}$

In [0]:

# Testing

S = 120
K = 100
r = 0.05
d = 0
sigma = 0.2
T = 1

greeks = Greeks(S, K, r, d, sigma, T)
greeks_formula = greeks.get_all()

call_greeks = greeks_formula['call']
put_greeks = greeks_formula['put']

call_greeks = pd.DataFrame.from_dict(call_greeks, orient='index', columns = ['formula'])
put_greeks = pd.DataFrame.from_dict(put_greeks, orient='index', columns = ['formula'])

In [4]:

call_greeks

Out[4]:

	formula
delta	0.896455
gamma	0.007500
vega	21.600708
rho	81.405559
theta	-6.230349

In [5]:

put_greeks

Out[5]:

	formula
delta	-0.103545
gamma	0.007500
vega	21.600708
rho	-13.717384
theta	-1.474202

Testing¶

Test them all by comparing with the finite differencing price. That is let $\epsilon$ be a small number, and compute the price change for bumping the parameter for $\epsilon$ and divide by $\epsilon$ . To approximate the delta, for example, take

$\frac{1}{\epsilon}\left( BS(S + \epsilon, T, \sigma, r, d) - BS(S, T, \sigma, r, d) \right)$ The Gamma can be approximated by finite differencing the Delta, or by taking the formula

$\frac{1}{\epsilon^2} \left( BS(S + \epsilon, T, \sigma, r, d) - 2BS(S, T, \sigma, r, d) + BS(S - \epsilon, T, \sigma, r, d) \right)$ (Why does this work?)

In [9]:

e = 0.001 # error in finite difference method

f_S_c = lambda x: Option(S + x, K, r, d, sigma, T).european_call()
f_vol_c = lambda x: Option(S, K, r, d, sigma + x, T).european_call()
f_r_c = lambda x: Option(S, K, r + x, d, sigma, T).european_call()
f_T_c = lambda x: -Option(S, K, r, d, sigma, T + x).european_call()


delta_fc_c = finite_difference(f_S_c, e)
gamma_fc_c = finite_difference(f_S_c, e, second_order = True)
vega_fc_c = finite_difference(f_vol_c, e)
rho_fc_c = finite_difference(f_r_c, e)
theta_fc_c = finite_difference(f_T_c, e)

greeks_approx_c = [delta_fc_c, gamma_fc_c, vega_fc_c, rho_fc_c, theta_fc_c]

call_greeks['Finite Difference'] = greeks_approx_c
call_greeks['FD Abs error (FD)'] = np.abs(call_greeks.iloc[:,0] - call_greeks.iloc[:,1])
call_greeks

Out[9]:

	formula	Finite Difference	FD Abs error (FD)
delta	0.896455	0.896459	3.750064e-06
gamma	0.007500	0.007500	1.578147e-08
vega	21.600708	21.672831	7.212281e-02
rho	81.405559	81.418740	1.318111e-02
theta	-6.230349	-6.229884	4.648117e-04

In [10]:

f_S_p = lambda x: Option(S + x, K, r, d, sigma, T).european_put()
f_vol_p = lambda x: Option(S, K, r, d, sigma + x, T).european_put()
f_r_p = lambda x: Option(S, K, r + x, d, sigma, T).european_put()
f_T_p = lambda x: -Option(S, K, r, d, sigma, T + x).european_put()

delta_fc_p = finite_difference(f_S_p, e)
gamma_fc_p = finite_difference(f_S_p, e, second_order = True)
vega_fc_p = finite_difference(f_vol_p, e)
rho_fc_p = finite_difference(f_r_p, e)
theta_fc_p = finite_difference(f_T_p, e)

greeks_approx_p = [delta_fc_p, gamma_fc_p, vega_fc_p, rho_fc_p, theta_fc_p]

put_greeks['Finite Difference'] = greeks_approx_p
put_greeks['FD Abs error (FD)'] = np.abs(put_greeks.iloc[:,0] - put_greeks.iloc[:,1])
put_greeks

Out[10]:

	formula	Finite Difference	FD Abs error (FD)
delta	-0.103545	-0.103541	3.750054e-06
gamma	0.007500	0.007500	1.045240e-08
vega	21.600708	21.672831	7.212281e-02
rho	-13.717384	-13.656657	6.072673e-02
theta	-1.474202	-1.473856	3.459100e-04

In [0]:

# Question: Why does it work?

Graphs¶

Once you have all the formulas working and tested. Plot the following graphs and try to interpret them.

(i) The Delta of a call option as a function of spot.

(ii) The Delta of a call option as a function of time for in-the-money, out-of-the-money and at-the-money options.

(iii) The Gamma of a call option as a function of spot,

(iv) The Vega of a call option as a function of volatility, as a function of spot and as a function of time.

In [0]:

# (i) 

S_upper = 100
graph_num = 5
x = np.linspace(1, S_upper, S_upper)

for j in range(graph_num):
    K = randint(1,100)
    r = random()
    T = random()
    d = random()
    sigma = random()
    y = [Greeks(S, K, r, d, sigma, T).delta() for S in range(1, S_upper+1)]
    plt.plot(x,y)
    plt.xlabel('Spot')
    plt.ylabel('Delta')

Out[0]:

From above, we observe that delta of a call option is between $0$ and $1$ because it is a cdf. It is convex at the beginning then it becomes concave.

In [0]:

# (ii) 

K = 100
r = 0
d = 0
sigma = 0.2

S_upper = 200
T_upper = 10
x = np.linspace(1, S_upper, S_upper)

for t in reversed(range(1,T_upper)):
    y = [Greeks(S, K, r, d, sigma, t).delta() for S in range(1, S_upper+1)]
    plt.plot(x,y, label = str(t))
    plt.xlabel('Spot')
    plt.ylabel('Delta')
    plt.vlines(K, 0, 1)
    plt.legend()

Out[0]:

From above, if time to maturity decreases, then delta of an ITM call option will increase, OTM will decrease and ATM will approach to $0.5$

In [0]:

# (iii) 

S_upper = 100
graph_num = 5
x = np.linspace(1, S_upper, S_upper)

for j in range(graph_num):
    K = randint(1,100)
    r = random()
    T = random()
    d = random()
    sigma = random()
    y = [Greeks(S, K, r, d, sigma, T).gamma() for S in range(1, S_upper+1)]
    plt.plot(x,y)
    plt.xlabel('Spot')
    plt.ylabel('Gamma')

Out[0]:

From above, gamma of a call option is the normal density function.

(iv) The Vega of a call option as a function of volatility, as a function of spot and as a function of time.

In [0]:

# (iv) function of volatility

K = 100
r = 0
d = 0
T = 1

S_upper = 200
sigma_upper = 10
graph_num = 5
x = np.linspace(1, S_upper, S_upper)

for vol in range(1,sigma_upper):
    y = [Greeks(S, K, r, d, vol / 10, T).vega() for S in range(1, S_upper+1)]
    plt.plot(x,y, label = str(vol))
    plt.xlabel('Volatility')
    plt.ylabel('Vega')
    plt.legend()

Out[0]:

Recall that time to maturity and volatility have the same effect to all option Greeks.

In [0]:

# (iv)   # function of spot

S_upper = 100
graph_num = 5
x = np.linspace(1, S_upper, S_upper)

for j in range(graph_num):
    K = randint(1,100)
    r = random()
    T = random()
    d = random()
    sigma = random()
    y = [Greeks(S, K, r, d, sigma, T).vega() for S in range(1, S_upper+1)]
    plt.plot(x,y)
    plt.xlabel('Spot')
    plt.ylabel('Vega')

Out[0]:

In [0]:

# (iv)  # function of time

K = 100
r = 0
d = 0
sigma = 0.2

S_upper = 200
T_upper = 10
x = np.linspace(1, S_upper, S_upper)

for t in range(1,T_upper):
    y = [Greeks(S, K, r, d, sigma, t).vega() for S in range(1, S_upper+1)]
    plt.plot(x,y, label = str(t))
    plt.xlabel('Time to maturity')
    plt.ylabel('Vega')
    plt.legend()

Out[0]:

Monte Carlo Greeks¶

We can now try various Monte Carlo methods for computing the Greeks of vanilla options and compare them to the analytical formulas. Implement the following methods and compare to the formulas. How do the convergence speeds compare?

(i) Run the Monte Carlo twice. The second time with the parameter slightly bumped, and finite difference to get the Greek. Use different random numbers for the two simulations.

(ii) Do the same again but using the same random numbers for the two simulations. (Depending upon the language you are using, it will either default to different random numbers or default to the same ones. Setting the random number seed is the way to achieve either.)

(iii) Implement the pathwise derivative estimate method for the Delta.

(iv) Implement the likelihood ratio method for the Delta.

For more information on pathwise derivative estimate and likelihood ratio methods, please refer to the Glasserman's Monte Carlo Methods in Financial Engineering.

Intuitively, the pathwise derivative estimate method allows one to calculate greeks by using the ideas of

(1) interchanging derivative and integration (expectation), if applicable and

(2) the Chain rule.

For example, if we want to calculate delta of an European call option, since $S_T$ is a function of $S_0,$ then pathwise method implies that

$\begin{align*} \frac{\partial c}{\partial S_0} & = \frac{\partial}{\partial S_0} \mathbb{E}[e^{-rT}(S_T-K)^+] \\ & = \mathbb{E}\left[ \frac{\partial}{\partial S_0} e^{-rT}(S_T-K)^+ \right] \\ & = \mathbb{E} \left[e^{-rT} 1_{\{S_T>K\}} \frac{\partial S_T}{\partial S_0} \right] \\ & = \mathbb{E} \left[ e^{-rT} 1_{\{S_T>K\}} e^{(r-d - \frac{1}{2} \sigma^2)T + \sigma \sqrt{T} Z} \right] \\ & = \mathbb{E} \left[ e^{-rT} 1_{\{S_T>K\}} \frac{S_T}{S_0}\right], \end{align*}$ where the expectation can be calculated using Monte Carlo method.

The pathwise derivative method can be modified to calculate other greeks as well.

However, one drawback of the pathwise derivative method is that the payoff function must be differentiable or differentiable almost everywhere. In the example above, the payoff function $(S_T-K)^+$ is differentiable everywhere except at $S_T = K$ .

The likelihood ratio method uses the fact that the random variable has density function. Assume that $X = (X_1,...,X_n)$ has a probability density function $g_\theta$ with parameter $\theta.$ (One can think of $X_1,...,X_n$ as underlying financial instruments.) Then

$\mathbb{E}(Y) = \mathbb{E}_\theta(f(X_1,...,X_n)) = \int_{\mathbb{R}^n }f(x)\cdot g_\theta(x) dx.$ Assume that we can interchange derivativ and integration. Then

$\frac{d}{d\theta}\mathbb{E}_\theta(f(X_1,...,X_n)) = \int_{\mathbb{R}^n }f(x)\cdot \frac{ \frac{d g_\theta(x)}{d\theta} }{g_\theta(x)} g_\theta(x) dx = \mathbb{E}_\theta \left[ f(X) \cdot \frac{ \frac{d g_\theta(X)}{d\theta} }{g_\theta(X)} \right].$

For example, if we want to calculate delta of a European call option, then

$\theta = S_0, \quad f(S_T) = e^{-rT} (S_T-K)^+, \quad g_\theta(x) = \frac{1}{x\sigma\sqrt{T}} \phi(\xi(x))$ where

$\xi(x) = \frac{ \ln \left( \frac{x}{S_0} \right) - (r - \frac{1}{2}\sigma^2)T }{\sigma\sqrt{T}}.$

Note that $g_\theta$ is obtained by using $S_T = S_0 e^{ (r - \frac{1}{2}\sigma^2)T + \sigma \sqrt{T} Z }.$ Since $\frac{\phi'(x)}{\phi(x)} = -x,$ so

$\frac{\frac{d g_\theta(x)}{d S_0}}{g_\theta(x)} = \frac{\phi'(\xi(x))\cdot \xi'(x)}{\phi(\xi(x))} = -\xi(x) \xi'(x) = \frac{ \ln \left( \frac{x}{S_0} \right) - (r - \frac{1}{2}\sigma^2)T }{ S_0\sigma^2 T} .$

It follows that

$\Delta = \frac{d}{d S(0)}\mathbb{E} [e^{-rT} (S_T-K)^+] = \mathbb{E} \left[ e^{-rT}(S_T-K)^+ \left( \frac{ \ln \left( \frac{S_T}{S_0} \right) - (r - \frac{1}{2}\sigma^2)T }{ S_0\sigma^2 T} \right) \right].$

In [0]:

# (ii)

e = 0.001

S0 = 120
K = 100
r = 0.05
d = 0
sigma = 0.2
T = 1

num_steps = 100
num_paths = 17

# apply Monte Carlo to approximate european call's greeks
BS_mc_S_c = lambda x: black_scholes_monte_carlo_pricer('European call', S0 + x, K, r, d, sigma, T, num_steps, num_paths, False, 0)
BS_mc_r_c = lambda x: black_scholes_monte_carlo_pricer('European call', S0, K, r + x, d, sigma, T, num_steps, num_paths, False, 0)
BS_mc_vol_c = lambda x: black_scholes_monte_carlo_pricer('European call', S0, K, r, d, sigma + x, T, num_steps, num_paths, False, 0)
BS_mc_T_c = lambda x: - black_scholes_monte_carlo_pricer('European call', S0, K, r, d, sigma, T + x, num_steps, num_paths, False, 0)

delta_mc_c = finite_difference(BS_mc_S_c, e)
gamma_mc_c = finite_difference(BS_mc_S_c, e, True)
vega_mc_c = finite_difference(BS_mc_vol_c, e)
rho_mc_c = finite_difference(BS_mc_r_c, e)
theta_mc_c = finite_difference(BS_mc_T_c, e)

greeks_mc_c = [delta_mc_c, gamma_mc_c, vega_mc_c, rho_mc_c, theta_mc_c]

call_greeks['MC Greeks'] = greeks_mc_c
call_greeks['MC abs error'] = np.abs(call_greeks.iloc[:, 0] - call_greeks.iloc[:,3])

# apply Monte Carlo to approximate european put's greeks
BS_mc_S_p = lambda x: black_scholes_monte_carlo_pricer('European put', S0 + x, K, r, d, sigma, T, num_steps, num_paths, False, 0)
BS_mc_r_p = lambda x: black_scholes_monte_carlo_pricer('European put', S0, K, r + x, d, sigma, T, num_steps, num_paths, False, 0)
BS_mc_vol_p = lambda x: black_scholes_monte_carlo_pricer('European put', S0, K, r, d, sigma + x, T, num_steps, num_paths, False, 0)
BS_mc_T_p = lambda x: - black_scholes_monte_carlo_pricer('European put', S0, K, r, d, sigma, T + x, num_steps, num_paths, False, 0)

delta_mc_p = finite_difference(BS_mc_S_p, e)
gamma_mc_p = finite_difference(BS_mc_S_p, e, True)
vega_mc_p = finite_difference(BS_mc_vol_p, e)
rho_mc_p = finite_difference(BS_mc_r_p, e)
theta_mc_p = finite_difference(BS_mc_T_p, e)

greeks_mc_p = [delta_mc_p, gamma_mc_p, vega_mc_p, rho_mc_p, theta_mc_p]

put_greeks['MC Greeks'] = greeks_mc_p
put_greeks['MC abs error'] = np.abs(put_greeks.iloc[:, 0] - put_greeks.iloc[:,3])

In [12]:

call_greeks

Out[12]:

	formula	Finite Difference	FD Abs error (FD)	MC Greeks	MC abs error
delta	0.896455	0.896459	3.750064e-06	0.896991	0.000536
gamma	0.007500	0.007500	1.578147e-08	0.002801	0.004699
vega	21.600708	21.672831	7.212281e-02	21.858899	0.258191
rho	81.405559	81.418740	1.318111e-02	81.346387	0.059172
theta	-6.230349	-6.229884	4.648117e-04	-6.244217	0.013868

In [13]:

put_greeks

Out[13]:

	formula	Finite Difference	FD Abs error (FD)	MC Greeks	MC abs error
delta	-0.103545	-0.103541	3.750054e-06	-0.103460	0.000085
gamma	0.007500	0.007500	1.045240e-08	0.002801	0.004699
vega	21.600708	21.672831	7.212281e-02	21.640966	0.040258
rho	-13.717384	-13.656657	6.072673e-02	-13.667977	0.049407
theta	-1.474202	-1.473856	3.459100e-04	-1.469397	0.004805

Note that the vega generated using Monte Carlo deviates from true value significantly compared to other Greeks.

Recall that the pathwise derivative esimate method propses that

$\begin{align*} \Delta = \frac{\partial c}{\partial S_0} & = \mathbb{E} \left[ e^{-rT} 1_{\{S_T>K\}} \frac{S_T}{S_0}\right], \\ \frac{\partial c}{\partial \sigma} & = \mathbb{E}\left[ e^{-rT} 1_{\{ S_T>K \}}\left( \frac{\log \left( \frac{S_T}{S_0}\right) - (r + 0.5\sigma^2)T}{\sigma} \right) S_T \right] . \end{align*}$

Note that pathwise derivative estimate method fails for $Y = e^{-rT} 1_{\{S_T>K\}}$ because

$0 = \mathbb{E} \left[ \frac{\partial Y}{\partial S_0} \right] \neq \frac{\partial}{\partial S_0} \mathbb{E}[Y].$ Therefore, the method cannot be used to calculate gamma of european call option.

On the other hand, pathwise derivative estimate method for European put option gives

$\begin{align*} \Delta = \frac{\partial P}{\partial S_0} & = \mathbb{E} \left[ -e^{-rT} 1_{\{S_T<K\}} \frac{S_T}{S_0}\right], \\ \frac{\partial c}{\partial \sigma} & = \mathbb{E}\left[ -e^{-rT} 1_{\{ S_T<K \}}\left( \frac{\log \left( \frac{S_T}{S_0}\right) - (r + 0.5\sigma^2)T}{\sigma} \right) S_T \right] . \end{align*}$

In [0]:

# (iii) pathwise derivative estimate method

S0 = 120
K = 100
r = 0.05
d = 0
sigma = 0.2
T = 1

num_steps = 100
num_paths = 17

call_greeks['Pathwise'] = 0
call_greeks['Pathwise Abs err'] = 0

put_greeks['Pathwise'] = 0
put_greeks['Pathwise Abs err'] = 0

# generate stock price paths
ST = GBM_fd(S0, r, d, sigma, T, num_steps, 2**num_paths, plot = False, seed = None)

# pathwise estimate for delta of European call option
dpayoff_dS0 = np.exp(-r * T)  * np.piecewise(ST, [ST > K, ST <= K], [1, 0])
payoff = dpayoff_dS0 * ST / S0
call_greeks.loc['delta', 'Pathwise'] = monte_carlo(payoff)
call_greeks.loc['delta', 'Pathwise Abs err'] = np.abs(call_greeks.loc['delta', 'formula'] - call_greeks.loc['delta', 'Pathwise'])

# print('Analytical delta: ', greeks.delta('call'))
# print('Pathwise method delta: ' , monte_carlo(payoff))
# print('abs difference:', np.abs(greeks.delta('call') - monte_carlo(payoff)))

# pathwise estimate for vega of European call option
d1_similar = (np.log( ST / S0 ) - (r + 0.5 * sigma**2) * T) / (sigma)

payoff = dpayoff_dS0 * d1_similar * ST
call_greeks.loc['vega', 'Pathwise'] = monte_carlo(payoff)
call_greeks.loc['vega', 'Pathwise Abs err'] = np.abs(call_greeks.loc['vega', 'formula'] - call_greeks.loc['vega', 'Pathwise'])

# print('Analytical vega: ', greeks.vega())
# print('Pathwise method vega: ' , monte_carlo(payoff))
# print('abs difference:', np.abs(greeks.vega() - monte_carlo(payoff)))

# pathwise estimate for delta of European put option
dpayoff_dS0_put = - np.exp(-r * T)  * np.piecewise(ST, [ST > K, ST <= K], [0, 1])
payoff = dpayoff_dS0_put * ST / S0

put_greeks.loc['delta', 'Pathwise'] = monte_carlo(payoff)
put_greeks.loc['delta', 'Pathwise Abs err'] = np.abs(put_greeks.loc['delta', 'formula'] - put_greeks.loc['delta', 'Pathwise'])

# print('Analytical delta: ', greeks.delta('put'))
# print('Pathwise method delta: ' , monte_carlo(payoff))
# print('abs difference:', np.abs(greeks.delta('put') - monte_carlo(payoff)))

# pathwise estimate for vega for European put option
payoff = dpayoff_dS0_put * d1_similar * ST
put_greeks.loc['vega', 'Pathwise'] = monte_carlo(payoff)
put_greeks.loc['vega', 'Pathwise Abs err'] = np.abs(put_greeks.loc['vega', 'formula'] - put_greeks.loc['vega', 'Pathwise'])

# print('Analytical vega: ', greeks.vega())
# print('Pathwise method vega: ' , monte_carlo(payoff))
# print('abs difference:', np.abs(greeks.vega() - monte_carlo(payoff)))

In [15]:

call_greeks

Out[15]:

	formula	Finite Difference	FD Abs error (FD)	MC Greeks	MC abs error	Pathwise	Pathwise Abs err
delta	0.896455	0.896459	3.750064e-06	0.896991	0.000536	0.896830	0.000375
gamma	0.007500	0.007500	1.578147e-08	0.002801	0.004699	0.000000	0.000000
vega	21.600708	21.672831	7.212281e-02	21.858899	0.258191	21.515335	0.085374
rho	81.405559	81.418740	1.318111e-02	81.346387	0.059172	0.000000	0.000000
theta	-6.230349	-6.229884	4.648117e-04	-6.244217	0.013868	0.000000	0.000000

In [16]:

put_greeks

Out[16]:

	formula	Finite Difference	FD Abs error (FD)	MC Greeks	MC abs error	Pathwise	Pathwise Abs err
delta	-0.103545	-0.103541	3.750054e-06	-0.103460	0.000085	-0.103114	0.000431
gamma	0.007500	0.007500	1.045240e-08	0.002801	0.004699	0.000000	0.000000
vega	21.600708	21.672831	7.212281e-02	21.640966	0.040258	21.579024	0.021685
rho	-13.717384	-13.656657	6.072673e-02	-13.667977	0.049407	0.000000	0.000000
theta	-1.474202	-1.473856	3.459100e-04	-1.469397	0.004805	0.000000	0.000000

Recall that likelihood ratio method proposes that

In [19]:

# (iv) likelihood ratio method

S0 = 120
K = 100
r = 0.05
d = 0
sigma = 0.2
T = 1

num_steps = 100
num_paths = 17

ST = GBM_fd(S0, r, d, sigma, T, num_steps, 2**num_paths, plot = False, seed = None)
xi = (np.log( ST / S0 ) - (r - 0.5 * sigma**2) * T) / (S0 * sigma**2 * T)
payoff = np.exp(-r * T) * np.maximum(ST - K,0) * xi
print('Analytical delta: ', greeks.delta('call'))
print('Likelihood ratio method delta: ' , monte_carlo(payoff))
print('abs Difference: ' , np.abs(greeks.delta('call') - monte_carlo(payoff)))

Analytical delta:  0.8964550230770805
Likelihood ratio method delta:  0.8933323733097936
abs Difference:  0.0031226497672868225