Modelo Heston sobre la moneda

Question

Podríamos tener la fórmula de Heston modelo de moneda como (bajo la Risk-neutral measure para $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$

Normalmente, estimamos los parámetros del modelo observando los precios de las opciones Call y Put con diferentes vencimientos.

Sin embargo, para el caso de la moneda, ¿dónde puedo ver esos contratos de opciones negociables en el mercado? Como en CME ( https://www.cmegroup.com/trading/fx/g10/euro-fx_quotes_globex_options.html?optionProductId=59#optionProductId=8117&strikeRange=ATM ), la mayoría de las opciones se negocian en el Futures .

Así que si quiero estimar los parámetros del modelo para EUR-USD spot process como la del terminal de Bloomberg https://www.bloomberg.com/quote/EURUSD:CUR ¿Cómo debo proceder para estimar los parámetros del modelo?

Cualquier indicación será muy apreciada.

Answer 1

Últimamente he estado trabajando un poco en este problema. Por desgracia, en el contexto de las divisas, no es tan sencillo como en el caso de la renta variable, por dos razones:

1. Las opciones FX se negocian OTC en lugar de en bolsa, por lo que se necesita acceso a las pantallas de los corredores para negociarlas (por ejemplo, en BBG)
2. Las opciones FX se cotizan por (delta, tenor, vol) en lugar de (strike, tenor, precio) por lo que tenemos que hacer un poco de trabajo previo para obtener los strikes correspondientes a las opciones para nuestra calibración Heston

Una pantalla de opciones de EURUSD de BBG se parece a esto:

Las operaciones se realizan OTC entre clientes, pero muchas de ellas aún deben notificarse a la DTCC, y BBG tiene una pantalla que muestra algunos ejemplos de opciones OTC recientes que se negociaron:

El procedimiento exacto para convertirlos en pares (strike, precio) depende del par de divisas considerado, una gran referencia sobre las convenciones se encuentra en este documento pero resulta ser relativamente sencillo para el EURUSD. Como se describe en el documento, se necesita una función que se parece a esto:

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)

Después de hacer eso, tengo dos tablas ( NB se trata de un conjunto de datos diferente al que se muestra en la imagen de pantalla anterior, porque lo transcribí y calculé antes) - la tabla original que muestra el vol para cada par (delta, tenor), y la nueva que muestra el strike para cada par. La nueva tabla tiene el siguiente aspecto:

Ahora tenemos lo suficiente para calibrar una superficie de vol de Heston utilizando los triples (tenor, strike, vol) de cada opción observada (nb. también tendrás que ajustar las curvas de los tipos nacionales y extranjeros, pero eso es otra historia) - para mis opciones de arriba, la superficie tiene este aspecto:

Aquí hay un ejemplo de código (los datos de arriba están codificados en la parte superior) que generará la superficie vol de arriba para usted:

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