Greeks: option sensitivies, formula proofs and Python scripts

Option greeks: formula proofs and python implementation – Part 2

This documents is the second part of a general overview of vanilla options partial sensitivities (option greeks). In a first article we had covered 1st generation greeks, their formula, mathematical proof, and suggested an implementation in Python.

In this post we add some second order greeks such as Vanna and Charm.

Vanna

Vanna is the option’s Delta sensitivity to small changes in the underlying volatility. This measure is actually tantamount to sensitivity of the option’s Vega to small changes in the underlying asset price.

Formula

Let’s remind the Black-Scholes-Merton formula for Vega:

(1) $\begin{equation*} \nu_{c} = \nu_{p} = e^{-qT} N(d_{1}) \end{equation*}$

The call/put option Vanna will be:

(2) $\begin{equation*} \frac{\partial^2 c}{\partial S \partial \sigma} = \frac{\partial^2 p}{\partial S \partial \sigma} = N'(d_{1}) \frac{-e^{-qT} d_{2}}{\sigma} \end{equation*}$

Proof

(3) $\begin{equation*} \begin{split} \begin{aligned} \frac{\partial^2 c}{\partial S \partial \sigma} = \frac{\partial\Delta}{\partial\sigma} = {} & \frac{\partial\nu}{\partial S} \\ & {=} \frac{\partial (e^{-qT} N(d_{1}))}{\partial \sigma} \\ & {=} e^{-qT} \frac{\partial N(d_{1})}{\partial d_{1}} \frac{\partial d_{1}}{\partial \sigma} \\ & {=} e^{-qT} N'(d_{1}) \frac{\partial d_{1}}{\partial \sigma} \\ & {=} e^{-qT} N'(d_{1}) \frac{\partial \frac{\ln \frac{S}{X} + (r-q + \frac{\sigma^2}{2})T}{\sigma \sqrt T}}{\partial \sigma} \\ & {=} e^{-qT} N'(d_{1}) \frac {\sigma \sqrt T - d_{1}}{\sigma} \\ & {=} N'(d_{1}) \frac{-e^{-qT} d_{2}}{\sigma} \end{aligned} \end{split} \end{equation*}$

Code

Below you have the python script for Vanna calculation for a 1% change in the unerlying asset volatility.

import numpy as np
from math import sqrt, pi,log, e
from enum import Enum
import scipy.stats as stat
from scipy.stats import norm
import time

class BSMerton:
    def __init__(self, args):
        self.Type = int(args[0])                # 1 for a Call, - 1 for a put
        self.S = float(args[1])                 # Underlying asset price
        self.K = float(args[2])                 # Option strike K
        self.r = float(args[3])                 # Continuous risk fee rate
        self.q = float(args[4])                 # Dividend continuous rate
        self.T = float(args[5]) / 365.0         # Compute time to expiry
        self.sigma = float(args[6])             # Underlying volatility
        self.sigmaT = self.sigma * self.T ** 0.5# sigma*T for reusability
        self.d1 = (log(self.S / self.K) + \
                   (self.r - self.q + 0.5 * (self.sigma ** 2)) \
                   * self.T) / self.sigmaT
        self.d2 = self.d1 - self.sigmaT
        [self.Premium] = self.premium()
        [self.Delta] = self.delta()
        [self.Theta] = self.theta()
        [self.Rho] = self.rho()
        [self.Vega] = self.vega()
        [self.Gamma] = self.gamma()
        [self.Phi] = self.phi()
        [self.Charm] = self.dDeltadTime()
        [self.Vanna] = self.dDeltadVol()

    def premium(self):
        tmpprem = self.Type * (self.S * e ** (-self.q * self.T) * norm.cdf(self.Type * self.d1) - \
				self.K * e ** (-self.r * self.T) * norm.cdf(self.Type * self.d2))
        return [tmpprem]

    ############################################
    ############ 1st order greeks ##############
    ############################################

    def delta(self):
        dfq = e ** (-self.q * self.T)
        if self.Type == 1:
            return [dfq * norm.cdf(self.d1)]
        else:
            return [dfq * (norm.cdf(self.d1) - 1)]

    # Vega for 1% change in vol
    def vega(self):
		return [0.01 * self.S * e ** (-self.q * self.T) * \
          norm.pdf(self.d1) * self.T ** 0.5]

    # Theta for 1 day change
    def theta(self):
        df = e ** -(self.r * self.T)
        dfq = e ** (-self.q * self.T)
        tmptheta = (1.0 / 365.0) \
            * (-0.5 * self.S * dfq * norm.pdf(self.d1) * \
               self.sigma / (self.T ** 0.5) + \
	        self.Type * (self.q * self.S * dfq * norm.cdf(self.Type * self.d1) \
            - self.r * self.K * df * norm.cdf(self.Type * self.d2)))
        return [tmptheta]

    def rho(self):
	    df = e ** -(self.r * self.T)
	    return [self.Type * self.K * self.T * df * 0.01 * norm.cdf(self.Type * self.d2)]

    def phi(self):
	    return [0.01* -self.Type * self.T * self.S * \
             e ** (-self.q * self.T) * norm.cdf(self.Type * self.d1)]

    ############################################
    ############ 2nd order greeks ##############
    ############################################

    def gamma(self):
	    return [e ** (-self.q * self.T) * norm.pdf(self.d1) / (self.S * self.sigmaT)]

    # Charm for 1 day change
    def dDeltadTime(self):
        dfq = e ** (-self.q * self.T)
        if self.Type == 1:
            return [(1.0 / 365.0) * -dfq * (norm.pdf(self.d1) * ((self.r-self.q) / (self.sigmaT) - self.d2 / (2 * self.T)) \
                            + (-self.q) * norm.cdf(self.d1))]
        else:
            return [(1.0 / 365.0) * -dfq * (norm.pdf(self.d1) * ((self.r-self.q) / (self.sigmaT) - self.d2 / (2 * self.T)) \
                            + self.q * norm.cdf(-self.d1))]

    # Vanna for 1% change in vol
    def dDeltadVol(self):
        return [0.01 * -e ** (-self.q * self.T) * self.d2 / self.sigma * norm.pdf(self.d1)]

    # Vomma
    def dVegadVol(self):
        return [0.01 * -e ** (-self.q * self.T) * self.d2 / self.sigma * norm.pdf(self.d1)]

