Keep the digital American at-expiry price accurate at low volatility

ql/pricingengines/americanpayoffatexpiry.cpp

@@ -19,9 +19,24 @@
 
 #include <ql/pricingengines/americanpayoffatexpiry.hpp>
 #include <ql/math/distributions/normaldistribution.hpp>
+#include <ql/mathconstants.hpp>
+#include <cmath>
 
 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<StrikedTypePayoff>& 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; 

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<SimpleQuote> spot(new SimpleQuote(0.0));
+    ext::shared_ptr<SimpleQuote> qRate(new SimpleQuote(0.0));
+    ext::shared_ptr<YieldTermStructure> qTS = flatRate(today, qRate, dc);
+    ext::shared_ptr<SimpleQuote> rRate(new SimpleQuote(0.0));
+    ext::shared_ptr<YieldTermStructure> rTS = flatRate(today, rRate, dc);
+    ext::shared_ptr<SimpleQuote> vol(new SimpleQuote(0.0));
+    ext::shared_ptr<BlackVolTermStructure> volTS = flatVol(today, vol, dc);
+
+    ext::shared_ptr<BlackScholesMertonProcess> stochProcess(new
+        BlackScholesMertonProcess(Handle<Quote>(spot),
+                                  Handle<YieldTermStructure>(qTS),
+                                  Handle<YieldTermStructure>(rTS),
+                                  Handle<BlackVolTermStructure>(volTS)));
+    ext::shared_ptr<PricingEngine> engine(
+                             new AnalyticDigitalAmericanEngine(stochProcess));
+
+    for (auto& value : values) {
+
+        ext::shared_ptr<StrikedTypePayoff> payoff(
+            new CashOrNothingPayoff(value.type, value.strike, 10.0));
+        Date exDate = today + timeToDays(value.t);
+        ext::shared_ptr<Exercise> 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 "