Modello Heston su valuta

Question

Potremmo avere la formula per il Hestonmodello per la valuta come (sotto il Risk-neutral measurefor$r_d$) -

$dS_t = \left( r_d - r_f \right) S_tdt+S_t \sqrt{V_t}dW^S$

$dV_t = a(\bar{V}- V_t)dt + \eta \sqrt{V_t}dW_t^V$

Tipicamente stimiamo i parametri del modello osservando i prezzi delle opzioni Call e Put con diverse scadenze.

Tuttavia, per il caso Valuta, dove posso vedere tali contratti di opzione negoziabili sul mercato? Come in ECM (https://www.cmegroup.com/trading/fx/g10/euro-fx_quotes_globex_options.html?optionProductId=59#optionProductId=8117&strikeRange=ATM), la maggior parte delle opzioni sono negoziate su Futures.

Quindi, se voglio stimare i parametri del modello per EUR-USD spot processcome quello nel terminale Bloomberghttps://www.bloomberg.com/quote/EURUSD:CUR, come devo procedere per stimare i parametri del modello?

Qualsiasi puntatore sarà molto apprezzato.

Answer 1

Ho lavorato un po' su questo problema ultimamente. Sfortunatamente nel contesto FX, non è così semplice come nel caso delle azioni, per due motivi:

1. Le opzioni FX negoziano OTC invece che in borsa, quindi è necessario accedere alle schermate del broker per negoziarle (ad es. su BBG)
2. Le opzioni FX sono quotate da (delta, tenor, vol) invece di (strike, tenor, price), quindi dobbiamo fare un po' di lavoro preliminare per ottenere le opzioni che corrispondono ai strike per la nostra calibrazione Heston

Una schermata delle opzioni EURUSD di BBG è simile a questa:

Gli scambi vengono effettuati OTC tra clienti, ma molti devono ancora essere segnalati al DTCC e BBG ha una schermata che mostra alcuni esempi di recenti opzioni OTC che sono state negoziate:

L'esatta procedura richiesta per trasformarle in coppie (strike, prezzo) dipende dalla coppia di valute in esame, un ottimo riferimento sulle convenzioni si trova in questo documento , ma risulta essere relativamente semplice per EURUSD. Come descritto nel documento, è necessaria una funzione simile a questa:

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)

Dopo averlo fatto, ho due tabelle ( NB questo è un set di dati diverso da quello mostrato nell'immagine dello schermo sopra, perché l'ho trascritto e calcolato in precedenza) - la tabella originale che mostra il vol per ciascuna coppia (delta, tenore), e quello nuovo che mostra lo sciopero per ogni coppia. La nuova tabella è simile a questa:

Ora abbiamo abbastanza per calibrare una superficie Heston vol usando le triple (tenore, strike, vol) di ciascuna opzione osservata (nb. dovrai anche adattare le curve dei tassi nazionali ed esteri, ma questa è un'altra storia) - per le mie opzioni sopra , la superficie si presenta così:

Ecco un esempio di codice (i dati sopra sono codificati in alto) che genererà la superficie vol sopra per te:

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))