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)
The call/put option Vanna will be:
(2)
Proof
(3)
Code
Below you have the python script for Vanna calculation for a 1% change in the unerlying asset volatility.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 |
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 |
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:
This is the shape of the Charm sensitivity (Delta variation for a 1-day change in the time to expiry):
Karan Pillai
Thanks for the proof and codes, super helpful! I was wondering if you’d ever consider doing the same for Vomma and DvegaDtime?
Tim
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
Hi Tim!
Thanks for your comment.
N() stands for the cumulative distribution function (CDF) of the standard normal distribution N(0,1).
(1)
N'() is the probability density function of the standard normal distribution N(0,1).
(2)
Hope this clarifies
Best Regards
Team Smile of Thales
Tim
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.