When you value an option, the variables in the model (e.g. stock price, time to expiration) are known except volatility, which is an estimate.

If models were completely correct, the volatility surface across strike prices and maturities would be flat. In practice.

This is not the case as you’ll see.

Let's go!

Build a volatility surface with Python

A volatility surface plots the level of implied volatility in 3D space. The days to expiration are on the X-axis, the strike price is on the Y-axis, and implied volatility is on the Z-axis.

Implied volatility is the market’s expectations of volatility over the life of an option. To find implied volatility, you need three things: the market price of the option, a pricing model, and a root finder. You can then find the volatility that sets the price from the model equal to the price of the market with the root finder. "The volatility implied by the market." The volatility surface is found by repeating this for all options and plotting the results.

Pricing models assume volatility is the same for all strike prices and maturities. Volatility surfaces show this assumption is not true. Volatility has skew (volatility is different across strike prices) and a term structure (volatility is different across maturities). Quants use volatility surfaces to help calibrate models and price OTC derivatives that don’t trade on exchanges.

Imports and setup

Start by importing the libraries you need. I use yfinance to get options data for free.

1import numpy as np
2import pandas as pd
3import yfinance as yf
4import datetime as dt
5
6import matplotlib.pyplot as plt

yfinance returns data for all strikes for a single expiration at a time. It’s easier to work with all strikes and expirations at the same time so write a function to combine the expirations.

1def option_chains(ticker):
2    """
3    """
4    asset = yf.Ticker(ticker)
5    expirations = asset.options
6    
7    chains = pd.DataFrame()
8    
9    for expiration in expirations:
10        # tuple of two dataframes
11        opt = asset.option_chain(expiration)
12        
13        calls = opt.calls
14        calls['optionType'] = "call"
15        
16        puts = opt.puts
17        puts['optionType'] = "put"
18        
19        chain = pd.concat([calls, puts])
20        chain['expiration'] = pd.to_datetime(expiration) + pd.DateOffset(hours=23, minutes=59, seconds=59)
21        
22        chains = pd.concat([chains, chain])
23    
24    chains["daysToExpiration"] = (chains.expiration - dt.datetime.today()).dt.days + 1
25    
26    return chains

This function first gets all the expirations. Then it loops through each expiration and gets the option chain. It adds a column for option type, changes the expiration date to be at the end of the day, the combines each option chain together in a DataFrame. Finally, it computes the number of days until expiration.

Analyze skew and term structure

yfinance provides an estimate of implied volatility so you don’t have to compute it. This is ok for a quick analysis. In practice, quants derive their own implied volatility using custom models.

Start by downloading the data and getting the call options. The options data are in a pandas DataFrame which makes it easy.

1options = option_chains("SPY")
2
3calls = options[options["optionType"] == "call"]

Next, pick an expiration so you can plot the volatility skew.

1# print the expirations
2set(calls.expiration)
3
4# select an expiration to plot
5calls_at_expiry = calls[calls["expiration"] == "2023-01-20 23:59:59"]
6
7# filter out low vols
8filtered_calls_at_expiry = calls_at_expiry[calls_at_expiry.impliedVolatility >= 0.001]
9
10# set the strike as the index so pandas plots nicely
11filtered_calls_at_expiry[["strike", "impliedVolatility"]].set_index("strike").plot(
12    title="Implied Volatility Skew", figsize=(7, 4)
13)

The result is an image that looks like this.

In practice, quants use their own models to calculate implied volatility. They also filter out outliers and use smoothing algorithms. Notice the implied volatility varies with each strike. In particular, it is lowest at the $750 strike.

This is known as a volatility smile.

Next, build a volatility term structure. Pick a strike price to plot by expiration.

1# select an expiration to plot
2calls_at_strike = options[options["strike"] == 400.0]
3
4# filter out low vols
5filtered_calls_at_strike = calls_at_strike[calls_at_strike.impliedVolatility >= 0.001]
6
7# set the strike as the index so pandas plots nicely
8filtered_calls_at_strike[["expiration", "impliedVolatility"]].set_index("expiration").plot(
9    title="Implied Volatility Term Structure", figsize=(7, 4)
10)

The result is an image that looks like this.

Implied volatility is decreasing as the expiration dates get further out. This tells you the market expectation of volatility is lower in the future than it is today. You'll often see spikes in the term structure when big economic news is scheduled. The effect is caused by traders bidding up the prices of options in expectation of market swings.

Plot a volatility surface

By putting both charts together, you get the volatility surface. In derivatives pricing and trading, volatility surfaces are very important. Quants use the surface to price and trade other more exotic derivatives and look for ​market mispricings​. Volatility surfaces are also used to determine profit and loss by “marking trades to model.”

# pivot the dataframe
surface = (
    calls[['daysToExpiration', 'strike', 'impliedVolatility']]
    .pivot_table(values='impliedVolatility', index='strike', columns='daysToExpiration')
    .dropna()
)

# create the figure object
fig = plt.figure(figsize=(10, 8))

# add the subplot with projection argument
ax = fig.add_subplot(111, projection='3d')

# get the 1d values from the pivoted dataframe
x, y, z = surface.columns.values, surface.index.values, surface.values

# return coordinate matrices from coordinate vectors
X, Y = np.meshgrid(x, y)

# set labels
ax.set_xlabel('Days to expiration')
ax.set_ylabel('Strike price')
ax.set_zlabel('Implied volatility')
ax.set_title('Call implied volatility surface')

# plot
ax.plot_surface(X, Y, z)

The result is the famous volatility surface.