Initially import all the modules we will be using for our notebook
import math
import numpy as np
import numpy.random as npr
# COMMAND PROMPT: pip install matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
import os
The payoff of a European call option on an index at maturity is given by $h(S_T)\equiv max\{S_T-K, 0\}$, where $S_T$ is the index level at maturity date $T$ and $K$ is the strike price. Given a risk-neutral measure for the relevant stochastic process (e.g., geometric Brownian motion), the price of such an option is given by the formula:
$$C_{0}=e^{-rT}\mathbb{E}_{0}^{Q}[h(S_{T})]=e^{-rT}\intop_{0}^{\infty}h(s)q(s)ds$$The equation below provides the respective Monte Carlo estimator for the European option, where $\tilde{S}_{T}^{i}$ is the $T$-th simulated index level at maturity.
$$\tilde{C}_{0}=e^{-rT}\frac{1}{I}\sum_{i=1}^{I}h(\tilde{S}_{T}^{i})$$Consider the following parameterization for the geometric Brownian motion and the valuation function gbm_mcs_stat()
, taking as a parameter only the strike price. Here, only the index level at maturity is simulated. As a reference, consider the case with a strike price of $K = 105$:
def bsm_call_value(S0, K, T, r, sigma):
''' Valuation of European call option in BSM model.
Analytical formula.
Parameters
==========
S0: float
initial stock/index level
K: float
strike price
T: float
maturity date (in year fractions)
r: float
constant risk-free short rate
sigma: float
volatility factor in diffusion term
Returns
=======
value: float
present value of the European call option
'''
from math import log, sqrt, exp
from scipy import stats
S0 = float(S0)
d1 = (log(S0 / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
d2 = (log(S0 / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sqrt(T))
# stats.norm.cdf --> cumulative distribution function
# for normal distribution
value = (S0 * stats.norm.cdf(d1, 0.0, 1.0) -
K * exp(-r * T) * stats.norm.cdf(d2, 0.0, 1.0))
return value
def gen_sn(M, I):
''' Function to generate random numbers for simulation.
Parameters
==========
M: int
number of time intervals for discretization
I: int
number of paths to be simulated
'''
sn = npr.standard_normal((M + 1, I))
return sn
S0 = 100.
r = 0.05
sigma = 0.25
T = 1.0
I = 50000
K = 105
def gbm_mcs_stat(K):
''' Valuation of European call option in Black-Scholes-Merton
by Monte Carlo simulation (of index level at maturity)
Parameters
==========
K: float
(positive) strike price of the option
Returns
=======
C0: float
estimated present value of European call option
'''
sn = gen_sn(1, I)
# Simulate Prices at Maturity
# calculate payoff at maturity
# calculate discounted price
return C0
# The Monte Carlo estimator value for the European call option.
gbm_mcs_stat(K)
9.962393419194543
# The analytical solution is
bsm_call_value(S0, K, T, r, sigma)
10.00220211715488
Next, consider the dynamic simulation approach and allow for European put options in addition to the call option. The function gbm_mcs_dyna()
implements the algorithm. The code also compares option price estimates for a call and a put stroke at the same level:
def gbm_mcs_dyna(K):
''' Valuation of European options in Black-Scholes-Merton
by Monte Carlo simulation (of index level paths)
Parameters
==========
K: float
(positive) strike price of the option
option : string
type of the option to be valued ('call', 'put')
Returns
=======
C0: float
estimated present value of European call option
'''
dt = T / M
# simulation of index level paths
S = np.zeros((M + 1, I))
S[0] = S0
sn = gen_sn(M, I)
for t in range(1, M + 1):
S[t] = S[t - 1] * np.exp((r - 0.5 * sigma ** 2) * dt
+ sigma * math.sqrt(dt) * sn[t])
# calculation of payoff
hT = np.maximum(S[-1] - K, 0)
# calculation of MCS estimator
C0 = math.exp(-r * T) * np.mean(hT)
return C0
# The Monte Carlo estimator value for the European call option (Dynamic Simulation Approach)
gbm_mcs_dyna(K)
10.088345691009375
Consider the stochastic differential equation $\mathrm{d} X_{t}=a\left(X_{t}\right) \mathrm{d} t+b\left(X_{t}\right) \mathrm{d} W_{t}$ with initial condition $X_{0}=x_{0},$ where $W_{t}$ stands for the Wiener process, and suppose that we wish to solve this SDE on some internal of time $[0, T]$. Then the Euler-Maruyama approximation to the true solution $X$ is the Markov chain $Y$ defined as follows:
# Calculate Price with these parameters again
S0 = 100.
r = 0.05
sigma = 0.25
T = 1.0
I = 5000
K = 105
10.012785056220235
(Next Time!)