[Quantlib-users] Numerical Gamma Problem Using FDDividendAmericanEngine

[Ney S. <ney...@gm...>]

 Hi QuantLib Community,

I am using QuantLib to develop a back-tester, primarily for listed US
options initially.  For now I am using FDDividendAmericanEngine to create
American options with a dividend yield, but I have had problems when
calculating a numerical gamma.  I know finite differences are not
guaranteed to be stable, but in the minimum it is important that I find a
solution before implementing this methodology.

With basic bump/shock methods using relative bumps of 1 bip, the gamma
values I determined were incredibly small and wrong, while using a 100 bip
bump was in an acceptable range when compared to the .gamma() method in the
code, and a European option's gamma when passed the same parameters.  The
rest of the numerical greeks so far have been fine.  I found a quantstack
post flagging the issue, but no work around was presented.  Can anyone
share any insights into the underlying issue and potential solutions,
whether an augmented bump method or way to optimize relative bump size
based on underlying spot price, implied volatiltiy, time to expiry, etc?
Any help is much appreciated.  Below is a script of the code that I am
using, with output for the QuantLib gamma method and my numerical gamma.

Thanks,

Ney

import QuantLib as ql

#option parameters
option_type = ql.Option.Call
expiry = ql.Date(15,1,2016)
strike = 65.0
rf_rate = ql.SimpleQuote(0.01)
spot = ql.SimpleQuote(42.41)
evaluation_date = ql.Date(11,2,2015)
mkt_px = 0.45
dividend_yield = ql.SimpleQuote(0.030417)
implied_vol = ql.SimpleQuote(0.50)

# set calendar, day count, and evaluation date
calendar = ql.UnitedStates()
day_counter = ql.ActualActual()
ql.Settings.instance().evaluationDate = evaluation_date

# set rate, div, and vol curves
rf_rate_curve = ql.FlatForward(0, calendar, ql.QuoteHandle(rf_rate),
day_counter)
dividend_yield_curve = ql.FlatForward(0, calendar,
ql.QuoteHandle(dividend_yield), day_counter)
vol_curve = ql.BlackConstantVol(0, calendar, ql.QuoteHandle(implied_vol),
day_counter)

# create process, set mesh size w/ time_steps & grid_points, and create
engine
process = ql.BlackScholesMertonProcess(ql.QuoteHandle(spot),
                                               ql.YieldTermStructureHandle(
dividend_yield_curve),
                                               ql.YieldTermStructureHandle(
rf_rate_curve),
                                               ql.
BlackVolTermStructureHandle(vol_curve))
time_steps = calendar.businessDaysBetween(evaluation_date, expiry)
grid_points = time_steps
engine = ql.FDDividendAmericanEngine(process, time_steps, grid_points)

# create option
payoff = ql.PlainVanillaPayoff(option_type, strike)
exercise = ql.AmericanExercise(evaluation_date, expiry)
amerOption = ql.VanillaOption(payoff, exercise)
amerOption.setPricingEngine(engine)

# unbumped analytical prices
analytical_price = amerOption.NPV()
spot_ref = spot.value()
analytical_delta = amerOption.delta()
analytical_gamma = amerOption.gamma()

# relative bumps
spot_bump_1_bip = spot_ref * 0.0001
spot_bump_100_bips = spot_ref * 0.01

# bumped up prices
spot.setValue(spot_ref + spot_bump_100_bips)
price_spot_up_100_bips = amerOption.NPV()
spot.setValue(spot_ref + spot_bump_1_bip)
price_spot_up_1_bip = amerOption.NPV()

# bumped down prices
spot.setValue((spot_ref - spot_bump_100_bips))
price_spot_down_100_bips = amerOption.NPV()
spot.setValue((spot_ref - spot_bump_1_bip))
price_spot_down_1_bip = amerOption.NPV()

# spot back to original price
spot.setValue(spot_ref)

# gamma bump functions
numerical_gamma_100_bips = (price_spot_up_100_bips - 2 * analytical_price +
price_spot_down_100_bips) / spot_bump_100_bips**2
numerical_gamma_1_bip = (price_spot_up_1_bip - 2 * analytical_price +
price_spot_down_1_bip) / spot_bump_1_bip**2

print("price_spot_up_100_bips:   ", price_spot_up_100_bips)
print("price_spot_up_1_bip:      ", price_spot_up_1_bip)
print("analytical price:         ", analytical_price)
print("price spot down 1 bip:    ", price_spot_down_1_bip)
print("price spot down 100 bips: ", price_spot_down_100_bips)

print("\nnumerical gamma 100 bips: ", numerical_gamma_100_bips)
print("numerical gamma 1 bip:    ", numerical_gamma_1_bip)