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 "