Monte Carlo Simulation and Option Pricing with Stochastic Processes¶
Monte Carlo simulation a computational method that utilies random numbers. It is among the most important numerical techniques in finance, this mainly stems from the fact it is the most flexible method when it comes to the evaluating integrals.
Stochatic proccesses are equations with integrals that have an element of randomized input that can be used to describe the random movement of certain phenomena over time, in this option and stock values!
Utilsing Monte Carlo simulation we can evaluate these stochastic procces to be able to value financial instruments!
step 1 - Random number generation¶
Stochastic proccess require random number input, these random numbers can be drawn from a variety of distributions dependent upon the process.
Put simplistically we can evaluate a stochastic process integral by inputting a vector of random numbers drawn from the descriptive distribution! Each random number drawn represents a simulation of how the stochastic process which describe our financial instrument of interest over one time interval.
Below I demonstrate how python can be used to draw random numbers from both a unifrom and standard distribution and how these numbers can be shaped into arrays to feed our multiple simulations and mutliple time intervals we may want to interrogate.
import math
import numpy as np
import numpy.random as random
import matplotlib.pyplot as plt
plt.style.use('seaborn')
#fixes the seed value
random.seed(1)
# fixes the number of decimal places returns by the numpy functions, suppress gets rid of scientific format
np.set_printoptions(precision=3, suppress=True)
#rand() by default returns uniformly distibuted numbers between 0 and 1, inputting a single number returns a 1d array of that size fo random numbers
random.rand(10)
#two digits returns a 2d array of the equivalent shape
random.rand(10,2)
#lower limit
a = -0.01
#upper limit
b=0.1
# we can translate the numbers to be between limits as follows
uniform_random_numbers = random.rand(1000)*(b-a)+a
#generate 10 standard normally distributed numbers
SD_random_numbers = random.standard_normal(1000)
fig, (ax1, ax2) = plt.subplots(1,2,figsize=(14,6))
ax1.hist(uniform_random_numbers , bins=30, color = 'r')
ax1.set_title()
ax2.hist(SD_random_numbers, bins=30)
Step 2 BLack Scholes Merton for simulating index movements¶
The Black Scholes Merton setup can be used to simulate the random movement of a stock over a defined time period the equation is as follows:
S_t = S_0 exp((r-sigma^2/2)T+ sigmasqrt(T)*Z)
where:
S_t, is the stock price at time t
S_0, is the stock price at time 0
r, is the constant riskless short rate
sigma, is the constant volatilty (standard deviaiton of returns of S)
z , standard normally distributed random variable
below we are going to project what the value of the APPL stock price will be in 4 years using monte carlo simulation on the above equation! The more simulations we run the more accurate our expectatoin for the outcome will be at the trade off of increased computing time. I will use 100,000 sims
import pandas as pd
#import the appl stock data over the last ten years to eastblish the volatility of the stock
APPL = pd.read_csv(r'C:\Users\Adam\Desktop\Personal\Python Syntax\Python Input Files\finance\AAPL.csv', parse_dates=True ,index_col=0)
APPL.sort_index(ascending=True, inplace=True)
#resample the data to find the annual returns
APPL_Yearly = APPL.resample('1Y').last()
APPL_Yearly['Annual_Return'] = np.log((APPL_Yearly['Price'])/APPL_Yearly['Price'].shift(1))
APPL_Yearly.dropna(inplace=True)
#calculate the stdev of the retunrns
APPL_std_y = APPL_Yearly['Annual_Return'].std()
""" APPL standard deviation of returns over the last 10 years has been 22.7% """
"""
As of writing this project
APPL stock price is at $311
APPL stock volatility = 22.7% (above)
USA 1 Yr treasury rate sits at 0.16%
inputting these paramters into the equation we get the follow: """
S_0 = 311
r = 0.0016
sigma = 0.227
T=4 #number of years
I= 100000 # number of simulations
""" we can create the output via vectorized expression - the output is an array of all the simulated stock price outcomes """
S_t = S_0 * np.exp((r-0.5*sigma**2)* T + sigma*math.sqrt(T)* random.standard_normal(I))
""" plotting the results"""
fig1, ax = plt.subplots(figsize=(16,8))
ax.hist(S_t,bins=50)
ax.set_title('Simulated outcomes for APPL sotck price in 4 years using simple BSM')
ax.set_xlabel('APPL Stock Price ($)')
ax.set_ylabel('frequency')
""" taking the mean of the outcomes essentially multiples each outcome by its probability
of occuring = 1/I to return to us the mean expectation for the outcome! In this case we get $592.6
The value of APPL is expected to be $312 in 4 years"""
print(" The expectation for the APPL stock price in 4 years is: {}".format(S_t.mean()))
Step 3 - Dynamic BLack Scholes Merton - Geometric Brownian Motion¶
A stochastic process, roughly speaking is a sequence of random variables.Stochastic processes used in finance exhibit the Markov property, which means t+1 is related to t, but not related to any prior time periods (the proces sis memoryless).
Incorporating this dynamic markov property into the Black Scholes Merton up we get the Geometric Brownian motion equation!
If we decrease the size of our time interval and have multiple sequential intervals we increase the simulation realism and thereforeincrease the accuracy of our evaluation again, yet again at the expense of computation time. We are going to use time intervals of 1 monthto evaluate again the apple stock price in 4 years!
Becuase we are doing 100,000 simulations this time over 48 months we need to set up a 2d array of the stock price movements to record the movementover each interval for each simulation.
r = 0.0016 #Current US short rate
sigma = 0.227 #AAPL annual vol
T=4 #years
I= 10000 # number of simulations
M=48 # number of months within overall time period
dt = T/M #time interval
S = np.zeros((M+1,I)) # establish a 2d array of the stock prices at each month for each simulation
S[0] = 311 # establish the first month stock price
""" loops through each time interval recording the updated share price"""
for t in range (1, M+1):
S[t] = S[t-1]*np.exp((r-0.5*sigma**2)* dt + sigma*math.sqrt(dt)* random.standard_normal(I))
""" plot the dynamically simulated geometric Brownaian motion paths"""
fig1, ax = plt.subplots(figsize=(16,8))
ax.plot(S[:,:10], lw=1.5) # only plot 10 simulations and set line width to 1.5
ax.set_title('Simulated outcomes for APPL stock price in 4 years using dynamic BSM')
ax.set_xlabel('time (months)')
ax.set_ylabel('APPL Stock Price ($)')
#it is expected that APPl stock will reach 592.6
print(" The expectation for the APPL stock price in 4 years is: {}".format(S[48].mean()))
step 4 - Heston Stocastic Volatility model¶
A simple sotchatsic process is the most basic way of simulating how a stock price may move. In reality we know that interest rates and stock volatility can change. Adding some level of fluctuation into their values would create a more realistic view of reality. Therefore in the 1990s we had the introduction of the Heston stocastic Volatility model, where a stochastic element was added in to describe the volatility paramter.
The original differental equation S_t = S_0 exp((r-sigma^2/2)T+ sigmasqrt(T)*Z)
now has sigma = sqrt(V[t]) where:
V_T is a square root diffusion process (also known as a mean reversion process as volatility levels like to revert back to a mean level)
where V_t = V[t-1] + kappa(theta-np.maximum(vh[t-1],0))dt + sigma np.sqrt(np.maximum(vh[t-1],0))math.sqrt(dt) vol_Z[1])
kappa, is the mean reversion value (how quickly the volatility returns to the long term value) theta , the long term value (mean reversion value) for the volatility.
IMPORTANT -- In the volatility process we have VOL_Z instead of Z to represent the fact that the random processes aren't the same but they are correlated!
vol_Z * Z = rho
To account for the correlation between the two processes, I need to calculate te Cholesky decomposition of the correlation matrix as follows:
""" firstly we define the level of correlation between the two processes"""
rho = 0.6
"""then the matrix is defined as follows: """
corr_mat = np.zeros((2,2))
corr_mat[0,:] = [1.0,rho]
corr_mat[1,:] = [rho, 1.0]
""" then we can calculated the Cholesky matrix from the correlation matrix"""
cho_mat = np.linalg.cholesky(corr_mat)
Basically we use the matrix so that when we draw the random numbers for the brownation process we can matrix multiple them by thecholesky matrix to get a correlated set of random numbers. As it is likely there is a relationship between volatility and stock price movementso we want to be creating linked scenarios not unrealistic ones.
- Now we can pull both stochastic processes together to get a new expectation for the APPL stock price in 4 years, but this time using the more realistic Heston Stochastic Volatility Model!
""" step 1 set up the paramters and the random variable matricies """
S0 = 311 #APPL stock price now
v0 = 0.228 #Appl annual vol
r = 0.0016 #current US short rate
kappa = 3 # mean reversion rate
theta = 0.15 #long term volatility rate
sigma = 0.1 #constant volatility of volatility processs
T=4 #number of years
# we need to set up the parameters of the siumulation again
I= 100000 # number of simulations
M=48 # number of months
dt = T/M #time interval
# we need to create a 3d array of standard normal random numbers to simulate each month for each simulation and each stocastic process
ran_num = random.standard_normal((2,M+1,I))
v = np.zeros_like(ran_num[0]) #select the random number for the first brownian process for volatility
vh = np.zeros_like(ran_num[0]) # we need to establish 2 matrics as the vh matrics will be uses to not allow any negative values the volatility is not allowed to go negative
#we need to eastblish the values for volatility as these feed in as a variable into the movements for the stock prices
v[0] = v0
vh[0] = v0
""" step 2 compute the volatility values """
for t in range(1, M +1):
ran = np.dot(cho_mat, ran_num[:,t,:]) #picks out the relevant random number subset and transforms it via the cholesky matrix
vh[t] = vh[t-1] + kappa*(theta - np.maximum(vh[t-1],0))*dt + sigma * np.sqrt(np.maximum(vh[t-1],0))* math.sqrt(dt) * ran[1]
""" step 3 remove any negative volatility values """
# any negative volatilities that are produced by the stochastic process are replaced with 0s you cant have negative volatility
v = np.maximum(vh,0)
""" step 4 simulate the stock prices """
#establish a matrix for the simulation outcomes
S = np.zeros_like(ran_num[0])
#set the AAPL stock price now
S[0] = S0
# loop through all the months and calulcate the stock prices for each simulation bringing in the correpsonding values of v at each time period for each simulation from above
for t in range(1, M+1):
ran = np.dot(cho_mat, ran_num[:,t,:]) #picks out the relevant random number subset and transforms it via the cholesky matrix
S[t] = S[t-1] * np.exp((r - 0.5 * v[t])* dt + np.sqrt(v[t]) * math.sqrt(dt)* ran[0]) # this time we use ran[1] not ran[2] the random numbers aren't the same but they are correlated
fig, (ax1 , ax2) =plt.subplots(1,2, figsize = (16,8))
ax1.hist(S[-1], bins=50)
ax1.set_title('Simulated APPL stock prices in 4 years')
ax1.set_xlabel('AAPL Stock Price ($)')
ax1.set_ylabel('frequency')
ax2.hist(v[-1],bins=50)
ax2.set_title('Simulated APPL stock Volatility in 4 years (reverts to long term mean)')
ax2.set_xlabel('volatility')
fig, (ax1, ax2) = plt.subplots(2,1, sharex=True, figsize=(16,8))
ax1.plot(S[:,:10], lw=1.5)
ax1.set_title('Simulated APPL stock prices over 4 years')
ax1.set_ylabel('AAPL Stock Price ($)')
ax2.plot(v[:,:10], lw=1.5)
ax2.set_title('Simulated APPL stock Volatility over 4 years (reverts to long term mean)')
ax2.set_xlabel('Time (Months)')
ax2.set_ylabel('volatility')
""" we find that the expectation for the appl stock after 4 years is $352 """
print(" The expectation for the APPL stock price in 4 years is: {}".format(S[-1].mean()))
Step 5 - Valuing Options - European Option¶
One of the most important applications of Monte Carlo simulation is the valuation of contingent claims ( options, derivatives, hybrids instruments etc..) """
The pay off of a euopean call option on an index or appl stock in our exmaple is give by:
- h(S_t) = max (S_t - K, 0)
where K, is the strike price and S_t is the price of the stock at maturity of the option (time t). under a risk neutral expectation we can calculate the value of the option as follows:
we can simulate S_t using stochastic processes as shown above and then we use the projected stock prices to calculate the projected payoffs.
We can then multiply each payoff by its probability of taking place which is 1/number of simulations ( which is essentially the mean of the payoffs). This is our expected pay off - we then discount that at the risk free rate to obtain the expected present value of the option.
Below I will code a function for calulcating european options then I will calculate the current value of a european call option on APPL stock with a maturity of 4 years and a strike price of $330!
def Eur_Option_price (Stock_Price, r, sigma, Maturity, strike_price, Simulations, Call=True, dt=1):
''' Valuation of american option using black-Scholes-Merton by Monte Carlo simulation (of index level at maturity)
Args:
Stock_Price (float): current stock price of index
r float): current short rate
sigma (float): annual volatilty of stock
Maturity (float): Time till maturity in year fractions
strike_price (float): the strike price of the option
Simulations(int): the number of simulations you want to run
Returns:
float: estimated present value of European call options
'''
import numpy.random as random
random.seed(1)
dt = dt #we want to work in year time frames
S = np.zeros((Maturity+1,Simulations)) # establish a zero 2d array for the expected stock prices at each month for each simulation
S[0] = Stock_Price # current stock price
# establish all the random numbers for each simulation at each year value
SN_numbers = random.standard_normal((Maturity, Simulations))
payoffs=[]
""" loops through each time interval recording the updated share price"""
for t in range (1,Maturity+1):
# every simulation for each year are carried out at once using vectorisation
S[t] = S[t-1]*np.exp((r-0.5*sigma**2)* dt + sigma*math.sqrt(dt)* SN_numbers[t-1])
#calculate payoff at maturity
if Call:
payoffs.append(np.maximum(S[-1]-strike_price,0))
else:
payoffs.append(np.maximum(strike_price - S[-1],0))
return math.exp(-r*Maturity) * np.mean(payoffs)
APPL_Call_Price_Eur = Eur_Option_price(311, 0.00016, 0.228, 4, 330, 10000,Call=True )
print(" The present value of a Euro call option on APPL stock price with maturity 4 years and a strike price of $330 is: ${}".format(APPL_Call_Price_Eur))
Step 6 - Valuing Options - American Option¶
American options differ from european options in that they can be exercised at any time. Therefor they need to be calclated differently.
If we consider the stock price movements in terms of a binomial tree payoffs no longer are just caluclated at maturity but need to be calculated at each branch of the tree with the maximum payoff at that time being calculated (the maximum between waiting till till next interval or exercising the option at that time.
The payoffs from each point of the binomial tree are then discounted to the present moment in time to calculate the present value of the american option!
def Ame_Option_price (Stock_Price, r, sigma, Maturity, strike_price, Simulations, Call=True, dt=1):
''' Valuation of american call option using black-Scholes-Merton by Monte Carlo simulation (of index level at maturity)
Args:
Stock_Price (float): current stock price of index
r float): current short rate
sigma (float): annual volatilty of stock
Maturity (float): Time till maturity in year fractions
strike_price (float): the strike price of the option
Simulations(int): the number of simulations you want to run
Returns:
float: estimated present value of american call options
'''
import numpy.random as random
random.seed(1)
dt = dt #we want to work in year time frames
# 1 branch dicount variable
df = math.exp(-r*dt)
S = np.zeros((Maturity+1,Simulations)) # establish a zero 2d array for the expected stock prices at each month for each simulation
S[0] = Stock_Price # current stock price
SN_numbers = random.standard_normal((Maturity, Simulations))
""" loops through each time interval recording the updated share price"""
for t in range (1,Maturity+1):
#stock price movements are calulcated just as for the european wit GBM
S[t] = S[t-1]*np.exp((r-0.5*sigma**2)* dt + sigma*math.sqrt(dt)* SN_numbers[t-1])
#calculate payoff at every point in the stock price matrix
if Call:
#this time we calculate the payoffs at every end of year (every branch)
us_payoffs = np.maximum(S-strike_price,0)
else:
us_payoffs = np.maximum(strike_price - S,0)
#now we carry out a least squares Regression fit to find the maximum payoff at eveyr point the max between (payoff at that point or holding)
#we first take a copy as the payoffs between an american and eropean would be the same if you held to maturity but we are going to iterate back through the rest and find out if not activating early would increase the payoff
V = np.copy(us_payoffs)
#We iterate back from the end of the binomal tree to value the american option at each brnach
for t in range(Maturity - 1, 0, -1):
#we start with the 2nd to last stock price branches and compare there value to the discounts payoff value to see if its worth more or less and iterarte back from there
#reg is the expected pay off next time interval / current stock price this is where we add in the probability around we dont know how the stock will move next
reg = np.polyfit(S[t], V[t + 1] * df,7)
#gives us the payoff next interval
wait_payoff = np.polyval(reg, S[t])
#we then test to see which payoff is greater the present vlaue of waiting for the payoff next interval or activiating the payoffs now. The value returned is the greater of the two
V[t] = np.where(wait_payoff > us_payoffs[t], V[t+1]*df, us_payoffs[t])
#the first row of branches are then averaged to factor in the probability element of each branch and discounted once more to the current positoin, the top of the tree
value = df * np.mean(V[1])
return value
APPL_Call_Price_Ame = Ame_Option_price(311, 0.00016, 0.228, 4, 330, 10000, True )
print(" The present value of a American call option on APPL stock price with maturity 4 years and a strike price of $330 is: ${}".format(APPL_Call_Price_Ame))
Step 7 - European vs American Option¶
As epected the US option has a higher value this is becuase the extra funcitnoality with an american option leads to greater value compared to the european option, the european option with the same strike price would act as a lower bound on the potential value of an american option """
We can see this demonstrate below across multiple strike prices:
euro_value = []
american_value = []
k_list = np.arange(200, 280, 5)
for k in k_list:
euro_value.append(Eur_Option_price(311, 0.00016, 0.228, 4, k, 100000, Call=True ))
american_value.append(Ame_Option_price(311, 0.00016, 0.228, 4, k, 100000, Call=True ))
euro_value = np.array(euro_value)
american_value = np.array(american_value)
fig ,(ax1 , ax2) = plt.subplots(2,1, sharex =True, figsize=(10,6))
ax1.plot(k_list, euro_value, 'b' ,label='Euro Call')
ax1.plot(k_list,american_value,'ro', label = 'US Call' )
ax1.set_ylabel('put option value')
ax1.legend(loc=0)
wi=1
ax2.bar(k_list, (american_value - euro_value)/euro_value*100)
ax2.set_xlabel('Strike')
ax2.set_ylabel('Early exercise premium in %')
