?Ā Part I: Risk-neutral valuation, Monte Carlo integration vs. the Black-Scholes formula
You can find the code in the GitHub repository for this article.
The Black-Scholes (BS) pricing model is still a de facto standard method of pricing financial options.
Even though there has been much research into improved and possibly more realistic models, the fact is that the BS model is implicitly assumed in the way most option prices are quoted in practice, in terms of the model parameter called implied volatility. Anyone working with options in the finance industry will have to know the basics of this valuation method. In this tutorial, we will guide you through the minimal amount of theory needed to implement the pricing method in Python and then look at a basic calculation example.
As you go over the article, feel free to watch the associated explainer video:
The Black-Scholes Model
Let StĀ be the price of a stock at time t. The Black-Scholes model is based on the Geometric Brownian Motion (GBM) model which implies that the logarithmic return of the stock price can be described by a normal distribution whose variance is proportional to the time step.
It is also a part of the definition of the GBM that these log-returns are statistically independent when measured over disjoint time intervals. We assume that time tĀ is measured in years, while StĀ is measured in the currency in which the stock price is denominated, such as USD or EUR. Ī¼Ā is called the expected return and ĻĀ is the volatility. Now letās assume that the initial price S0 = 1 and letās say expected return is 5%Ā and the volatility is 20%, in other words, Ī¼ = 0.05 Ā and Ļ = 0.2.
The reason for the term
in the expression for the location parameter of the normal distribution is that the convexity of the exponential introduces a bias which has to be compensated for; this phenomenon is also known as Jensenās inequality.
The distribution of the future stock price one year from now, given that the price right now is equal to 1, is known as the lognormal distribution and can be plotted in Python using the scipy.stats submodule and the matplotlib package to create the actual plot:
from scipy.stats import lognorm import matplotlib.pyplot as plt import numpy as np import math mu = 0.05 sigma = 0.2 S0 = 1 x = np.linspace(0,2,1000) y = lognorm.pdf(x, sigma,0,math.exp(mu)) plt.plot(x,y) plt.show()
We can see that the most likely outcome, aka the āmodeā or maximum of the density function, is located slightly to the left of 1, corresponding with the mathematical fact that the mode of the lognormal distribution can be shown in this case to equal:
This might disappoint some people given that the expected return was supposed to be 5%, however the distribution is skewed and so the mode is actually not equal to the mean which is:
In other words, the most likely outcome is that we lose 1% on the investment, however the skewed probability of large positive returns implies that on average we will still make a 5%Ā return!
If we want to see this stock price process from a different perspective we can do a path simulation, i.e. randomly generate hypothetical trajectories of the stock price according to the probability law of the GBM. We can do this in Python just using the numpy package. In the example below we have simulated 50 realizations of the stock price path over 1 year, divided into 100 uniform time increments:
import numpy as np
import matplotlib.pyplot as plt
Nsim = 30
t0 = 0
t1 = 1
Nt = 100
mu=0.05
sigma=0.2
S0 = 1
t = np.linspace(t0,t1,Nt)
dt = (t1-t0)/Nt
S = np.zeros([Nsim,Nt])
S[:,0] = S0
for j in range(0, Nt-1):
S[:,j+1] = S[:,j]*np.exp((mu-sigma**2/2)*dt + sigma*np.sqrt(dt)*np.random.normal(0,1, Nsim))
for j in range(0,Nsim):
plt.plot(t,S[j,:])
plt.show()
The results are seen in the next figure, it may look vaguely like real stock price charts, at least if we ignore the fact that real stock prices sometimes make very sudden, sharp downward or upward movements, as the result of some crisis or other random event or revelation which suddenly changes the marketās perception of the value of a stock.
Now letās consider a European call option on the stock initiated at time t = 0, with a notional amount of 100 shares, expiry date t = 1Ā and strike price $1.1. The option contract gives the holder the right but not the obligation to buy 100 shares of the stock at the expiry date one year from now, for the price per share of 1 dollar and 10 cents. The fact that the option has to be exercised at the specific date of expiry is what is meant by a European option, there are also American options which can be exercised any time until expiry and there are even other types which have different conventions with regard to exercise.
What would be the realized profit/loss from this contract, as a function of S1, the stock price at expiry? If I could buy the 100 shares for $110, and the market price has increased to above $1.1, of course, I could immediately turn around and sell the share again for the higher price. So the payoff would be 100 * S1 – 1.1. However if the stock price fell or didnāt increase by more than 10%, my resell value would not exceed what I had paid for the shares, so I would not have any interest in exercising the option. The payoff would then be zero. So the payoff would in any event be given by the random variable:
We can plot this function to visualize the asymmetric dependence of payoff on the final outcome of the stock price:
import numpy as np
import matplotlib.pyplot as plt
k = 1.1
def payoff(x):
return 100*np.maximum(0,x-k)
x=np.linspace(0,2, 100)
y=payoff(x)
plt.plot(x,y)
plt.xlabel('Stock price at expiry')
plt.ylabel('Payoff')
plt.show()
The Risk-Neutral World
The Black-Scholes model implies that the value of an option at a time tĀ prior to expiry should be equal to the expected present value of its future payoff, with just a little twist: the expectation is not computed using the actual probability distribution of the stock price, even if we did actually believe in the modelās statistical assumptions about the world. Instead, the expectation should be taken according to a risk-neutral probability distribution, which means that the expected return Ī¼ is replaced by the risk-free interest rate rĀ while volatility is unchanged. The risk-free interest rate is the rate of return that an investor could expect to receive by lending money without running the risk of the borrower defaulting; usually, short-term government bond rates are used as a proxy for risk-free interest rates but even this assumption may be debatable these days. In our imaginary risk-neutral world, log-returns would have the distribution given by
The option price at time 0 would then be obtained by computing the expectation
where EQ denotes the risk-neutral expectation. Now letās set this up in Python and compute the price, we can use Monte Carlo integration to compute the expectation, meaning that we draw a large number of random samples from this probability distribution (which would correspond to the end values of simulation paths like those we showed earlier),Ā and compute the mean to give our estimate of the expectation. According to the law of large numbers this estimate approximates the true expectation to arbitrary precision if we just make the sample size large enough.
import numpy as np
import math
Nsim = 10000
amount_underlying = 100
strike = 1.1
sigma = 0.2
mu = 0.06
r = 0.015
def payoff(x):
return amount_underlying * np.maximum(0, x-strike)
num0 = np.random.normal(0,1,Nsim)
S0 = 15
S1 = np.exp(r-sigma**2/2+sigma*num0)
C0 = math.exp(-r)*np.mean(payoff(S1))
print(C0)
Now executing the code in Python we get the following result:
D:\Finxter\Tutorials\Black-Scholes-1>python riskneutral.py 4.555089461101134
What this means in practical terms is that with a share price of $1, an implied volatility level of 20%, and a risk-free interest rate of 1.5%, we should expect to pay $4.555 today (plus some transaction fee) for an option to buy the 100 shares in one year at $1.1 per share.
Exact Computation via the Black-Scholes Formula
The method of Monte Carlo integration to obtain the risk-neutral expectation of the discounted payoff function is a very general method in the sense that it would work for all European options no matter what payoff function we assumed, and we could even have experimented with other assumptions on the stock price process. However, in the simplified lognormal world of the BS model, it turns out that thereās actually a closed-form equation that describes the call-option price, the so-called Black-Scholes formula:
where C0 Ā is the price of the option at the start of the contract, KĀ is the strike price, tĀ is time to expiry,Ā N Ā is the cumulative distribution function of the standard normal distribution, while d1Ā and d2Ā are given by
Letās also transcribe these functions into Python code:
import numpy as np
from scipy.stats import norm
amount_underlying = 100
strike = 1.1
sigma = 0.2
mu = 0.06
r = 0.015
S0 = 1
t = 1
def fun_d1(sigma,k,t,r,x):
return (np.log(x/k) + (r+sigma**2/2)*t)/(sigma*np.sqrt(t))
def fun_d2(sigma,k,t,r,x):
return fun_d1(sigma,k,t,r,x) - sigma*np.sqrt(t)
def call_value(amount_underlying, sigma,k,t,r,x):
d1 = fun_d1(sigma,k,t,r,x)
d2 = fun_d2(sigma,k,t,r,x)
temp = norm.cdf(d1)*x-norm.cdf(d2)*k*np.exp(-r*t)
return amount_underlying * temp
C0 = call_value(amount_underlying, sigma,strike,t,r,S0)
print(C0)Now if we run this with Python, we get the following result:
D:\Finxter\Tutorials\Black-Scholes-1>python bsformula.py 4.775025500484964
Compared to our Monte Carlo result above, we see that thereās a difference in the third decimal, 4.775 vs 4.777. We used 10000 samples for our simulation, letās run it again with 1000 times the sample size, changing the Nsim parameter to 10,000,000:
D:\Finxter\Tutorials\Black-Scholes-1>python riskneutral.py 4.774596150369479
Now we are getting closer to the formula-based calculation, indicating that the two different methods are indeed equivalent; to get exactly identical answers we would have to go to the limit of an infinitely large sample size.