Skip to content

Instantly share code, notes, and snippets.

@huangt

huangt/option_greeks.py

Created August 8, 2013 23:51
Show Gist options
  • Save huangt/6189923 to your computer and use it in GitHub Desktop.
Save huangt/6189923 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
option_greeks.py
import sys
import os
import math
"""
Black-Scholes-Merton method to calculate option premium for forward or futures contract
@param callPutFlag
@param F future price
@param X strike price
@param T time to expiration
@param r risk-free interest rate
@param v implied volatility
@return Double option premium
"""
def BlackScholes(callPutFlag, f, x, t, r, v):
d1 = (math.log(f / x) + (math.pow(v,2) / 2) * t) / (v * math.sqrt(t))
d2 = d1 - v * math.sqrt(t)
premium = 0
if (callPutFlag == "c"):
premium = math.exp(-1 * r * t) * (f * CND(d1) - x * CND(d2))
elif (callPutFlag == "p"):
premium = math.exp(-1 * r * t) * (x * CND(-1 * d2) - f * CND(-1 * d1))
return premium
"""
Cumulative Normal Distribution Function using Hart Algorithm
Hart, J.F. et al, 'Computer Approximations', Wiley 1968
@param X
@return Double cumulative normal distribution
"""
def CND(x):
y = abs(x)
if (x > 37):
cnd = 1
return cnd
elif (x < -37):
cnd = 0
return cnd
else:
exponential = math.exp(-1 * y * y / 2)
if (y < 7.07106781186547):
sumA = 0.0352624965998911 * y + 0.700383064443688
sumA = sumA * y + 6.373962203531650
sumA = sumA * y + 33.91286607838300
sumA = sumA * y + 112.0792914978709
sumA = sumA * y + 221.2135961699311
sumA = sumA * y + 220.2068679123761
sumB = 0.0883883476483184 * y + 1.75566716318264
sumB = sumB * y + 16.06417757920695
sumB = sumB * y + 86.78073220294608
sumB = sumB * y + 296.5642487796737
sumB = sumB * y + 637.3336333788311
sumB = sumB * y + 793.8265125199484
sumB = sumB * y + 440.4137358247522
cnd = exponential * sumA / sumB
else:
sumA = y + 0.65
sumA = y + 4 / sumA
sumA = y + 3 / sumA
sumA = y + 2 / sumA
sumA = y + 1 / sumA
cnd = exponential / (sumA * 2.506628274631001)
if (x > 0):
cnd = 1 - cnd
return cnd
def phi(x):
'Cumulative distribution function for the standard normal distribution'
return (1.0 + math.erf(x / math.sqrt(2.0))) / 2.0
"""
calculates delta, b = 0 for futures option
Delta is the option's sensitivity to small changes in the underlying asset price
@param callPutFlag
@param F future price
@param X strike price
@param T time to expiration
@param r risk-free interest rate
@param b cost of carry
@param v implied volatility
@return Double Delta
"""
def Delta(callPutFlag, f, x, t, r, b, v):
delta = 0
d1 = (math.log(f / x) + (b + v * v / 2) * t) / (v * math.sqrt(t))
if (callPutFlag == "c"):
delta = math.exp(-1 * r * t) * CND(d1)
else:
delta = -1 * math.exp(-1 * r * t) * CND(-1 * d1)
return delta
"""
calculates vega, b = 0 for futures option
Vega is the option's sensitivity to a small change in the volatility of the underlying asset.
@param F future price
@param X strike price
@param T time to expiration
@param r risk-free interest rate
@param b cost of carry
@param v implied volatility
@return Double Vega
"""
def Vega(f, x, t, r, b, v):
d1 = (math.log(f / x) + (b + v * v / 2) * t) / (v * math.sqrt(t))
nd1 = (1 / math.sqrt(2 * math.pi)) * math.exp(-1 * d1 * d1 / 2)
vega = f * math.exp((b - r) * t) * (nd1 * math.sqrt(t))
return vega
def onePercentVega(f, x, t, r, b, v):
vega = Vega(f, x, t, r, b, v)
return vega / 100
"""
Gamma is the delta's sensitivity to small changes in the underlying asset price
@param F future price
@param X strike price
@param T time to expiration
@param r risk-free interest rate
@param b cost of carry
@param v implied volatility
@return Double Gamma
"""
def Gamma(f, x, t, r, b, v):
d1 = (math.log(f / x) + (b + v * v / 2) * t) / (v * math.sqrt(t))
nd1 = (1 / math.sqrt(2 * math.pi)) * math.exp(-1 * d1 * d1 / 2)
gamma = (nd1 * math.exp((b - r) * t)) / (f * v * math.sqrt(t))
return gamma
def Theta(callPutFlag, f, x, t, r, b, v):
d1 = (math.log(f / x) + (b + v * v / 2) * t) / (v * math.sqrt(t))
d2 = d1 - v * math.sqrt(t)
nd1 = (1 / math.sqrt(2 * math.pi)) * math.exp(-1 * d1 * d1 / 2)
theta = 0
if (callPutFlag == "p"):
nnd1 = CND(-1 * d1)
nnd2 = CND(-1 * d2)
theta = (-1 * f * math.exp((b - r) * t) * nd1 * v) / (2 * math.sqrt(t)) + (b - r) * f * math.exp((b - r) * t) * nnd1 + r * x * math.exp(-1 * r * t) * nnd2
elif (callPutFlag == "c"):
nd1 = CND(d1)
nd2 = CND(d2)
theta = (-1 * f * math.exp((b - r) * t) * nd1 * v) / (2 * math.sqrt(t)) - (b - r) * f * math.exp((b - r) * t) * nd1 - r * x * math.exp(-1 * r * t) * nd2
# divide by 365 to get theta for a one-day time decay
return theta / 365
"""
Newton-Raphson method to calculate implied volatility for future
b = 0 gives Black(1976) futures option model.
@param callPutFlag
@param F future price
@param X strike price
@param T time to expiration
@param r risk-free interest rate
@param b cost of carry
@param cm market price for option
@param epsilon desired degree of accuracy
@return Double implied volatility
"""
def ImpliedVolatilityNR(callPutFlag, f, x, t, r, b, cm, epsilon):
# Manaster and Koehler seed value (vi)
vi = math.sqrt(abs(math.log(f / x)) * 2 / t)
ci = BlackScholes(callPutFlag, f, x, t, r, vi)
vegai = Vega(f, x, t, r, 0.0, vi)
minDiff = abs(cm - ci)
while (abs(cm - ci) >= epsilon and abs(cm - ci) <= minDiff):
vi = vi - (ci - cm) / vegai
ci = BlackScholes(callPutFlag, f, x, t, r, vi)
vegai = Vega(f, x, t, r, 0.0, vi)
minDiff = abs(cm - ci)
if (abs(cm - ci) < epsilon):
impliedVolatility = vi
else:
impliedVolatility = None
return impliedVolatility
print "CND with Hart Algorithm\n"
print CND(-0.54)
print "CND with python erf()\n"
print phi(-0.54)
premium = BlackScholes("c", 187.93, 195.0, 15.0/365.0, 0, 0.15525)
print "call, f=5569.97, x=6000, t=139, r=0.0231, v=0.18611 \n"
print "premium: " + str(premium)
delta = Delta("c", 187.93, 195.0, 15.0/365.0, 0, 0, 0.15525)
print "delta: " + str(delta)
impliedVol = ImpliedVolatilityNR("c", 187.93, 195.0, 15.0/365.0, 0, 0, 0.36, 0.0001)
print "implied vol: " + str(impliedVol)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment