As we all know, the prices of options are highly correlated with stock prices, while the stock prices follow a geometric Brownian motion (GBM). In this group project, we investigate the pricing of different options (e.g. European Option, Asian Option, and American Option) under different Monte Carlo Methods. We first outline the prior knowledge of the Monte Carlo method and its improvement (e.g. naive Monte Carlo, antithetic Monte Carlo and control variate), and then proceed to price the three options.
For European option, since it has an analytical solution, we simply use a naive Monte Carlo to simulate its price and then compare it with the Black-Scholes formula.
For Asian option, since it is a strong path-dependent option, it can be classified into two types: discrete arithmetic average Asian option and discrete geometric average Asian option. For the latter, if the underlying asset price follows a log-normal distribution, because the geometric mean of a series of log-normal distribution remains log-normal distribution, which can be quickly drawn from the analytical solution of the geometric average Asian option pricing formula, so there is no need to use Monte Carlo to simulate. For the former, since the pricing formula can not be directly deduced using the Black-Scholes formula, we use three Monte Carlo simulation methods and their Quasi-versions (a total of six methods) to simulate.
For American option, since its prices do not have a general analytical form because of its early exercise feature, we must use some kind of numerical method. Unfortunately, general Monte Carlo method cannot be applied either because of its feature. So a workaround emerged: the least square Monte Carlo method (almost the current standard method for valuing American options using Monte Carlo simulation). We use three least square Monte Carlo simulation methods and their Quasi-versions (a total of six methods) to simulate.
Monte Carlo is an extremely important type of numerical method in finance. It can also be called a computer simulation method. Its birth is based on the theory of probability, mathematical statistics, the central limit theorem and the law of large numbers as the main theoretical basis.
The theory of Monte Carlo can be roughly expressed as: Establish a probability model related to solving the problem, so that the value of the required solution or its function can be expressed as the mathematical expectation of the built model. Carry out a large number of sampling data simulations on the model, and use the arithmetic mean of random variables generated by sampling as the approximate estimated value of the obtained solution. Therefore, it can be simply considered that the Monte Carlo algorithm calculates integrals by means of random experiments, that is, the mathematical expectation that the integral to be calculated is regarded as a random variable subject to a special distribution density function.
In this article, we use Monte Carlo to estimate the price of three kinds of options.
We assume that the quantity $\mu$ we need to calculated is the mathematical expectation of a random variable $\eta$, $\mu = \mathbb{E}[\eta]$. Then the Monte Carlo for estimating $\mu$ is to perform repeated sampling on to generate independent and identically distributed random variable sequences $\eta_1, \eta_2, ……, \eta_n$, and then calculate the sample mean $\bar{\eta} = \frac{1}{n} \sum_{k=1}^n \eta_k$ Accroding to the law of large numbers, we know that $\text{Pr}(\lim_{n\to\infty} \bar{\eta_n} = \mu) = 1$.
Assume that the variance of $\eta$, $\text{Var}[\eta] = \sigma^2$, then $$ \text{Pr}\left(|\eta_n-\mu| < Z_{\frac{\sigma}{2}} \frac{\sigma}{\sqrt n}\right) \approx \frac{1}{\sqrt{2\pi}} \int^{Z_{\frac{\sigma}{2}}}_{-Z_{\frac{\sigma}{2}}} e^{- \frac{t^2}{2}} dt = 1-\sigma. $$
So the error of Monte Carlo is $Z_{\frac{\sigma}{2}} \frac{\sigma}{\sqrt n}$, the error convergence rate is $\mathcal{O}(n^{-\frac{1}{2}})$.
From 1.2. we can know the error is determined by $\sigma$ and $n$. Under the premise of sampling the same $\eta$, $\sigma$ is fixed. So if we want to improve accuracy, we have to increse $n$, but this comes at the expense of speed. We want to explore ways that can improve accuracy without sacrificing speed.
Antithetic Variates Method uses the negative correlation between variables. For example, for $Y = \int^1_0 f(y)dy$, we can use Monte Carlo to get two estimated values of $Y$ : $$ Y_a = \frac{1}{N} \sum^N_{i=1} f(u_i) = \frac{1}{N} \sum^N_{i=1} Y_i, $$ $$ Y_b = \frac{1}{N} \sum^N_{i=1} f(1-u_i) = \frac{1}{N} \sum^N_{i=1} Y_i', $$ where $u_i \in [0,1)$ is a random number generated uniformly on the interval $[0,1)$.
Since the two estimated values $Y_a$, $Y_b$ are related, the two sets of observations obtained should also be related theoretically. So we call them Antithetic Variates. The average value of them, $Y_{AV}$, can be given as $$ Y_{AV} = \frac{Y_a +Y_b}{2} = \frac{1}{N} \sum^n_{i=1} \frac{f(u_i)+f(1-u_i)}{2}= \frac{1}{N} \sum^n_{i=1} \frac{Y_i+ Y_i'}{2}, $$
$$\renewcommand{\Var}{\text{Var}} \Var(Y_{AV})= \Var\left(\frac{Y_i+ {Y_i}'}{2}\right) = \sigma_{AV}^2.$$Acccroding to the central limit theorem, as $N$ goes large enough, $$ \frac{Y_{AV}-\mathbb{E}[Y]}{\frac{\sigma_{AV}}{\sqrt N}} \backsim \mathcal{N}(0,1), $$ where $\sigma_{AV} = \sqrt{\frac{1}{N} \sum^n_{i=1}(\bar{Y_i} - Y_{AV})^2}$.
So the ratio of variance under Antithetic Variates Method and variance under the basic Monte Carlo is $$ \frac{\frac{1}{2} \Var(Y_i)}{\Var(\bar{Y_i})} = \frac{\frac{1}{2} \Var(Y_i)}{\Var\left(\frac{Y_i + {Y_i}'}{2}\right)} = \frac{\frac{1}{2} \Var(Y_i)}{\frac{1}{2} \Var(Y_i)(1+\rho)} = \frac{1}{1+\rho}. $$
It's easy to find that, only when $\rho <0$ can the variance be reduced compared to basic Monte Carlo.
The concept of control variates is that in order to estimate a quantity from data, the information about another can be used to adjust the original estimator.
Suppose we want to estimate some mean value $\mathbb{E}(h(X)) = \mu_h$ and this has been done by using a basic Monte Carlo estimator $\bar{X}_h$. Suppose that the same data set, {$X_i, 1 ≤ i ≤ n$}, can be used to estimate a known quantity, so that this sample can be used to adjust the naive estimator because the actual estimation error can be computed. If we let $\mathbb{E}(g(X)) = \mu_g$ which is known, at the same time, we can also use the basic Monte Carlo estimator $\bar{F_g}$. Then $\mathbb{E}(\bar{X_g} - \mathbb{E}(g(X))) =0$. If $\bar{X_g}$ is correlated with $\bar{X}_h$, we can improve the estimator of $\mu_h$ by using $$ \bar{X_h} + \hat{a} (\bar{X_g} - \mathbb{E}(g(X))), $$
where $\hat{a}$ should be chosen to minimize the mean square error:
$$\renewcommand{\Cov}{\text{Cov}} \Var[\bar{X_h} + \hat{a} (\bar{X_g} - \mathbb{E}(g(X)))] = \Var(\bar{X_h}) + \hat a ^2 \Var(\bar{X_g}) + 2\hat a \Cov(\bar{X_h}, \bar{X_g}). $$In order to minimize the variance $$ \left.\frac{d \Var[\bar{X_h} + \hat{a} (\bar{X_g} - \mathbb{E}(g(X)))]}{d \hat a} \right|_{a^*}= 2a^* \Var(\bar{X_g}) + 2\Cov(\bar{X_h}, \bar{X_g})=0, $$ $$ a^*= - \frac{\Cov(\bar{X_h}, \bar{X_g})}{\Var(\bar{X_g})}. $$
When the classic Monte Carlo algorithm takes random numbers, most of them are implemented using the randn function. Although the random numbers given in this function meet most random requirements, there are still some generated random numbers that cannot be ignored. Large deviations may have strong randomness and uneven distribution of numbers. Therefore, in order to avoid this phenomenon, it is necessary to improve the Monte Carlo algorithm. We introduce a new concept called low deviation sequence. The main function of the sequence is to make the generated random numbers uniform. The higher the uniformity of the random numbers is, the smaller the deviation will be, and the more obvious the effect of the improvement of the classic Monte Carlo algorithm. To put it simply, the quasi-Monte Carlo algorithm is to simulate the price of derivative securities with the method of low discrepancy sequence, which guarantees a small error from the source.
In this article, two low discrepancy sequences we will use is called the Faure sequence[1] and Halton sequence[2], which seem to have become more and more popular in mathematical finance simulations. To further deal with the these sequences (the numbers in which lie in the range $(0,1)$ ), we will use Moro's inverse transtormation [3] to convert it into a normally distributed sequence. The basic idea of this transformation is to separately convert the head and tail of the input sequence, and its implementation requires only a few constants. The code can be seen below:
def MoroInvCND(probs):
"""Function to convert a low discrepancy sequence into a normal distributed sequence.
input : probs , 1 * n , dtype=np.float64
output : zs , 1 * n , dtype=np.float64
"""
zs = np.zeros_like(probs, dtype=np.float64)
# Define the constants in the transformation process.
a = [2.50662823884, -18.61500062529, 41.39119773534, -25.44106049637]
b = [1, -8.47351093090, 23.08336743743, -21.06224101826, 3.13082909833]
c = [7.7108870705487895, 2.7772013533685169, 0.3614964129261002, 0.0373418233434554,
0.0028297143036967, 0.0001625716917922, 0.0000080173304740, 0.0000003840919865, 0.0000000129707170]
k1 = 0.4179886424926431
k2 = 4.2454686881376569
s = len(probs)
probs = probs-0.5
for i in range(s):
if probs[i] <= 0.42 and probs[i] >= -0.42:
zs[i] = probs[i]*((a[0]+a[1]*(probs[i]**2)+a[2]*(probs[i]**4)+a[3]*(probs[i]**6))/(
b[0]+b[1]*(probs[i]**2)+b[2]*(probs[i]**4)+b[3]*(probs[i]**6)+b[4]*(probs[i]**8)))
else:
z = k1*(2*np.log(-np.log(0.5-np.abs(probs[i])))-k2)
d = [0]*11
for j in range(8, 0, -1):
d[j] = 2*z*d[j+1]-d[j+2]+c[j]
f = z*d[1]-d[2]+c[0]/2
if probs[i] < 0.5 and probs[i] > 0.41:
zs[i] = f
else:
zs[i] = -f
return zs
Here, we import all the packages used in this group project. Meanwhile, we use the seaborn package to set up our plot style.
import numpy as np
from statsmodels.tools import sequences
import openturns as ot
import matplotlib.pyplot as plt
import seaborn as sns
import math
from scipy import stats
sns.set(style='ticks', palette='Set2')
Among these packages, OpenTURNS may not be widely used, so you can consider installing it through
%pip install openturns
A European option is a version of an options contract that limits execution to its expiration date. In other words, if the underlying security such as a stock has moved in price, an investor would not be able to exercise the option early and take delivery of or sell the shares. Instead, the call or put action will only take place on the date of expiry.
We learned it in class, so here we just give out the final formula
For call
$$ C = S e^{-d\tau} N(d_1) - K e^{-r\tau} N(d_2) $$For put
$$ P = K e^{-r\tau} N(-d_2) - S e^{-d\tau} N(-d_1). $$where $\tau$ is time to expiry, $N(\cdot)$ denotes the cdf for standard normal, and
$$ d_1 = \frac{\log\left(\frac{Se^{-d\tau}}{Ke^{-r\tau}}\right)}{\sigma\sqrt\tau}+ \frac{\sigma\sqrt\tau}2, \qquad d_2 = d_1 - \sigma \sqrt\tau $$class BSM:
def __init__(self, type, S0, K, T, r, sigma):
try:
assert isinstance(type, str)
self.type = type
self.S0 = float(S0)
self.K = float(K)
assert T > 0
self.T = float(T)
assert r >= 0
self.r = float(r)
assert sigma >= 0
self.sigma = float(sigma)
except ValueError:
print("Please pass correct parameters")
if type != 'call' and type != 'put':
raise ValueError("Please use 'call' or 'put' as option type.")
@property
def price(self):
d1 = (np.log(self.S0 / self.K) + (self.r + 0.5 * self.sigma**2)
* self.T)/(self.sigma * np.sqrt(self.T))
d2 = (np.log(self.S0 / self.K) + (self.r - 0.5 * self.sigma**2)
* self.T)/(self.sigma * np.sqrt(self.T))
def N(d): return stats.norm.cdf(d, 0, 1)
if self.type == 'call':
value = self.S0 * N(d1) - self.K * np.exp(-self.r * self.T) * N(d2)
else:
value = self.S0 * N(d1) - self.K * np.exp(-self.r * self.T) * \
N(d2) - self.S0+self.K*np.exp(-self.r*self.T)
return value
Based on the risk-neutral option pricing principle, any asset is risk-neutral to the holder under the risk-neutral probability measurement. The risk-neutral probability can be used to obtain the expected return of the option and then it will be risk-free Discount, so that the value of the option at the initial moment is obtained. It can be seen that calculating the price of an option is to calculate an expected value, and the Monte Carlo algorithm is used to estimate the expected value. Generally, the Monte Carlo algorithm for option pricing is divided into the following steps:
Under the premise of the risk-neutral measurement, simulate the path of the underlying asset price of the option, divide the time interval $[0, T]$ into n sub-intervals $0=t_0<t_1<……<t_n=T$, and the discrete form of the underlying asset price process is $$ S^j (t_{i+1}) = S^j (t_i) e^{(r- \frac{1}{2}\sigma^2)(t_{i+1}-t_i) + \sigma \sqrt{t_{i+1}-t_i} Z_i},\quad Z_i\sim \mathcal N(0,1). $$
Calculate the return to maturity of the option under this path, and try to obtain the discount of the return to maturity of the option based on the risk-free interest rate $$ C^j = e^{-rT} \max(0, S^j_T -K). $$
Repeating the previous two steps, a large number of sample samples of the discount value of option returns can be obtained.
Calculate the sample average value, and then obtain the Monte Carlo simulation value of the option price $$ C_{MC} = \frac{1}{m} \sum_{j=1}^m C^j. $$
class EuropeanOptionsPricing:
""" Class for American Options pricing.
type : str , 'put' or 'call'
S0 : float , initial stock price
K : float, strike price
T : float , time to expire (in years)
num_T : int , grid for time (number of total grid points)
r : float, the riskless discount rate
sigma : float , volatility
num_sim : int, number of simulations
MCtype : str , here we achieved antithetic variates method and control variates method (plus the Normal Monte Carlo method)
isquasi : bool , True for quasi-Monte Carlo
"""
def __init__(self, type, S0, K, T, num_T, r, sigma, num_sim, MCtype='Antithetic', isquasi=False):
try:
assert isinstance(type, str)
self.type = type
self.S0 = float(S0)
self.K = float(K)
assert T > 0
self.T = float(T)
assert isinstance(num_T, int)
assert num_T > 0
self.num_T = num_T
assert r >= 0
self.r = float(r)
assert sigma >= 0
self.sigma = float(sigma)
assert isinstance(num_sim, int)
assert num_sim > 0
self.num_sim = num_sim
assert MCtype in ['Antithetic', 'Control', 'Normal']
self.MCtype = MCtype
assert isinstance(isquasi, bool)
self.isquasi = isquasi
except ValueError:
print("Please pass correct parameters")
if type != 'call' and type != 'put':
raise ValueError("Please use 'call' or 'put' as option type.")
self.dT = self.T/self.num_T
self.discount = np.exp(-self.r*self.dT)
if self.MCtype == 'Antithetic':
self.num_sim = self.num_sim*2
@property
def MC_price_matrix(self):
""" Return (quasi) Monte Carlo price matrix
Rows: time
Columns: price
"""
np.random.seed(123)
mtx = np.zeros((self.num_T+1, self.num_sim), dtype=np.float64)
mtx[0, :] = self.S0
if not self.isquasi:
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
brownian = np.random.standard_normal(self.num_sim)
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
brownian = np.random.standard_normal(self.num_sim//2)
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
else:
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
faure_seqs = np.array(ot.HaltonSequence(
self.num_T).generate(self.num_sim)).T
faure_seq = faure_seqs[t-1]
brownian = MoroInvCND(faure_seq)
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
faure_seqs = np.array(ot.HaltonSequence(
self.num_T).generate(self.num_sim//2)).T
faure_seq = faure_seqs[t-1]
brownian = MoroInvCND(faure_seq)
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
return mtx
@property
def MC_payoff(self):
"""Return the inner value of the option"""
if self.type == 'call':
payoff = np.maximum(self.MC_price_matrix-self.K,
np.zeros((self.num_T+1, self.num_sim), dtype=np.float64))
else:
payoff = np.maximum(self.K-self.MC_price_matrix,
np.zeros((self.num_T+1, self.num_sim), dtype=np.float64))
return payoff
@property
def value_vector(self):
value_mtx = self.MC_payoff
return value_mtx[-1, :]
@property
def controlled_value_vector(self):
control_variables = self.MC_price_matrix[self.num_T, :]
beta = np.cov(control_variables, self.value_vector)[
0][1]/np.var(self.value_vector)
return self.value_vector - beta * (control_variables - self.S0*np.exp(self.r * self.T))
@property
def price(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return np.sum(self.value_vector)/self.num_sim
else:
return np.sum(self.controlled_value_vector)/self.num_sim
@property
def stddev(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return np.std(self.value_vector)
else:
return np.std(self.controlled_value_vector)
@property
def delta(self):
dS0 = self.S0/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0+dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0-dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dS0)
@property
def gamma(self):
dS0 = self.S0/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0+dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0-dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_1.delta-approx_option_2.delta)/(2*dS0)
@property
def vega(self):
dsigma = self.sigma/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r, self.sigma+dsigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r, self.sigma-dsigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dsigma)
@property
def rho(self):
dr = self.r/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r+dr, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r-dr, self.sigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dr)
@property
def theta(self):
dtime = 1/252
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T+dtime, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T-dtime, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_2.price-approx_option_1.price)/(2*dtime)
simulation_times = np.arange(500, 10000, 500)
instances = {'Normal': [], 'ControlVariate': [], 'AntitheticVariate': [
], 'QuasiMC&Normal': [], 'QuasiMC&ControlVariate': [], 'QuasiMC&AntitheticVariate': []}
for s in simulation_times:
instances['Normal'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Normal'))
instances['ControlVariate'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Control'))
instances['AntitheticVariate'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s)))
instances['QuasiMC&Normal'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Normal', isquasi=True))
instances['QuasiMC&ControlVariate'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Control', isquasi=True))
instances['QuasiMC&AntitheticVariate'].append(EuropeanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), isquasi=True))
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
plt.plot(simulation_times, [BSM('call', 36, 30, 1, 0.06,
0.2).price]*len(simulation_times), '--', color='grey')
for method in instances:
plt.plot(simulation_times, [
instance.price for instance in instances[method]], label=method)
plt.legend()
plt.xlabel('Times of simulation')
plt.ylabel('Price')
plt.show()
From the figure above, it's clear that all methods tend to converge as $N$ goes large enough.
QuasiMC&ControlVariate method, QuasiMC&AntitheticVariate method, ControlVariate methord and AntitheticVariate methord are all more effective than normal Monte Carlo.
Moreover, we can examine the standard deviation of the samples:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
stddevs = []
for method in instances:
stddevs.append(instances[method][-1].stddev)
plt.bar(instances.keys(), stddevs, color=[
'#66c2a5', '#fc8d62', '#8da0cb', '#e78ac3', '#a6d854', '#ffd92f'])
plt.ylabel('Stddev')
plt.xticks(rotation=30)
plt.show()
From this figure, we can see that ControlVariate method andQuasiMC&ControlVariate method have the smallest Standard deviation.
In conclusion, the Monte Carlo method is a powerful pricing method for discrete arithmetic average European options, and the effect of ControlVariate (QuasiMc) method is best.
Asian option is a special type of option whose payoff is determined by the average underlying price over some pre-set period of time. This is different from the case of the usual European option and American option, where the payoff of the option contract depends on the price of the underlying asset at exercise. There are two ways to calculate the average price: geometric average and arithmetic average. One advantage of Asian options is that these reduce the risk of market manipulation of the underlying instrument at maturity. Another advantage of Asian options involves the relative cost of Asian options compared to European or American options. Because of the averaging feature, Asian options reduce the volatility inherent in the option; therefore, Asian options are typically cheaper than European or American options.
If the assets price obeys the lognormal distribution, the geometric mean of multiple lognormal distributions also obeys the lognormal distribution.
Divide the time interval $[0, T]$ into n sub-intervals $0=t_0<t_1<\cdots<t_n=T$, assume $S_1, S_2,\cdots, S_n$ are price of assets at $t_1, t_2,\cdots, t_n$, then the geometric mean of assets price is $$ I^G = \left(\prod_{i=1}^n S_i \right)^{\frac{1}{n}}. $$
Let $Y = \frac{\sum_{i=1}^n \ln S_i}{n}$， we can know that $Y$ obeys normal distribution:
$$\begin{align*} \mathbb{E}(Y) &= \mathbb{E}\left(\frac{\sum_{i=1}^n \ln S_i}{n}\right) \\ &= \frac{n\ln S_0+(1+2+\cdots+n)\left(r-\frac{\sigma^2}{2}\right)t}{n}\\ &= \ln S_0 + \frac{n+1}{2n}\left(r-\frac{\sigma^2}{2}\right)T,\\ \Var(Y) &= \frac{1}{n^2}\Var(\ln S_1+\ln S_2+\cdots+\ln S_n)\\ &= \frac{1}{n^2}\Var(\ln(S_1 S_2……S_n))\\ &= \frac{1}{n^2}\Var\left(\ln \frac{S_n}{S_{n-1}}+2\ln \frac{S_{n-1}}{S_{n-2}}+3\ln \frac{S_{n-2}}{S_{n-3}}+\cdots+n\ln \frac{S_1}{S_0}\right)+n\ln S_0\\ &=\frac{(n+1)(2n+1)}{6n^2} \sigma^2 T. \end{align*} $$As a result $$ \mathbb{E}(I^G) = e^{\ln S_0+(\mu_G - \frac{\sigma^2}{2})T + \frac{\sigma^2 T}{2}} = S_0e^{\mu_G T}. $$ Hence, the price of Geometric Average Asian Call Option is $$ \begin{aligned} C(S_0, K, T) &= e^{-rT} \mathbb{E}(I^G - K)^+\\ &=S_0 e^{(\mu_G -r)T}N(d_1) - Ke^{-rT}N(d_2), \end{aligned} $$ where $d_1 = \frac{\ln \frac{S_0}{K}+ \left(\mu_G+\frac{\sigma^2}{2}\right)T}{\sigma \sqrt T}$, $d_2 = d_1 - \sigma\sqrt T$.
Because Geometric Average Asian Option have a definite mathematical solution, there is no need to use Monte Carlo to simulate. Unfortunately, Geometric Average Asian Option is rarely used in real life, and we use Arithmetic Average Asian Options more often. Now let's turn to Arithmetic Average Asian Options.
class AsianOptionsPricing:
""" Class for American Options pricing.
S0 : float , initial stock price
K : float, strike price
T : float , time to expire (in years)
num_T : int , grid for time (number of total grid points)
r : float, the riskless discount rate
sigma : float , volatility
num_sim : int, number of simulations
MCtype : str , here we achieved antithetic variates method and control variates method (plus the Normal Monte Carlo method)
isquasi : bool , True for quasi-Monte Carlo
"""
def __init__(self, type, S0, K, T, num_T, r, sigma, num_sim, MCtype='Antithetic', isquasi=False):
try:
assert isinstance(type, str)
self.type = type
self.S0 = float(S0)
self.K = float(K)
assert T > 0
self.T = float(T)
assert isinstance(num_T, int)
assert num_T > 0
self.num_T = num_T
assert r >= 0
self.r = float(r)
assert sigma >= 0
self.sigma = float(sigma)
assert isinstance(num_sim, int)
assert num_sim > 0
self.num_sim = num_sim
assert MCtype in ['Antithetic', 'Control', 'Normal']
self.MCtype = MCtype
assert isinstance(isquasi, bool)
self.isquasi = isquasi
except ValueError:
print("Please pass correct parameters")
if type != 'call' and type != 'put':
raise ValueError("Please use 'call' or 'put' as option type.")
self.dT = self.T/self.num_T
self.discount = np.exp(-self.r*self.T)
if self.MCtype == 'Antithetic':
self.num_sim = self.num_sim*2
@property
def MC_price_matrix(self):
""" Return (quasi) Monte Carlo price matrix
Rows: time
Columns: price
"""
np.random.seed(1234)
mtx = np.zeros((self.num_T+1, self.num_sim), dtype=np.float64)
mtx[0, :] = self.S0
if not self.isquasi:
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
brownian = np.random.standard_normal(self.num_sim)
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
brownian = np.random.standard_normal(self.num_sim//2)
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
else:
num_sim_ = self.num_sim//(2 if self.MCtype == 'Antithetic' else 1)
faure_seqs = np.array(ot.FaureSequence(
self.num_T).generate(num_sim_)).T
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
brownian = MoroInvCND(faure_seqs[t-1])
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
brownian = MoroInvCND(faure_seqs[t-1])
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
return mtx
@property
def value_vector(self):
"""Return the C_i of i_th simulated option"""
payoff = self.MC_price_matrix[1:, :].sum(axis=0)
payoff = payoff / self.num_T
if self.type == 'call':
payoff = np.maximum(payoff-self.K,
np.zeros(self.num_sim, dtype=np.float64))
else:
payoff = np.maximum(self.K-payoff,
np.zeros(self.num_sim, dtype=np.float64))
return payoff
# (num_sum)-dimensional vector
@property
def E_CX(self): # expectation of the control variable CX=\sum{S_i}
E_CX_ = 0.0
for i in range(self.num_T):
E_CX_ = E_CX_ + self.S0 * \
math.exp((self.r-self.sigma**2/2)*(i+1)*self.dT)
return E_CX_
@property
def controlled_value_vector(self):
"""Return the CX (the control variable) of i_th simulated option"""
control_variables = self.MC_price_matrix[1:, :].sum(axis=0)
beta = np.cov(control_variables, self.value_vector)[
0][1]/np.cov(control_variables)
return self.value_vector - beta * (control_variables - self.E_CX)
# (num_sum)-dimensional vector
@property
def price(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return self.value_vector.mean() * self.discount
else:
return self.controlled_value_vector.mean() * self.discount
@property
def stddev(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return np.std(self.value_vector*self.discount)
else:
return np.std(self.controlled_value_vector*self.discount)
simulation_times = np.arange(100, 5000, 500)
instances = {'Normal': [], 'ControlVariate': [], 'AntitheticVariate': [
], 'QuasiMC&Normal': [], 'QuasiMC&ControlVariate': [], 'QuasiMC&AntitheticVariate': []}
for s in simulation_times:
instances['Normal'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Normal'))
instances['ControlVariate'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Control'))
instances['AntitheticVariate'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Antithetic'))
instances['QuasiMC&Normal'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Normal', isquasi=True))
instances['QuasiMC&ControlVariate'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Control', isquasi=True))
instances['QuasiMC&AntitheticVariate'].append(AsianOptionsPricing(type='call', S0=100, K=105, T=1, num_T=16,
r=0.15, sigma=0.05, num_sim=int(s), MCtype='Antithetic', isquasi=True))
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
for method in instances:
plt.plot(simulation_times, [
instance.price for instance in instances[method]], label=method)
plt.legend()
plt.xlabel('Times of simulation')
plt.ylabel('Price')
plt.show()
From the figure above, it's clear that all methods tend to converge as $N$ goes large enough.
QuasiMC&ControlVariate methods, ControlVariate methord and AntitheticVariate methord are all more effective than normal Monte Carlo.
Moreover, we can examine the standard deviation of the samples:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
stddevs = []
for method in instances:
stddevs.append(instances[method][-1].stddev)
plt.bar(instances.keys(), stddevs, color=[
'#66c2a5', '#fc8d62', '#8da0cb', '#e78ac3', '#a6d854', '#ffd92f'])
plt.ylabel('Stddev')
plt.xticks(rotation=30)
plt.show()
From this figure, we can see that ControlVariate method andQuasiMC&ControlVariate method have the smallest Standard deviation.
In conclusion, the Monte Carlo method is a powerful pricing method for discrete arithmetic average Asian options, and the effect of ControlVariate (QuasiMc) method is best.
An American option is an important kind of option that can be exercised on any day during its lifetime. That is, the option holder can choose to exercise or not to exercise the option on any business day prior to the option expiration date. American option is usually divided into American call option and American put option. During the exercise period of American call option, if the underlying asset $S$ is larger than the strike price $K$, the holder can buy the underlying asset at the strike price K at any time. At this time, the holder buys the underlying asset at price $K$ and sells it in the market at the instant price $S$, gaining profit $S-K$; if the price of the underlying asset is always lower than the strike price $K$ during the exercise period, the holder can choose not to buy any underlying asset, and the profit gained at this time is zero. Similar logic can be applied to the American put option. From the definition, we can see that American options are similar to European options in all aspects except for the execution time.
Unfortunately, unlike the European Options whose could be easily determined through the classical Black-Scholes equation, the prices of the more widely used American options do not have a general analytical form because of its early exercise feature. The American option pricing problem is actually an optimal stopping time problem, that's to say, the key is to find the moment when the return from exercising the option is greater than the expected return from continuing to hold the option, Therefore, some kind of numerical method must be applied. However, since the Monte Carlo method uses a forward solving approach, we cannot calculate the expected return of continuing to hold the option at each moment, and thus cannot decide whether to execute the option immediately or to continue to hold the option. So a workaround emerged: the least square Monte Carlo method (almost the current standard method for valuing American options using Monte Carlo simulation).
The basic idea of the Least-Square Monte Carlo algorithm (LSMC), first published in 2001[4], as its name implies, is to use least-square regression on a set of basis functions to estimate the price of the American Option. As mentioned before, at any time before the expiry date, the holder of the American Option compares the return of immediate exercise to the continuation value, and if the former one exceeds the latter one, the option will be exercised. Hence, the optimal strategy depends on the conditional expectation of continuing to hold the option, and such conditional expectation can be estimated from cross-sectional information through least-square algorithm. Technically, we can do the regression of the return obtained ex-post from the continuation to the function of state variable values, through which we can further get the estimation of the conditional expectation for each exercise time.
To set up the basic valuation framework, we assume a complete probability space $(\Omega, \mathcal{F}_T, P)$ and time horizon $[0,T]$, where $\Omega$ is the set of all possible states with sample path denoted by $\omega$, $\mathcal{F}_T$ is the sigma field of distinguishable events at time $T$, and $P$ is the probability measure defined on the elements of $\mathcal F_T$. Considering the no-arbitage paradigm, we can further introduce the risk-neutral pricing measure $Q$. Then. we define the path of cash flow $C(\omega,s;t,T)$ generated by the option, where the option is not exercised at or prior to time $t$ and $t < s \le T$.
In order to make the algorithm practically available, we have to discretize the time axis: i.e., exercise is only permitted at discrete times $0 < t_1 < t_2 < \cdots < t_K =T$, and at every allowed date, the holder should reconsider the stratagy. Obviously, as $K$ goes larger, the simulation will be more effective. With risk-neutral measure $Q$, the value of continuation $F(\omega;t_k)$ at time $t_k$ can be expressed as
$$F(\omega; t_k)=\mathbb{E}_Q\left[\sum_{j=k+1}^K \exp\left(-\int_{t_k}^{t_j}r(\omega,s) ds\right)C(\omega,t_j;t_k,T)\middle| \mathcal{F}_{t_k}\right].$$Hence, the pricing problem is converted into comparing the immediate exercise value and the conditional expectation of the value of continuation and such a conditional expectation can be approximated by a linear combination of some basis functions.
The cash flow at the expiry date $t_K=T$, is easy to get, since there are no more continuation values, or in other words, it is just the payoff of an option in the terminal value. Then, we consider the case at time $t_{K-1}$, the exercise value can be simply calculated, and from the assumption above, the continuation value can be expressed as $$F(\omega; t_{K-1})=\sum_{j=0}^{\infty} a_j (t_K)B_j (S_{t_K}),$$ where $S$ is the spot price, $a_j$'s are coefficients and $B_j$ are basis functions.
Similar logic can be applied to other $t_k$'s, and whenever the exercise value exceeds the condition value, the option will be exercised. Following this progress, we can step backward until we reach $t_1$.
Here we propose a Python Class for Amecican Options pricing following the ideas above, where the riskless discount rate is assumed to be constant out of simplicity. Practically, we use the polynomial fit method for linear square regression, and the stock prices at the time of expiry date $S_T$ as the control variates (obviously, its expectation equals to $S_0 \mathrm{e}^{r T}$). The code can be seen below:
class AmericanOptionsPricing:
""" Class for American Options pricing.
type : str , 'put' or 'call'
S0 : float , initial stock price
K : float, strike price
T : float , time to expire (in years)
num_T : int , grid for time (number of total grid points)
r : float, the riskless discount rate
sigma : float , volatility
num_sim : int, number of simulations
MCtype : str , here we achieved antithetic variates method and control variates method (plus the Normal Monte Carlo method)
isquasi : bool , True for quasi-Monte Carlo
"""
def __init__(self, type, S0, K, T, num_T, r, sigma, num_sim, MCtype='Antithetic', isquasi=False):
try:
assert isinstance(type, str)
self.type = type
self.S0 = float(S0)
self.K = float(K)
assert T > 0
self.T = float(T)
assert isinstance(num_T, int)
assert num_T > 0
self.num_T = num_T
assert r >= 0
self.r = float(r)
assert sigma >= 0
self.sigma = float(sigma)
assert isinstance(num_sim, int)
assert num_sim > 0
self.num_sim = num_sim
assert MCtype in ['Antithetic', 'Control', 'Normal']
self.MCtype = MCtype
assert isinstance(isquasi, bool)
self.isquasi = isquasi
except ValueError:
print("Please pass correct parameters")
if type != 'call' and type != 'put':
raise ValueError("Please use 'call' or 'put' as option type.")
self.dT = self.T/self.num_T
# different from the Asian Option, here self.discount denotes the discount between t and t+dT.
self.discount = np.exp(-self.r*self.dT)
if self.MCtype == 'Antithetic':
self.num_sim = self.num_sim*2
@property
def MC_price_matrix(self):
""" Return (quasi) Monte Carlo price matrix
Rows: time
Columns: price
"""
np.random.seed(123)
mtx = np.zeros((self.num_T+1, self.num_sim), dtype=np.float64)
mtx[0, :] = self.S0
if not self.isquasi:
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
brownian = np.random.standard_normal(self.num_sim)
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
brownian = np.random.standard_normal(self.num_sim//2)
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
else:
if self.MCtype == 'Normal' or self.MCtype == 'Control':
for t in range(1, self.num_T + 1):
faure_seqs = np.array(ot.FaureSequence(
self.num_T).generate(self.num_sim)).T
faure_seq = faure_seqs[t-1]
brownian = MoroInvCND(faure_seq)
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
elif self.MCtype == 'Antithetic':
for t in range(1, self.num_T+1):
faure_seqs = np.array(ot.FaureSequence(
self.num_T).generate(self.num_sim//2)).T
faure_seq = faure_seqs[t-1]
brownian = MoroInvCND(faure_seq)
brownian = np.concatenate((brownian, -brownian))
mtx[t, :] = (mtx[t-1, :]*np.exp((self.r-self.sigma**2/2)
* self.dT+self.sigma*brownian*np.sqrt(self.dT)))
return mtx
@property
def MC_payoff(self):
"""Return the inner value of the option"""
if self.type == 'call':
payoff = np.maximum(self.MC_price_matrix-self.K,
np.zeros((self.num_T+1, self.num_sim), dtype=np.float64))
else:
payoff = np.maximum(self.K-self.MC_price_matrix,
np.zeros((self.num_T+1, self.num_sim), dtype=np.float64))
return payoff
@property
def value_vector(self):
value_mtx = np.zeros_like(self.MC_payoff)
value_mtx[-1, :] = self.MC_payoff[-1, :]
for t in range(self.num_T-1, 0, -1):
reg = np.polyfit(
self.MC_price_matrix[t, :], value_mtx[t+1, :]*self.discount, 5)
continuation_values = np.polyval(reg, self.MC_price_matrix[t, :])
value_mtx[t, :] = np.where(
self.MC_payoff[t, :] > continuation_values, self.MC_payoff[t, :], value_mtx[t+1, :]*self.discount)
return value_mtx[1, :]*self.discount
@property
def controlled_value_vector(self):
control_variables = self.MC_price_matrix[self.num_T, :]
beta = np.cov(control_variables, self.value_vector)[
0][1]/np.var(self.value_vector)
return self.value_vector - beta * (control_variables - self.S0*np.exp(self.r * self.T))
@property
def price(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return np.sum(self.value_vector)/self.num_sim
else:
return np.sum(self.controlled_value_vector)/self.num_sim
@property
def stddev(self):
if self.MCtype == 'Antithetic' or self.MCtype == 'Normal':
return np.std(self.value_vector)
else:
return np.std(self.controlled_value_vector)
@property
def delta(self):
dS0 = self.S0/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0+dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0-dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dS0)
@property
def gamma(self):
dS0 = self.S0/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0+dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0-dS0, self.K, self.T, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_1.delta-approx_option_2.delta)/(2*dS0)
@property
def vega(self):
dsigma = self.sigma/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r, self.sigma+dsigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r, self.sigma-dsigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dsigma)
@property
def rho(self):
dr = self.r/100
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r+dr, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T, self.num_T, self.r-dr, self.sigma, self.num_sim)
return (approx_option_1.price-approx_option_2.price)/(2*dr)
@property
def theta(self):
dtime = 1/252
approx_option_1 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T+dtime, self.num_T, self.r, self.sigma, self.num_sim)
approx_option_2 = AmericanOptionsPricing(
self.type, self.S0, self.K, self.T-dtime, self.num_T, self.r, self.sigma, self.num_sim)
return (approx_option_2.price-approx_option_1.price)/(2*dtime)
To have a better understanding of the differences between the methods above, we first compare the option prices predicted by them. Here we choose a hypothetical call option with an initial price of 36, an strike price of 30, a risk-free rate of 0.06, an expiry time of 1 year (divided into 16 parts), and a volatility of 0.05 for the underlying asset. We start with generating some instances of our class AmericanOptionsPricing using six different methods:
simulation_times = np.arange(100, 500, 100)
instances = {'Normal': [], 'ControlVariate': [], 'AntitheticVariate': [
], 'QuasiMC&Normal': [], 'QuasiMC&ControlVariate': [], 'QuasiMC&AntitheticVariate': []}
for s in simulation_times:
instances['Normal'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Normal'))
instances['ControlVariate'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Control'))
instances['AntitheticVariate'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s)))
instances['QuasiMC&Normal'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Normal', isquasi=True))
instances['QuasiMC&ControlVariate'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), MCtype='Control', isquasi=True))
instances['QuasiMC&AntitheticVariate'].append(AmericanOptionsPricing('call', 36, 30, 1, 16,
0.06, 0.2, int(s), isquasi=True))
Then we can plot the price of the American option obtained from the simulation versus the simulation times:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
for method in instances:
plt.plot(simulation_times, [
instance.price for instance in instances[method]], label=method)
plt.legend()
plt.xlabel('Times of simulation')
plt.ylabel('Price')
plt.show()
From the figure above, it's clear that all methods tend to converge as $N$ goes large enough.
ControlVariate methord and AntitheticVariate methord are all more effective than normal Monte Carlo.
Moreover, we can examine the standard deviation of the samples:
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
stddevs = []
for method in instances:
stddevs.append(instances[method][-1].stddev)
plt.bar(instances.keys(), stddevs, color=[
'#66c2a5', '#fc8d62', '#8da0cb', '#e78ac3', '#a6d854', '#ffd92f'])
plt.ylabel('Stddev')
plt.xticks(rotation=30)
plt.show()
From this figure, we can see that ControlVariate method has the smallest Standard deviation.
In conclusion, the Monte Carlo method is a powerful pricing method for discrete arithmetic average European options, and the effect of ControlVariate method is best.
[1] Faure, H. (1982). Discrépance de suites associées à un système de numération (en dimension s). Acta arithmetica, 41(4), 337-351.
[2] Halton, J. H. (1964). Algorithm 247: Radical-inverse quasi-random point sequence. Communications of the ACM, 7(12), 701-702.
[3] Moro, B. (1995). The Full Monte. Risk, 8(2).
[4] Longstaff, F. A., & Schwartz, E. S. (2001). Valuing American options by simulation: a simple least-squares approach. The review of financial studies, 14(1), 113-147.