In the end I found that fitting a SABR smile to each tenor (borrowing a result from this answer) was sufficient to build a local vol surface that was smooth and well-behaved enough to build a variance surface worked nicely. I also fitted a Heston model to it, and the two surfaces do look fairly similar. Here is the final code and the fits generated (the long snippet at the very bottom is required to generate these plots, and also contains the raw data required)
Firstly, looping over each tenor and fitting a SABR smile:
# This is the 'SABR-solution'... fit a SABR smile to each tenor, and let the vol surface interpolate
# between them. Below, we're using the python minimizer to do a fit to the provided smiles
calibrated_params = {}
# params are sigma_0, beta, vol_vol, rho
params = [0.4, 0.6, 0.1, 0.2]
fig, i = plt.figure(figsize=(6, 42)), 1
for tte, group in full_df.groupby('tte'):
fwd = group.iloc[0]['fwd']
expiry = group.iloc[0]['expiry']
strikes = group.sort_values('strike')['strike'].values
vols = group.sort_values('strike')['vol'].values
def f(params):
params[0] = max(params[0], 1e-8) # Avoid alpha going negative
params[1] = max(params[1], 1e-8) # Avoid beta going negative
params[2] = max(params[2], 1e-8) # Avoid nu going negative
params[3] = max(params[3], -0.999) # Avoid nu going negative
params[3] = min(params[3], 0.999) # Avoid nu going negative
calc_vols = np.array([
ql.sabrVolatility(strike, fwd, tte, *params)
for strike in strikes
])
error = ((calc_vols - np.array(vols))**2 ).mean() **.5
return error
cons = (
{'type': 'ineq', 'fun': lambda x: x[0]},
{'type': 'ineq', 'fun': lambda x: 0.99 - x[1]},
{'type': 'ineq', 'fun': lambda x: x[1]},
{'type': 'ineq', 'fun': lambda x: x[2]},
{'type': 'ineq', 'fun': lambda x: 1. - x[3]**2}
)
result = optimize.minimize(f, params, constraints=cons, options={'eps': 1e-5})
new_params = result['x']
calibrated_params[tte] = {'v0': new_params[0], 'beta': new_params[1], 'alpha': new_params[2], 'rho': new_params[3], 'fwd': fwd}
newVols = [ql.sabrVolatility(strike, fwd, tte, *new_params) for strike in strikes]
# Start next round of optimisation with this round's parameters, they're probably quite close!
params = new_params
plt.subplot(len(tenors), 1, i)
i = i+1
plt.plot(strikes, vols, marker='o', linestyle='none', label='market {}'.format(expiry))
plt.plot(strikes, newVols, label='SABR {0:1.2f}'.format(tte))
plt.title("Smile {0:1.3f}".format(tte))
plt.grid()
plt.legend()
plt.show()
generates a sequence of plots like this, all of which mostly fit quite well:
which generates SABR params at each tenor looking like this (for this example I've set foreign and domestic discount curves to be flat):
Then I calibrated a local vol model and a Heston vol model, which actually both look quite close together:
# Fit a local vol surface to a strike-tenor grid extrapolated according to SABR
strikes = np.linspace(1.0, 1.5, 21)
expiration_dates = [calc_date + ql.Period(int(365 * x), ql.Days) for x in params.index]
implied_vols = []
for tte, row in params.iterrows():
fwd, v0, beta, alpha, rho = row['fwd'], row['v0'], row['beta'], row['alpha'], row['rho']
vols = [ql.sabrVolatility(strike, fwd, tte, v0, beta, alpha, rho) for strike in strikes]
implied_vols.append(vols)
implied_vols = ql.Matrix(np.matrix(implied_vols).transpose().tolist())
local_vol_surface = ql.BlackVarianceSurface(calc_date, calendar, expiration_dates, strikes, implied_vols, day_count)
# Fit a Heston model to the data as well
v0 = 0.005; kappa = 0.01; theta = 0.0064; rho = 0.0; sigma = 0.01
heston_process = ql.HestonProcess(dom_dcf_curve, for_dcf_curve, ql.QuoteHandle(ql.SimpleQuote(spot)), v0, kappa, theta, sigma, rho)
heston_model = ql.HestonModel(heston_process)
heston_engine = ql.AnalyticHestonEngine(heston_model)
# Set up Heston 'helpers' to calibrate to
heston_helpers = []
for idx, row in full_df.iterrows():
vol = row['vol']
strike = row['strike']
tenor = ql.Period(row['expiry'])
helper = ql.HestonModelHelper(tenor, calendar, spot, strike, ql.QuoteHandle(ql.SimpleQuote(vol)), dom_dcf_curve, for_dcf_curve)
helper.setPricingEngine(heston_engine)
heston_helpers.append(helper)
lm = ql.LevenbergMarquardt(1e-8, 1e-8, 1e-8)
heston_model.calibrate(heston_helpers, lm, ql.EndCriteria(5000, 100, 1.0e-8, 1.0e-8, 1.0e-8))
theta, kappa, sigma, rho, v0 = heston_model.params()
feller = 2 * kappa * theta - sigma ** 2
print(f"theta = {theta:.4f}, kappa = {kappa:.4f}, sigma = {sigma:.4f}, rho = {rho:.4f}, v0 = {v0:.4f}, spot = {spot:.4f}, feller = {feller:.4f}")
heston_handle = ql.HestonModelHandle(heston_model)
heston_vol_surface = ql.HestonBlackVolSurface(heston_handle)
# Plot the two vol surfaces ...
plot_vol_surface([local_vol_surface, heston_vol_surface], plot_years=np.arange(0.1, 1.0, 0.1), plot_strikes=np.linspace(1.05, 1.45, 20))
We expect the local vol model to price vanillas correctly but give unrelistic vol dynamics, while we expect Heston to give better vol dynamics but not price vanillas so well, but by calibrating a leverage function and using a Heston stochastic local vol model we can possibly get the best of both worlds - and this is also a good test that the local vol surface we've created is well behaved
# Calculate the Dupire instantaneous vol surface
local_vol_surface.setInterpolation('bicubic')
local_vol_handle = ql.BlackVolTermStructureHandle(local_vol_surface)
local_vol = ql.LocalVolSurface(local_vol_handle, dom_dcf_curve, for_dcf_curve, ql.QuoteHandle(ql.SimpleQuote(spot)))
# Calibrating a leverage function
end_date = ql.Date(21, 9, 2021)
generator_factory = ql.MTBrownianGeneratorFactory(43)
timeStepsPerYear = 182
nBins = 101
calibrationPaths = 2**19
stoch_local_mc_model = ql.HestonSLVMCModel(local_vol, heston_model, generator_factory, end_date, timeStepsPerYear, nBins, calibrationPaths)
leverage_functon = stoch_local_mc_model.leverageFunction()
plot_vol_surface(leverage_functon, funct='localVol', plot_years=np.arange(0.5, 0.98, 0.1), plot_strikes=np.linspace(1.05, 1.35, 20))
which produces a nice looking leverage function, which is close to 1 everywhere (indicating that the raw Heston fit was already quite good)
Boilerplate code to generate above images (including the FX delta-to-strike conversion):
import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
from scipy.stats import norm
from scipy import optimize, stats
import QuantLib as ql
calc_date = ql.Date(1, 9, 2020)
def plot_vol_surface(vol_surface, plot_years=np.arange(0.1, 3, 0.1), plot_strikes=np.arange(70, 130, 1), funct='blackVol'):
if type(vol_surface) != list:
surfaces = [vol_surface]
else:
surfaces = vol_surface
fig = plt.figure(figsize=(10, 6))
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(plot_strikes, plot_years)
Z_array, Z_min, Z_max = [], 100, 0
for surface in surfaces:
method_to_call = getattr(surface, funct)
Z = np.array([method_to_call(float(y), float(x))
for xr, yr in zip(X, Y)
for x, y in zip(xr, yr)]
).reshape(len(X), len(X[0]))
Z_array.append(Z)
Z_min, Z_max = min(Z_min, Z.min()), max(Z_max, Z.max())
# In case of multiple surfaces, need to find universal max and min first for colourmap
for Z in Z_array:
N = (Z - Z_min) / (Z_max - Z_min) # normalize 0 -> 1 for the colormap
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0.1, facecolors=cm.coolwarm(N))
m = cm.ScalarMappable(cmap=cm.coolwarm)
m.set_array(Z)
plt.colorbar(m, shrink=0.8, aspect=20)
ax.view_init(30, 300)
def generate_multi_paths_df(process, num_paths=1000, timestep=24, length=2):
"""Generates multiple paths from an n-factor process, each factor is returned in a seperate df"""
times = ql.TimeGrid(length, timestep)
dimension = process.factors()
rng = ql.GaussianRandomSequenceGenerator(ql.UniformRandomSequenceGenerator(dimension * timestep, ql.UniformRandomGenerator()))
seq = ql.GaussianMultiPathGenerator(process, list(times), rng, False)
paths = [[] for i in range(dimension)]
for i in range(num_paths):
sample_path = seq.next()
values = sample_path.value()
spot = values[0]
for j in range(dimension):
paths[j].append([x for x in values[j]])
df_paths = [pd.DataFrame(path, columns=[spot.time(x) for x in range(len(spot))]) for path in paths]
return df_paths
# Define functions to map from delta to strike
def strike_from_spot_delta(tte, fwd, vol, delta, dcf_for, put_call):
sigma_root_t = vol * np.sqrt(tte)
inv_norm = norm.ppf(delta * put_call * dcf_for)
return fwd * np.exp(-sigma_root_t * put_call * inv_norm + 0.5 * sigma_root_t * sigma_root_t)
def strike_from_fwd_delta(tte, fwd, vol, delta, put_call):
sigma_root_t = vol * np.sqrt(tte)
inv_norm = norm.ppf(delta * put_call)
return fwd * np.exp(-sigma_root_t * put_call * inv_norm + 0.5 * sigma_root_t * sigma_root_t)
# World State for Vanilla Pricing
spot = 1.17858
rateDom = 0.0
rateFor = 0.0
calendar = ql.NullCalendar()
day_count = ql.Actual365Fixed()
# Set up the flat risk-free curves
riskFreeCurveDom = ql.FlatForward(calc_date, rateDom, ql.Actual365Fixed())
riskFreeCurveFor = ql.FlatForward(calc_date, rateFor, ql.Actual365Fixed())
dom_dcf_curve = ql.YieldTermStructureHandle(riskFreeCurveDom)
for_dcf_curve = ql.YieldTermStructureHandle(riskFreeCurveFor)
tenors = ['1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y']
deltas = ['ATM', '35D Call EUR', '35D Put EUR', '25D Call EUR', '25D Put EUR', '15D Call EUR', '15D Put EUR', '10D Call EUR', '10D Put EUR', '5D Call EUR', '5D Put EUR']
vols = [[7.255, 7.428, 7.193, 7.61, 7.205, 7.864, 7.261, 8.033, 7.318, 8.299, 7.426],
[7.14, 7.335, 7.07, 7.54, 7.08, 7.836, 7.149, 8.032, 7.217, 8.34, 7.344],
[7.195, 7.4, 7.13, 7.637, 7.167, 7.984, 7.286, 8.226, 7.394, 8.597, 7.58],
[7.17, 7.39, 7.11, 7.645, 7.155, 8.031, 7.304, 8.303, 7.438, 8.715, 7.661],
[7.6, 7.827, 7.547, 8.105, 7.615, 8.539, 7.796, 8.847, 7.952, 9.308, 8.222],
[7.285, 7.54, 7.26, 7.878, 7.383, 8.434, 7.671, 8.845, 7.925, 9.439, 8.344],
[7.27, 7.537, 7.262, 7.915, 7.425, 8.576, 7.819, 9.078, 8.162, 9.77, 8.713],
[7.275, 7.54, 7.275, 7.935, 7.455, 8.644, 7.891, 9.188, 8.283, 9.922, 8.898],
[7.487, 7.724, 7.521, 8.089, 7.731, 8.742, 8.197, 9.242, 8.592, 9.943, 9.232],
[7.59, 7.81, 7.645, 8.166, 7.874, 8.837, 8.382, 9.354, 8.816, 10.065, 9.51]]
# Convert vol surface to strike surface (we need both)
full_option_surface = []
for i, name in enumerate(deltas):
delta = 0.5 if name == "ATM" else int(name.split(" ")[0].replace("D", "")) / 100.
put_call = 1 if name == "ATM" else -1 if name.split(" ")[1] == "Put" else 1
for j, tenor in enumerate(tenors):
expiry = calc_date + ql.Period(tenor)
tte = day_count.yearFraction(calc_date, expiry)
fwd = spot * for_dcf_curve.discount(expiry) / dom_dcf_curve.discount(expiry)
for_dcf = for_dcf_curve.discount(expiry)
vol = vols[j][i] / 100.
# Assume that spot delta used out to 1Y (used to be this way...)
if tte < 1.:
strike = strike_from_spot_delta(tte, fwd, vol, put_call*delta, for_dcf, put_call)
else:
strike = strike_from_fwd_delta(tte, fwd, vol, put_call*delta, put_call)
full_option_surface.append({"vol": vol, "fwd": fwd, "expiry": tenor, "tte": tte, "delta": put_call*delta, "strike": strike, "put_call": put_call, "for_dcf": for_dcf, "name": name})
full_df = pd.DataFrame(full_option_surface)
display_df = full_df.copy()
display_df['call_delta'] = 1 - (display_df['put_call'].clip(0) - display_df['delta'])
df = display_df.set_index(['tte', 'call_delta']).sort_index()[['strike']].unstack()
df = df.reindex(sorted(df.columns, reverse=True), axis=1)
fig = plt.figure(figsize=(12,9))
plt.subplot(2,1,1)
plt.plot(full_df['tte'], full_df['strike'], marker='o', linestyle='none', label='strike grid')
plt.title("Option Strike Grid, tte vs. K")
plt.grid()
plt.xlim(0, 2.1)
df
