Zeliade's whitepaper describes an efficient way to calibrate Gatheral's SVI model. This code is largely based on bramjochems's F# implementation in MyExcelLib.
In [1]:
%pylab inline
import numpy as np
import pandas as pd
today = pd.to_datetime('2015-05-08')
futures = pd.read_csv('FTSE100-Futures-20150508.csv', dtype='U,U,f8,f8,f8,f8,f8,f8,i8,i8,i8,i8,i8', sep='\t')
options = pd.read_csv('FTSE100-Options-20150508.csv', dtype='U,U,i8,U,f8,f8,f8,f8,f8,f8,f8,i8,i8,i8,i8', sep='\t')
In [2]:
options.MONTH = pd.DatetimeIndex(pd.to_datetime(options.MONTH, format='%b%y')) + pd.datetools.WeekOfMonth(week=3,weekday=4)
futures.CONTRACTMONTH = pd.DatetimeIndex(pd.to_datetime(futures.CONTRACTMONTH, format='%y-%b')) + pd.datetools.WeekOfMonth(week=3,weekday=4)
Build the forward curve as a step function between the points:
In [3]:
forward_curve = pd.Series(futures.PRICE.values, futures.CONTRACTMONTH.values, name='FUTUREPRICE')
forward_curve = forward_curve.resample('D', fill_method='bfill')
forward_curve = forward_curve.reindex(index=pd.date_range(today, forward_curve.index.max()), method='bfill', copy=False);
Use a flat yield curve for simplicity
In [4]:
yield_curve = pd.Series([0.01], index=forward_curve.index, name='RISKFREERATE')
Join the forward curve and yield curve onto the options, and calculate the moneyness and time to expiry in years (t).
In [5]:
options = pd.merge(options, forward_curve.to_frame(), left_on='MONTH', right_index=True, how='left')
options = pd.merge(options, yield_curve.to_frame(), left_on='MONTH', right_index=True, how='left')
options['MONEYNESS'] = np.log(options.STRIKE / options.FUTUREPRICE)
options['t'] = options.apply(lambda row: (row.MONTH - today).days / 365.25, axis=1)
Only use out-of-the-money options to calibrate the surface
In [6]:
options = options[
(((options.PC == 'C') & (options.STRIKE > options.FUTUREPRICE)) |
((options.PC == 'P') & (options.STRIKE < options.FUTUREPRICE)))]
We'll convert the option prices into implied volatities using the Black's model implementation in vollib, dropping any prices that fail:
In [7]:
from vollib.black import implied_volatility
def calc_iv(row):
return implied_volatility.implied_volatility_of_discounted_option_price(
row.PRICE,
row.FUTUREPRICE,
row.STRIKE,
row.RISKFREERATE,
row.t,
row.PC.lower())
options['IMPLIEDVOL'] = options.apply(calc_iv, axis=1)
options = options[options.IMPLIEDVOL > 0]
Do the parameters fall inside Zeniade's compact and convex domain D?
In [8]:
import sys
def acceptable(S, a, d, c, vT):
eps = sys.float_info.epsilon
return -eps <= c and c <= 4 * S + eps and abs(d) <= min(c, 4 * S - c) + eps and -eps <= a and a <= vT.max() + eps
The error function from the optimizer, in terms of a, d and c so we don't need to keep dividing by T:
In [9]:
def sum_of_squares(x, a, d, c, vT):
diff = a + d * x + c * sqrt(x * x + 1) - vT
return (diff * diff).sum()
The key part of the Zeliade paper - identify the values of a, d and c that produce the lowest sum-of-squares error when the gradient of the f function is zero, ensuring that the parameters remain in the domain D identified in the paper.
In [10]:
def solve_grad((S, M), x, vT):
ys = (x - M) / S
y = ys.sum()
y2 = (ys * ys).sum()
y2one = (ys * ys + 1).sum()
ysqrt = sqrt(ys * ys + 1).sum()
y2sqrt = (ys * sqrt(ys * ys + 1)).sum()
v = vT.sum()
vy = (vT * ys).sum()
vsqrt = (vT * sqrt(ys * ys + 1)).sum()
matrix = [
[1, y, ysqrt],
[y, y2, y2sqrt],
[ysqrt, y2sqrt, y2one]
]
vector = [v, vy, vsqrt]
_a, _d, _c = np.linalg.solve(np.array(matrix), np.array(vector))
if acceptable(S, _a, _d, _c, vT):
return _a, _d, _c, sum_of_squares(ys, _a, _d, _c, vT)
_a, _d, _c, _cost = None, None, None, None
for matrix, vector, clamp_params in [
([[1,0,0],[y,y2,y2sqrt],[ysqrt,y2sqrt,y2one]], [0, vy, vsqrt], False), # a = 0
([[1,0,0],[y,y2,y2sqrt],[ysqrt,y2sqrt,y2one]], [vT.max(),vy,vsqrt], False), # a = _vT.max()
([[1,y,ysqrt],[0,-1,1],[ysqrt,y2sqrt,y2one]], [v,0,vsqrt], False), # d = c
([[1,y,ysqrt],[0,1,1],[ysqrt,y2sqrt,y2one]], [v,0,vsqrt], False), # d = -c
([[1,y,ysqrt],[0,1,1],[ysqrt,y2sqrt,y2one]], [v,4*S,vsqrt], False), # d <= 4*s-c
([[1,y,ysqrt],[0,-1,1],[ysqrt,y2sqrt,y2one]], [v,4*S,vsqrt], False), # -d <= 4*s-c
([[1,y,ysqrt],[y,y2,y2sqrt],[0,0,1]], [v,vy,0], False), # c = 0
([[1,y,ysqrt],[y,y2,y2sqrt],[0,0,1]], [v,vy,4*S], False), # c = 4*S
([[1,y,ysqrt],[0,1,0],[0,0,1]], [v,0,0], True), # c = 0, implies d = 0, find optimal a
([[1,y,ysqrt],[0,1,0],[0,0,1]], [v,0,4*S ], True), # c = 4s, implied d = 0, find optimal a
([[1,0,0],[0,-1,1],[ysqrt,y2sqrt,y2one]], [0,0,vsqrt], True), # a = 0, d = c, find optimal c
([[1,0,0],[0,1,1],[ysqrt,y2sqrt,y2one]], [0,0,vsqrt], True), # a = 0, d = -c, find optimal c
([[1,0,0],[0,1,1],[ysqrt,y2sqrt,y2one]], [0,4*S,vsqrt], True), # a = 0, d = 4s-c, find optimal c
([[1,0,0],[0,-1,1],[ysqrt,y2sqrt,y2one]], [0,4*S,vsqrt], True) # a = 0, d = c-4s, find optimal c
]:
a, d, c = np.linalg.solve(np.array(matrix), np.array(vector))
if clamp_params:
dmax = min(c, 4 * S - c)
a = min(max(a, 0), vT.max())
d = min(max(d, -dmax), dmax)
c = min(max(c, 0), 4 * S)
cost = sum_of_squares(ys, a, d, c, vT)
if acceptable(S, a, d, c, vT) and (_cost is None or cost < _cost):
_a, _d, _c, _cost = a, d, c, cost
assert _cost is not None, "S=%s, M=%s" % (S, M)
return _a, _d, _c, _cost
Use scipy's optimizer to minimize the error function
In [11]:
from scipy import optimize
def solve_grad_get_score((S, M), x, vT):
return solve_grad((S, M), x, vT)[3]
def calibrate(df):
vT = df.IMPLIEDVOL * df.IMPLIEDVOL * df.t
res = optimize.minimize(solve_grad_get_score, [.1, .0], args=[df.MONEYNESS, vT], bounds=[(0.001, None), (None, None)])
assert res.success
S, M = res.x
a, d, c, _ = solve_grad((S, M), df.MONEYNESS, vT)
T = df.t.max() # should be the same for all rows
A, P, B = a / T, d / c, c / (S * T)
assert T >= 0 and S >= 0 and abs(P) <= 1
return A, P, B, S, M
SVI formula for the volatility (not variance):
In [12]:
def svi(A, P, B, S, M, x):
return sqrt(A + B * (P * (x - M) + sqrt((x - M) * (x - M) + S * S)))
Example calibration of a single slice:
In [163]:
curve = options[options.MONTH == pd.to_datetime('2015-12-25')].sort('STRIKE')
A, P, B, S, M = calibrate(curve)
def calculated_vol(row):
return svi(A, P, B, S, M, row.MONEYNESS)
curve['CALCULATEDVOL'] = curve.apply(calculated_vol, axis=1)
curve.plot(x='MONEYNESS', y=['IMPLIEDVOL', 'CALCULATEDVOL'])
Out[163]:
To build a surface, linearly interpolate between slices, or use the first or last slice only if the given date is outside the curve:
In [180]:
from functools import partial
calibrated = [(month, partial(svi, *calibrate(df))) for month, df in options.groupby('t')]
slices = pd.Series([x[1] for x in calibrated], [x[0] for x in calibrated])
def volsurface(slices, t, moneyness):
ts = []
if t >= slices.index.min():
ts.append(slices.index.get_loc(t, method='ffill'))
if t <= slices.index.max():
ts.append(slices.index.get_loc(t, method='bfill'))
if len(ts) == 1:
return slices.iat[ts[0]](moneyness)
period = slices.index[ts[1]] - slices.index[ts[0]]
fraction = (t - slices.index[ts[0]]) / period
return (1 - fraction) * slices.iat[ts[0]](moneyness) + fraction * slices.iat[ts[1]](moneyness)
Plot a slice across expiries for a single moneyness:
In [178]:
r = frange(0, 1, 0.01)
pd.Series([volsurface(slices, t, -.5) for t in r], r).plot()
Out[178]:
Finally draw the surface:
In [179]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
moneynesses = linspace(-0.7, 0.7, 100)
expiries = linspace(0, 1, 50)
X, Y = meshgrid(moneynesses, expiries)
Z = np.array([volsurface(slices, t, moneyness) for t, moneyness in zip(np.ravel(X), np.ravel(Y))]).reshape(X.shape)
fig = plt.figure(figsize=(20,15))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.set_xlabel('Moneyness')
ax.set_ylabel('Expiry')
ax.set_zlabel('Vol')
p = ax.plot_wireframe(X, Y, Z, rstride=1, cstride=1)