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 [6]:
import numpy as np
import scipy.stats as si
import sympy as sy
import sympy.statistics as systats
In [7]:
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 [8]:
euro_vanilla_call(50, 100, 1, 0.05, 0.25)
Out[8]:
0.027352509369436617
In [20]:
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 [21]:
euro_vanilla_put(50, 100, 1, 0.05, 0.25)
Out[21]:
45.150294959440842

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

In [9]:
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 [10]:
euro_vanilla(50, 100, 1, 0.05, 0.25, option = 'put')
Out[10]:
45.150294959440842

Sympy implementation for Exact Results

In [28]:
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 = 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))
    
    call = (S * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2))
    
    return call
In [29]:
euro_call_sym(50, 100, 1, 0.05, 0.25)
Out[29]:
-25.0*erf(1.22379436111989*sqrt(2)) - 22.5614712250357 + 47.5614712250357*erf(1.34879436111989*sqrt(2))
In [16]:
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 [30]:
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 = 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))
    
    if option == 'call':
        result = (S * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2))
    if option == 'put':
        result = (K * sy.exp(-r * T) * N.cdf(-d2) - S * N.cdf(-d1))
        
    return result
In [31]:
sym_euro_vanilla(50, 100, 1, 0.05, 0.25, option = 'put')
Out[31]:
-25.0*erf(1.22379436111989*sqrt(2)) + 22.5614712250357 + 47.5614712250357*erf(1.34879436111989*sqrt(2))

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 [32]:
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 [33]:
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 [34]:
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 [35]:
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 = systats.Normal(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) * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2)
    
    return call
In [37]:
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 = systats.Normal(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) * N.cdf(-d2) - S * sy.exp(-q * T) * N.cdf(-d1)
    
    return put

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

In [38]:
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 = systats.Normal(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) * N.cdf(d1) - K * sy.exp(-r * T) * N.cdf(d2)
    if option == 'put':
        result = K * sy.exp(-r * T) * N.cdf(-d2) - S * sy.exp(-q * T) * N.cdf(-d1)
        
    return result