Interpolación de la superficie FX Vol a partir de una retícula no uniforme de huelga contra tenor

Question

TL;DR

Estoy tratando de ajustar una superficie de vol a las comillas de las opciones FX del mercado con el fin de construir un modelo de vol local para fijar el precio. A diferencia de las opciones que cotizan en bolsa, que suelen tener una bonita cuadrícula rectangular de strikes y tenores, las opciones de divisas tienden a negociarse OTC y las comillas disponibles no proporcionan una cuadrícula uniforme.

¿Cuál es el enfoque sensato para la interpolación 2D en mallas no uniformes? Las ideas que tenía eran:

• Crear una cuadrícula más fina de puntos e interpolar los valores de los mismos (por ejemplo, utilizando scipy.interpolate.griddata que se muestra a continuación), y construir la superficie superficie de vol para eso (aunque esto parece un desperdicio)
• Aplique un poco de transformación a los strikes de las opciones para uniformemente (estirando más los primeros tenores que los últimos) que los posteriores) y luego utilizar un interpolador de cuadrícula 2D estándar

Con el tiempo me gustaría construir un modelo en QuantLib utilizando ql.BlackVarianceSurface que actualmente requiere una cuadrícula rectangular de vols.

Me gustaría saber qué enfoques ha adoptado la gente, incluidos los peligros de la interpolación en 2D y los problemas de extrapolación.

Más detalles sobre el problema

Este es un ejemplo de una superficie de voltaje de divisas cotizada por el mercado:

Una vez que esto es convertidos en triples (strike, tenor, vol) las huelgas son algo así:

Esto nos da una malla no uniforme de vols, graficada en una superficie 2D se ven así (en tte y en root tte):

Fundido a una rejilla cuadrada utilizando scipy.interpolate.griddata y bi-interpolado:

Answer 1

Al final descubrí que adaptar una sonrisa SABR a cada tenor (tomando prestado un resultado de esta respuesta ) fue suficiente para construir una superficie local de vol que fuera lo suficientemente suave y bien comportada como para construir una superficie de varianza que funcionara bien. También ajusté un modelo de Heston, y las dos superficies se parecen bastante. Aquí está el código final y los ajustes generados (el fragmento largo en la parte inferior se requiere para generar estos gráficos, y también contiene los datos en bruto necesarios)

En primer lugar, se hace un bucle sobre cada tenor y se coloca una sonrisa SABR:

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

genera una secuencia de tramas como ésta, que en su mayoría encajan bastante bien:

que genera parámetros SABR en cada tenor con el siguiente aspecto (para este ejemplo he establecido que las curvas de descuento extranjeras y nacionales sean planas):

Luego calibré un modelo de vol local y un modelo de vol de Heston, que en realidad ambos se parecen bastante:

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

Esperamos que el modelo de vol. local valore correctamente las vainillas pero dé una dinámica de vol. poco realista, mientras que esperamos que Heston dé una mejor dinámica de vol. pero no valore tan bien las vainillas, pero calibrando una función de apalancamiento y utilizando un Modelo estocástico de vol local de Heston posiblemente podamos obtener lo mejor de ambos mundos - y esto es también una buena prueba de que la superficie local de vol que hemos creado se comporta bien

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

que produce una función de apalancamiento de buen aspecto, que se acerca a 1 en todas partes (lo que indica que el ajuste de Heston en bruto ya era bastante bueno)

Código de base para generar las imágenes anteriores (incluyendo la conversión de delta a golpe de FX):

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

Answer 2

Hace unas semanas intenté algo parecido en Quantlib python. Un poco más simple en comparación con su enfoque, creo:

1. empezar con una convención de comilla delta estándar para vols de FX (10D puts, 25D puts,ATM,25D call, 10D call)
2. Calcular el valor monetario de las opciones para obtener el conjunto de precios de ejercicio (será un conjunto de precios de ejercicio grande, ya que cada vencimiento de la opción tendrá precios de ejercicio únicos correspondientes a las comillas del valor monetario de la fuente original).
3. interpolar los vols que faltan para el conjunto completo de strikes para cada vencimiento - lo hice usando la función BlackVarianceSurface en Quantlib. De este modo, tenía una cuadrícula completa de vencimientos/contratiempos
4. Finalmente tomé estos datos y probé una calibración Heston y conecté la salida a una función HestonBlackVolSurface

Los resultados no fueron muy buenos, ya que los vols implícitos de Heston no reprodujeron realmente los vols de mi fuente de entrada con exactitud, pero eso probablemente tenga más que ver con mi mala calibración y los valores de la fuente de entrada ficticia que utilicé. Sin embargo, fue un ejercicio que valió la pena.