import numpy as np

from math import sqrt, pi,log, e

from enum import Enum

import scipy.stats as stat

from scipy.stats import norm

import time

class BSMerton:

def __init__(self, args):

self.Type = int(args[0]) # 1 for a Call, - 1 for a put

self.S = float(args[1]) # Underlying asset price

self.K = float(args[2]) # Option strike K

self.r = float(args[3]) # Continuous risk fee rate

self.q = float(args[4]) # Dividend continuous rate

self.T = float(args[5]) / 365.0 # Compute time to expiry

self.sigma = float(args[6]) # Underlying volatility

self.sigmaT = self.sigma * self.T ** 0.5# sigma*T for reusability

self.d1 = (log(self.S / self.K) + \

(self.r - self.q + 0.5 * (self.sigma ** 2)) \

* self.T) / self.sigmaT

self.d2 = self.d1 - self.sigmaT

[self.Premium] = self.premium()

[self.Delta] = self.delta()

[self.Theta] = self.theta()

[self.Rho] = self.rho()

[self.Vega] = self.vega()

[self.Gamma] = self.gamma()

[self.Phi] = self.phi()

[self.Charm] = self.dDeltadTime()

[self.Vanna] = self.dDeltadVol()

def premium(self):

tmpprem = self.Type * (self.S * e ** (-self.q * self.T) * norm.cdf(self.Type * self.d1) - \

self.K * e ** (-self.r * self.T) * norm.cdf(self.Type * self.d2))

return [tmpprem]

############################################

############ 1st order greeks ##############

############################################

def delta(self):

dfq = e ** (-self.q * self.T)

if self.Type == 1:

return [dfq * norm.cdf(self.d1)]

else:

return [dfq * (norm.cdf(self.d1) - 1)]

# Vega for 1% change in vol

def vega(self):

return [0.01 * self.S * e ** (-self.q * self.T) * \

norm.pdf(self.d1) * self.T ** 0.5]

# Theta for 1 day change

def theta(self):

df = e ** -(self.r * self.T)

dfq = e ** (-self.q * self.T)

tmptheta = (1.0 / 365.0) \

* (-0.5 * self.S * dfq * norm.pdf(self.d1) * \

self.sigma / (self.T ** 0.5) + \

self.Type * (self.q * self.S * dfq * norm.cdf(self.Type * self.d1) \

- self.r * self.K * df * norm.cdf(self.Type * self.d2)))

return [tmptheta]

def rho(self):

df = e ** -(self.r * self.T)

return [self.Type * self.K * self.T * df * 0.01 * norm.cdf(self.Type * self.d2)]

def phi(self):

return [0.01* -self.Type * self.T * self.S * \

e ** (-self.q * self.T) * norm.cdf(self.Type * self.d1)]

############################################

############ 2nd order greeks ##############

############################################

def gamma(self):

return [e ** (-self.q * self.T) * norm.pdf(self.d1) / (self.S * self.sigmaT)]

