# Quantlib issue with BlackVarianceSurface diffusing with the wrong vol when there are either holes or arbitrages in early maturities

I'm new to Quantlib in Python and I'm running into a quite awkward situation.

I have a vol surface from the market on the SPX index. Not all strikes/maturities are populated. In addition, there could potentially be arbitrages in short maturities since we have lots of daily options and I am not correcting yet for arbitrages.

Yet when building a BlackVarianceSurface with this data and price a vanilla option, I expect that option to be priced with a constant vol interpolated from the grid.

When I pass the surface object to the pricer I am getting a very different price from simply passing the constant vol taken from the surface. The price with the surface seems to take into account very short term vols and arbitrages. I have played with the surface deleting very short term vols and/or far strikes and the price of the option actually changes (?!?). When I say the price is different, I mean instead of getting the price of an option with 22 implied vol, I'm getting that of an option with a 38 implied vol instead.

The price with the constant vol is exactly what I expect and it gives me the proper break even vol. So I could clearly interpolated from the grid for every price and get the proper price.

Note that the option (strike,maturity) exists in the grid given in input so it's neither an interpolation nor an extrapolation issue (although playing with these settings slightly changes the price too but to a much lesser extent and compatible with interpolation noise).

It seems very counterproductive to have to interpolate manually from the surface before every single pricing as opposed to passing the surface object and let QLIB get the proper vol under the cover.

Could you help me understand exactly how QLIB handles the surface and what diffusion process is built from that surface ? We know it's not a local vol process that is built. Consequently why is it not a deterministic process based on the interpolated vol ?

Thank you very very much.

Here's some sample code that compares a constant vol with a surface vol pricing. Notice how the theta changes between the 2 cases leading to a very different break even vol.

import pandas as pd, QuantLib as ql, math, numpy as np
from datetime import datetime

def from_YYYYMMDD_to_ql_date(ymd):
ymd = int(ymd)
year = math.floor(ymd / 10000)
month = int(ymd / 100) % 100
day = ymd % 100

ql_date = ql.Date(day, month, year)
return ql_date

def run_test():

####################################################################
####            Global Variables
####################################################################
pricing_date = datetime.strptime('20220916', '%Y%m%d')
ql_pricing_date = ql.Date(16, 9, 2022)
calendar = ql.UnitedStates(ql.UnitedStates.Settlement)
day_counter = ql.ActualActual()
spot = 3876.00
spot_handle = ql.QuoteHandle(ql.SimpleQuote(spot))
flat_ts = ql.YieldTermStructureHandle(ql.FlatForward(ql_pricing_date, 0, day_counter))
dividend_yield = ql.YieldTermStructureHandle(ql.FlatForward(ql_pricing_date, 0, day_counter))
K = 3800
ql_date_expiry = ql.Date(16, 12, 2022)
exercise = ql.EuropeanExercise(ql_date_expiry)
payoff = ql.PlainVanillaPayoff(ql.Option.Put, strike=K)

####################################################################
####            Build Vol surface
####################################################################
list_opts = []
expirations = [datetime.strptime('20220921', '%Y%m%d'), datetime.strptime('20221118', '%Y%m%d'), datetime.strptime('20221216', '%Y%m%d'), ]
list_opts.append([pricing_date, expirations[0], 3800., 45.])
list_opts.append([pricing_date, expirations[0], 3900., 40.])

list_opts.append([pricing_date, expirations[1], 3800., 27.])
list_opts.append([pricing_date, expirations[1], 3900., 25.])

list_opts.append([pricing_date, expirations[2], 3800., 37.])
list_opts.append([pricing_date, expirations[2], 3900., 31.])

df = pd.DataFrame(list_opts, columns=['date', 'maturity', 'strike', 'vol'])

df['ymd'] = df['maturity'].dt.strftime('%Y%m%d')
df['q_exp_dates'] = df['ymd'].apply(from_YYYYMMDD_to_ql_date)

strikes = sorted(df.strike.unique())
expirations = sorted(df['q_exp_dates'].unique())

df['strike_idx'] = df['strike'].apply(lambda x: strikes.index(x))
df['expiry_idx'] = df['q_exp_dates'].apply(lambda x: expirations.index(x))

volMat_np = np.full((len(strikes), len(expirations)), np.nan)

for strk, exp, vol in zip(df['strike_idx'], df['expiry_idx'], df['vol']):
volMat_np[strk][exp] = vol / 100

df_vals = pd.DataFrame(volMat_np, index=strikes, columns=expirations)
print(df_vals)

ql_Matrix = ql.Matrix(len(strikes), len(expirations), np.nan)
for i in range(len(strikes)):
for j in range(len(expirations)):
ql_Matrix[i][j] = volMat_np[i, j]

####################################################################
####            We build 1 vol surface and a constant vol = BlackVol(surface)
####################################################################

vol_process = ql.BlackVarianceSurface(ql_pricing_date, calendar, expirations, strikes,
ql_Matrix, day_counter)

