Plotting Volatility Surface for Options
by Mary Lin, Tom Starke and Michelle Lin
This blog post is a revised edition of Tom’s original blog post with a newer data set. More information, source code & inspiration can be found here. Code for this blog post is in our Github repository.
Options are complex instruments with many moving parts. Specifically, options are contracts that grant the right, but not the obligation to buy or sell an underlying asset at a set price on or before a certain date. The right to buy is called a call option and the right to sell is a put option.
We use the Black Scholes formula for pricing an option. Its inputs include
- current underlying price
- options strike price
- time until expiration as a percent of a year
- implied volatility
- risk free interest rate
Each input affects the the other inputs and consequently the option price some way. In this blog post we provide an exploratory analysis for those interested in visualising the 3D relationship between implied volatility, strike prices and time to expiry.
Implied volatility versus strike price: the volatility smile shows implied volatility increases when option strike prices are further away.
Implied volatility versus time to expiration: The volatility cone shows implied volatility is higher when the option is close to expiry, holding the strike constant.
Below is Python code that shows how to plot the implied volatility surface with both time to expiration and strike price as features.
from pandas_datareader.data import Options from dateutil.parser import parse from datetime import datetime from numpy import * import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm from matplotlib.colors import LogNorm #from implied_vol import BlackScholes from functools import partial from scipy import optimize import numpy as np from scipy.interpolate import griddata
After importing in the libraries, we create a function for Black Scholes.
def CND(X): return float(norm.cdf(X))
def BlackScholes(v,CallPutFlag = 'c',S = 100.,X = 100.,T = 1.,r = 0.01): try: d1 = (log(S/X)+(r+v*v/2.)*T)/(v*sqrt(T)) d2 = d1-v*sqrt(T) if CallPutFlag=='c': return S*CND(d1)-X*exp(-r*T)*CND(d2) else: return X*exp(-r*T)*CND(-d2)-S*CND(-d1) except: return 0
Initially I plotted the features as data points using a generic 3D plot. Note this block of code will plot only the data points as opposed to a surface plot.
def plot3D(X,Y,Z): fig = plt.figure() ax = Axes3D(fig, azim = -29, elev = 50) ax.plot(X,Y,Z,'o') plt.xlabel("expiry") plt.ylabel("strike") plt.show()
The code for calculating implied volatility:
def calc_impl_vol(price = 5., right = 'c', underlying = 100., strike = 100., time = 1., rf = 0.01, inc = 0.001): f = lambda x: BlackScholes(x,CallPutFlag=right,S=underlying,X=strike,T=time,r=rf)-price return optimize.brentq(f,0.,5.)
Next we want to retrieve the numbers for the option’s features or inputs.
We append to a list called “vals” the time to expiry, strike price and implied volatility for an option.
These can now be referenced as vals, vals and vals respectively.
def get_surf(ticker): q = Options(ticker, 'yahoo').get_all_data() q.reset_index(inplace=True) q.set_index('Symbol',inplace=True) vals =  print(q.head()) for index, row in q.iterrows(): if row['Type'] == 'put': underlying = float(row['Underlying_Price']) price = (float(row['Ask'])+float(row['Bid']))/2.0 expd = (row['Expiry'] - datetime.now()).days exps = (row['Expiry'] - datetime.now()).seconds exp = (expd*24.*3600. + exps) / (365.*24.*3600.) try: impl = calc_impl_vol(price, 'p', underlying, float(row['Strike']), exp) vals.append([exp,float(row['Strike']),impl]) except: pass vals = array(vals).T mesh_plot2(vals,vals,vals) #if you want to call the 3d plot above use this code instead: #plot3D(vals,vals,vals) #if you want to call both plots use this code instead: #combine_plots(vals,vals,vals)
meshgrid() is part of the Numpy library. It creates a rectangular grid out of an array of x and y values. Imagine it as overlaying a grid on top of the data points. We are able to change the resolution using the linspace() command, which returns evenly spaced numbers over a specified interval. The higher the number in linspace, the higher the resolution.
To add the third dimension z, we use griddata(), which is part of the SciPy library. griddata() interpolates between the mesh data points. That means, it will determine the value of z using interpolation between the data points and the overlaying ‘grid’.
To optimise the volatility surface visualisation, we can do two things: 1) smooth the volatility surface, and 2) add the data points on top of the surface plot.
To smooth the surface, I re-adjusted the resolution and applied a ‘linear’ interpolation method in griddata(). What that does is use piecewise linear interpolation to connect the data points and determine the value of the third dimension. Another option is ‘cubic’, which uses cubic splines to determine the value of the third dimension. The function for make_surf() now looks like:
def make_surf(X,Y,Z): XX,YY = meshgrid(linspace(min(X),max(X),230),linspace(min(Y),max(Y),230)) ZZ = griddata(array([X,Y]).T,array(Z),(XX,YY), method='linear') return XX,YY,ZZ
To plot the data points from the generic 3D plot on top of the surface plot, I created a new function combine_plots() that plots both the data points and the surface plot.
def plot3D(X,Y,Z,fig,ax): ax.plot(X,Y,Z,'o', color = 'pink') plt.xlabel("expiry") plt.ylabel("strike") def mesh_plot2(X,Y,Z,fig,ax): XX,YY,ZZ = make_surf(X,Y,Z) ax.plot_surface(XX,YY,ZZ, color = 'white') ax.contour(XX,YY,ZZ) plt.xlabel("expiry") plt.ylabel("strike") def combine_plots(X,Y,Z): fig = plt.figure() ax = Axes3D(fig, azim = -29, elev = 50) mesh_plot2(X,Y,Z,fig,ax) plot3D(X,Y,Z,fig,ax) plt.show()
In the example below, we create the surface plot for Apple (AAPL).