I've been working on this problem a little bit lately. Unfortunately in the FX context, it's not quite as straight-forward as in the equities case, for two reasons:
- FX options trade OTC instead of on exchange, so you need access to broker screens to trade them (eg. on BBG)
- FX Options are quoted by (delta, tenor, vol) instead of (strike, tenor, price) so we have to do a bit of pre-work to get the options corresponding strikes for our Heston calibration
A EURUSD options screen from BBG looks something like this:
Trades are done OTC between clients, but many still need to be reported to the DTCC, and BBG has a screen showing some examples of recent OTC options that were traded:
The exact procedure required to turn these into (strike, price) pairs depends on the currency pair under consideration, a great reference on the conventions is found in this paper, but it turns out to be relatively simple for EURUSD. As described in the paper, you need a function that looks like this:
import numpy as np
from scipy.stats import norm
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)
strike = strike_from_fwd_delta(tte, fwd, vol, put_call*delta, put_call)
After doing that, I've got two tables (NB this is a different dataset to that shown in the screen image above, because I transcribed and calculated it earlier) - the original table showing the vol for each (delta, tenor) pair, and the new one showing the strike for each pair. The new table looks something like this:
Now we have enough to calibrate a Heston vol surface using the (tenor, strike, vol) triples from each observed option (nb. you'll also have to fit domestic and foreign rates curves, but that's another story) - for my options above, the surface looks like this:
Here is a sample of code (the data above is hard-coded at the top) that will generate the vol surface above for you:
import numpy as np
from matplotlib import pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.mplot3d import Axes3D
import QuantLib as ql
strikes = [1.1787, 1.1788, 1.1794, 1.1804, 1.1815, 1.1846, 1.1873, 1.1909, 1.1978, 1.2046, 1.1833, 1.1854, 1.1891, 1.1942, 1.1995, 1.2092, 1.2178, 1.2263, 1.2426, 1.2574, 1.1741, 1.1725, 1.1702, 1.1673, 1.1646, 1.1619, 1.1598, 1.158, 1.1561, 1.1556, 1.1871, 1.1906, 1.197, 1.2056, 1.2143, 1.2301, 1.2441, 1.2571, 1.2814, 1.3034, 1.1708, 1.1678, 1.1632, 1.1574, 1.1517, 1.1442, 1.1379, 1.1327, 1.1241, 1.1179, 1.192, 1.1977, 1.2078, 1.2214, 1.2351, 1.2605, 1.2834, 1.304, 1.3402, 1.374, 1.1664, 1.1618, 1.1542, 1.1445, 1.1349, 1.1206, 1.1081, 1.0979, 1.0805, 1.0667, 1.1956, 1.2028, 1.2157, 1.233, 1.2506, 1.2839, 1.3147, 1.3419, 1.3876, 1.4314, 1.1635, 1.1577, 1.1479, 1.1354, 1.1231, 1.1035, 1.0859, 1.0718, 1.0483, 1.0288, 1.2012, 1.211, 1.2284, 1.2519, 1.2758, 1.3228, 1.3668, 1.4053, 1.4677, 1.5291, 1.1589, 1.1513, 1.1381, 1.1212, 1.1046, 1.0763, 1.0505, 1.0301, 0.997, 0.9687]
vols = [0.0726, 0.0714, 0.072, 0.0717, 0.076, 0.0728, 0.0727, 0.0728, 0.0749, 0.0759, 0.0743, 0.0733, 0.074, 0.0739, 0.0783, 0.0754, 0.0754, 0.0754, 0.0772, 0.0781, 0.0719, 0.0707, 0.0713, 0.0711, 0.0755, 0.0726, 0.0726, 0.0728, 0.0752, 0.0764, 0.0761, 0.0754, 0.0764, 0.0764, 0.0811, 0.0788, 0.0791, 0.0793, 0.0809, 0.0817, 0.0721, 0.0708, 0.0717, 0.0716, 0.0761, 0.0738, 0.0742, 0.0746, 0.0773, 0.0787, 0.0786, 0.0784, 0.0798, 0.0803, 0.0854, 0.0843, 0.0858, 0.0864, 0.0874, 0.0884, 0.0726, 0.0715, 0.0729, 0.073, 0.078, 0.0767, 0.0782, 0.0789, 0.082, 0.0838, 0.0803, 0.0803, 0.0823, 0.083, 0.0885, 0.0885, 0.0908, 0.0919, 0.0924, 0.0935, 0.0732, 0.0722, 0.0739, 0.0744, 0.0795, 0.0793, 0.0816, 0.0828, 0.0859, 0.0882, 0.083, 0.0834, 0.086, 0.0872, 0.0931, 0.0944, 0.0977, 0.0992, 0.0994, 0.1006, 0.0743, 0.0734, 0.0758, 0.0766, 0.0822, 0.0834, 0.0871, 0.089, 0.0923, 0.0951]
expiries = ['1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y', '1W', '2W', '1M', '2M', '3M', '6M', '9M', '1Y', '18M', '2Y']
rate = 0.0
today = ql.Date(1, 9, 2020)
spot = 1.1786
usd_calendar = ql.NullCalendar()
# Set up the flat risk-free curves
usd_curve = ql.FlatForward(today, 0.0, ql.Actual365Fixed())
eur_curve = ql.FlatForward(today, 0.0, ql.Actual365Fixed())
usd_rates_ts = ql.YieldTermStructureHandle(usd_curve)
eur_rates_ts = ql.YieldTermStructureHandle(eur_curve)
v0 = 0.005; kappa = 0.01; theta = 0.0064; rho = 0.0; sigma = 0.01
heston_process = ql.HestonProcess(usd_rates_ts, eur_rates_ts, 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 strike, vol, expiry in zip(strikes, vols, expiries):
tenor = ql.Period(expiry)
helper = ql.HestonModelHelper(tenor, usd_calendar, spot, strike, ql.QuoteHandle(ql.SimpleQuote(vol)), usd_rates_ts, eur_rates_ts)
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}")
# Plot the vol surface ...
heston_handle = ql.HestonModelHandle(heston_model)
heston_vol_surface = ql.HestonBlackVolSurface(heston_handle)
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)
plot_vol_surface(heston_vol_surface, plot_years=np.arange(0.1, 2.0, 0.1), plot_strikes=np.linspace(1.0, 1.5, 30))
