Interpolación de FX Vol Surface de rumbo no uniforme vs cuadrícula tenor

Aug 16 2020

TL; DR

Estoy tratando de ajustar una superficie de volumen a las cotizaciones de opciones de divisas del mercado para construir un modelo de volumen local con el que fijar el precio. A diferencia de las opciones enumeradas que normalmente tienen una cuadrícula rectangular agradable de strikes y vencimientos, las opciones FX tienden a negociarse OTC y las cotizaciones disponibles no proporcionan una cuadrícula uniforme.

¿Cuál es un enfoque sensato para la interpolación 2D en cuadrículas no uniformes? Las ideas que tuve fueron:

  • Cree una cuadrícula cuadrada más fina de puntos e interpole los valores para esos (por ejemplo, usando que se scipy.interpolate.griddatamuestra a continuación), y construya la superficie de volumen para eso (aunque esto parece un desperdicio)
  • Aplique algo de transformación a los golpes de opción para distribuirlos uniformemente (estirando los tenores anteriores más que los últimos) y luego utilizando un interpolador de cuadrícula 2D estándar

Eventualmente me gustaría construir un modelo en el QuantLibuso ql.BlackVarianceSurface, que actualmente requiere una cuadrícula rectangular de vols.

Me encantaría saber qué enfoques han adoptado las personas, incluidos los peligros de la interpolación 2D y los problemas de extrapolación.

Más detalles sobre el problema

A continuación se muestra un ejemplo de una superficie FX vol cotizada por el mercado:

Una vez que esto se convierte en (strike, tenor, vol) triples, los strikes se verán así:

Esto nos da una cuadrícula no uniforme de vols, trazada en una superficie 2D y se ven así (en tte y en raíz tte):

Transmitir a una cuadrícula usando scipy.interpolate.griddatay bi-interpolado:

Respuestas

3 user35980 Aug 16 2020 at 17:54

Probé algo en este sentido en Quantlib python hace unas semanas. Un poco más simple en comparación con su enfoque, creo:

  1. Comience con una convención de cotización delta estándar para FX vols (put 10D, put 25D, cajero automático, llamada 25D, llamada 10D)
  2. calcular el valor monetario de las opciones para obtener el conjunto de acciones (este será un conjunto de acciones grande ya que cada vencimiento de la opción tendrá acciones únicas correspondientes a las cotizaciones monetarias de la fuente original)
  3. interpolar los vols faltantes para el conjunto completo de strikes para cada vencimiento; hice esto usando la función BlackVarianceSurface en Quantlib. Por lo tanto, tenía una cuadrícula completa de vencimientos / huelgas.
  4. Finalmente tomé estos datos y probé una calibración de Heston y conecté la salida a una función HestonBlackVolSurface

Los resultados no fueron excelentes ya que Heston implicaba que los vols no reproducían realmente los vols de mi fuente de entrada con precisión, pero eso probablemente tenga más que ver con mi mala calibración y los valores de fuente de entrada ficticios que usé. Sin embargo, fue un ejercicio que valió la pena.

En caso de que pueda ser útil, mi código Quantlib está a continuación:

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
4 StackG Sep 30 2020 at 12:59

Al final, descubrí que adaptar una sonrisa SABR a cada tenor (tomando prestado un resultado de esta respuesta ) era suficiente para construir una superficie de volumen local que era lo suficientemente suave y de buen comportamiento para construir una superficie de variación que funcionaba bien. También le instalé un modelo Heston, y las dos superficies se ven bastante similares. Aquí está el código final y los ajustes generados (el fragmento largo en la parte inferior es necesario para generar estos gráficos y también contiene los datos sin procesar necesarios)

En primer lugar, recorriendo cada tenor y colocando 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 gráficos como este, que en su mayoría se ajustan bastante bien:

que genera parámetros SABR en cada tenor con este aspecto (para este ejemplo, establecí curvas de descuento nacionales y extranjeras planas):

Luego calibré un modelo de volumen local y un modelo de volumen de Heston, que en realidad se ven bastante juntos:

# 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 volumen local fije el precio de las vainillas correctamente pero proporcione una dinámica de volumen poco realista, mientras que esperamos que Heston proporcione una mejor dinámica de volumen pero que no fije tan bien las vainillas, pero calibrando una función de apalancamiento y utilizando un modelo de volumen local estocástico de Heston , posiblemente podamos obtener lo mejor de ambos mundos, y esta también es una buena prueba de que la superficie de volumen local 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 aspecto agradable, que está cerca de 1 en todas partes (lo que indica que el ajuste de Heston en bruto ya era bastante bueno)

Código repetitivo para generar las imágenes anteriores (incluida la conversión FX delta-to-strike):

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