En caso de que pueda ser útil mi código Quantlib está abajo:

def deltavolquotes(ccypair,fxcurve):

from market import curveinfo

sheetname = ccypair + '_fx_volcurve'
df = pd.read_excel('~/iCloud/python_stuff/finance/marketdata.xlsx', sheet_name=sheetname)
curveinfo = curveinfo(ccypair, 'fxvols')
calendar = curveinfo.loc['calendar', 'fxvols']
daycount = curveinfo.loc['curve_daycount', 'fxvols']
settlement = curveinfo.loc['curve_sett', 'fxvols']
flat_vol = ql.SimpleQuote(curveinfo.loc['flat_vol', 'fxvols'])
flat_vol_shift = ql.SimpleQuote(0)
used_flat_vol = ql.CompositeQuote(ql.QuoteHandle(flat_vol_shift), ql.QuoteHandle(flat_vol), f)
vol_shift = ql.SimpleQuote(0)
calculation_date = fxcurve.referenceDate()
settdate = calendar.advance(calculation_date, settlement, ql.Days)

date_periods = df[ccypair].tolist()
atm = [ql.CompositeQuote(ql.QuoteHandle(vol_shift), ql.QuoteHandle(ql.SimpleQuote(i)), f) for i in
       df['ATM'].tolist()]
C25 = [ql.CompositeQuote(ql.QuoteHandle(vol_shift), ql.QuoteHandle(ql.SimpleQuote(i)), f) for i in
       df['25C'].tolist()]
P25 = [ql.CompositeQuote(ql.QuoteHandle(vol_shift), ql.QuoteHandle(ql.SimpleQuote(i)), f) for i in
       df['25P'].tolist()]
C10 = [ql.CompositeQuote(ql.QuoteHandle(vol_shift), ql.QuoteHandle(ql.SimpleQuote(i)), f) for i in
       df['10C'].tolist()]
P10 = [ql.CompositeQuote(ql.QuoteHandle(vol_shift), ql.QuoteHandle(ql.SimpleQuote(i)), f) for i in
       df['10P'].tolist()]
dates = [calendar.advance(settdate, ql.Period(i)) for i in date_periods]
yearfracs = [daycount.yearFraction(settdate, i) for i in dates]
dvq_C25 = [ql.DeltaVolQuote(0.25, ql.QuoteHandle(i), j, 0) for i, j in zip(C25, yearfracs)]
dvq_P25 = [ql.DeltaVolQuote(-0.25, ql.QuoteHandle(i), j, 0) for i, j in zip(P25, yearfracs)]
dvq_C10 = [ql.DeltaVolQuote(0.10, ql.QuoteHandle(i), j, 0) for i, j in zip(C10, yearfracs)]
dvq_P10 = [ql.DeltaVolQuote(-0.10, ql.QuoteHandle(i), j, 0) for i, j in zip(P10, yearfracs)]

info=[settdate,calendar,daycount,df,used_flat_vol,vol_shift,flat_vol_shift,date_periods]

return atm,dvq_C25,dvq_P25,dvq_C10,dvq_P10,dates,yearfracs,info

def fxvolsurface(ccypair,FX,fxcurve,curve):

atm,dvq_C25,dvq_P25,dvq_C10,dvq_P10,dates,yearfracs,info = deltavolquotes(ccypair,fxcurve)
settdate = info[0]
calendar=info[1]
daycount=info[2]
df=info[3]
used_flat_vol=info[4]
vol_shift=info[5]
flat_vol_shift=info[6]
date_periods=info[7]

blackdc_C25=[ql.BlackDeltaCalculator(ql.Option.Call,j.Spot,FX.value(),
                                   fxcurve.discount(i)/fxcurve.discount(settdate),
                                   curve.discount(i)/curve.discount(settdate),
                                   j.value()*(k**0.5))
                                   for i,j,k in zip(dates,dvq_C25,yearfracs)]
blackdc_C10=[ql.BlackDeltaCalculator(ql.Option.Call,j.Spot,FX.value(),
                                   fxcurve.discount(i)/fxcurve.discount(settdate),
                                   curve.discount(i)/curve.discount(settdate),
                                   j.value()*(k**0.5))
                                   for i,j,k in zip(dates,dvq_C10,yearfracs)]
