## Black-Scholes Formula and Python Implementation

The Black-Scholes model was first introduced by Fischer Black and Myron Scholes in 1973 in the paper "The Pricing of Options and Corporate Liabilities". Since being published, the model has become a widely used tool by investors and is still regarded as one of the best ways to determine fair prices of options.

The purpose of the model is to determine the price of a vanilla European call and put options (option that can only be exercised at the end of its maturity) based on price variation over time and assuming the asset has a lognormal distribution.

## Assumptions¶

To determine the price of vanilla European options, several assumptions are made:

• European options can only be exercised at expiration
• No dividends are paid during the option's life
• Market movements cannot be predicted
• The risk-free rate and volatility are constant
• Follows a lognormal distribution

## Non-Dividend Paying Black-Scholes Formula¶

In Black-Scholes formulas, the following parameters are defined.

• $S$, the spot price of the asset at time $t$
• $T$, the maturity of the option. Time to maturity is defined as $T - t$
• $K$, strike price of the option
• $r$, the risk-free interest rate, assumed to be constant between $t$ and $T$
• $\sigma$, volatility of underlying asset, the standard deviation of the asset returns

#### $N(d)$ is the cumulative distribution of the standard normal variable Z¶

$$N(d) = \frac{1}{\sqrt{2\pi}} \int_{-\infty}^d e^{-\frac{1}{2}x^2} dx$$

$C(S,t)$ is the value at time $t$ of a call option and $P(S,t)$ is the value at time $t$ of a put option.

The Black-Scholes call formula is given as:

$$C(S,t) = SN(d_1) - Ke^{-r(T - t)} N(d_2)$$

The put formula is given:

$$P(S,t) = Ke^{-r(T - t)}N(-d_2) - SN(-d_1)$$

Where:

$$d_1 = \frac{\ln \left(\frac{S}{K} \right) + \left(r + \frac{\sigma^2}{2} \right)(T - t)}{\sigma \sqrt{T - t}}$$
$$d_2 = d_1 - \sigma \sqrt{T - t} = \frac{\ln \left(\frac{S}{K} \right) + \left(r - \frac{\sigma^2}{2}\right)(T - t)}{\sigma \sqrt{T - t}}$$

## Python Implementation of Black-Scholes formula for non-dividend paying options¶

In :
import numpy as np
import scipy.stats as si
import sympy as sy
from sympy.stats import Normal, cdf
from sympy import init_printing
init_printing()

In :
def euro_vanilla_call(S, K, T, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

call = (S * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))

return call

In :
euro_vanilla_call(50, 100, 1, 0.05, 0.25)

Out:
$\displaystyle 0.027352509369436617$
In :
def euro_vanilla_put(S, K, T, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

put = (K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * si.norm.cdf(-d1, 0.0, 1.0))

return put

In :
euro_vanilla_put(50, 100, 1, 0.05, 0.25)

Out:
$\displaystyle 45.15029495944084$

The next function can be called with 'call' or 'put' for the option parameter to calculate the desired option

In :
def euro_vanilla(S, K, T, r, sigma, option = 'call'):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

if option == 'call':
result = (S * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))
if option == 'put':
result = (K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * si.norm.cdf(-d1, 0.0, 1.0))

return result

In :
euro_vanilla(50, 100, 1, 0.05, 0.25, option = 'put')

Out:
$\displaystyle 45.15029495944084$

Sympy implementation for Exact Results

