गैर-वर्दी हड़ताल बनाम टेनर ग्रिड से एफएक्स वॉल्यूम सरफेस का इंटरपोलेशन

Aug 16 2020

टीएल, डॉ

मैं मूल्य के साथ एक स्थानीय वॉल्यूम मॉडल बनाने के लिए FX विकल्प उद्धरण को बाजार में उतारने के लिए एक वॉल्यूम सतह फिट करने की कोशिश कर रहा हूं। सूचीबद्ध विकल्पों के विपरीत, जिनमें आम तौर पर स्ट्राइक्स और टेनर्स का एक अच्छा आयताकार ग्रिड होता है, एफएक्स विकल्प ओटीसी का व्यापार करते हैं और उपलब्ध उद्धरण एक समान ग्रिड प्रदान नहीं करते हैं।

गैर-समान ग्रिड पर 2D-प्रक्षेप के लिए लेने के लिए एक समझदार दृष्टिकोण क्या है? विचार मैं था:

  • बिंदुओं का महीन वर्ग ग्रिड बनाएं और उन लोगों के लिए मानों को प्रक्षेपित करें (जैसे scipy.interpolate.griddataनीचे दिखाया गया है), और उस के लिए वॉल्यूम सतह का निर्माण करें (हालांकि यह बेकार लगता है)
  • विकल्प स्ट्राइक में कुछ परिवर्तन लागू करें ताकि उन्हें समान रूप से फैलाने के लिए (बाद के लोगों की तुलना में पहले के किरायेदारों को खींचकर) फिर एक मानक 2D ग्रिड इंटरपोलर का उपयोग किया जा सके

अंततः मैं एक मॉडल का QuantLibउपयोग करना चाहूंगा ql.BlackVarianceSurface, जिसके लिए वर्तमान में आयताकार ग्रिड की आवश्यकता है।

मैं यह सुनना पसंद करूंगा कि लोगों ने किसी भी 2 डी-इंटरपोलेशन खतरों और एक्सट्रपलेशन के मुद्दों सहित क्या दृष्टिकोण अपनाए हैं।

समस्या पर आगे विस्तार से

यहां बाजार द्वारा उद्धृत एफएक्स वॉल्यूम सतह का एक उदाहरण दिया गया है:

एक बार जब इसे स्ट्राइक, टेनॉर, वॉल्यूम) में बदल दिया जाता है, तो स्ट्राइक इस तरह दिखती है:

यह हमें एक गैर-समान ग्रिड प्रदान करता है, जिसे एक 2 डी सतह पर प्लॉट किया जाता है, वे इस तरह दिखते हैं (टीटीई और रूट में):

एक वर्ग ग्रिड का उपयोग scipy.interpolate.griddataऔर द्वि-प्रक्षेपित कास्ट :

जवाब

3 user35980 Aug 16 2020 at 17:54

मैंने कुछ सप्ताह पहले क्वांटलिब अजगर में इन पंक्तियों के साथ कुछ करने की कोशिश की। मुझे लगता है कि आपके दृष्टिकोण की तुलना में थोड़ा अधिक सरल:

  1. एफएक्स वोल्ट (10 डी पुट, 25 डी पुट, एटीएम, 25 डी कॉल, 10 डी नोट) के लिए एक मानक डेल्टा उद्धरण सम्मेलन के साथ शुरू करें
  2. स्ट्राइक सेट प्राप्त करने के लिए विकल्पों की मात्रा की गणना करें (यह एक बड़ी स्ट्राइक सेट होगी क्योंकि प्रत्येक विकल्प की परिपक्वता में मूल स्रोत के मनीनेस कोट्स के अनुरूप अद्वितीय स्ट्राइक होंगे)
  3. प्रत्येक परिपक्वता के लिए हड़ताल के पूर्ण सेट के लिए लापता वॉल्यूम को प्रक्षेपित करें - मैंने क्वांटलिब में BlackVarianceSurface फ़ंक्शन का उपयोग करके ऐसा किया। इस प्रकार मेरे पास परिपक्वता / हड़ताल का पूरा ग्रिड था
  4. मैं अंत में इस डेटा को ले लिया और एक Heston अंशांकन की कोशिश की और एक HestonBlackVolSurface फ़ंक्शन में आउटपुट प्लग किया

परिणाम बहुत अच्छे नहीं थे, क्योंकि हेस्टन के आरोपित खंडों ने वास्तव में सटीकता के साथ मेरे इनपुट स्रोत खंडों को पुन: पेश नहीं किया, लेकिन शायद मेरे खराब अंशांकन और मेरे द्वारा उपयोग किए गए डमी इनपुट स्रोत मूल्यों के साथ ऐसा करना अधिक है। फिर भी यह एक सार्थक अभ्यास था।

मामले में यह उपयोगी हो सकता है मेरा क्वांटलिब कोड नीचे है:

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

अंत में मैंने पाया कि प्रत्येक किराएदार को एक SABR स्माइल फिट करना ( इस उत्तर से एक परिणाम उधार लेना ) एक स्थानीय वॉल्यूम सतह बनाने के लिए पर्याप्त था जो सुचारू रूप से काम करने के लिए पर्याप्त चिकनी और अच्छी तरह से व्यवहार किया गया था। मैंने इसके लिए एक हेस्टन मॉडल भी फिट किया है, और दो सतहें काफी समान दिखती हैं। यहाँ अंतिम कोड है और उत्पन्न फिट (इन भूखंडों को उत्पन्न करने के लिए बहुत नीचे की ओर लंबा स्निपेट आवश्यक है, और इसमें आवश्यक कच्चा डेटा भी शामिल है)

सबसे पहले, प्रत्येक टेनर पर लूपिंग और एक 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()

इस तरह के भूखंडों का एक क्रम उत्पन्न करता है, जिनमें से ज्यादातर ज्यादातर अच्छी तरह से फिट होते हैं:

जो इस तरह दिखने वाले प्रत्येक टेनर में SABR पैरामेट्स बनाता है (इस उदाहरण के लिए मैंने फ्लैट होने के लिए विदेशी और घरेलू छूट घटता है):

फिर मैंने एक स्थानीय वॉल्यूम मॉडल और एक हेस्टन वॉल्यूम मॉडल को कैलिब्रेट किया, जो वास्तव में दोनों एक साथ काफी करीब दिखते हैं:

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

हम उम्मीद करते हैं कि स्थानीय वॉल्यूम मॉडल वैनिलस की कीमत सही ढंग से तय करेंगे, लेकिन अविश्वसनीय वोल्ट डायनामिक्स देते हैं, जबकि हम हेस्टन से बेहतर वॉल्यूम डायनामिक्स देने की अपेक्षा करते हैं, लेकिन वेनिला को इतनी अच्छी कीमत नहीं देते हैं, लेकिन लीवरेज फ़ंक्शन को कैलिब्रेट करके और हेस्टन एनोचैस्टिक स्थानीय वॉल्यूम मॉडल का उपयोग करके हम संभवतः प्राप्त कर सकते हैं। दोनों दुनिया के सर्वश्रेष्ठ - और यह भी एक अच्छी परीक्षा है कि हमने जो स्थानीय वॉल्यूम सतह बनाई है वह अच्छी तरह से व्यवहार की है

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

जो एक अच्छा दिखने वाला उत्तोलन कार्य करता है, जो हर जगह 1 के करीब है (यह दर्शाता है कि कच्चा हेस्टन फिट पहले से ही काफी अच्छा था)

बॉयलरप्लेट कोड ऊपर छवियों को उत्पन्न करने के लिए (एफएक्स डेल्टा-टू-स्ट्राइक रूपांतरण सहित):

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