Pricing options by Monte Carlo simulation is amongst the most popular ways to price certain types of financial options. This article will give a brief overview of the mathematics involved in simulating option prices using Monte Carlo methods, Python code snippets and a few examples. Monte Carlo methods according to Wikipedia:
"Monte Carlo methods, or Monte Carlo experiments, are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. The underlying concept is to use randomness to solve problems that might be deterministic in principle. They are often used in physical and mathematical problems and are most useful when it is difficult or impossible to use other approaches. Monte Carlo methods are mainly used in three problem classes: optimization, numerical integration, and generating draws from a probability distribution."
In order to simulate the price of a European call option, first we must decide on the process that the stock price follows throughout the life of the option \((T- t)\). In the financial literature stocks are said to follow geometric brownian motion. Assume that the stock price \(S\) , in questions pays annual dividend \(q\) and has an expected return \(\mu\) equal to the risk-free rate \(r\) - \(q\) , the volatility \(\sigma\) is assumed to be constant.
The stock price can be modeled by a stochastic differential equation.
Essentially this is a differential equation in which at least one of the terms is a random process. First it may be useful to consider an ordinary differential equation in the context of our problem. Let's consider the case when volatility is 0 i.e. the stock price can be described like a deposit in a savings account paying \(\mu\) per annum The change in any given time increment is then given by
\(dS = \mu S dt\)
Given the price of the stock now \(S_{0}\) we then know with certainty the price \(S_{T}\) at given time \(T\) by separating and intergrating as follows:
\(\displaystyle \int_{0}^{T}\frac{dS}{S} = \displaystyle \int_{0}^T \mu dt \)
Which gives:
\(S_{T} = S_{0}e^{\mu T}\)
It may be useful to notice now that we can write the result above as \(ln(S_{T}) = ln(S_{0}) + \displaystyle \int_{0}^T\mu dt\)
However, since stock prices do exhibit randomness we need to include a stochastic term in the equation above. We can't simply integrate to get a nice result as we have in the equation above, in order to capture the randomness inherent in stock markets we add another term and are SDE is defined as follows:
\(dS= S\mu dt + S\sigma dW(t)\)
Where \(W_{t}\) is a Wiener Process. The equation above is now in the form of an Ito process.
In order to proceed a short word on Ito's Lemma:
Ito's Lemma shown below, states if a random variable follows an Ito Process (example above) then another twice differentiable function \(G\) described by the stock price \(S\) and time \(t\) also follows at Ito process:
(the notation below has been changed from here to keep it consistent with the equations above for the purposes of stock options)
\(dG = ( \frac{\partial G}{\partial S}S\mu + \frac{\partial G}{\partial t} + \frac{1}{2}\frac{\partial^2 G}{\partial S^2}S^2\sigma^2 )dt + \frac{\partial G}{\partial S}S\sigma dW(t)\)
We could apply Ito's lemma to \(G = S\) in order to obtain arithmetic Brownian motion, however using \(G = ln(S)\) which gives a nice property that the stock price is strictly greater than 0. So applying Ito's lemma to \(ln(S)\) first we calculate the partial derivatives with respect to \(t\) and \(S\) as follows:
\(G = ln(S)\)
\( \frac{\partial G}{\partial S} = \frac{1}{S}\) , \(\frac{\partial G}{\partial t} = 0\) , \(\frac{\partial^2 G}{\partial S^2} = -\frac{1}{S^2}\)
Plugging the partial derivatives into Ito's lemma gives:
\(\begin{align*}dG & =( \frac{1}{S}S\mu + 0 - \frac{1}{2}\frac{1}{S^2}S^2\sigma^2 )dt + \frac{1}{ S}S\sigma dW(t) \\ & = (\mu - \frac{\sigma^2}{2})dt\ + \sigma dW(t)\end{align*}\)
Therefore the distrubiton of \(ln(S_{T}) - ln(S_{0}) \) = \((\mu - \frac{\sigma^2}{2})T\ + \sigma \sqrt T\)
The distibution of the stock price at expiration is given by rearraging the equation above an taking the exponential of both sides:
\(S_{T} =S_{0}e^{(\mu - \frac{\sigma^2}{2})dt + \sigma dW(t) }\)
The above can also be written as:
\(ln(S_{T}) = ln(S_{0}) + \displaystyle \int_{0}^t(\mu - \frac{\sigma^2}{2})dt + \displaystyle \int_{0}^t \sigma dW(t)\), \(for \hspace{0.25cm} t \in[0,\cdots,T]\)
Which makes it easier to work with in Python.
Python Example
import numpy as np
import matplotlib.pyplot as plt
def geo_paths(S, T, r, q, sigma, steps, N):
"""
Inputs
#S = Current stock Price
#K = Strike Price
#T = Time to maturity 1 year = 1, 1 months = 1/12
#r = risk free interest rate
#q = dividend yield
# sigma = volatility
Output
# [steps,N] Matrix of asset paths
"""
dt = T/steps
#S_{T} = ln(S_{0})+\int_{0}^T(\mu-\frac{\sigma^2}{2})dt+\int_{0}^T \sigma dW(t)
ST = np.log(S) + np.cumsum(((r - q - sigma**2/2)*dt +\
sigma*np.sqrt(dt) * \
np.random.normal(size=(steps,N))),axis=0)
return np.exp(ST)
S = 100 #stock price S_{0}
K = 110 # strike
T = 1/2 # time to maturity
r = 0.05 # risk free risk in annual %
q = 0.02 # annual dividend rate
sigma = 0.25 # annual volatility in %
steps = 100 # time steps
N = 1000 # number of trials
paths= geo_paths(S,T,r, q,sigma,steps,N)
plt.plot(paths);
plt.xlabel("Time Increments")
plt.ylabel("Stock Price")
plt.title("Geometric Brownian Motion")
Pricing a European call option with a strike of 110 and comparing to Black-Scholes Price.
Update(31/7/2023) An earlier version of this article had the numerator in the python script wrong due to missing parenthesis a thanks to Reda for spotting the error.
def black_scholes_call(S,K,T,r,q,sigma):
"""
Inputs
#S = Current stock Price
#K = Strike Price
#T = Time to maturity 1 year = 1, 1 months = 1/12
#r = risk free interest rate
#q = dividend yield
# sigma = volatility
Output
# call_price = value of the option
"""
d1 = (np.log(S/K) + (r - q + sigma**2/2)*T) / (sigma*np.sqrt(T))
d2 = d1 - sigma* np.sqrt(T)
call = S * np.exp(-q*T)* norm.cdf(d1) - K * np.exp(-r*T)*norm.cdf(d2)
return call
payoffs = np.maximum(paths[-1]-K, 0)
option_price = np.mean(payoffs)*np.exp(-r*T) #discounting back to present value
bs_price = black_scholes_call(S,K,T,r,q,sigma)
print(f"Black Scholes Price is {bs_price}")
print(f"Simulated price is {option_price}")
Black Scholes Price is 3.7451887662410783
Simulated price is 4.1795126304375065
There is considerable difference between the two prices due to the low sample size chosen. Let's try changing N to 100000 and running the script again.
N=100000
paths= geo_paths(S, T, r, q,sigma, steps, N)
payoffs = np.maximum(paths[-1]-K, 0)
option_price = np.exp(-r*T)*np.mean(payoffs)
bs_price = black_scholes_call(S,K,T,r,q,sigma)
print(f"Black Scholes Price is {bs_price}")
print(f"Simulated price is {option_price}")
Black Scholes Price is 3.7451887662410783
Simulated price is 3.8268221823897663
As we increase N towards infinity the price approaches the Black-Scholes price, due to Central Limit Theorem.
A visual representation of what is happening above
n, bins, patches = plt.hist(paths[-1],bins=250);
for c, p in zip(bins, patches):
if c > K:
plt.setp(p, 'facecolor', 'green')
else:
plt.setp(p, 'facecolor', 'red')
plt.axvline(K, color='black', linestyle='dashed', linewidth=2,label="Strike")
plt.title("Distribution of $S_{T}$")
plt.xlabel("$S_{T}$")
plt.ylabel('Count')
plt.legend()
The Monte Carlo Algorithm prices the option as \(call =e^{-rT}[\frac{1}{N} \sum\limits_{i=1}^{N } (S_{T}-K)^+]\) consider the \(^+\) in the previous equation to be only the green values from the plot above.
Path Dependent Options
It may seem like the above was largely unnecessary since we have the Black-Scholes equation, since it takes longer and is less accurate. However, there are a number of cases where a closed form solution is not readily available. Consider again the plot of paths at the beginning of the document. Let's say for some reason someone wants to buy an option that allows the holder to exercise at the most favorable price throughout the specified time interval. Take for example, if the stock in question follows the path below, the holder of this option would be able to choose \(S_{max}\) (dashed red line below).
A visual for comparison
max_=np.max(paths,axis=0)
n, bins, patches = plt.hist(max_,bins=250);
for c, p in zip(bins, patches):
if c > K:
plt.setp(p, 'facecolor', 'green')
else:
plt.setp(p, 'facecolor', 'red')
The price is calculated similarly to the vanilla option: \(lookback =e^{-rT}[\frac{1}{N} \sum\limits_{i=1}^{N } (S_{max}-K)^+]\)
Pricing a lookback with fixed strike of 110.
max_=np.max(paths,axis=0)
payoffs = np.maximum(max_-K, 0)
lookback_price = np.mean(payoffs)*np.exp(-r*T)
print(f"lookback price is {lookback_price}")
Out:
lookback price is 7.076210421600172
The answer above is almost surely underestimating the true value of this option, try setting the steps to a much larger value and notice that the price increases dramatically. This makes sense as Geometric Brownian Motion assumes infinitely divisible time throughout the life of the option, and if we sample at 100 increments over a 6 month period, approx once every 1.2 days, we will miss a lot the highs and therefore undervalue the option. For this reason the steps parameter at the beginning of the document should be adjusted accordingly. Perhaps Python isn't the best tool for this type of calculation. However, it serves to illustrate the concept.
There are many more applications of Monte Carlo methods for option pricing. Links will be posted below to future articles on this topic.