Interpolating FX Volatility Surfaces from Irregular Market Quotes

# Interpolation of FX Vol Surface from non-uniform strike vs tenor grid


# Interpolation of FX Vol Surface from non-uniform strike vs tenor grid












TL;DR

I'm trying to fit a vol surface to market FX options quotes in order to build a local vol model to price with. Unlike listed options that typically have a nice rectangular grid of strikes and tenors, FX options tend to trade OTC and the quotes available don't provide a uniform grid.

What is a sensible approach to take for 2D-interpolation on non-uniform grids? Ideas I had were:

- Create a finer square grid of points and interpolate values for those (eg. using `scipy.interpolate.griddata` shown below), and build the vol surface for that (although this seems wasteful)

- Apply some transform to the option strikes to spread them out uniformly (stretching out the earlier tenors more than the later ones) then using a standard 2D grid interpolator

Eventually I'd like to build a model in `QuantLib` using `ql.BlackVarianceSurface`, which currently requires a rectangular grid of vols.

I'd love to hear what approaches people have taken, including any 2D-interpolation dangers, and issues of extrapolation.

Further detail on the problem

Here is an example of an FX vol surface quoted by the market:

Once this is converted into (strike, tenor, vol) triples the strikes look something like this:

This gives us a non-uniform grid of vols, plotted on a 2D surface they look like this (in tte and in root tte):

Cast to a square grid using scipy.interpolate.griddata and bi-interpolated:

## Answer by user35980 (score 7, accepted)

https://quant.stackexchange.com/a/57383

I tried something along these lines in Quantlib python a few weeks ago. Slightly more simple compared to your approach I think:

- start with a standard delta quote convention for FX vols (10D puts, 25D puts,ATM,25D call, 10D call)

- calculate the moneyness of the options to obtain the strike set (this will be a large strike set since each option maturity will have unique strikes corresponding to the moneyness quotes of the original source)

- interpolate the missing vols for the full set of strikes for each maturity - I did this using the BlackVarianceSurface function in Quantlib. Thus I had a full grid of maturities/strikes

- I finally took this data and tried a Heston calibration and plugged the output into a HestonBlackVolSurface function

The results weren't great since the Heston implied vols didn't really reproduce my input source vols with accuracy but that's probably more to do with my poor calibration and the dummy input source values I used. Nevertheless it was a worthwhile exercise.

In case it may be helpful my Quantlib code is below:

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

## Answer by StackG (score 8)

https://quant.stackexchange.com/a/58349

In the end I found that fitting a SABR smile to each tenor (borrowing a result from this answer) was sufficient to build a local vol surface that was smooth and well-behaved enough to build a variance surface worked nicely. I also fitted a Heston model to it, and the two surfaces do look fairly similar. Here is the final code and the fits generated (the long snippet at the very bottom is required to generate these plots, and also contains the raw data required)

Firstly, looping over each tenor and fitting a SABR smile:

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

generates a sequence of plots like this, all of which mostly fit quite well:

which generates SABR params at each tenor looking like this (for this example I've set foreign and domestic discount curves to be flat):

Then I calibrated a local vol model and a Heston vol model, which actually both look quite close together:

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

We expect the local vol model to price vanillas correctly but give unrelistic vol dynamics, while we expect Heston to give better vol dynamics but not price vanillas so well, but by calibrating a leverage function and using a Heston stochastic local vol model we can possibly get the best of both worlds - and this is also a good test that the local vol surface we've created is well behaved

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

which produces a nice looking leverage function, which is close to 1 everywhere (indicating that the raw Heston fit was already quite good)

Boilerplate code to generate above images (including the FX delta-to-strike conversion):

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