Monte Carlo Simulation in Python; Predict European Call Option Prices
Simulation is a procedure in which random numbers are generated according to probabilities assumed to be associated with a source of uncertainty, such as a new product’s sales or, more appropriately for our purposes, stock prices, interest rates, exchange rates or commodity prices. Outcomes associated with these random drawings are then analyzed to determine the likely results and the associated risk. Oftentimes this technique is called Monte Carlo simulation, being named for the city of Monte Carlo, which is noted for its casinos.
The gambling analogy notwithstanding, Monte Carlo simulation is a legitimate and widely used technique for dealing with uncertainty in many aspects of business operations. For our purposes, it has been shown to be an accurate method of pricing options and particularly useful for path-dependent options and others for which no known formula exists.
To facilitate an understanding of the technique, we shall look at how Monte Carlo simulation has been used to price standard European options. Of course, we know that the Black-Scholes model is the correct method of pricing these options so Monte Carlo simulation is not really needed. It is useful, however, to conduct this experiment because it demonstrates the accuracy of the technique for a simple option of which the exact price is easily obtained from a known formula.
The assumptions of the Black-Scholes model imply that for a given stock price at time t, simulated changes in the stock price at a future time t+Δt can be generated by the following formula:
where S is the current stock price, ΔS is the change in the stock price, r is the continuously compounded risk-free rate, σ is the volatility of the stock and Δt is the length of the time interval over which the stock price change occurs. The variable ε, is a random number generated from a standard normal probability distribution. Recall that the standard normal random variable has a mean of 0.0, a standard deviation of 1.0 and occurs with a frequency corresponding to that associated with the famous bell shaped curve.
Generating future stock prices according to the above formula is actually quite easy. A standard normal random variable can be approximated with a slight adjustment to Rand() function.
import numpy as np
import pandas as pd
from scipy.stats import norm
from numpy.random import randn
from numpy import random as rn
import scipy.stats as si
from matplotlib import pyplot as plt
from IPython.display import Image
%matplotlib inlineThe Rand() function produces a uniform random number between 0 and 1, meaning that it generates numbers between 0 and 1 with equal probability.
def Rand():
d = rn.uniform(0, 1, 1)[0]
return (d)Rand()
0.05886276967175863
A good approximation for a standard normal variable is obtained by the following formula:
Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() + Rand() - 6
or simply 12 uniform random numbers minus 6.
This approximation is based on the fact that the distribution of the sum of twelve uniformly distributed random numbers between 0 and 1 will have a mean of six and a standard deviation of 1. By subtracting 6, we adjust the mean to zero without changing the standard deviation. What we obtain is technically not normally distributed but is symmetric with mean zero and standard deviation of 1, which are three properties associated with the normal distribution. The procedure is widely accepted as a quick and reasonable approximation but might not pass the most demanding tests for normality. On the other hand, many other random number generators would also not pass the most demanding tests either. If you generate a large enough number of random draws, my experience has been that the procedure converges nicely to the Black-Scholes model, which gives a great deal of confidence in it.
After generating one standard normal random variable, you then simply insert it into the right hand side of the above formula for ΔS. This gives the price change over the life of the option, which is then added to the current price to obtain the price of the asset at expiration. We then compute the price of the option at expiration according to the standard formulas, Max(0;ST - X) for a call or Max(0;X - ST) for a put, where X is the exercise price and ST is the asset price at expiration. This produces one possible option value at expiration. You then repeat this procedure many thousands of times, take the average value of the call at expiration and discount that value at the risk-free rate. Some users compute the standard deviation of the call prices in order to obtain a feel for the possible error in estimating the price.
Let us price a call option. The stock price is 164, the exercise price is 165, the risk-free rate is 0.0521, the volatility is 0.29 and the time to expiration is 0.0959. Inserting the above approximation formula for a standard normal random variable in any cell in an Excel spreadsheet produces a random number. Suppose that number is 0.733449. Inserting into the formula for ΔS gives
164⋅(0.0521)⋅(0.0959)+164⋅(0.29)⋅(0.733449)⋅√0.0959=11.62
Thus, the simulated value of the stock at expiration is 164 + 11.62 = 175.62. At that price, the option will be worth Max(0;175.62 - 165) = 10.62 at expiration. We then draw another random number. Suppose we get -0.18985. Inserting into the formula for ΔS, we obtain
164⋅(0.0521)⋅(0.0959)+164⋅(0.29)⋅(-0.18985)⋅√0.0959=-1.98
which gives us a stock price at expiration of 164 -1.98 = 162.02, leading to an option price of Max(0;162.02 - 165) = 0. We repeat this procedure several thousand times, after which we take an average of the simulated option prices and then discount that average to the present using the present value formula exp(-(0.0521)⋅(0.0959)).
Naturally every simulation is different because each set of random numbers is different.
Let’s write a code for the black & scholes model’s price for a European call option:
def PolanitzerNormsdist(x):
PolanitzerNormsdist = si.norm.cdf(x,0.0,1.0)
return (PolanitzerNormsdist)def PolanitzerBlackScholesCall(StockPrice, StrikePrice, Maturity, RiskFreeRate, Volatility):
d1 = (np.log(StockPrice/StrikePrice)+(RiskFreeRate+0.5*Volatility**2)*Maturity)/(Volatility*np.sqrt(Maturity))
d2 = (np.log(StockPrice/StrikePrice)+(RiskFreeRate-0.5*Volatility**2)*Maturity)/(Volatility*np.sqrt(Maturity))
PolanitzerBlackScholesCall = StockPrice*PolanitzerNormsdist(d1)-StrikePrice*np.exp(-RiskFreeRate*Maturity)*PolanitzerNormsdist(d2)
return(PolanitzerBlackScholesCall)
A Monte Carlo procedure written in python produced the following values for this call, whose actual Black-Scholes price is 5.79.
# Assumptions:
StockPrice = 164
StrikePrice = 165
Maturity = 0.0959
RiskFreeRate = 0.0521
Volatility = 0.29a = PolanitzerBlackScholesCall(StockPrice, StrikePrice, Maturity, RiskFreeRate, Volatility)
print("\033[1m The price of the European call option based on the Black-Scholes Model is {:.3}".format(a))
The price of the European call option based on the Black-Scholes Model is 5.79
Let’s write a code for the Monte Carlo model’s price for a European call option:
def PolanitzerNormsinv(x):
PolanitzerNormsinv = si.norm.ppf(x)
return (PolanitzerNormsinv)def PolanitzerMonteCarloSimulationCall(StockPrice, StrikePrice, Maturity, RiskFreeRate, DividendYield, Volatility, Simulations):
drift = (RiskFreeRate - DividendYield - 0.5*Volatility**2) * Maturity
shock = Volatility * np.sqrt(Maturity)
SUM = 0
for i in range(1,Simulations):
St = StockPrice * np.exp(drift + shock * PolanitzerNormsinv(Rand()))
SUM = SUM + np.max((St-StrikePrice,0))
PolanitzerMonteCarloSimulationCall = np.exp(-RiskFreeRate*Maturity) * SUM / Simulations
return(PolanitzerMonteCarloSimulationCall)
Let’s find the Monte Carlo model’s price for different sample sizes, n:
Simulations = [1000, 10000, 50000, 100000]
DividendYield = 0
MonteCarloSimulationV = []
for i in range(0,len(n)):
MonteCarloSimulationV.append(round(PolanitzerMonteCarloSimulationCall(StockPrice, StrikePrice, Maturity, RiskFreeRate, DividendYield, Volatility, Simulations[i]),2))
# initialize list of lists
data = [Simulations, MonteCarloSimulationV]
INDEX = [“n”,”Call Price”]# Create the pandas DataFrame
df = pd.DataFrame(data, index = INDEX)
# print dataframe.
df.T
It would appear that a sample of at least 50,000 is required for the simplest case of a standard European option.
The option price obtained from a Monte Carlo simulation is a sample average. Thus, its standard deviation is the standard deviation of the sample divided by the square root of the sample size. Although this is a point from the most elementary principles of statistics, it is worthwhile to distinguish between the standard deviation of a sample and the standard deviation of the sample mean. The latter is the former divided by the square root of the sample size. Hence, the sample mean is much less volatile than the sample values themselves.
Consequently, the error reduces at the rate of 1 over the square root of the sample size. Notice what this means for increasing the sample accuracy by increasing the sample size. Suppose σ is the standard deviation of the sample. We first conduct a Monte Carlo simulation using n1 random drawings. Since the option value is a sample mean, the standard deviation of our estimate of the option value:
Now suppose we wanted to reduce that standard deviation in half. How much larger must the sample be? Let this new sample size be n2. Then its standard deviation of the estimate of the option price is
Now note that,
if and only if
Thus, to achieve a 50% reduction in error, i.e., a 50% increase in accuracy, we must quadruple the number of random drawings. That is, the standard error reduces only at the rate of the square root of the sample size, not at the rate of the sample size itself.
About the Author
Roi Polanitzer, CFV, QFV, FEM, F.IL.A.V.F.A., FRM, CRM, PDS, is a well-known authority in Israel the field of business valuation and has written hundreds of papers that articulate many of the concepts used in modern business valuation around the world. Mr. Polanitzer is the Owner and Chief Appraiser of Intrinsic Value — Independent Business Appraisers, a business valuation firm headquartered in Rishon LeZion, Israel. He is also the Owner and Chief Data Scientist of Prediction Consultants, a consulting firm that specializes in advanced analysis and model development.
Over more than 17 years, he has performed valuation engagements for mergers and acquisitions, purchase price allocation (PPA) valuations, goodwill impairment test valuations, embedded option and real option valuations, employee stock option (ESOP) valuations, common stock valuations (409A), splitting equity components and complicated equity/liability instrument valuations (PWERM / CCM / OPM), contingent liability, guarantees and loan valuations, independent expert opinions for litigation purposes, damage quantifications, balancing resources between spouses due to divorce proceedings and many other kinds of business valuations. Mr. Polanitzer has testified in courts and tribunals across the country and from time to time participates in mediation proceedings between spouses.
Mr. Polanitzer holds an undergraduate degree in economics and a graduate degree in business administration, majoring in finance, both from the Ben-Gurion University of the Negev. He is a Full Actuary (Fellow), a Corporate Finance Valuator (CFV), a Quantitative Finance Valuator (QFV) and a Financial and Economic Modeler (FEM) from the Israel Association of Valuators and Financial Actuaries (IAVFA). Mr. Polanitzer is the Founder of the IAVFA and currently serves as its chairman.
Mr. Polanitzer’s professional recognitions include being designated a Financial Risk Manager (FRM) by the Global Association of Risk Professionals (GARP), a Certified Risk Manager (CRM) by the Israel Association of Risk Managers (IARM), as well as being designated a Python Data Analyst (PDA), a Machine Learning Specialist (MLS), an Accredited in Deep Learning (ADL) and a Professional Data Scientist (PDS) by the Professional Data Scientists’ Israel Association (PDSIA). Mr. Polanitzer is the Founder of the PDSIA and currently serves as its CEO.
He is the editor of IAVFA’s weekly newsletter since its inception (primarily for the professional appraisal community in Israel).
Mr. Polanitzer develops and teaches business valuation professional trainings and courses for the Israel Association of Valuators and Financial Actuaries, and frequently speaks on business valuation at professional meetings and conferences in Israel. He also developed IAVFA’s certification programs in the field of valuation and he is responsible for writing the IAVFA’s statement of financial valuation standards.