option_vol = vol_process.blackVol(ql_date_expiry,K)
constant_vol = ql.BlackConstantVol(ql_pricing_date, calendar, option_vol, day_counter)

####################################################################
####            Build Process
####################################################################

process_surface = ql.BlackScholesMertonProcess(spot_handle, dividendTS=dividend_yield, riskFreeTS=flat_ts,
volTS=ql.BlackVolTermStructureHandle(vol_process))

process_constant = ql.BlackScholesMertonProcess(spot_handle, dividendTS=dividend_yield, riskFreeTS=flat_ts,
volTS=ql.BlackVolTermStructureHandle(constant_vol))

####################################################################
####            Build Option to price
####################################################################
option_surface = ql.DividendVanillaOption(payoff, exercise, [], [])
engine_surface = ql.FdBlackScholesVanillaEngine(process_surface, 200, 200)
option_surface.setPricingEngine(engine_surface)

option_constant = ql.DividendVanillaOption(payoff, exercise, [], [])
engine_constant = ql.FdBlackScholesVanillaEngine(process_constant, 200, 200)
option_constant.setPricingEngine(engine_constant)

prc_cst = option_constant.NPV()
gamma_cst = option_constant.gamma()
theta_cst = option_constant.theta() / 252. # biz day theta
vol_be_sqrt252_cst = round(math.sqrt(-2. * theta_cst / (gamma_cst * spot ** 2)) * math.sqrt(252),4)      # calculate theta-gamma break even vol to estimate what vol is being used by quantlib
print(f'Option details')
print(f'-------------------------------------------------')
print(f'Maturity: {ql_date_expiry}')
print(f'Strike: {K}')
print(f'Vol from surface: {option_vol}')
print(f'')
print(f'-------------------------------------------------')
print(f'Constant vol case')
print(f'-----------------')
print(f'Option price: {prc_cst}')
print(f'Option gamma: {gamma_cst}')
print(f'Option theta: {theta_cst}')
print(f'Break even Vol = {vol_be_sqrt252_cst} should be {option_vol}')

prc_surface = option_surface.NPV()
gamma_surface = option_surface.gamma()
theta_surface = option_surface.theta() / 252. # biz day theta
vol_be_sqrt252_surface = round(math.sqrt(-2. * theta_surface / (gamma_surface * spot ** 2)) * math.sqrt(252),4)      # calculate theta-gamma break even vol to estimate what vol is being used by quantlib
print(f'-------------------------------------------------')
print(f'Vol Surface case')
print(f'----------------')
print(f'Option price: {prc_surface}')
print(f'Option gamma: {gamma_surface}')
print(f'Option theta: {theta_surface}')
print(f'Break even Vol = {vol_be_sqrt252_surface} should be {option_vol}')

if __name__ == '__main__':
run_test()


We get the following results:

            September 21st, 2022  November 18th, 2022  December 16th, 2022
3800.0                  0.45                 0.27                 0.37
3900.0                  0.40                 0.25                 0.31
Option details
-------------------------------------------------
Maturity: December 16th, 2022
Strike: 3800
Vol from surface: 0.37

-------------------------------------------------
Constant vol case
-----------------
Option price: 246.10555275972283
Option gamma: 0.0005462073912642451
Option theta: -2.2348817578128877
Break even Vol = 0.3705 should be 0.37
-------------------------------------------------
Vol Surface case
----------------
Option price: 246.10595637993455
Option gamma: 0.0005462035022605052
Option theta: -3.3101198905503284
Break even Vol = 0.4509 should be 0.37


If you change the vol surface to something which is not arbitrable but still has high very short term vol, the problem persists and even accentuates. In all cases the difference appears on the theta. The price of the option does not change.

list_opts.append([pricing_date, expirations[0], 3800., 45.])
list_opts.append([pricing_date, expirations[0], 3900., 40.])

list_opts.append([pricing_date, expirations[1], 3800., 27.])
list_opts.append([pricing_date, expirations[1], 3900., 25.])

list_opts.append([pricing_date, expirations[2], 3800., 27.])
list_opts.append([pricing_date, expirations[2], 3900., 25.])


We get the following:

Option details
-------------------------------------------------
Maturity: December 16th, 2022
Strike: 3800
Vol from surface: 0.27

-------------------------------------------------
Constant vol case
-----------------
Option price: 170.495477037466
Option gamma: 0.0007462620578669159
Option theta: -1.6258437988242298
Break even Vol = 0.2703 should be 0.27
-------------------------------------------------
Vol Surface case
----------------
Option price: 170.49560335893932
Option gamma: 0.0007462597949764497
Option theta: -4.538109308076177
Break even Vol = 0.4517 should be 0.27

