From ef6e214690cae8a559909415a1f2fdff36cc1427 Mon Sep 17 00:00:00 2001 From: klaus spanderen Date: Mon, 23 Dec 2024 23:39:26 +0100 Subject: [PATCH 1/4] first light --- .../basket/operatorsplittingspreadengine.cpp | 44 ++++++++++---- test-suite/basketoption.cpp | 59 +++++++++++++++++++ 2 files changed, 91 insertions(+), 12 deletions(-) diff --git a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp index a7bfc9dea1..e73b304ba5 100644 --- a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp +++ b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp @@ -17,10 +17,10 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ -#include - #include #include +#include + #include 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 < 1000*QL_EPSILON) { + 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); } } diff --git a/test-suite/basketoption.cpp b/test-suite/basketoption.cpp index f7d1fb9b7b..44043f1f59 100644 --- a/test-suite/basketoption.cpp +++ b/test-suite/basketoption.cpp @@ -2512,6 +2512,65 @@ 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 zeroFlat(flatRate(today, Real(0), dc)); + + const auto processGen = [&](Real spot, Volatility vol) { + return ext::make_shared( + Handle(ext::make_shared(spot)), + zeroFlat, zeroFlat, + Handle(flatVol(today, vol, dc)) + ); + }; + + BasketOption basketOption( + ext::make_shared( + ext::make_shared(Option::Put, Real(50)) + ), + ext::make_shared(maturity) + ); + + const Real eps = 1e-5; + + + const auto engineGen = [&](Volatility vol) { + return ext::make_shared( + 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 = 1e-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); + + + //std::cout << std::setprecision(16) << 0.5*(lNpv + rNpv) << " " << calculated - 0.5*(lNpv + rNpv) << std::endl; +} BOOST_AUTO_TEST_SUITE_END() From a5b388ff4d7161b94861a3ca567a1e7072839a87 Mon Sep 17 00:00:00 2001 From: klaus spanderen Date: Mon, 30 Dec 2024 11:38:56 +0100 Subject: [PATCH 2/4] formatting --- .../basket/operatorsplittingspreadengine.cpp | 40 +++++++++---------- test-suite/basketoption.cpp | 5 +-- 2 files changed, 21 insertions(+), 24 deletions(-) diff --git a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp index e73b304ba5..3c449492ec 100644 --- a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp +++ b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp @@ -60,9 +60,9 @@ namespace QuantLib { 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 + 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)); + *(rs*d1*d2 + (1-rho_*rho_)*variance1)); if (order_ == First) return callPutParityPrice(kirkCallNPV + 0.5*oPlt); @@ -92,24 +92,24 @@ namespace QuantLib { const Real vol22 = vol2*vol2; const Real vol23 = vol22*vol2; - if (rs < 1000*QL_EPSILON) { - 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); - } + 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; diff --git a/test-suite/basketoption.cpp b/test-suite/basketoption.cpp index 44043f1f59..a755ec9842 100644 --- a/test-suite/basketoption.cpp +++ b/test-suite/basketoption.cpp @@ -2559,7 +2559,7 @@ BOOST_AUTO_TEST_CASE(testNoDivByZeroOperatorSplitting) { const Real calculated = basketOption.NPV(); const Real diff =std::abs(calculated - expected); - const Real tol = 1e-8; + 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) @@ -2567,9 +2567,6 @@ BOOST_AUTO_TEST_CASE(testNoDivByZeroOperatorSplitting) { << "\n expected: " << expected << "\n diff: " << diff << "\n tolerance: " << tol); - - - //std::cout << std::setprecision(16) << 0.5*(lNpv + rNpv) << " " << calculated - 0.5*(lNpv + rNpv) << std::endl; } From cf28d744943ffa49ddbe21022dce76f31364715c Mon Sep 17 00:00:00 2001 From: klaus spanderen Date: Mon, 30 Dec 2024 11:39:48 +0100 Subject: [PATCH 3/4] formatting --- ql/pricingengines/basket/operatorsplittingspreadengine.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp index 3c449492ec..5270f63c70 100644 --- a/ql/pricingengines/basket/operatorsplittingspreadengine.cpp +++ b/ql/pricingengines/basket/operatorsplittingspreadengine.cpp @@ -103,7 +103,7 @@ namespace QuantLib { 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)* + (-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)); From c76169b6e7f1e411593f471c3756a5efb2179f1b Mon Sep 17 00:00:00 2001 From: klaus spanderen Date: Mon, 30 Dec 2024 13:22:30 +0100 Subject: [PATCH 4/4] formatting --- test-suite/basketoption.cpp | 48 ++++++++++++++++++------------------- 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/test-suite/basketoption.cpp b/test-suite/basketoption.cpp index a755ec9842..7c30586679 100644 --- a/test-suite/basketoption.cpp +++ b/test-suite/basketoption.cpp @@ -2531,36 +2531,36 @@ BOOST_AUTO_TEST_CASE(testNoDivByZeroOperatorSplitting) { BasketOption basketOption( ext::make_shared( - ext::make_shared(Option::Put, Real(50)) - ), - ext::make_shared(maturity) + ext::make_shared(Option::Put, Real(50)) + ), + ext::make_shared(maturity) ); - const Real eps = 1e-5; + const Real eps = 1e-5; - + const auto engineGen = [&](Volatility vol) { - return ext::make_shared( - processGen(160, 0.25*3), - processGen(100, vol*150.0/100.0), - 1/3.0, OperatorSplittingSpreadEngine::Second - ); + return ext::make_shared( + 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) + 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