# Charm for 1 day change

def dDeltadTime(self):

dfq = e ** (-self.q * self.T)

if self.Type == 1:

return [(1.0 / 365.0) * -dfq * (norm.pdf(self.d1) * ((self.r-self.q) / (self.sigmaT) - self.d2 / (2 * self.T)) \

+ (-self.q) * norm.cdf(self.d1))]

else:

return [(1.0 / 365.0) * -dfq * (norm.pdf(self.d1) * ((self.r-self.q) / (self.sigmaT) - self.d2 / (2 * self.T)) \

+ self.q * norm.cdf(-self.d1))]

# Vanna for 1% change in vol

def dDeltadVol(self):

return [0.01 * -e ** (-self.q * self.T) * self.d2 / self.sigma * norm.pdf(self.d1)]

# Vomma

def dVegadVol(self):

return [0.01 * -e ** (-self.q * self.T) * self.d2 / self.sigma * norm.pdf(self.d1)]

Here is the piece of code that you can use to calculate and chart the Vanna surface displayed above (the python file that contains the Delta calculation above is called “OptionsAnalytics.py”). This sample was already available in the first part of this overview of BS greeks.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math
from matplotlib import cm
import OptionsAnalytics
from OptionsAnalytics import BSMerton

# Option parameters
sigma = 0.12        # Flat volatility
strike = 105.0      # Fixed strike
epsilon = 0.4       # The % on the left/right of Strike. 
                    # Asset prices are centered around Spot ("ATM Spot")
shortexpiry = 30    # Shortest expiry in days
longexpiry = 720    # Longest expiry in days
riskfree = 0.00     # Continuous risk free rate
divrate = 0.00      # Continuous div rate

# Grid definition
dx, dy = 40, 40     # Steps throughout asset price and expiries axis

# xx: Asset price axis, yy: expiry axis, zz: greek axis
xx, yy = np.meshgrid(np.linspace(strike*(1-epsilon), (1+epsilon)*strike, dx), \
    np.linspace(shortexpiry, longexpiry, dy))
print "Calculating greeks ..."
zz = np.array([BSMerton([1,x,strike,riskfree,divrate,y,sigma]).Vanna for
               x,y in zip(np.ravel(xx), np.ravel(yy))])
zz = zz.reshape(xx.shape)

# Plot greek surface
print "Plotting surface ..."
fig = plt.figure()
fig.suptitle('Call Delta',fontsize=20)
ax = fig.gca(projection='3d')
surf = ax.plot_surface(xx, yy, zz,rstride=1, cstride=1,alpha=0.75,cmap=cm.RdYlBu)
ax.set_xlabel('Asset price')
ax.set_ylabel('Expiry')
ax.set_zlabel('Delta')

# Plot 3D contour
zzlevels = np.linspace(zz.min(),zz.max(),num=8,endpoint=True)
xxlevels = np.linspace(xx.min(),xx.max(),num=8,endpoint=True)
yylevels = np.linspace(yy.min(),yy.max(),num=8,endpoint=True)
cset = ax.contourf(xx, yy, zz, zzlevels, zdir='z',offset=zz.min(),
                   cmap=cm.RdYlBu,linestyles='dashed')
cset = ax.contourf(xx, yy, zz, xxlevels, zdir='x',offset=xx.min(),
                   cmap=cm.RdYlBu,linestyles='dashed')
cset = ax.contourf(xx, yy, zz, yylevels, zdir='y',offset=yy.max(),
                   cmap=cm.RdYlBu,linestyles='dashed')

for c in cset.collections:
    c.set_dashes([(0, (2.0, 2.0))]) # Dash contours

plt.clabel(cset,fontsize=10, inline=1)

ax.set_xlim(xx.min(),xx.max())
ax.set_ylim(yy.min(),yy.max())
ax.set_zlim(zz.min(),zz.max())

#ax.relim()
#ax.autoscale_view(True,True,True)

# Colorbar
colbar = plt.colorbar(surf, shrink=1.0, extend='both', aspect = 10)
l,b,w,h = plt.gca().get_position().bounds
ll,bb,ww,hh = colbar.ax.get_position().bounds
colbar.ax.set_position([ll, b+0.1*h, ww, h*0.8])

# Show chart
plt.show()

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

import math

from matplotlib import cm

import OptionsAnalytics

from OptionsAnalytics import BSMerton

# Option parameters

sigma = 0.12 # Flat volatility

strike = 105.0 # Fixed strike

epsilon = 0.4 # The % on the left/right of Strike.

# Asset prices are centered around Spot ("ATM Spot")

shortexpiry = 30 # Shortest expiry in days