• Can you share some sample code to reproduce the problem? Sep 22, 2022 at 8:36
• Added some sample code for your reference Sep 27, 2022 at 23:55
• Does anybody have the same issue ? Oct 3, 2022 at 12:22
• It wasn't reported before as far as I know. You can try asking on the QuantLib mailing list, where you'll reach a lot more QuantLib users. Thanks for the code, by the way, but I haven't had time to try it out yet. Oct 4, 2022 at 21:28

I have run your code above and can confirm same results. I have also reproduced in C++ and seeing same thing. Initially I thought maybe it's daycount issue (365 v 252) but don't believe that to be the case, even though that ratio would come close to resolving the discrepancy in this specific case - I think that's just coincidental. I also thought there may be an issue with the number of FD grid steps, but the result is relatively insensitive to the # of points used. If you ask the variance surface for the variance at a given strike, it's giving you back the right value. I attached a debugger and think the problem lies somewhere with the interaction between the BlackScholesMertonProcess, the FdmBlackScholesMesher, and the solver. It's deeper into QuantLib than my knowledge extends.

If there is a QuantLib expert willing to take a look, I can clean up the C++ test case and post here if helpful.

@volPMNYC -- if your Oct 4 post is a duplicate of this, can you delete it so we focus the attention/conversation here.

• Thank you @tangent360 for taking the time. Oct 13, 2022 at 19:32
• Thank you. I also thought of a daycount issue but on other examples the difference is too large (see the last example above where there's a ratio of 3 on the theta for a 2m option). I see 2 coincidences too important to be ignored. #1 it happens exclusively on theta and theta is computed from the first node. #2 this happens even more when very short term vols are much higher than the vol used to diffuse in the grid. Consequently I would think FdBlackScholesVanillaEngine probably has some kind of bug on the theta calculation where some short term vol from the curve is used too. Thoughts ? Oct 13, 2022 at 19:44

I believe what's happening is that the difference operator, FdmBlackScholesOp, is calling blackForwardVariance(t1, t2, strike) on your BlackVarianceSurface to get the variance to apply to each time step as it evolves the solution back to t=0 on the FD grid. Because your surface has such pronounced term skew, with Sep expiry trading 18 vols over Nov, the value decay is much higher in the early days of the option life until after the 21sep22 expiry. So even though the localVol flag on the FdBlackScholesVanillaEngine is set to false, it's still aware of the term structure of vol for a given strike, and it's allocating the variance over the life of the option in a non-uniform manner, which sums to the same total variance, and hence the correct option price, but shifts the greeks so that they are consistent with the market term structure.

Put it another way, based on your surface, the 16dec22 p3800 is ~245 and the 20sep22 p3800 is ~40, the calendar price is telling you that your theta will be a lot higher than the analytic BSM theta and the FD model is reflecting this.

An interesting observation, that stumped me for some time, is that advancing the evaluationDate by 1 day and re-running the FD model, will give you the same value change as predicted by the analytic BSM theta, leading one to think that the FD theta is incorrect. However, what should probably happen is for the Dec vols to come down to reflect the loss of value as a result of dropping a high variance day at the start of the option life.

Maybe a QuantLib expert can corroborate this thesis.

Below is the theta evolution in the early stages of the option life taken from the FD grid. Within 9 days the theta has fallen below the analytic BSM theta.

• We get the same price with a constant vol as we get with the surface. It refutes the idea that variance(t1,t2,strike) is used at each node. If it did we’d have a significantly different price as the total variance (the sum of all the “instantaneous” variances) would be different from the constant vol case. But I agree it may happen exclusively for the theta (other Greeks are the same as in the constant vol case). It would be great if a Quantlib contributor could confirm if this is happening exclusively for the theta. ? And if there’s a way to shut down this effect while not in local vol. Oct 16, 2022 at 11:06
• "... the total variance (the sum of all the “instantaneous” variances) would be different from the constant vol case" ... ...this is not correct. The forward variance function is going to return interval variances that, when summed over adjacent steps on [t0, T], will sum to the total variance computed from the IV you provided for T. So the prices will match with no divs and no rates. You can use the FD engine with a flat vol term structure if you want a uniform FD grid, or use the vol surface with the AnalyticEuropeanEngine or Bjerksund-Stensland for an American closed form approximation Oct 17, 2022 at 7:52
• True. You're right (in order to achieve that, it can only do it for 1 strike, the strike of the option meaning there's 1 variance and 1 variance only per column of the grid to guarantee it sums up to the vol(K,T)). Isn't there a way to speed up the process by telling it to NOT compute these semi-local variances and instead use exclusively vol(K,T) ? Oct 17, 2022 at 15:10
• Didn't fully follow what you meant but, yes, it is only using one variance per time step and that is the variance computed from the termstructure you provided for a specific strike ... not a Local Vol approach ... neighboring strikes don't affect it. As for the speed, I wouldn't let it concern me because a) looking up the variance for each time step is a tiny % of the total execution time involved and b) you're using python -- the overhead from that is going to be multiples of the execution time in C++ Oct 17, 2022 at 20:06