Deriving the Black-Schole Equation

Solving the Black-Schole Equation

Recap

European Vanilla

$\frac{\partial V}{\partial t} + \frac{1}{2} σ^{2} S \frac{\partial^{2} V}{\partial S^{2}} + r S \frac{\partial V}{\partial S} - r V = 0$

V is the option price at the time $t$ . So $V = V (S, t)$
S is the asset spot price
t is the time to expiry (in years)
$σ$ is the asset diffusion term (its stochastic element)
$r$ is the annualized continuously compounded risk-free rate (imaginary friend)

In the case of a European Call Option with no-dividend, the BSE has solution:

$C = S N (d_{1}) - K e^{- r (T - t)} N (d_{2})$

And in the case of a European Put Option with no-dividend, the BSE has solution: $P = K e^{- r (T - t)} N (- d_{2}) - S N (- d_{1})$

where, $d_{1} = \frac{1}{σ \sqrt{T - t}} [l n (\frac{S}{K}) + (r + \frac{σ^{2}}{2}) (T - t)]$

$d_{2} = \frac{1}{σ \sqrt{T - t}} [l n (\frac{S}{K}) + (r - \frac{σ^{2}}{2}) (T - t)] = d 1 - σ \sqrt{T - t}$

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

K is the strike price

European Vanilla with Dividends

And in the case of a dividend (ok … assuming continuous dividend yield):

$C = S e^{- D (T - t)} N (d_{1}) - K e^{- r (T - t)} N (d_{2})$

$P = K e^{- r (T - t)} N (- d_{2}) - S e^{- D (T - t)} N (- d_{1})$

$d_{1} = \frac{1}{σ \sqrt{T - t}} [l n (\frac{S}{K}) + (r - D + \frac{σ^{2}}{2}) (T - t)]$

$d_{2} = \frac{1}{σ \sqrt{T - t}} [l n (\frac{S}{K}) + (r - D - \frac{σ^{2}}{2}) (T - t)] = d 1 - σ \sqrt{T - t}$

The Greeks

Greek	Description	Formula	Call Option	Put Option
Delta	Sensitivity of option value to changes in asset price	$\frac{\partial V}{\partial S}$	$N (d_{1})$	$- N (- d_{!})$
Gamma	Sensitivity of Delta to changes in asset price	$\frac{\partial^{2} V}{\partial S^{2}}$	$\frac{N (d_{1})}{S σ \sqrt{t}}$
Vega	Sensitivity of option value to changes in volatility	$\frac{\partial V}{\partial σ}$	$S N (d_{1}) \sqrt{t}$
Theta	Sensitivity of option value to changes in time	$\frac{\partial V}{\partial t}$
Rho	Sensitivity of option value to change in risk-free rate	$\frac{\partial V}{\partial r}$	$K t e^{- r t} N (d_{2})$	$- K t e^{- r t} N (- d_{2})$

Create a function for numerical computation

Using Python

import numpy as np
import pandas as pd
from scipy.stats import norm
import matplotlib.pyplot as plt 


class option_pricing: 
  
  """
  To price European Style options without dividends
  """
  
  def __init__(self, spot, strike, rate, dte, sigma): 
    
    # assign our variables
    self.spot = spot
    self.strike = strike
    self.rate = rate
    self.dte = dte    # days to expiration (in years)
    self.sigma = sigma
    
    # to avoid zero division, let not allow strike of 0
    if self.strike == 0: 
      raise ZeroDivisionError('The Strike Price cannot be 0')
    else: 
      self._d1_ = (np.log(self.spot / self.strike) + (self.rate + (self.sigma**2 / 2)) * dte) / (self.sigma * self.dte**0.5)
      self._d2_ = self._d1_ - (self.sigma * self.dte**0.5) 
    
    for i in ['callPrice', 'putPrice', 'callDelta', 'putDelta', 'gamma']: 
      self.__dict__[i] = None
      
    [self.callPrice, self.putPrice] = self._price() 
    [self.callDelta, self.putDelta] = self._delta()
    self.gamma = self._gamma()
    
  def _price(self): 
      if self.sigma == 0 or self.dte == 0: 
        call = maximum(0.0, self.spot - self.strike)
        put = maximum(0.0, self.strike - self.spot) 
      else: 
        call = (self.spot * norm.cdf(self._d1_)) - (self.strike * np.e**(- self.rate * self.dte) * norm.cdf(self._d2_))
        put = (self.strike * np.e**(- self.rate * self.dte) * norm.cdf(-self._d2_)) - (self.spot * norm.cdf(-self._d1_))
      return [call, put] 
    
  def _delta(self): 
    if self.sigma == 0 or self.dte == 0: 
      call = 1.0 if self.spot > self.strike else 0.0
      put = -1.0 if self.spot < self.strike else 0.0
    else: 
      call = norm.cdf(self._d1_)
      put = -norm.cdf(-self._d1_)
    return [call, put]
  
  def _gamma(self): 
    return norm.cdf(self._d1_) / (self.spot * self.sigma * self.dte**0.5)

from tabulate import tabulate

option = option_pricing(100, 100, 0, 1, 0.2)

header = ['Call Price', 'Put Price', 'Call Delta', 'Gamma']
table = [[option.callPrice, option.putPrice, option.callDelta, option.gamma]]
print(tabulate(table, header))

  Call Price    Put Price    Call Delta      Gamma
------------  -----------  ------------  ---------
     7.96557      7.96557      0.539828  0.0269914

Retrieving option data using Yahoo finance

import yfinance as yf 

amd = yf.Ticker('AMD')
amd_hist = amd.history(start = '2022-01-01')
options = amd.option_chain('2023-12-15')

from datetime import datetime
dte = (datetime(2023, 12, 15) - datetime.today()).days 

log_returns = np.log(amd_hist['Close'] / amd_hist['Close'].shift(1)).dropna()
historical_vol = log_returns.std() * np.sqrt(dte)

spot = 116; strike = 120; rate = 0.05

amd_opt = option_pricing(spot=spot, strike=strike, rate=rate, dte=dte/365, sigma=historical_vol)

print(f'The BS model for AMD 147 days ahead is {amd_opt.callPrice:0.4f}')

The BS model for AMD 147 days ahead is 2.7995

df = options.calls[(options.calls['strike'] >= 90) & (options.calls['strike'] < 150)]
df.reset_index(drop = True, inplace = True)

df1 = pd.DataFrame({'strike': df['strike'], 'price': df['lastPrice'], 'impl_vol': df['impliedVolatility']})
df1['delta'] = df1['gamma'] = 0.

for i in range(len(df1)): 
  df['delta'].iloc[i] = option_pricing(spot=spot, strike=df['strike'].iloc[i], rate, dte=dte, sigma = df1['impl_vol'].iloc[i]).callDelta

for i in range(len(df1)):
    
    df1['Delta'].iloc[i] = option_pricing(spot,df1['strike'].iloc[i],rate,dte,df1['impl_vol'].iloc[i]).callDelta
    df1['Gamma'].iloc[i] = option_pricing(spot,df1['strike'].iloc[i],rate,dte,df1['impl_vol'].iloc[i]).gamma