QuantLib: Analytical Greeks and Numerical Greeks do not match?

Question

I use the Black Scholes Merton (BSM) model from QuantLib to calculate Call options price and its analytical Greeks. I also manually calculate its Numerical Greeks (Theta, Vega), but the results do not match. Does anyone know what went wrong with my implementation? Thanks.

import numpy as np
import pandas as pd
import QuantLib as ql

    
def bsm_quantlib(numerical=False):  
    spot = 100          # spot price
    strike = 120        # strike price
    rf_rate = 0.035     # risk-free annual interest rate
    vol = 0.16          # annual volatility or sigma
    div = 0.01          # annual dividend rate

    eval_date = ql.Date(7, 1, 2023)
    expiry_date = ql.Date(7, 1, 2024)    
    ql.Settings.instance().evaluationDate = eval_date

    calendar = ql.UnitedStates(ql.UnitedStates.NYSE)
    day_counter = ql.Actual365Fixed() 
    payoff = ql.PlainVanillaPayoff(ql.Option.Call, strike)

    spot_quote = ql.SimpleQuote(spot)
    rf_quote = ql.SimpleQuote(rf_rate)
    vol_quote = ql.SimpleQuote(vol)

    spot_handle = ql.QuoteHandle(spot_quote)
    vol_handle = ql.QuoteHandle(vol_quote)
    rf_handle = ql.QuoteHandle(rf_quote)
    div_handle = ql.QuoteHandle(ql.SimpleQuote(div))

    dividend_yield = ql.YieldTermStructureHandle(
        ql.FlatForward(0, calendar, div_handle , day_counter))

    risk_free_curve = ql.YieldTermStructureHandle(
        ql.FlatForward(0, calendar, rf_handle, day_counter))

    volatility = ql.BlackVolTermStructureHandle(
        ql.BlackConstantVol(0, calendar, vol_handle, day_counter))

    engine = ql.AnalyticEuropeanEngine(
        ql.BlackScholesMertonProcess(
            spot_handle, dividend_yield, risk_free_curve, volatility))
    
    exercise = ql.EuropeanExercise(expiry_date)
    option = ql.VanillaOption(payoff, exercise)
    option.setPricingEngine(engine)

    greeks = (
        numerical_greeks(option, spot_quote, vol_quote, eval_date) 
        if numerical 
        else analytical_greeks(option))
    
    return greeks | dict(
        call_price=option.NPV(), 
        spot=spot, 
        strike=strike, 
        tau=(expiry_date - eval_date) / 365.0, 
        riskfree_rate=rf_rate, 
        volatility=vol,
        dividend=div)

where analytical_greeks() and numerical_greeks() are

def analytical_greeks(option):
    return dict(
        Greeks='QuantLib Analytical',
        delta=option.delta(),
        theta=option.thetaPerDay(),
        vega=option.vega()/100)
    
            
def numerical_greeks(option, spot_quote, vol_quote, eval_date):
    # delta
    p0 = option.NPV()
    s0 = spot_quote.value()
    v0 = vol_quote.value()
    
    h = 0.01
    spot_quote.setValue(s0 + h)
    pplus = option.NPV()
    spot_quote.setValue(s0 - h)
    pminus = option.NPV()
    spot_quote.setValue(s0)    
    delta = (pplus - pminus) / (2*h)

    # vega
    vol_quote.setValue(v0 + h)
    pplus = option.NPV()
    vol_quote.setValue(v0)
    vega = (pplus - p0) / h

    # theta
    ql.Settings.instance().evaluationDate = eval_date + 365
    pplus = option.NPV()
    ql.Settings.instance().evaluationDate = eval_date
    theta = (pplus - p0)
    
    return dict(
        Greeks='QuantLib Numerical',
        delta=delta,
        theta=theta/365,
        vega=vega/100)

For comparison, I also include the calculation using the py_vollib package. FYI, I had to modify the date slightly so that alltau=1.0

import numpy as np
import pandas as pd
import datetime as dt

from py_vollib.black_scholes_merton import black_scholes_merton as bsm

from py_vollib.black_scholes_merton.greeks.numerical import (
    delta as delta_bsm_n, 
    theta as theta_bsm_n, 
    vega as vega_bsm_n)

from py_vollib.black_scholes_merton.greeks.analytical import (
    delta as delta_bsm_a, 
    theta as theta_bsm_a, 
    vega as vega_bsm_a)


def bsm_vollib(numerical=False): 
    flag = 'c'          # call options
    spot = 100          # spot price
    strike = 120        # strike price
    rf_rate = 0.035     # risk-free annual interest rate
    vol = 0.16          # annual volatility or sigma
    div = 0.01          # annual dividend rate

    eval_date = dt.datetime(2023, 7, 2)
    expiry_date = dt.datetime(2024, 7, 1)    
    tau = (expiry_date - eval_date).days / 365.0
    
    price = bsm(flag, spot, strike, tau, rf_rate, vol, div)
    if numerical:     
        greeks = dict(
            Greeks='Vollib Analytical',
            delta=delta_bsm_n(flag, spot, strike, tau, rf_rate, vol, div),
            theta=theta_bsm_n(flag, spot, strike, tau, rf_rate, vol, div),
            vega=vega_bsm_n(flag, spot, strike, tau, rf_rate, vol, div))           
    else:
        greeks = dict(
            Greeks='Vollib Numerical',
            delta=delta_bsm_a(flag, spot, strike, tau, rf_rate, vol, div),
            theta=theta_bsm_a(flag, spot, strike, tau, rf_rate, vol, div),
            vega=vega_bsm_a(flag, spot, strike, tau, rf_rate, vol, div))           

    return greeks | dict(
        call_price=price,
        spot=spot,
        strike=strike,
        tau=tau,
        riskfree_rate=rf_rate,
        volatility=vol,
        dividend=div)

Here is the result. The theta and vega from py_vollib Numerical and Analytical greek are almost identical. But with QuantLib, they are off quite a bit. Any idea why?

pd.DataFrame([
    bsm_quantlib(numerical=False),
    bsm_quantlib(numerical=True),
    bsm_vollib(numerical=False),
    bsm_vollib(numerical=True),
]).round(4)

Answer 1

How close the numerical Greeks are to the analytic Greeks depends on how one calculates them. The analytic Greeks correspond to the mathematical derivatives, so in general smaller increments should give you a better approximation; for instance, if I set h to 0.001 in your code (i.e., 10 times smaller) I get a the numerical vega which is much closer to the analytic one, and if I set h to 0.0001 (100 times smaller) the numerical vega is the same as the analytic one within the 4 decimals you're using for output.

For theta, you can't go smaller than 1 day, so you can try using 1 day instead of 365 when calculating it. In this particular case, however, you would get a numerical theta of 0, because of a snag in your setup. You're passing to the curves 0 settlement days and the NYSE calendar, and you're setting the evaluation date to January 7th 2023, which happens to be a Saturday. This means that the reference date of the curve is rolled to the next Monday (and, therefore, that you're actually calculating the formula with tau < 1). If you increase the evaluation date by 1 day to calculate the theta, you'll set it to the 8th, Sunday, which also gets rolled to the next Monday, so there's no change in price.

If you want to reproduce the analytic results, you can either (a) change your evaluation date to a business date or (b) use calendar = ql.NullCalendar() instead, which considers all days as business days. I tried (b) in your code, together with shifting the evaluation date by 1 day, and I got a numerical theta which is the same as the analytic one within the 4 decimals in the output.