Black-Scholes formula and greeks: Python implementation

• Black-Scholes formula and greeks;
• implementation of option strategies;
• Examples of study of portfolio greeks;

Hello, all!

In this post, I´ll share with you a plug and play implementation in Python for the Black-Scholes formula in addition to the greeks for vanilla European calls and puts. The implementation is simple, given that all results are obtained directly from the implementation of closed analytical formulas.

Let´s first remind the Black-Scholes formula, for calls:

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

and for puts:

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

with:
$d_1=\frac{ln(S/K) + (r + \sigma^2 /2)(T-t)}{\sigma \sqrt{T-t}}$
$d_2 = d_1-\sigma \sqrt{T-t}$
and: $N(x)=\frac{1}{\sqrt{2 \pi}} \int_{- \infty}^{x} exp(\frac{z^2}{2})dz$
where $(T-t)$ is the time to maturity, $K$ is the strike price, $S$ is the asset price at time $t$, $r$ is the risk free interest rate, and $\sigma$ is the volatility.
The greeks (discussed in this post) are:
Delta:
$\Delta = \frac{\partial C}{\partial S}$
Call:
$\Delta = N(d_1)$
Put:
$\Delta = -N(-d_1)$
Gamma:
$\Gamma = \frac{\partial^2 C}{\partial S^2}$
Call & Put:
$\Gamma = \frac{N'(d_1)}{S \sigma \sqrt{T-t} }$
Vega:
$Vega = \frac{\partial C}{\partial \sigma}$
Call & Put:
$Vega = S N'(d_1) \sqrt{T-t}$
Rho:
$\rho = \frac{\partial C}{\partial r}$
Call:
$\rho = K(T-t)e^{-r(T-t)}N(d_2)$
Put:
$\rho = -K(T-t)e^{-r(T-t)}N(-d_2)$
Theta:
$\theta = \frac{\partial C}{\partial t}$
Call:
$\theta = \frac{S \sigma N'(d_1)}{\sigma \sqrt{T-t} } - rKe^{-r(T-t)} N(d_2)$
Put:
$\theta = \frac{S \sigma N'(d_1)}{\sigma \sqrt{T-t} } + rKe^{-r(T-t)} N(-d_2)$

Implementation:

First, we load the packages needed for this implementation:

	import scipy.stats as scistat
	import numpy as np

view raw load.py hosted with ❤ by GitHub

the values of $d_1$ and $d_2$ are common for both calls and puts, so we define:

	def d1(S, K, r, sigma, T):
	return (np.log(S/K) + (r+sigma*sigma/2)*T)/(sigma*np.sqrt(T))
	def d2(S, K, r, sigma, T):
	return d1(S, K, r, sigma, T) - sigma*np.sqrt(T)

view raw BS_d1_d2.py hosted with ❤ by GitHub

Now, we define a class for Calls including functions for the prices and the greeks:

	'''
	Input parameters:
	S -> asset price
	K -> strike price
	r -> interest rate
	sigma -> volatility
	T -> time to maturity
	'''
	class Call:
	def Price(S, K, r, sigma, T):
	return np.maximum(S - K, 0) if T==0 else S*scistat.norm.cdf(d1(S, K, r, sigma, T)) - K*np.exp(-r*T)*scistat.norm.cdf(d2(S, K, r, sigma, T))
	def Delta(S, K, r, sigma, T):
	return scistat.norm.cdf(d1(S, K, r, sigma, T))
	def Gamma(S, K, r, sigma, T):
	return scistat.norm.pdf(d1(S, K, r, sigma, T))/(S*sigma*np.sqrt(T))
	def Vega(S, K, r, sigma, T):
	return S*scistat.norm.pdf(d1(S, K, r, sigma, T))*np.sqrt(T)
	def Theta(S, K, r, sigma, T):
	aux1 = -S*scistat.norm.pdf(d1(S, K, r, sigma, T))*sigma/(2*np.sqrt(T))
	aux2 = -r*K*np.exp(-r*T)*scistat.norm.cdf(d2(S, K, r, sigma, T))
	return aux1+aux2
	def Rho(S, K, r, sigma, T):
	return K*T*np.exp(-r*T)*scistat.norm.cdf(d2(S, K, r, sigma, T))
	'''
	Range calculations
	'''
	def Get_range_value(Smin, Smax, Sstep, K, r, sigma, T, num_curves, value="Price"):
	vec = np.linspace(Smin, Smax, (Smax - Smin) / Sstep)
	vecT = np.linspace(0,T,num_curves, endpoint=True)
	if value=="Price":
	return vec,vecT, [[Call.Price(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Delta":
	return vec, vecT, [[Call.Delta(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Gamma":
	return vec, vecT, [[Call.Gamma(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Vega":
	return vec, vecT, [[Call.Vega(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Theta":
	return vec, vecT, [[Call.Theta(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Rho":
	return vec, vecT, [[Call.Rho(S, K, r, sigma, t) for S in vec] for t in vecT]

view raw BS_Calls.py hosted with ❤ by GitHub

another class for the puts and corresponding greeks:

	'''
	Input parameters:
	S -> asset price
	K -> strike price
	r -> interest rate
	sigma -> volatility
	T -> time to maturity
	'''
	class Put:
	def Price(S, K, r, sigma, T):
	return np.maximum(K-S,0) if T==0 else K*np.exp(-r*T)*scistat.norm.cdf(-1*d2(S, K, r, sigma, T)) - S*scistat.norm.cdf(-1*d1(S, K, r, sigma, T))
	def Delta(S, K, r, sigma, T):
	return scistat.norm.cdf(d1(S, K, r, sigma, T)) - 1
	def Gamma(S, K, r, sigma, T):
	return scistat.norm.pdf(d1(S, K, r, sigma, T))/(S*sigma*np.sqrt(T))
	def Vega(S, K, r, sigma, T):
	return S*scistat.norm.pdf(d1(S, K, r, sigma, T))*np.sqrt(T)
	def Theta(S, K, r, sigma, T):
	aux1 = -S*scistat.norm.pdf(d1(S, K, r, sigma, T))*sigma/(2*np.sqrt(T))
	aux2 = r*K*np.exp(-r*T)*scistat.norm.cdf(-1*d2(S, K, r, sigma, T))
	return aux1+aux2
	def Rho(S, K, r, sigma, T):
	return -K*T*np.exp(-r*T)*scistat.norm.cdf(-1*d2(S, K, r, sigma, T))
	def Get_range_value(Smin, Smax, Sstep, K, r, sigma, T, num_curves, value="Price"):
	vec = np.linspace(Smin, Smax, (Smax - Smin) / Sstep)
	vecT = np.linspace(0,T,num_curves, endpoint=True)
	if value=="Price":
	return vec,vecT, [[Put.Price(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Delta":
	return vec, vecT, [[Put.Delta(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Gamma":
	return vec, vecT, [[Put.Gamma(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Vega":
	return vec, vecT, [[Put.Vega(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Theta":
	return vec, vecT, [[Put.Theta(S, K, r, sigma, t) for S in vec] for t in vecT]
	elif value=="Rho":
	return vec, vecT, [[Put.Rho(S, K, r, sigma, t) for S in vec] for t in vecT]

view raw BS_Puts.py hosted with ❤ by GitHub

The next function simply receives information of the option type (call/put), the parameter to calculate (price, delta, gamma, and so on) in addition to the numerical input data:

	'''
	Input parameters:
	num_curves -> number of different "times to maturity" to be considered
	opt -> string "Call" / "Put"
	Smin/Smax/Sstep -> range of spot price
	T -> maturity
	value -> string ["Price", "Delta", "Gamma", "Vega", "Rho", "Theta"]
	'''
	def GetValue(opt, Smin, Smax, Sstep, K, r, sigma, T, num_curves=1, value="Price"):
	if opt == "Call":
	obj_opt = Call
	elif opt == "Put":
	obj_opt= Put
	return obj_opt.Get_range_value(Smin, Smax, Sstep, K, r, sigma, T, num_curves, value)

view raw BS_getValue.py hosted with ❤ by GitHub

OK, now we can test the implementation by using, for instance:

	def Plot_Data(x, vec_y, vec_label, ylabel="Option price", xlabel="Spot price"):
	for t in range(len(vec_label)):
	plt.plot(x, vec_y[t], label="T-t="+str(vec_label[t]))
	plt.xlabel(xlabel)
	plt.ylabel(ylabel)
	plt.legend(loc="best")
	plt.show()
	Smin = 10
	Smax = 120
	Sstep = 1
	K = 60
	r = 0.05
	sigma = 0.5
	T = 1
	num_curves = 3
	get = ["Price", "Delta", "Gamma", "Theta", "Rho", "Vega"]
	for g in get:
	x, vecT, vec_y = BS.GetValue("Call",Smin, Smax, Sstep, K, r, sigma, T, num_curves,g)
	Plot_Data(x, vec_y,vecT, g)

view raw BS_test_visualize.py hosted with ❤ by GitHub

which produces, for example:

By using this implementation, one can easily study strategies, such as butterflies and strangles:

	K1 = 50
	K2 = 60
	K3 = 70
	def Butterfly (opt, Smin, Smax, Sstep, K1, K2, K3, r, sigma, T, num_curves, value="Price"):
	x, vecT, opt1 = BS.GetValue(opt, Smin, Smax, Sstep, K1, r, sigma, T, num_curves, value)
	x, vecT, opt2 = BS.GetValue(opt, Smin, Smax, Sstep, K2, r, sigma, T, num_curves, value)
	x, vecT, opt3 = BS.GetValue(opt, Smin, Smax, Sstep, K3, r, sigma, T, num_curves, value)
	butterfly = [np.array(opt1[t])- 2*np.array(opt2[t])+np.array(opt3[t]) for t in range(len(opt1))]
	return x, butterfly, vecT
	def Strangle (Smin, Smax, Sstep, K1, K2, r, sigma, T, num_curves):
	x, vecT, call = BS.GetValue("Call",Smin, Smax, Sstep, K1, r, sigma, T, num_curves, value)
	x, vecT, put = BS.GetValue("Put",Smin, Smax, Sstep, K2, r, sigma, T, num_curves, value)
	strangle = [np.array(call[t]) + np.array(put[t]) for t in range(len(call))]
	return x, strangle, vecT
	x, butterfly_put, vecT = Butterfly ("Put", Smin, Smax, Sstep, K1, K2, K3, r, sigma, T, num_curves)
	Plot_Data(x, butterfly_put, vecT, "Put (Butterfly)")
	K1 = 50
	K2 = 90
	x, strangle, vecT = Strangle (Smin, Smax, Sstep, K1, K2, r, sigma, T, num_curves)
	Plot_Data(x, strangle, vecT, "Strangle")

view raw BS_strategies.py hosted with ❤ by GitHub

which produces:

I hope you enjoyed this post! Good calculations for everyone!
May the Force be with you.

#=================================

Diogo de Moura Pedroso

LinkedIn: www.linkedin.com/in/diogomourapedroso