diff --git a/ql/experimental/finitedifferences/fdextoujumpvanillaengine.cpp b/ql/experimental/finitedifferences/fdextoujumpvanillaengine.cpp index 186e77a6940..7200d212c11 100644 --- a/ql/experimental/finitedifferences/fdextoujumpvanillaengine.cpp +++ b/ql/experimental/finitedifferences/fdextoujumpvanillaengine.cpp @@ -43,9 +43,11 @@ namespace QuantLib { const ext::shared_ptr& process, const ext::shared_ptr& rTS, Size tGrid, Size xGrid, Size yGrid, + const ext::shared_ptr& shape, const FdmSchemeDesc& schemeDesc) : process_(process), rTS_(rTS), + shape_(shape), tGrid_(tGrid), xGrid_(xGrid), yGrid_(yGrid), @@ -58,9 +60,9 @@ namespace QuantLib { = rTS_->dayCounter().yearFraction(rTS_->referenceDate(), arguments_.exercise->lastDate()); const ext::shared_ptr ouProcess( - process_->getExtendedOrnsteinUhlenbeckProcess()); + process_->getExtendedOrnsteinUhlenbeckProcess()); const ext::shared_ptr xMesher( - new FdmSimpleProcess1dMesher(xGrid_, ouProcess,maturity)); + new FdmSimpleProcess1dMesher(xGrid_, ouProcess,maturity)); const ext::shared_ptr yMesher( new ExponentialJump1dMesher(yGrid_, @@ -73,7 +75,7 @@ namespace QuantLib { // 2. Calculator const ext::shared_ptr calculator( - new FdmExtOUJumpModelInnerValue(arguments_.payoff, mesher)); + new FdmExtOUJumpModelInnerValue(arguments_.payoff, mesher, shape_)); // 3. Step conditions const ext::shared_ptr conditions = diff --git a/ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp b/ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp index 3d29f302b03..3da791aa7d4 100644 --- a/ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp +++ b/ql/experimental/finitedifferences/fdextoujumpvanillaengine.hpp @@ -28,6 +28,7 @@ #include #include #include +#include namespace QuantLib { @@ -38,17 +39,20 @@ namespace QuantLib { : public GenericEngine { public: - FdExtOUJumpVanillaEngine( - const ext::shared_ptr& p, - const ext::shared_ptr& rTS, - Size tGrid = 50, Size xGrid = 200, Size yGrid = 50, - const FdmSchemeDesc& schemeDesc=FdmSchemeDesc::Hundsdorfer()); + typedef FdmExtOUJumpModelInnerValue::Shape Shape; + FdExtOUJumpVanillaEngine( + const ext::shared_ptr& p, + const ext::shared_ptr& rTS, + Size tGrid = 50, Size xGrid = 200, Size yGrid = 50, + const ext::shared_ptr& shape = ext::shared_ptr(), + const FdmSchemeDesc& schemeDesc=FdmSchemeDesc::Hundsdorfer()); void calculate() const; private: const ext::shared_ptr process_; const ext::shared_ptr rTS_; + const ext::shared_ptr shape_; const Size tGrid_, xGrid_, yGrid_; const FdmSchemeDesc schemeDesc_; }; diff --git a/ql/methods/finitedifferences/meshers/exponentialjump1dmesher.cpp b/ql/methods/finitedifferences/meshers/exponentialjump1dmesher.cpp index 945f98f1e6a..ba3f3415f7c 100644 --- a/ql/methods/finitedifferences/meshers/exponentialjump1dmesher.cpp +++ b/ql/methods/finitedifferences/meshers/exponentialjump1dmesher.cpp @@ -38,13 +38,15 @@ namespace QuantLib { QL_REQUIRE(steps > 1, "minimum number of steps is two"); const Real start = 0.0; - const Real end = 1.0-eps; + const Real end = 1.0-eps; const Real dx = (end-start)/(steps-1); - + const Real scale = 1/(1-std::exp(-beta/jumpIntensity)); + for (Size i=0; i < steps; ++i) { const Real p = start + i*dx; - locations_[i] = -1.0/eta*std::log(1.0-p); + locations_[i] = scale*(-1.0/eta*std::log(1.0-p)); } + for (Size i=0; i < steps-1; ++i) { dminus_[i+1] = dplus_[i] = locations_[i+1]-locations_[i]; } diff --git a/test-suite/swingoption.cpp b/test-suite/swingoption.cpp index 3727bbd9042..0a5458b2ed8 100644 --- a/test-suite/swingoption.cpp +++ b/test-suite/swingoption.cpp @@ -20,7 +20,10 @@ #include "swingoption.hpp" #include "utilities.hpp" +#include #include +#include +#include #include #include #include @@ -34,6 +37,7 @@ #include #include #include +#include #include #include #include @@ -60,8 +64,7 @@ namespace { new ExtendedOrnsteinUhlenbeckProcess(speed, volatility, x0[0], constant(x0[0]))); return ext::make_shared( - ouProcess, x0[1], beta, - jumpIntensity, eta); + ouProcess, x0[1], beta, jumpIntensity, eta); } } @@ -434,6 +437,162 @@ void SwingOptionTest::testExtOUJumpSwingOption() { } } +namespace { + class SwingPdePricing { + public: + typedef FdSimpleExtOUJumpSwingEngine::Shape Shape; + + SwingPdePricing(const ext::shared_ptr& process, + const ext::shared_ptr& option, + const ext::shared_ptr& shape) + : process_(process), option_(option), shape_(shape) {} + + Real operator()(Real x) const { + const ext::shared_ptr rTS( + flatRate(0.0, Actual365Fixed())); + + const Size gridX = 200; + const Size gridY = 100; + const Size gridT = 100; + + option_->setPricingEngine( + ext::make_shared( + process_, rTS, + Size(gridT)/x, Size(gridX/x), Size(gridY/x), shape_)); + + return option_->NPV(); + } + + private: + const ext::shared_ptr process_; + const ext::shared_ptr option_; + const ext::shared_ptr shape_; + }; +} + +void SwingOptionTest::testKlugeChFVanillaPricing() { + BOOST_TEST_MESSAGE("Testing Kluge PDE Vanilla Pricing in" + " comparison to moment matching..."); + + SavedSettings backup; + + Date settlementDate = Date(22, November, 2019); + Settings::instance().evaluationDate() = settlementDate; + DayCounter dayCounter = Actual365Fixed(); + Date maturityDate = settlementDate + Period(6, Months); + const Time t = dayCounter.yearFraction(settlementDate, maturityDate); + + const Real f0 = 30; + + const Real x0 = 0.0; + const Real y0 = 0.0; + + const Real beta = 5.0; + const Real eta = 5.0; + const Real lambda = 4.0; + const Real alpha = 4.0; + const Real sig = 1.0; + + const ext::shared_ptr klugeProcess = + ext::make_shared( + ext::make_shared( + alpha, sig, x0, constant(0.0)), + y0, beta, lambda, eta); + + const Real strike = f0; + + const ext::shared_ptr option = + ext::make_shared( + ext::make_shared(Option::Call, strike), + ext::make_shared(maturityDate)); + + typedef FdSimpleExtOUJumpSwingEngine::Shape Shape; + const ext::shared_ptr shape(ext::make_shared()); + + const Real ps = std::log(f0) + - sig*sig/(4*alpha)*(1-std::exp(-2*alpha*t)) + - lambda/beta*std::log((eta-std::exp(-beta*t))/(eta-1.0)); + + shape->push_back(Shape::value_type(t, ps)); + + const Real expected = + RichardsonExtrapolation( + SwingPdePricing(klugeProcess, option, shape), 4.0)(2.0, 1.5); + + const Real stdDev = std::sqrt((((2 - 2*std::exp(-2*beta*t))*lambda) + /(beta*eta*eta) + ((1 - std::exp(-2*alpha*t))*sig*sig)/alpha)/2.); + + const Real bsNPV = blackFormula(Option::Call, strike, f0, stdDev); + + const Real g1 = ((2 - 2*std::exp(-3*beta*t))*lambda)/(beta*eta*eta*eta) + / (stdDev*stdDev*stdDev); + + const Real g2 = 3*(std::exp((alpha + beta)*t) + * square()(2*alpha*std::exp(2*alpha*t)*(-1 + std::exp(2*beta*t)) + *lambda + beta*std::exp(2*beta*t)*(-1 + std::exp(2*alpha*t)) + *eta*eta*sig*sig) + + 16*alpha*alpha*beta*std::exp((5*alpha + 3*beta)*t)*lambda + *std::sinh(2*beta*t)) + / (4.*alpha*alpha*beta*beta + *std::exp(5*(alpha + beta)*t)*eta*eta*eta*eta) + / (stdDev*stdDev*stdDev*stdDev) - 3.0; + + const Real d = (std::log(f0/strike) + 0.5*stdDev*stdDev)/stdDev; + + // Jurczenko E., Maillet B. and Negrea B., + // Multi-Moment Approximate Option Pricing Models: + // A General Comparison (Part 1) + // https://papers.ssrn.com/sol3/papers.cfm?abstract_id=300922 + const NormalDistribution n; + const Real q3 = 1/Factorial::get(3)*f0*stdDev*(2*stdDev - d)*n(d); + const Real q4 = 1/Factorial::get(4)*f0*stdDev*(d*d - 3*d*stdDev - 1)*n(d); + const Real q5 = 10/Factorial::get(6)*f0*stdDev*( + d*d*d*d - 5*d*d*d*stdDev - 6*d*d + 15*d*stdDev + 3)*n(d); + + // Corrado C. and T. Su, (1996-b), + // “Skewness and Kurtosis in S&P 500 IndexReturns Implied by Option Prices”, + // Journal of Financial Research 19 (2), 175-192. + const Real ccs3 = bsNPV + g1*q3; + const Real ccs4 = ccs3 + g2*q4; + + // Rubinstein M., (1998), “Edgeworth Binomial Trees”, + // Journal of Derivatives 5 (3), 20-27. + const Real cr = ccs4 + g1*g1*q5; + + const Volatility expectedImplVol = blackFormulaImpliedStdDevLiRS( + Option::Call, strike, f0, expected, 1.0)/std::sqrt(t); + + const Volatility bsImplVol = blackFormulaImpliedStdDevLiRS( + Option::Call, strike, f0, bsNPV, 1.0)/std::sqrt(t); + + const Volatility ccs3ImplVol = blackFormulaImpliedStdDevLiRS( + Option::Call, strike, f0, ccs3, 1.0)/std::sqrt(t); + + const Volatility ccs4ImplVol = blackFormulaImpliedStdDevLiRS( + Option::Call, strike, f0, ccs4, 1.0)/std::sqrt(t); + + const Volatility crImplVol = blackFormulaImpliedStdDevLiRS( + Option::Call, strike, f0, cr, 1.0)/std::sqrt(t); + + const Real tol[] = {0.01, 0.0075, 0.005, 0.004}; + const std::string methods[] = { + "Second Order", "Third Order", "Fourth Order", "Rubinstein"}; + + const Real calculated[] = {bsImplVol, ccs3ImplVol, ccs4ImplVol, crImplVol}; + + for (Size i=0; i < 4; ++i) { + const Real diff = std::fabs(calculated[i] - expectedImplVol); + if (diff > tol[i]) { + BOOST_ERROR("failed to reproduce vanilla option implied volatility " + "with moment matching" + << "\n calculated: " << calculated[i] + << "\n expected: " << expectedImplVol + << "\n difference: " << diff + << "\n tolerance: " << tol[i] + << "\n method: " << methods[i]); + } + } +} test_suite* SwingOptionTest::suite(SpeedLevel speed) { test_suite* suite = BOOST_TEST_SUITE("Swing-Option Test"); @@ -443,6 +602,8 @@ test_suite* SwingOptionTest::suite(SpeedLevel speed) { suite->add(QUANTLIB_TEST_CASE(&SwingOptionTest::testFdBSSwingOption)); suite->add(QUANTLIB_TEST_CASE( &SwingOptionTest::testFdmExponentialJump1dMesher)); + suite->add(QUANTLIB_TEST_CASE( + &SwingOptionTest::testKlugeChFVanillaPricing)); if (speed <= Fast) { suite->add(QUANTLIB_TEST_CASE( diff --git a/test-suite/swingoption.hpp b/test-suite/swingoption.hpp index fd63378270a..59e96a62bb1 100644 --- a/test-suite/swingoption.hpp +++ b/test-suite/swingoption.hpp @@ -33,6 +33,7 @@ class SwingOptionTest { static void testExtOUJumpSwingOption(); static void testFdmExponentialJump1dMesher(); static void testExtOUJumpVanillaEngine(); + static void testKlugeChFVanillaPricing(); static boost::unit_test_framework::test_suite* suite(SpeedLevel); };