QuantLib: ¿Las griegas analíticas y las griegas numéricas no coinciden?

Question

Utilizo el modelo Black Scholes Merton (BSM) de QuantLib para calcular el precio de las opciones Call y sus griegas analíticas. También calculo manualmente sus griegas numéricas ( Theta, Vega ), pero los resultados no coinciden. ¿Alguien sabe qué ha fallado en mi aplicación? Gracias.

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)

donde analytical_greeks() y numerical_greeks() son

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)

Para comparar, también incluyo el cálculo utilizando el py_vollib paquete. Para su información, tuve que modificar ligeramente la fecha para que todos los tau=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)

He aquí el resultado. El theta y vega de py_vollib El griego numérico y el griego analítico son casi idénticos. Pero con QuantLib, se desvían bastante. ¿Alguna idea de por qué?

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

Answer 1

Lo cerca que estén las griegas numéricas de las griegas analíticas depende de cómo se calculen. Las griegas analíticas corresponden a las derivadas matemáticas, por lo que, en general, los incrementos más pequeños deberían proporcionar una mejor aproximación; por ejemplo, si establezco h a 0.001 en su código (es decir, 10 veces más pequeño) obtengo una vega numérica mucho más cercana a la analítica, y si establezco h a 0.0001 (100 veces menor) la vega numérica es la misma que la analítica dentro de los 4 decimales que estás utilizando para la salida.

En el caso de theta, no puede ser inferior a 1 día, por lo que puede intentar utilizar 1 día en lugar de 365 al calcularlo. En este caso concreto, sin embargo, obtendrías un theta numérico de 0, debido a un inconveniente en tu configuración. Estás pasando a las curvas 0 días de liquidación y el calendario NYSE, y estás fijando la fecha de evaluación en el 7 de enero de 2023, que resulta ser un sábado. Esto significa que la fecha de referencia de la curva se pasa al lunes siguiente (y, por tanto, que en realidad estás calculando la fórmula con tau < 1). Si aumenta la fecha de evaluación en 1 día para calcular theta, la fijará en el día 8, domingo, que también se traslada al lunes siguiente, por lo que no hay cambios en el precio.

Si desea reproducir los resultados analíticos, puede (a) cambiar su fecha de evaluación por una fecha comercial o (b) utilizar calendar = ql.NullCalendar() que considera todos los días como hábiles. He probado (b) en tu código, junto con el desplazamiento de la fecha de evaluación en 1 día, y he obtenido un theta numérico que es el mismo que el analítico dentro de los 4 decimales en la salida.