#!/usr/bin/env python
# coding: utf-8

# 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)


# In[ ]:


# 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[ ]:


# 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


# In[5]:


put_greeks


# ## 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


# 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


# In[ ]:


# 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[ ]:


# (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')


# 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[ ]:


# (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()


# 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[ ]:


# (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')


# 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[ ]:


# (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()


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

# In[ ]:


# (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')


# In[ ]:


# (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()


# ## 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[ ]:


# (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


# In[13]:


put_greeks


# 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[ ]:


# (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


# In[16]:


put_greeks


# Recall that likelihood ratio method proposes 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[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)))