In :
def euro_call_sym(S, K, T, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

N = Normal('x', 0.0, 1.0)

d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

call = (S * cdf(N)(d1) - K * sy.exp(-r * T) * cdf(N)(d2))

return call

In :
euro_call_sym(50, 100, 1, 0.05, 0.25)

Out:
$\displaystyle - 25 \operatorname{erf}{\left(1.22379436111989 \sqrt{2} \right)} - 22.5614712250357 + 47.5614712250357 \operatorname{erf}{\left(1.34879436111989 \sqrt{2} \right)}$
In :
def euro_put_sym(S, K, T, r, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

N = systats.Normal(0.0, 1.0)

d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

put = (K * sy.exp(-r * T) * N.cdf(-d2) - S * N.cdf(-d1))

return put


Sympy implementation of the above function that enables one to specify a call or put result.

In :
def sym_euro_vanilla(S, K, T, r, sigma, option = 'call'):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#sigma: volatility of underlying asset

N = Normal('x', 0.0, 1.0)

d1 = (sy.ln(S / K) + (r + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

if option == 'call':
result = (S * cdf(N)(d1) - K * sy.exp(-r * T) * cdf(N)(d2))
if option == 'put':
result = (K * sy.exp(-r * T) * cdf(N)(-d2) - S * cdf(N)(-d1))

return result

In :
sym_euro_vanilla(50, 100, 1, 0.05, 0.25, option = 'put')

Out:
$\displaystyle - 25 \operatorname{erf}{\left(1.22379436111989 \sqrt{2} \right)} + 22.5614712250357 + 47.5614712250357 \operatorname{erf}{\left(1.34879436111989 \sqrt{2} \right)}$

## Dividend Paying Black-Scholes Formula¶

For assets that pay dividends, the Black-Scholes formula is rather similar to the non-dividend paying asset formula; however, a new parameter $q$, is added.

• $S$, the spot price of the asset at time $t$
• $T$, the maturity of the option. Time to maturity is defined as $T - t$
• $K$, strike price of the option
• $r$, the risk-free interest rate, assumed to be constant between $t$ and $T$
• $\sigma$, volatility of underlying asset, the standard deviation of the asset returns
• $q$, the dividend rate of the asset. This is assumed to pay dividends at a continuous rate

In this case, the $q$ parameter is now included in $C(S,t)$ and $P(S,t)$.

$$C(S,t) = Se^{-q(T - t)} N(d_1) - Ke^{-r(T - t)} N(d_2)$$
$$P(S,t) = Ke^{-r(T - t)} N(-d_2) - Se^{-q(T - t)} N(-d_1)$$

Then, $d_1$ and $d_2$ are slightly modified to include the continuous dividends

$$d_1 = \frac{ln \left(\frac{S}{K} \right) + \left(r - q + \frac{\sigma^2}{2} \right)(T - t)}{\sigma \sqrt{T - t}}$$
$$d_2 = d_1 - \sigma \sqrt{T - t} = \frac{ln (\frac{S}{K}) + (r - q - \frac{\sigma^2}{2})(T - t)}{\sigma \sqrt{T - t}}$$

## Python Implementation¶

In :
def black_scholes_call_div(S, K, T, r, q, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

call = (S * np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))

return call

In :
def black_scholes_put_div(S, K, T, r, q, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

put = (K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0))

return put


Implementation that can be used to determine the put or call option price depending on specification

In :
def euro_vanilla_dividend(S, K, T, r, q, sigma, option = 'call'):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

d1 = (np.log(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))
d2 = (np.log(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * np.sqrt(T))

if option == 'call':
result = (S * np.exp(-q * T) * si.norm.cdf(d1, 0.0, 1.0) - K * np.exp(-r * T) * si.norm.cdf(d2, 0.0, 1.0))
if option == 'put':
result = (K * np.exp(-r * T) * si.norm.cdf(-d2, 0.0, 1.0) - S * np.exp(-q * T) * si.norm.cdf(-d1, 0.0, 1.0))

return result


Sympy Implementation of Black-Scholes with Dividend-paying asset

In :
def black_scholes_call_div_sym(S, K, T, r, q, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

N = Normal('x', 0.0, 1.0)

d1 = (sy.ln(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

call = S * sy.exp(-q * T) * cdf(N)(d1) - K * sy.exp(-r * T) * cdf(N)(d2)

return call

In :
def black_scholes_call_put_sym(S, K, T, r, q, sigma):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

N = Normal('x', 0.0, 1.0)

d1 = (sy.ln(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

put = K * sy.exp(-r * T) * cdf(N)(-d2) - S * sy.exp(-q * T) * cdf(N)(-d1)

return put


Sympy implementation of pricing a European put or call option depending on specification

In :
def sym_euro_vanilla_dividend(S, K, T, r, q, sigma, option = 'call'):

#S: spot price
#K: strike price
#T: time to maturity
#r: interest rate
#q: rate of continuous dividend paying asset
#sigma: volatility of underlying asset

N = Normal('x', 0.0, 1.0)

d1 = (sy.ln(S / K) + (r - q + 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))
d2 = (sy.ln(S / K) + (r - q - 0.5 * sigma ** 2) * T) / (sigma * sy.sqrt(T))

if option == 'call':
result = S * sy.exp(-q * T) * cdf(N)(d1) - K * sy.exp(-r * T) * cdf(N)(d2)
if option == 'put':
result = K * sy.exp(-r * T) * cdf(N)(-d2) - S * sy.exp(-q * T) * cdf(N)(-d1)

return result