Optiver lecture: Trading options as a market maker¶

# QuantBook Analysis Tool 

# For more information see [https://www.quantconnect.com/docs/research/overview]

qb = QuantBook()

​

# In this notebook I will look at the SPY

# Initialize the securities objects 

spy = qb.AddEquity("SPY")

spy_options = qb.AddOption("SPY")

​

# We need to specify the stikes that are to our interest

# Looking at 10 stikes under up to 15 stikes over the ATM point

# With timedelt we filter for options expiring in els than 30 days

spy_options.SetFilter(-10,+15, timedelta(0), timedelta(30))

# Define the period we're interested in.

# Looking at 30 minutes here

start_time =  datetime(2021, 11, 30, 9, 30) 

end_time = datetime(2021, 11, 30, 10, 0)

​

# Looking at 20 minutes of minute-resolution data

underlying_history = qb.History(spy.Symbol, start_time, end_time, Resolution.Minute)

option_history = qb.GetOptionHistory(spy.Symbol, start_time, end_time, Resolution.Minute)

​

# Get the 1 year risk-free rate from the yield curve.

# The only available resolution is daily. We'll deal with this later

yield_curve = qb.AddData(USTreasuryYieldCurveRate, "USTYCR")

oneyear_rate=qb.History(yield_curve.Symbol, start_time-timedelta(1), end_time, Resolution.Daily).oneyear.droplevel(0)

data = option_history.GetAllData()

data.head()

print(f"Option expiries in this dataset: \n {option_history.GetExpiryDates()}\n")

print(f"Selected strike prices: \n {option_history.GetStrikes()}\n")

print(f"Row indexers: \n {data.index.names}\n")

print(f"Columns: \n {data.columns}\n")

# Only keep the bid and ask columns and drop the expiry indexes since we work with only one expiry and one underlying

reordered_data = data.unstack(level=4)[['bidclose', 'askclose']].droplevel([0,3])

​

# Only keep call options and drop the index afterwards

call_data = reordered_data[reordered_data.index.isin(["Call"], level=1)].droplevel(1).T

​

# split up the bid and ask prices in separate DataFrames

call_bid = call_data[call_data.index.isin(['bidclose'], level=0)].droplevel(0)

call_ask = call_data[call_data.index.isin(['askclose'], level=0)].droplevel(0)

​

call_bid.head()

# We can also we can simplify the dataframe of the underlying time series

underlying_data = underlying_history[['bidclose', 'askclose']].droplevel(0)

underlying_data.head()

import numpy as np

​

# Begin by listing all the variables we have to our disposal

​

C = call_bid.iloc[:, 0]

S = underlying_data['bidclose']

X = pd.Series(np.ones(len(underlying_data)) * call_bid.columns[0], index=underlying_data.index)

t = pd.Series(np.ones(len(underlying_data)) * ((data.index[0][0] - call_bid.index[0]).days / 252),index=underlying_data.index)  # Assuming 252 days in a year

r = pd.Series(np.ones(len(underlying_data)) * oneyear_rate.values[0], index=underlying_data.index)

​

BS = pd.concat([C, S, X, t, r], axis=1)

BS.columns = ['C', 'S', 'X', 't', 'r']

BS.head()

​

# Applying the formula above to calculate Implied Volatility (IV)

IV = ((((2*np.pi)**0.5)/(BS.S+BS.X)) * (BS.C-((BS.S-BS.X)/2) 

    + ((((BS.C-((BS.S-BS.X)/2)))**2)-((BS.S-BS.X)/np.pi))**0.5)) / (t**0.5)

BS['IV'] = IV

BS.head()

IVs = list()

​

for strike_index in range(len(call_bid.columns)):

    C = call_bid.iloc[:, strike_index]

    S = underlying_data['bidclose']

    X = pd.Series(np.ones(len(underlying_data)) * call_bid.columns[strike_index], index=underlying_data.index)

    t = pd.Series(np.ones(len(underlying_data)) * ((data.index[0][0] - call_bid.index[0]).days / 252),index=underlying_data.index)  # Assuming 252 days in a year

    r = pd.Series(np.ones(len(underlying_data)) * oneyear_rate.values[0], index=underlying_data.index)

    BS = pd.concat([C, S, X, t, r], axis=1)

    BS.columns = ['C', 'S', 'X', 't', 'r']

​

    IV = ((((2*np.pi)**0.5)/(BS.S+BS.X)) * (BS.C-((BS.S-BS.X)/2) 

    + ((((BS.C-((BS.S-BS.X)/2)))**2)-((BS.S-BS.X)/np.pi))**0.5)) / (t**0.5)

​

    IVs.append(IV)

​

IV_series = pd.concat(IVs, axis=1)

IV_series.columns = call_bid.columns

IV_series.head()

import matplotlib.pyplot as plt

# from mpl_toolkits import mplot3d

fig = plt.figure(figsize=(15,10))

​

plt.plot(IV_series.T*100)