longexpiry = 720 # Longest expiry in days

riskfree = 0.00 # Continuous risk free rate

divrate = 0.00 # Continuous div rate

# Grid definition

dx, dy = 40, 40 # Steps throughout asset price and expiries axis

# xx: Asset price axis, yy: expiry axis, zz: greek axis

xx, yy = np.meshgrid(np.linspace(strike*(1-epsilon), (1+epsilon)*strike, dx), \

np.linspace(shortexpiry, longexpiry, dy))

print "Calculating greeks ..."

zz = np.array([BSMerton([1,x,strike,riskfree,divrate,y,sigma]).Vanna for

x,y in zip(np.ravel(xx), np.ravel(yy))])

zz = zz.reshape(xx.shape)

# Plot greek surface

print "Plotting surface ..."

fig = plt.figure()

fig.suptitle('Call Delta',fontsize=20)

ax = fig.gca(projection='3d')

surf = ax.plot_surface(xx, yy, zz,rstride=1, cstride=1,alpha=0.75,cmap=cm.RdYlBu)

ax.set_xlabel('Asset price')

ax.set_ylabel('Expiry')

ax.set_zlabel('Delta')

# Plot 3D contour

zzlevels = np.linspace(zz.min(),zz.max(),num=8,endpoint=True)

xxlevels = np.linspace(xx.min(),xx.max(),num=8,endpoint=True)

yylevels = np.linspace(yy.min(),yy.max(),num=8,endpoint=True)

cset = ax.contourf(xx, yy, zz, zzlevels, zdir='z',offset=zz.min(),

cmap=cm.RdYlBu,linestyles='dashed')

cset = ax.contourf(xx, yy, zz, xxlevels, zdir='x',offset=xx.min(),

cmap=cm.RdYlBu,linestyles='dashed')

cset = ax.contourf(xx, yy, zz, yylevels, zdir='y',offset=yy.max(),

cmap=cm.RdYlBu,linestyles='dashed')

for c in cset.collections:

c.set_dashes([(0, (2.0, 2.0))]) # Dash contours

plt.clabel(cset,fontsize=10, inline=1)

ax.set_xlim(xx.min(),xx.max())

ax.set_ylim(yy.min(),yy.max())

ax.set_zlim(zz.min(),zz.max())

#ax.relim()

#ax.autoscale_view(True,True,True)

# Colorbar

colbar = plt.colorbar(surf, shrink=1.0, extend='both', aspect = 10)

l,b,w,h = plt.gca().get_position().bounds

ll,bb,ww,hh = colbar.ax.get_position().bounds

colbar.ax.set_position([ll, b+0.1*h, ww, h*0.8])

# Show chart

plt.show()

Shape

When we run our script, this is the Vanna shape that we obtain:

Options Greeks Python Vanna

This is the shape of the Charm sensitivity (Delta variation for a 1-day change in the time to expiry):

Options Greeks Python Charm

4 Comments

Karan Pillai

August 31, 2016

Thanks for the proof and codes, super helpful! I was wondering if you’d ever consider doing the same for Vomma and DvegaDtime?

Tim

September 1, 2015

Hi – nice website! Just trying to follow some of your derivations and interested to know what you mean by the N'() operator (i.e. what does the hash mark ‘ mean) and how does this differ from the regular pmf of a standard normal distribution, that appears as just N() in your notation (no dash).
thanks
Tim

- Team Smile of Thales
  
  September 2, 2015
  
  Hi Tim!
  Thanks for your comment.
  N() stands for the cumulative distribution function (CDF) of the standard normal distribution N(0,1).
  
  (1) $\begin{equation*} N(x)= \frac {1}{\sqrt{2\pi}} \int_{-\propto}^x \! {e ^{\frac{-t^2}{2}}} \, \mathrm{d}t \end{equation*}$
  
  N'() is the probability density function of the standard normal distribution N(0,1).
  
  (2) $\begin{equation*} N'(x)= \frac {e ^{\frac{-1}{2}x^2}}{\sqrt{2\pi}} \end{equation*}$
  
  Hope this clarifies
  Best Regards
  Team Smile of Thales
  
  - Tim
    
    September 8, 2015
    
    Thanks, that’s great. (I worked through the derivations by hand, taking the derivatives myself and eventually realised that’s what they must have been). Good luck with the site.

The smile of Thales

A journey into Finance and Computation

Recent posts

Archives

Greeks: option sensitivies, formula proofs and Python scripts – Part 2

Option greeks: formula proofs and python implementation – Part 2

Vanna

Formula

Proof

Code

Shape

4 Comments

Karan Pillai

Tim

Team Smile of Thales

Tim

Leave a Reply Cancel reply