**Part I:** Risk-neutral pricing, Monte Carlo integration vs. the Black-Scholes formula

The Black-Scholes pricing model is the de facto standard way of pricing options in today’s financial markets.

Notwithstanding the many reasons for distrust of its underlying stochastic model for stock price dynamics, the fact is that market quotes of options prices are in fact most commonly quoted in terms of the value of a free parameter in the BS model, known as the implied volatility, that forces the model price in a specific instance to conform to the observed market price.

Since the BS model is actually not a faithful representation of reality, this creates an inconsistency between different calculation instances with respect to the fitting of this parameter, the upshot being that implied volatility varies with different strike prices and maturities, creating the so-called volatility surface, also known as the “volatility smile” when focusing on the dependence on strike.

Contents

## The Black-Scholes Model

Let be the price of a stock at time *t*, the Black-Scholes model assumes the Geometric Brownian Motion 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:

We assume that time *t* is measured in years, while 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 and consider a European call option contract 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 in, for the price of 1.1. The fact that the option has to be exercised at the specific date of expiry is what is meant by European; and American option can be exercised at any time until expiry, and a Bermudan option as a hybrid between the two, allowing exercise at a number of specific dates, all before or coinciding with the expiry date.

What would be the realized profit/loss from this contract, as a function of , 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 .

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

.

## 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 a little twist: the expectation is not computed according to the actual probability distribution of the stock, even assuming the model had the correct statistical assumptions.

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 while volatility is unchanged. The risk-free interest rate is the rate that an investor could expect to receive by lending money into the market 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 the soundness of this assumption may be debatable these days. In our imaginary risk-neutral world, logreturns would have the distribution given by

The option price at time 0 would then be obtained by computing the expectation

,

where 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 and compute the mean to give our estimate of the expectation. According to the strong law of large numbers, this estimate has probability 1 of converging to the true expectation 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:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.7771459645311785

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.777 today (plus some transaction fee) for an option to buy the 100 shares in one year at a strike price of $1.1.

## 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 works for all European options no matter what payoff function we stipulate, and we could even experiment with other assumptions on the stock price process. However, in the restricted 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 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

and

## Python Code

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:FinxterTutorialsBlack-Scholes-1>python bsformula.py

4.775025500484964

Comparing 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. We run it three times to illustrate the fact that we get closer to the formula, but there’s still some fluctuations due to the randomness in our sample:

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.774596150369479

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.773673080552248

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.7739917521923525

With 100 million samples we seem to have converged at least to the third decimal, and fluctuations are a bit smaller but still measurable:

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.775559370501232

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.775186078778582

D:FinxterTutorialsBlack-Scholes-1>python riskneutral.py

4.775651011152207

The post Black-Scholes Option Pricing in Python first appeared on Finxter.