​

plt.suptitle("Evolution of implied volatility skew", size = 20, y=1)

plt.title("SPY options on 2021-11-30 from 9:30 to 10:00")

​

plt.ylabel("Implied Volatility (%)")

plt.xlabel("Strike price")

​

lbl = [i.strftime("%H-%M") for i in IV_series.index]

plt.legend(lbl, loc=(1.01,0.1), title="Time")

​

plt.show()

# In-the-money option seem to have lower volatility of the IV itself

fig = plt.figure(figsize=(15,10))

​

plt.plot(IV_series*100)

​

plt.suptitle("Evolution of implied volatility over time, per strike", size = 20, y=1)

plt.title("SPY options on 2021-11-30 from 9:30 to 10:00")

​

plt.ylabel("Implied Volatility (%)")

plt.xlabel("Time (day-hour-minute)")

​

lbl = ["$"+str(int(i)) for i in IV_series.columns]

plt.legend(lbl, loc=(1.01,0.2), title="Strike price")

plt.show()

# Corners on the curved (bottom) parts of the lines seem to be pricing inefficiencies that get corrected

​

fig = plt.figure(figsize=(15,10))

ax = fig.add_subplot(111, projection='3d')

​

nrows = IV_series.shape[0]

ncols = IV_series.shape[1]

​

x = np.array([[i.minute] * ncols for i in (IV_series.index - timedelta(minutes=1))]).ravel() # x coordinates 

y = np.array([i for i in IV_series.columns] * nrows) # y coordinates 

z = np.array(IV_series*100).ravel() # z coordinates

dx = np.ones(nrows*ncols) # length along x-axis 

dy = np.ones(nrows*ncols) # length along y-axis 

dz = np.array(IV_series*100).ravel() # length along z-axis 

​

ax.set_xlabel("Time (Minutes)")

ax.set_ylabel("Strike price")

ax.set_zlabel("Implied volatility (%)")

​

plt.suptitle("Evolution of implied volatility skew", size = 20, y=1)

plt.title("SPY options on 2021-11-30 from 9:30 to 10:00")

​

ax.text2D(-0.05,-0.025,s="Inefficiencies ->", color='green', size=13)

ax.text2D(-0.007,-0.04,s="Inefficiencies ->", color='green', size=13)

​

ax.plot(x, y, z, dx, dy, dz)

plt.show()

# Let's go back to a singe option to simulate hedging it

# DataFrame from step 3 with strike of 453.0

BS['IV'] = ((((2*np.pi)**0.5)/(BS.S+BS.X)) * (BS.C-((BS.S-BS.X)/2) 

    + ((((BS.C-((BS.S-BS.X)/2)))**2)-((BS.S-BS.X)/np.pi))**0.5)) / (t**0.5)

BS.head()

# Source: https://aaronschlegel.me/measure-sensitivity-derivatives-greeks-python.html

import scipy.stats as si

def delta_call(S, X, t, r, IV):

    d1 = (np.log(S / X) + (r + 0.5 * IV ** 2) * t) / (IV * np.sqrt(t))

    delta_call = si.norm.cdf(d1, 0.0, 1.0)

    return delta_call

# With a delta of 0.74 we need to buy 74 option for over option contract bought

# (given that the lot size of option contracts is 100)

BS['delta'] = BS.apply(lambda x: delta_call(x.S, x.X, x.t, x.r, x.IV), axis=1)

BS.head()

fig, (ax1, ax2) = plt.subplots(2,1, sharex=True)

​

ax1.plot(BS.delta)

ax1.set_title("Evolution of an option's delta", size=10)

ax1.set_ylabel("Delta")

ax1.set_xlabel("Time (Day, Hour, Mintue)")

​

ax2.plot(BS.S)

ax2.set_title("Price of the underlying (SPY)", size=10)

ax2.set_ylabel("Price")

ax2.set_xlabel("Time (Day, Hour, Mintue)")

​

fig.subplots_adjust(hspace=0.4)

plt.show()

# If the option is bought at 09:31, then we would want to sell 74.38 shares

pf = BS.delta *(-1)

​

# At 09:32 the stock price seems to have gone up, as did the delta

# This means we would want to sell an additional 0.02 shares

# Etc. for every minute that goes by

trades = pf.diff()

trades[0] =  -0.01  # BS.delta[0] *(-1) For visual purposes

trades = pd.concat([trades, BS.S], axis=1)

trades.columns = ['trades', 'price'] 

​

print(pf[0])

trades.head()

# Liquidate the portfolio in the last trade

trades.iloc[0,0] = BS.delta[0] *(-1)

trades.iloc[-1,0] = sum(trades.trades) *(-1)

trades.head()

# Simulate PnL of hedging 

# Watch out: Here we assumed no transaction costs, no spread costs and perfect liquidity for our trades

PnL = sum(trades.trades*(trades.price/trades.price[0]))

print(f"{round(PnL*100, 2)}%")