From 10fc864471501ca2a2400ff93471cb5c51f76201 Mon Sep 17 00:00:00 2001 From: mkzung <103102868+mkzung@users.noreply.github.com> Date: Wed, 23 Sep 2026 01:55:48 +0100 Subject: [PATCH] Keep the digital American at-expiry price accurate at low volatility AmericanPayoffAtExpiry drops the (H/S)^(2 mu) N(D2) term whenever N(D2) underflows to zero. At low volatility the power it multiplies is large enough that the product is not negligible: at vol 0.1%, r -2% and q 3%, a down-and-in put on a barrier at 95 paying 10 at expiry priced at 1.0002 against 1.0176. Where N(D2) underflows or the power overflows, the product is now formed in logs, with the asymptotic tail of log N. Over 1152 options with volatility from 10% down to 1e-6, master had 4 values more than 1e-6 away from the closed form in 50-digit arithmetic, up to 1.7%; with this change the largest difference is 1.6e-8. --- ql/pricingengines/americanpayoffatexpiry.cpp | 28 ++++++++-- test-suite/digitaloption.cpp | 59 ++++++++++++++++++++ 2 files changed, 83 insertions(+), 4 deletions(-) diff --git a/ql/pricingengines/americanpayoffatexpiry.cpp b/ql/pricingengines/americanpayoffatexpiry.cpp index b89f0ff62da..36294806f4e 100644 --- a/ql/pricingengines/americanpayoffatexpiry.cpp +++ b/ql/pricingengines/americanpayoffatexpiry.cpp @@ -19,9 +19,24 @@ #include #include +#include +#include namespace QuantLib { + namespace { + + // log N(x), with the asymptotic tail where N(x) underflows + Real logCumNormal(Real x) { + if (x > -30.0) + return std::log(CumulativeNormalDistribution()(x)); + Real z2 = 1.0/(x*x); + return -0.5*x*x - std::log(-x) - 0.5*std::log(M_TWOPI) + + std::log1p(-z2*(1.0 - z2*(3.0 - 15.0*z2))); + } + + } + AmericanPayoffAtExpiry::AmericanPayoffAtExpiry( Real spot, DiscountFactor discount, DiscountFactor dividendDiscount, Real variance, const ext::shared_ptr& payoff, @@ -170,10 +185,15 @@ namespace QuantLib { Y_ = 1.0; } else { X_ = 1.0; - if (cum_d2_ == 0.0) - Y_ = 0.0; // check needed on some extreme cases - else - Y_ = std::pow(Real(strike_/spot_), Real(2.0*mu_)); + Y_ = std::pow(Real(strike_/spot_), Real(2.0*mu_)); + if (variance_ >= QL_EPSILON && (cum_d2_ == 0.0 || !std::isfinite(Y_))) { + // at small variance the power can overflow where N(D2) underflows; + // the product is finite, so it is taken in logs + cum_d2_ = std::exp(2.0*mu_*log_H_S_ + logCumNormal(D2_)); + Y_ = 1.0; + } else if (cum_d2_ == 0.0) { + Y_ = 0.0; + } } if (!knock_in_) Y_ *= -1.0; diff --git a/test-suite/digitaloption.cpp b/test-suite/digitaloption.cpp index 9a6db45c6a9..2cb7cdf2408 100644 --- a/test-suite/digitaloption.cpp +++ b/test-suite/digitaloption.cpp @@ -435,6 +435,65 @@ BOOST_AUTO_TEST_CASE(testCashAtExpiryOrNothingAmericanValues) { } } +BOOST_AUTO_TEST_CASE(testCashAtExpiryOrNothingAmericanLowVolatility) { + + BOOST_TEST_MESSAGE("Testing American cash-(at-expiry)-or-nothing " + "digital option at low volatility..."); + + // expected values are the closed form evaluated in 50-digit arithmetic + DigitalOptionData values[] = { + // type, strike, spot, q, r, t, vol, value, tol + { Option::Call, 105.00, 100.00, 0.03, 0.08, 1.0, 0.001, 8.2035184454191468, 1e-8, true }, + { Option::Put, 95.00, 100.00, 0.03, -0.02, 1.0, 0.001, 1.0176365525852901, 1e-8, true }, + { Option::Call, 115.00, 100.00, 0.00, 0.03, 5.0, 0.003, 8.0820190229249992, 1e-8, true }, + { Option::Put, 85.00, 100.00, 0.03, 0.00, 5.0, 0.003, 0.32750723018931787, 1e-8, true } + }; + + DayCounter dc = Actual360(); + Date today = Date::todaysDate(); + + ext::shared_ptr spot(new SimpleQuote(0.0)); + ext::shared_ptr qRate(new SimpleQuote(0.0)); + ext::shared_ptr qTS = flatRate(today, qRate, dc); + ext::shared_ptr rRate(new SimpleQuote(0.0)); + ext::shared_ptr rTS = flatRate(today, rRate, dc); + ext::shared_ptr vol(new SimpleQuote(0.0)); + ext::shared_ptr volTS = flatVol(today, vol, dc); + + ext::shared_ptr stochProcess(new + BlackScholesMertonProcess(Handle(spot), + Handle(qTS), + Handle(rTS), + Handle(volTS))); + ext::shared_ptr engine( + new AnalyticDigitalAmericanEngine(stochProcess)); + + for (auto& value : values) { + + ext::shared_ptr payoff( + new CashOrNothingPayoff(value.type, value.strike, 10.0)); + Date exDate = today + timeToDays(value.t); + ext::shared_ptr amExercise(new AmericanExercise(today, + exDate, + true)); + + spot->setValue(value.s); + qRate->setValue(value.q); + rRate->setValue(value.r); + vol->setValue(value.v); + + VanillaOption opt(payoff, amExercise); + opt.setPricingEngine(engine); + + Real calculated = opt.NPV(); + Real error = std::fabs(calculated - value.result); + if (!(error <= value.tol)) { + REPORT_FAILURE("value", payoff, amExercise, value.s, value.q, value.r, today, value.v, + value.result, calculated, error, value.tol, value.knockin); + } + } +} + BOOST_AUTO_TEST_CASE(testAssetAtExpiryOrNothingAmericanValues) { BOOST_TEST_MESSAGE("Testing American asset-(at-expiry)-or-nothing "