deal with removable singularity

ql/pricingengines/basket/operatorsplittingspreadengine.cpp

@@ -17,10 +17,10 @@
  FOR A PARTICULAR PURPOSE.  See the license for more details.
 */
 
-#include <ql/pricingengines/basket/operatorsplittingspreadengine.hpp>
-
 #include <ql/math/functional.hpp>
 #include <ql/math/distributions/normaldistribution.hpp>
+#include <ql/pricingengines/basket/operatorsplittingspreadengine.hpp>
+
 #include <utility>
 
 namespace QuantLib {
@@ -57,17 +57,18 @@ namespace QuantLib {
 
         const Real kirkCallNPV = df*(f1*N(d1) - (f2 + k)*N(d2));
 
-        const Real vv = (rho_*vol1 - sig2)*vol2/(sig_m*sig_m);
-        const Real oPlt = -sig2*sig2 * k * df * NormalDistribution()(d2) * vv
-                *( d2*(1 - rho_*vol1/sig2)
-                   - 0.5*sig_m * vv * k / (f2+k)
-                     * ( d1*d2 + (1-rho_*rho_)*squared(vol1/(rho_*vol1-sig2))));
+        const Real vs = vol2/(sig_m*sig_m);
+        const Real rs = squared(rho_*vol1 - sig2);
+
+        const Real oPlt = -sig2*sig2*k*df*NormalDistribution()(d2)*vs
+            *( -d2*rs/sig2 - 0.5*vs*sig_m * k / (f2+k)
+                *(rs*d1*d2  + (1-rho_*rho_)*variance1));
 
         if (order_ == First)
             return callPutParityPrice(kirkCallNPV + 0.5*oPlt);
 
         QL_REQUIRE(order_ == Second, "unknown approximation type");
-
+        
         /*
            In the original paper, the second-order approximation was computed using numerical differentiation.
            The following Mathematica scripts calculates the approximation to the n'th order.
@@ -87,6 +88,29 @@ namespace QuantLib {
 
         const Real R2 = f2+k;
         const Real R1 = f1/R2;
+        const Real vol12 = vol1*vol1;
+        const Real vol22 = vol2*vol2;
+        const Real vol23 = vol22*vol2;
+
+        if (rs < std::pow(QL_EPSILON, 0.625)) {
+            const Real vol24 = vol22*vol22;
+            const Real vol26 = vol22*vol24;
+            const Real k2 = k*k;
+            const Real R22 = R2*R2;
+            const Real R24 = R22*R22;
+            const Real kmR22 = squared(k-R2);
+            const Real lnR1 = std::log(R1);
+
+            const Real ooPlt = -0.0625*(k2*kmR22*vol26*(-8*R22*R24*(7*k2 - 7*k*R2 + R22)*vol12*vol12 
+                + kmR22*R24*vol12*(-112*k*R2 + 16*R22 + k2*(124 + 3*vol12))*vol22 - 2*kmR22*kmR22*R22*
+                (-28*k*R2 + 4*R22 + k2*(34 + 3*vol12))*vol24 + 3*k2*kmR22*kmR22*kmR22*vol26 - 4*k*(k -R2)*
+                R24*lnR1*(-4*R22*vol12 + 4*kmR22*vol22 + 3*k*(k - R2)*vol22*lnR1)))
+                /(std::exp(squared(-(R22*vol12) + kmR22*vol22 +    2*R22*lnR1)/(8*R24*vol12 - 8*kmR22*R22*vol22))*
+                M_SQRTPI*M_SQRT2*R22*R24*R2*squared(R22*vol12 - kmR22*vol22)*std::sqrt(vol12 - (kmR22*vol22)/R22));
+
+            return callPutParityPrice(kirkCallNPV + 0.5*oPlt + 0.125*ooPlt);        
+        } 
+
         const Real F2 = f2;
 
         const Real F22 = F2*F2;
@@ -97,9 +121,6 @@ namespace QuantLib {
         const Real iR22 = iR2*iR2;
         const Real iR23 = iR22*iR2;
         const Real iR24 = iR22*iR22;
-        const Real vol12 = vol1*vol1;
-        const Real vol22 = vol2*vol2;
-        const Real vol23 = vol22*vol2;
         const Real a = vol12 - 2*F2*iR2*rho_*vol1*vol2 + F22*iR22*vol22;
         const Real a2 = a*a;
         const Real b = a/2+std::log(R1);
@@ -158,7 +179,6 @@ namespace QuantLib {
                 4*f*F23*g*iR22*m*vol2 + f*F24*vol2*(-4*g*iR23*m + 2*g*iR22*x +
                 iR22*m*z))))))/(16.*a2*a2*c*e2*F23*M_SQRT2*M_SQRTPI*vol2);
 
-
         return callPutParityPrice(kirkCallNPV + 0.5*oPlt + 0.125*ooPlt);
     }
 }

test-suite/basketoption.cpp

@@ -2512,6 +2512,62 @@ BOOST_AUTO_TEST_CASE(testAccurateAmericanBasketOptions) {
                << "\n    tolerance:   " << tol);
 }
 
+BOOST_AUTO_TEST_CASE(testNoDivByZeroOperatorSplitting) {
+    BOOST_TEST_MESSAGE("Testing division by zero issue for the Operator Splitting engine...");
+
+    const DayCounter dc = Actual365Fixed();
+    const Date today = Date(5, December, 2024);
+    const Date maturity = today + Period(18, Months);
+
+    const Handle<YieldTermStructure> zeroFlat(flatRate(today, Real(0), dc)); 
+
+    const auto processGen = [&](Real spot, Volatility vol) {
+        return ext::make_shared<GeneralizedBlackScholesProcess>(
+            Handle<Quote>(ext::make_shared<SimpleQuote>(spot)),
+            zeroFlat, zeroFlat,
+            Handle<BlackVolTermStructure>(flatVol(today, vol, dc))
+        );
+    };
+
+    BasketOption basketOption(
+        ext::make_shared<SpreadBasketPayoff>(
+            ext::make_shared<PlainVanillaPayoff>(Option::Put, Real(50))
+        ),
+        ext::make_shared<EuropeanExercise>(maturity)
+    );
+
+    const Real eps = 1e-5;
+
+
+    const auto engineGen  = [&](Volatility vol) {
+        return ext::make_shared<OperatorSplittingSpreadEngine>(
+            processGen(160, 0.25*3),
+            processGen(100, vol*150.0/100.0),
+            1/3.0, OperatorSplittingSpreadEngine::Second
+        );
+    };
+
+    basketOption.setPricingEngine(engineGen(0.25 - eps));
+    const Real lNpv = basketOption.NPV();
+
+    basketOption.setPricingEngine(engineGen(0.25 + eps));
+    const Real rNpv = basketOption.NPV();
+
+    const Real expected = 0.5*(rNpv + lNpv);
+
+    basketOption.setPricingEngine(engineGen(0.25));
+    const Real calculated = basketOption.NPV();
+
+    const Real diff =std::abs(calculated - expected);
+    const Real tol = 5e-8;
+    if (std::isnan(calculated) || diff > tol)
+        BOOST_FAIL("failed to reproduce spread option price with OperatorSplittingEngine"
+           << std::fixed << std::setprecision(8)
+           << "\n    calculated:  " << calculated
+           << "\n    expected:    " << expected
+           << "\n    diff:        " << diff
+           << "\n    tolerance:   " << tol);
+}
 
 
 BOOST_AUTO_TEST_SUITE_END()