Keep the RS implied-standard-deviation approximation finite

ql/pricingengines/blackformula.cpp

@@ -290,18 +290,23 @@ namespace QuantLib {
         const Real b = std::exp(M_2_PI*y);
         const Real B = 4.0*(b + 1/b)
             - 2*K/F*(a + 1.0/a)*(ey2 + 1 - R2);
-        const Real C = (R2-squared(ey-1))*(squared(ey+1)-R2)/ey2;
+        const Real C = std::max(0.0, (R2-squared(ey-1))*(squared(ey+1)-R2)/ey2);
 
         const Real beta = 2*C/(B+std::sqrt(B*B+4*A*C));
         const Real gamma = -M_PI_2*std::log(beta);
 
+        const auto sqrtDiff = [](Real g, Real yy) {
+            return std::isfinite(g) ? Real(std::sqrt(g + yy) - std::sqrt(g - yy))
+                                    : Real(0.0);
+        };
+
         if (y >= 0.0) {
             const Real M0 = K*df*(
                 (type == Option::Call) ? Real(ey*Af(std::sqrt(2*y)) - 0.5)
                                        : 0.5-ey*Af(-std::sqrt(2*y)));
 
             if (marketValue <= M0)
-                return std::sqrt(gamma+y)-std::sqrt(gamma-y);
+                return sqrtDiff(gamma, y);
             else
                 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
         }
@@ -311,7 +316,7 @@ namespace QuantLib {
                                        : Af(std::sqrt(-2*y)) - 0.5*ey);
 
             if (marketValue <= M0)
-                return std::sqrt(gamma-y)-std::sqrt(gamma+y);
+                return sqrtDiff(gamma, -y);
             else
                 return std::sqrt(gamma+y)+std::sqrt(gamma-y);
         }

test-suite/blackformula.cpp

@@ -25,6 +25,7 @@
 #include "utilities.hpp"
 #include <ql/pricingengines/blackformula.hpp>
 #include <cmath>
+#include <iomanip>
 
 using namespace QuantLib;
 using namespace boost::unit_test_framework;
@@ -556,6 +557,43 @@ BOOST_AUTO_TEST_CASE(testDeepInTheMoneyPricesNotBelowIntrinsic) {
     }
 }
 
+
+BOOST_AUTO_TEST_CASE(testImpliedStdDevApproximationRSStaysFinite) {
+    BOOST_TEST_MESSAGE(
+        "Testing RS implied-standard-deviation approximation away from the money...");
+
+    const Real forward = 1.0;
+    const Real discount = 1.0;
+
+    for (Real ratio: {1.01, 1.1, 1.5, 2.0, 3.0, 5.0, 10.0}) {
+        const Real strike = forward / ratio;
+        for (Real stdDev: {0.01, 0.05, 0.1, 0.2, 0.4, 0.8, 1.5}) {
+            for (auto type: {Option::Call, Option::Put}) {
+                const Real price = blackFormula(type, strike, forward, stdDev, discount);
+                const Real approx = blackFormulaImpliedStdDevApproximationRS(
+                    type, strike, forward, price, discount);
+
+                if (!std::isfinite(approx))
+                    BOOST_ERROR("RS approximation is not finite"
+                                << std::setprecision(16)
+                                << "\n    forward:    " << forward
+                                << "\n    strike:     " << strike
+                                << "\n    std dev:    " << stdDev
+                                << "\n    price:      " << price
+                                << "\n    calculated: " << approx);
+
+                if (approx < 0.0)
+                    BOOST_ERROR("RS approximation is negative"
+                                << std::setprecision(16)
+                                << "\n    forward:    " << forward
+                                << "\n    strike:     " << strike
+                                << "\n    std dev:    " << stdDev
+                                << "\n    calculated: " << approx);
+            }
+        }
+    }
+}
+
 BOOST_AUTO_TEST_SUITE_END()
 
 BOOST_AUTO_TEST_SUITE_END()