blackdc_P25=[ql.BlackDeltaCalculator(ql.Option.Put,j.Spot,FX.value(),
                                   fxcurve.discount(i)/fxcurve.discount(settdate),
                                   curve.discount(i)/curve.discount(settdate),
                                   j.value()*(k**0.5))
                                   for i,j,k in zip(dates,dvq_P25,yearfracs)]
blackdc_P10=[ql.BlackDeltaCalculator(ql.Option.Put,j.Spot,FX.value(),
                                   fxcurve.discount(i)/fxcurve.discount(settdate),
                                   curve.discount(i)/curve.discount(settdate),
                                   j.value()*(k**0.5))
                                   for i,j,k in zip(dates,dvq_P10,yearfracs)]
C25_strikes=[i.strikeFromDelta(0.25) for i in blackdc_C25]
C10_strikes=[i.strikeFromDelta(0.10) for i in blackdc_C10]
P25_strikes=[i.strikeFromDelta(-0.25) for i in blackdc_P25]
P10_strikes=[i.strikeFromDelta(-0.10) for i in blackdc_P10]
ATM_strikes=[i.atmStrike(j.AtmFwd) for i,j in zip(blackdc_C25,dvq_C25)]
strikeset=ATM_strikes+C25_strikes+C10_strikes+P25_strikes+P10_strikes
strikeset.sort()
hestonstrikes=[P10_strikes,P25_strikes,ATM_strikes,C25_strikes,C10_strikes]
hestonvoldata=[df['10P'].tolist(),df['25P'].tolist(),df['ATM'].tolist(),df['25C'].tolist(),df['10C'].tolist()]

volmatrix=[]
for i in range(0,len(atm)):
    volsurface=ql.BlackVolTermStructureHandle(ql.BlackVarianceSurface(settdate,calendar,[dates[i]],
                                [P10_strikes[i],P25_strikes[i],ATM_strikes[i],C25_strikes[i],C10_strikes[i]],
                                [[dvq_P10[i].value()],[dvq_P25[i].value()],[atm[i].value()],[dvq_C25[i].value()],
                                 [dvq_C10[i].value()]],
                                daycount))
    volmatrix.append([volsurface.blackVol(dates[i],j,True) for j in strikeset])
volarray=np.array(volmatrix).transpose()
matrix = []
for i in range(0, volarray.shape[0]):
    matrix.append(volarray[i].tolist())
fxvolsurface=ql.BlackVolTermStructureHandle(
    ql.BlackVarianceSurface(settdate,calendar,dates,strikeset,matrix,daycount))

'''
process = ql.HestonProcess(fxcurve, curve, ql.QuoteHandle(FX), 0.01, 0.5, 0.01, 0.1, 0)
model = ql.HestonModel(process)
engine = ql.AnalyticHestonEngine(model)
print(model.params())
hmh = []
for i in range(0,len(date_periods)):
    for j in range(0,len(hestonstrikes)):
        helper=ql.HestonModelHelper(ql.Period(date_periods[i]), calendar, FX.value(),hestonstrikes[j][i],
                                    ql.QuoteHandle(ql.SimpleQuote(hestonvoldata[j][i])),fxcurve,curve)
        helper.setPricingEngine(engine)
        hmh.append(helper)
lm = ql.LevenbergMarquardt()
model.calibrate(hmh, lm,ql.EndCriteria(500, 10, 1.0e-8, 1.0e-8, 1.0e-8))
vs = ql.BlackVolTermStructureHandle(ql.HestonBlackVolSurface(ql.HestonModelHandle(model)))
vs.enableExtrapolation()'''

flatfxvolsurface = ql.BlackVolTermStructureHandle(
    ql.BlackConstantVol(settdate, calendar, ql.QuoteHandle(used_flat_vol), daycount))

fxvoldata=pd.DataFrame({'10P strike':P10_strikes,'25P strike':P25_strikes,'ATM strike':ATM_strikes,
                        '25C strike':C25_strikes,'10C strike':C10_strikes,'10P vol':df['10P'].tolist(),
                        '25P vol':df['25P'].tolist(),'ATM vol':df['ATM'].tolist(),
                        '25C vol':df['25C'].tolist(),'10C vol':df['10C'].tolist()})
fxvoldata.index=date_periods

fxvolsdf=pd.DataFrame({'fxvolsurface':[fxvolsurface,flatfxvolsurface],'fxvoldata':[fxvoldata,None]})
fxvolsdf.index=['surface','flat']
fxvolshiftsdf=pd.DataFrame({'fxvolshifts':[vol_shift,flat_vol_shift]})
fxvolshiftsdf.index=['surface','flat']

return fxvolshiftsdf,fxvolsdf