diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj index 3ac7544e02c..81efc70968c 100644 --- a/QuantLib.vcxproj +++ b/QuantLib.vcxproj @@ -1011,6 +1011,7 @@ + @@ -1019,6 +1020,7 @@ + @@ -2260,6 +2262,7 @@ + diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters index cba8a6501da..2a38c4cbf63 100644 --- a/QuantLib.vcxproj.filters +++ b/QuantLib.vcxproj.filters @@ -987,6 +987,9 @@ math + + math + math @@ -1164,6 +1167,9 @@ math\integrals + + math\integrals + math\integrals @@ -4850,6 +4856,9 @@ math + + math + math diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt index da79a0b8790..9e578e0d265 100644 --- a/ql/CMakeLists.txt +++ b/ql/CMakeLists.txt @@ -364,6 +364,7 @@ set(QL_SOURCES math/distributions/normaldistribution.cpp math/distributions/studenttdistribution.cpp math/errorfunction.cpp + math/expm1.cpp math/factorial.cpp math/incompletegamma.cpp math/integrals/discreteintegrals.cpp @@ -1426,6 +1427,7 @@ set(QL_HEADERS math/distributions/poissondistribution.hpp math/distributions/studenttdistribution.hpp math/errorfunction.hpp + math/expm1.hpp math/factorial.hpp math/fastfouriertransform.hpp math/functional.hpp @@ -1433,6 +1435,7 @@ set(QL_HEADERS math/incompletegamma.hpp math/integrals/discreteintegrals.hpp math/integrals/exponentialintegrals.hpp + math/integrals/expsinhintegral.hpp math/integrals/filonintegral.hpp math/integrals/gaussianorthogonalpolynomial.hpp math/integrals/gaussianquadratures.hpp diff --git a/ql/math/Makefile.am b/ql/math/Makefile.am index f7c8778f873..535de3132a5 100644 --- a/ql/math/Makefile.am +++ b/ql/math/Makefile.am @@ -17,6 +17,7 @@ this_include_HEADERS = \ comparison.hpp \ curve.hpp \ errorfunction.hpp \ + expm1.hpp \ factorial.hpp \ fastfouriertransform.hpp \ functional.hpp \ @@ -44,6 +45,7 @@ cpp_files = \ beta.cpp \ bspline.cpp \ errorfunction.cpp \ + expm1.cpp \ factorial.cpp \ incompletegamma.cpp \ matrix.cpp \ diff --git a/ql/math/all.hpp b/ql/math/all.hpp index e724b24b3ce..7df092f2b55 100644 --- a/ql/math/all.hpp +++ b/ql/math/all.hpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include diff --git a/ql/math/expm1.cpp b/ql/math/expm1.cpp new file mode 100644 index 00000000000..29f595bfc85 --- /dev/null +++ b/ql/math/expm1.cpp @@ -0,0 +1,54 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2023 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +#include +#include + +#include + +namespace QuantLib { + std::complex expm1(const std::complex& z) { + if (std::abs(z) < 1.0) { + const Real a = z.real(), b = z.imag(); + const Real exp_1 = std::expm1(a); + const Real cos_1 = -2*squared(std::sin(0.5*b)); + + return std::complex( + exp_1*cos_1 + exp_1 + cos_1, + std::sin(b)*std::exp(a) + ); + } + else { + return std::exp(z)-1.0; + } + } + + std::complex log1p(const std::complex& z) { + const Real a = z.real(), b = z.imag(); + if (std::abs(a) < 0.5 && std::abs(b) < 0.5) { + return std::complex( + 0.5*std::log1p(a*a + 2*a + b*b), + std::arg(1.0 + z) + ); + } + else { + return std::log(1.0+z); + } + } +} diff --git a/ql/math/expm1.hpp b/ql/math/expm1.hpp new file mode 100644 index 00000000000..c143a6f70a2 --- /dev/null +++ b/ql/math/expm1.hpp @@ -0,0 +1,35 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2023 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + + /*! \file expm1.hpp + \brief complex versions of expm1 and logp1 +*/ + +#ifndef quantlib_expm1_hpp +#define quantlib_expm1_hpp + +#include +#include + +namespace QuantLib { + std::complex expm1(const std::complex& z); + std::complex log1p(const std::complex& z); +} +#endif + diff --git a/ql/math/integrals/Makefile.am b/ql/math/integrals/Makefile.am index fcbcf543f5e..59bfbfe6fc6 100644 --- a/ql/math/integrals/Makefile.am +++ b/ql/math/integrals/Makefile.am @@ -5,6 +5,7 @@ this_include_HEADERS = \ all.hpp \ discreteintegrals.hpp \ exponentialintegrals.hpp \ + expsinhintegral.hpp \ filonintegral.hpp \ gausslaguerrecosinepolynomial.hpp \ gausslobattointegral.hpp \ diff --git a/ql/math/integrals/all.hpp b/ql/math/integrals/all.hpp index 5ee9cedebd1..ee03ff1a828 100644 --- a/ql/math/integrals/all.hpp +++ b/ql/math/integrals/all.hpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include diff --git a/ql/math/integrals/exponentialintegrals.cpp b/ql/math/integrals/exponentialintegrals.cpp index 8140542315b..b87d3836367 100644 --- a/ql/math/integrals/exponentialintegrals.cpp +++ b/ql/math/integrals/exponentialintegrals.cpp @@ -23,14 +23,12 @@ #include #include +#include #include +#include #include -#ifndef M_EULER_MASCHERONI - #define M_EULER_MASCHERONI 0.5772156649015328606065121 -#endif - namespace QuantLib { namespace exponential_integrals_helper { @@ -122,28 +120,61 @@ namespace QuantLib { } } + std::complex _Ei( + const std::complex& z, const std::complex& acc) { - std::complex E1(std::complex z) { - QL_REQUIRE(std::abs(z) <= 25.0, "Insufficient precision for |z| > 25.0"); + if (z.real() == 0.0 && z.imag() == 0.0 + && std::numeric_limits::has_infinity) { + return std::complex(-std::numeric_limits::infinity()); + } - std::complex s(0.0), sn(-z); - Size n; - for (n=2; n < 1000 && s + sn/Real(n-1) != s; ++n) { - s+=sn/Real(n-1); - sn *= -z/Real(n); - } + constexpr Real DIST = 4.5; + constexpr Real MAX_ERROR = 5*QL_EPSILON; - QL_REQUIRE(n < 1000, "series conversion issue"); + const Real z_inf = std::log(0.01*QL_MAX_REAL) + std::log(100.0); + QL_REQUIRE(z.real() < z_inf, "argument error " << z); - return -M_EULER_MASCHERONI - std::log(z) -s; - } + const Real z_asym = 2.0 - 1.035*std::log(MAX_ERROR); - std::complex Ei(std::complex z) { - QL_REQUIRE(std::abs(z) <= 25.0, "Insufficient precision for |z| > 25.0"); + using boost::math::sign; + const Real abs_z = std::abs(z); - std::complex s(0.0), sn=z; + const auto match = [=]( + const std::complex& z1, const std::complex& z2) + -> bool { + const std::complex d = z1 - z2; + return std::abs(d.real()) <= MAX_ERROR*std::abs(z1.real()) + && std::abs(d.imag()) <= MAX_ERROR*std::abs(z1.imag()); + }; + + if (z.real() > z_inf) + return std::complex(std::exp(z)/z) + acc; + if (abs_z > 1.1*z_asym) { + std::complex ei = acc + std::complex(Real(0.0), sign(z.imag())*M_PI); + std::complex s(std::exp(z)/z); + for (Size i=1; i <= std::floor(abs_z)+1; ++i) { + if (match(ei+s, ei)) + return ei+s; + ei += s; + s *= Real(i)/z; + } + QL_FAIL("series conversion issue for Ei(" << z << ")"); + } + + if (abs_z > DIST && (z.real() < 0 || std::abs(z.imag()) > DIST)) { + std::complex ei(0.0); + for (Size k = 47; k >=1; --k) { + ei = - Real(k*k)/(2.0*k + 1.0 - z + ei); + } + return (acc + std::complex(0.0, sign(z.imag())*M_PI)) + - std::exp(z)/ (1.0 - z + ei); + + QL_FAIL("series conversion issue for Ei(" << z << ")"); + } + + std::complex s(0.0), sn=z; Real nn=1.0; Size n; @@ -156,27 +187,65 @@ namespace QuantLib { sn *= -z / Real(2*n); } - QL_REQUIRE(n < 1000, "series conversion issue"); + QL_REQUIRE(n < 1000, "series conversion issue for Ei(" << z << ")"); + + const std::complex r + = (M_EULER_MASCHERONI + acc) + std::log(z) + std::exp(0.5*z)*s; - return M_EULER_MASCHERONI + std::log(z) + std::exp(0.5*z)*s; + if (z.imag() != Real(0.0)) + return r; + else + return std::complex(r.real(), acc.imag()); + } + + std::complex Ei(const std::complex&z) { + return _Ei(z, std::complex(0.0)); + } + + std::complex E1(const std::complex& z) { + if (z.imag() < 0.0) { + return -_Ei(-z, std::complex(0.0, -M_PI)); + } + else if (z.imag() > 0.0 || z.real() < 0.0) { + return -_Ei(-z, std::complex(0.0, M_PI)); + } + else { + return -Ei(-z); + } } // Reference: // https://functions.wolfram.com/GammaBetaErf/ExpIntegralEi/introductions/ExpIntegrals/ShowAll.html - std::complex Si(std::complex z) { - const std::complex i(0.0, 1.0); - - return 0.25*i*(2.0*(Ei(-i*z) - Ei(i*z)) - + std::log(i/z) - std::log(-i/z) - std::log(-i*z) - + std::log(i*z)); + std::complex Si(const std::complex& z) { + if (std::abs(z) <= 0.2) { + std::complex s(0), nn(z); + Size k; + for (k=2; k < 100 && s != s+nn; ++k) { + s += nn; + nn *= -z*z/((2.0*k-2)*(2*k-1)*(2*k-1))*(2.0*k-3); + } + QL_REQUIRE(k < 100, "series conversion issue for Si(" << z << ")"); + + return s; + } + else { + const std::complex i(0.0, 1.0); + return 0.5*i*(E1(-i*z) - E1(i*z) + - std::complex(0.0, ((z.real() >= 0 && z.imag() >= 0) + || (z.real() > 0 && z.imag() < 0) )? M_PI : -M_PI)); + } } - std::complex Ci(std::complex z) { + std::complex Ci(const std::complex& z) { const std::complex i(0.0, 1.0); - return 0.25*(2.0*(Ei(-i*z) + Ei(i*z)) - + std::log(i/z) + std::log(-i/z) - std::log(-i*z) - - std::log(i*z)) + std::log(z); + std::complex acc(0.0); + if (z.real() < 0.0 && z.imag() >= 0.0) + acc.imag(M_PI); + else if (z.real() <= 0.0 && z.imag() <= 0.0) + acc.imag(-M_PI); + + return -0.5*(E1(-i*z) + E1(i*z)) + acc; } } } diff --git a/ql/math/integrals/exponentialintegrals.hpp b/ql/math/integrals/exponentialintegrals.hpp index 957ca25d002..eff7d4fbbe3 100644 --- a/ql/math/integrals/exponentialintegrals.hpp +++ b/ql/math/integrals/exponentialintegrals.hpp @@ -20,23 +20,37 @@ /*! \file exponentialintegrals.hpp */ -#include +#ifndef quantlib_exponentail_integral_hpp +#define quantlib_exponentail_integral_hpp +#include #include +#ifndef M_EULER_MASCHERONI +#define M_EULER_MASCHERONI 0.5772156649015328606065121 +#endif + namespace QuantLib { /*! References: B. Rowe et al: GALSIM: The modular galaxy image simulation toolkit https://arxiv.org/abs/1407.7676 + + V. Pegoraro, P. Slusallek: + On the Evaluation of the Complex-Valued Exponential Integral + https://www.sci.utah.edu/~vpegorar/research/2011_JGT.pdf + */ + namespace ExponentialIntegral { Real Si(Real x); Real Ci(Real x); - std::complex Ci(std::complex z); - std::complex Si(std::complex z); - std::complex E1(std::complex z); - std::complex Ei(std::complex z); + std::complex Ci(const std::complex& z); + std::complex Si(const std::complex& z); + std::complex E1(const std::complex& z); + std::complex Ei(const std::complex& z); } } + +#endif diff --git a/ql/math/integrals/expsinhintegral.hpp b/ql/math/integrals/expsinhintegral.hpp new file mode 100644 index 00000000000..b62c711fdf6 --- /dev/null +++ b/ql/math/integrals/expsinhintegral.hpp @@ -0,0 +1,104 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2023 Klaus Spanderen + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +/*! \file expsinhintegral.hpp +*/ + +#ifndef quantlib_exp_sinh_integral_hpp +#define quantlib_exp_sinh_integral_hpp + +#include +#include +#include + +#if BOOST_VERSION >= 106900 +#define QL_BOOST_HAS_EXP_SINH + +#include +#include + +namespace QuantLib { + /*! The exp-sinh quadrature routine provided by boost is a rapidly convergent + numerical integration scheme for holomorphic integrands. The tolerance + is used against the error estimate for the L1 norm of the integral. + */ + class ExpSinhIntegral : public Integrator { + public: + explicit ExpSinhIntegral( + Real relTolerance = std::sqrt(std::numeric_limits::epsilon()), + Size maxRefinements = 9) + : Integrator(QL_MAX_REAL, Null()), + relTolerance_(relTolerance), + exp_sinh_(maxRefinements) {} + + // For half-infinite interval [0, \inf), the exp-sinh quadrature is provided by + Real integrate(const ext::function& f) const { + setNumberOfEvaluations(0); + Real error; + Real value = exp_sinh_.integrate( + [this, &f](Real x) -> Real { increaseNumberOfEvaluations(1); return f(x); }, + relTolerance_, &error); + setAbsoluteError(error); + + return value; + } + + protected: + Real integrate(const ext::function& f, Real a, Real b) + const override { + Real error; + Real value = exp_sinh_.integrate( + [this, &f](Real x) -> Real { + increaseNumberOfEvaluations(1); + return f(x); + }, + a, b, relTolerance_, &error); + setAbsoluteError(error); + + return value; + } + + private: + const Real relTolerance_; + mutable boost::math::quadrature::exp_sinh exp_sinh_; + }; +} + +#else + +namespace QuantLib { + + class ExpSinhIntegral : public Integrator { + public: + explicit ExpSinhIntegral( Real relTolerance = 0, Size maxRefinements = 0) + : Integrator(Null(), Null()) { + QL_FAIL("boost version 1.69 or higher is required in order to use ExpSinhIntegral"); + } + + protected: + Real integrate(const ext::function& f, Real a, Real b) const override { + QL_FAIL("boost version 1.69 or higher is required in order to use ExpSinhIntegral"); + } + }; + +} + +#endif + +#endif diff --git a/ql/math/integrals/gausslobattointegral.cpp b/ql/math/integrals/gausslobattointegral.cpp index 27187a44453..df2e2971786 100644 --- a/ql/math/integrals/gausslobattointegral.cpp +++ b/ql/math/integrals/gausslobattointegral.cpp @@ -138,8 +138,8 @@ namespace QuantLib { +432*(fmll+fmrr)+625*(fml+fmr)+672*fm); // avoid 80 bit logic on x86 cpu - volatile Real dist = acc + (integral1-integral2); - if(const_cast(dist)==acc || mll<=a || b<=mrr) { + const Real dist = acc + (integral1-integral2); + if(dist==acc || mll<=a || b<=mrr) { QL_REQUIRE(m>a && b>m,"Interval contains no more machine number"); return integral1; } diff --git a/ql/math/integrals/simpsonintegral.hpp b/ql/math/integrals/simpsonintegral.hpp index 69443578147..346cb6f2b3b 100644 --- a/ql/math/integrals/simpsonintegral.hpp +++ b/ql/math/integrals/simpsonintegral.hpp @@ -44,11 +44,14 @@ namespace QuantLib { // start from the coarsest trapezoid... Size N = 1; Real I = (f(a)+f(b))*(b-a)/2.0, newI; + increaseNumberOfEvaluations(2); + Real adjI = I, newAdjI; // ...and refine it Size i = 1; do { newI = Default::integrate(f,a,b,I,N); + increaseNumberOfEvaluations(N); N *= 2; newAdjI = (4.0*newI-I)/3.0; // good enough? Also, don't run away immediately diff --git a/ql/math/integrals/tanhsinhintegral.hpp b/ql/math/integrals/tanhsinhintegral.hpp index 31def3d45e7..61d915fb740 100644 --- a/ql/math/integrals/tanhsinhintegral.hpp +++ b/ql/math/integrals/tanhsinhintegral.hpp @@ -34,8 +34,6 @@ #include namespace QuantLib { - using tanh_sinh = boost::math::quadrature::tanh_sinh; - /*! The tanh-sinh quadrature routine provided by boost is a rapidly convergent numerical integration scheme for holomorphic integrands. The tolerance is used against the error estimate for the L1 norm of the integral. @@ -55,7 +53,12 @@ namespace QuantLib { Real integrate(const ext::function& f, Real a, Real b) const override { Real error; - Real value = tanh_sinh_.integrate(f, a, b, relTolerance_, &error); + Real value = tanh_sinh_.integrate( + [this, f](Real x)-> Real { + increaseNumberOfEvaluations(1); + return f(x); + }, + a, b, relTolerance_, &error); setAbsoluteError(error); return value; @@ -63,7 +66,7 @@ namespace QuantLib { private: const Real relTolerance_; - mutable tanh_sinh tanh_sinh_; + mutable boost::math::quadrature::tanh_sinh tanh_sinh_; }; } @@ -79,8 +82,9 @@ namespace QuantLib { } protected: - Real integrate(const ext::function& f, Real a, Real b) - const override { return 0.0; } + Real integrate(const ext::function& f, Real a, Real b) const override { + QL_FAIL("boost version 1.69 or higher is required in order to use TanhSinhIntegral"); + } }; } diff --git a/ql/math/integrals/trapezoidintegral.hpp b/ql/math/integrals/trapezoidintegral.hpp index 074499d8e8e..68954e217f6 100644 --- a/ql/math/integrals/trapezoidintegral.hpp +++ b/ql/math/integrals/trapezoidintegral.hpp @@ -61,10 +61,12 @@ namespace QuantLib { // start from the coarsest trapezoid... Size N = 1; Real I = (f(a)+f(b))*(b-a)/2.0, newI; + increaseNumberOfEvaluations(2); // ...and refine it Size i = 1; do { newI = IntegrationPolicy::integrate(f,a,b,I,N); + increaseNumberOfEvaluations(N*(IntegrationPolicy::nbEvalutions()-1)); N *= IntegrationPolicy::nbEvalutions(); // good enough? Also, don't run away immediately if (std::fabs(I-newI) <= absoluteAccuracy() && i > 5) diff --git a/ql/pricingengines/vanilla/analytichestonengine.cpp b/ql/pricingengines/vanilla/analytichestonengine.cpp index 8c2a9017ee0..369bb9fe41b 100644 --- a/ql/pricingengines/vanilla/analytichestonengine.cpp +++ b/ql/pricingengines/vanilla/analytichestonengine.cpp @@ -31,10 +31,18 @@ #include #include #include +#include #include +#include #include #include #include + +#include +#include + +#include +#include #include #if defined(QL_PATCH_MSVC) @@ -121,38 +129,11 @@ namespace QuantLib { // helper class for integration class AnalyticHestonEngine::Fj_Helper { public: - Fj_Helper(const VanillaOption::arguments& arguments, - const ext::shared_ptr& model, + Fj_Helper(Real kappa, Real theta, Real sigma, Real v0, + Real s0, Real rho, const AnalyticHestonEngine* engine, ComplexLogFormula cpxLog, - Time term, - Real ratio, - Size j); - - Fj_Helper(Real kappa, - Real theta, - Real sigma, - Real v0, - Real s0, - Real rho, - const AnalyticHestonEngine* engine, - ComplexLogFormula cpxLog, - Time term, - Real strike, - Real ratio, - Size j); - - Fj_Helper(Real kappa, - Real theta, - Real sigma, - Real v0, - Real s0, - Real rho, - ComplexLogFormula cpxLog, - Time term, - Real strike, - Real ratio, - Size j); + Time term, Real strike, Real ratio, Size j); Real operator()(Real phi) const; @@ -175,37 +156,11 @@ namespace QuantLib { const AnalyticHestonEngine* const engine_; }; - - AnalyticHestonEngine::Fj_Helper::Fj_Helper( - const VanillaOption::arguments& arguments, - const ext::shared_ptr& model, - const AnalyticHestonEngine* const engine, - ComplexLogFormula cpxLog, - Time term, Real ratio, Size j) - : j_ (j), //arg_(arguments), - kappa_(model->kappa()), theta_(model->theta()), - sigma_(model->sigma()), v0_(model->v0()), - cpxLog_(cpxLog), term_(term), - x_(std::log(model->process()->s0()->value())), - sx_(std::log(ext::dynamic_pointer_cast - (arguments.payoff)->strike())), - dd_(x_-std::log(ratio)), - sigma2_(sigma_*sigma_), - rsigma_(model->rho()*sigma_), - t0_(kappa_ - ((j_== 1)? model->rho()*sigma_ : Real(0))), - b_(0), g_km1_(0), - engine_(engine) - { - } - AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, Real theta, Real sigma, Real v0, Real s0, Real rho, const AnalyticHestonEngine* const engine, ComplexLogFormula cpxLog, - Time term, - Real strike, - Real ratio, - Size j) + Time term, Real strike, Real ratio, Size j) : j_(j), kappa_(kappa), @@ -226,23 +181,6 @@ namespace QuantLib { { } - AnalyticHestonEngine::Fj_Helper::Fj_Helper(Real kappa, - Real theta, - Real sigma, - Real v0, - Real s0, - Real rho, - ComplexLogFormula cpxLog, - Time term, - Real strike, - Real ratio, - Size j) - : j_(j), kappa_(kappa), theta_(theta), sigma_(sigma), v0_(v0), cpxLog_(cpxLog), term_(term), - x_(std::log(s0)), sx_(std::log(strike)), dd_(x_ - std::log(ratio)), sigma2_(sigma_ * sigma_), - rsigma_(rho * sigma_), t0_(kappa - ((j == 1) ? rho * sigma : Real(0))), b_(0), g_km1_(0), - engine_(nullptr) {} - - Real AnalyticHestonEngine::Fj_Helper::operator()(Real phi) const { const Real rpsig(rsigma_*phi); @@ -360,16 +298,169 @@ namespace QuantLib { } } + AnalyticHestonEngine::OptimalAlpha::OptimalAlpha( + const Time t, + const AnalyticHestonEngine* const enginePtr) + : t_(t), + fwd_(enginePtr->model_->process()->s0()->value() + * enginePtr->model_->process()->dividendYield()->discount(t) + / enginePtr->model_->process()->riskFreeRate()->discount(t)), + kappa_(enginePtr->model_->kappa()), + theta_(enginePtr->model_->theta()), + sigma_(enginePtr->model_->sigma()), + rho_(enginePtr->model_->rho()), + eps_(std::pow(2, -int(0.5*std::numeric_limits::digits))), + enginePtr_(enginePtr), + evaluations_(0) { + km_ = k(0.0, -1); + kp_ = k(0.0, 1); + } + + Real AnalyticHestonEngine::OptimalAlpha::alphaMax(Real strike) const { + const Real eps = 1e-8; + const auto cm = [this](Real k) -> Real { return M(k) - t_; }; + + Real alpha_max; + const Real adx = kappa_ - sigma_*rho_; + if (adx > 0.0) { + const Real kp_2pi = k(2*M_PI, 1); + + alpha_max = Brent().solve( + cm, eps_, 0.5*(kp_ + kp_2pi), (1+eps)*kp_, (1-eps)*kp_2pi + ) - 1.0; + } + else if (adx < 0.0) { + const Time tCut = -2/(kappa_ - sigma_*rho_*kp_); + if (t_ < tCut) { + const Real kp_pi = k(M_PI, 1); + alpha_max = Brent().solve( + cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi + ) - 1.0; + } + else { + alpha_max = Brent().solve( + cm, eps_, 0.5*(1.0 + kp_), + 1 + eps, (1-eps)*kp_ + ) - 1.0; + } + } + else { // adx == 0.0 + const Real kp_pi = k(M_PI, 1); + alpha_max = Brent().solve( + cm, eps_, 0.5*(kp_ + kp_pi), (1+eps)*kp_, (1-eps)*kp_pi + ) - 1.0; + } + + QL_REQUIRE(alpha_max >= 0.0, "alpha max must be larger than zero"); + + return alpha_max; + } + + std::pair + AnalyticHestonEngine::OptimalAlpha::alphaGreaterZero(Real strike) const { + const Real alpha_max = alphaMax(strike); + + return findMinima(eps_, std::max(2*eps_, (1.0-1e-6)*alpha_max), strike); + } + + Real AnalyticHestonEngine::OptimalAlpha::alphaMin(Real strike) const { + const auto cm = [this](Real k) -> Real { return M(k) - t_; }; + + const Real km_2pi = k(2*M_PI, -1); + + const Real alpha_min = Brent().solve( + cm, eps_, 0.5*(km_2pi + km_), (1-1e-8)*km_2pi, (1+1e-8)*km_ + ) - 1.0; + + QL_REQUIRE(alpha_min <= -1.0, "alpha min must be smaller than minus one"); + + return alpha_min; + } + + std::pair + AnalyticHestonEngine::OptimalAlpha::alphaSmallerMinusOne(Real strike) const { + const Real alpha_min = alphaMin(strike); + + return findMinima( + std::min(-1.0-1e-6, -1.0 + (1.0-1e-6)*(alpha_min + 1.0)), + -1.0 - eps_, strike + ); + } + + Real AnalyticHestonEngine::OptimalAlpha::operator()(Real strike) const { + try { + const std::pair minusOne = alphaSmallerMinusOne(strike); + const std::pair greaterZero = alphaGreaterZero(strike); + + if (minusOne.second < greaterZero.second) { + return minusOne.first; + } + else { + return greaterZero.first; + } + } + catch (const Error&) { + return -0.5; + } + + } + + Size AnalyticHestonEngine::OptimalAlpha::numberOfEvaluations() const { + return evaluations_; + } + + std::pair AnalyticHestonEngine::OptimalAlpha::findMinima( + Real lower, Real upper, Real strike) const { + const Real freq = std::log(fwd_/strike); + + return boost::math::tools::brent_find_minima( + [freq, this](Real alpha) -> Real { + ++evaluations_; + const std::complex z(0, -(alpha+1)); + return enginePtr_->lnChF(z, t_).real() + - std::log(alpha*(alpha+1)) + alpha*freq; + }, + lower, upper, int(0.5*std::numeric_limits::digits) + ); + } + + Real AnalyticHestonEngine::OptimalAlpha::M(Real k) const { + const Real beta = kappa_ - sigma_*rho_*k; + + if (k >= km_ && k <= kp_) { + const Real D = std::sqrt(beta*beta - sigma_*sigma_*k*(k-1)); + return std::log((beta-D)/(beta+D)) / D; + } + else { + const Real D_imag = + std::sqrt(-(beta*beta - sigma_*sigma_*k*(k-1))); + + return 2/D_imag + * ( ((beta>0.0)? M_PI : 0.0) - std::atan(D_imag/beta) ); + } + } + + Real AnalyticHestonEngine::OptimalAlpha::k(Real x, Integer sgn) const { + return ( (sigma_ - 2*rho_*kappa_) + + sgn*std::sqrt( + squared(sigma_-2*rho_*kappa_) + + 4*(kappa_*kappa_ + x*x/(t_*t_))*(1-rho_*rho_))) + /(2*sigma_*(1-rho_*rho_)); + } AnalyticHestonEngine::AP_Helper::AP_Helper( Time term, Real fwd, Real strike, ComplexLogFormula cpxLog, - const AnalyticHestonEngine* const enginePtr) + const AnalyticHestonEngine* const enginePtr, + const Real alpha) : term_(term), fwd_(fwd), strike_(strike), freq_(std::log(fwd/strike)), cpxLog_(cpxLog), - enginePtr_(enginePtr) { + enginePtr_(enginePtr), + alpha_(alpha), + s_alpha_(std::exp(alpha*freq_)) + { QL_REQUIRE(enginePtr != nullptr, "pricing engine required"); const Real v0 = enginePtr->model_->v0(); @@ -385,7 +476,7 @@ namespace QuantLib { break; case AndersenPiterbargOptCV: vAvg_ = -8.0*std::log(enginePtr->chF( - std::complex(0, -0.5), term).real())/term; + std::complex(0, alpha_), term).real())/term; break; case AsymptoticChF: phi_ = -(v0+term*kappa*theta)/sigma @@ -398,6 +489,16 @@ namespace QuantLib { *(v0 + kappa*theta*term) - 2*kappa*theta*std::atan(rho/std::sqrt(1-rho*rho)))) /(sigma*sigma); + case AngledContour: + vAvg_ = (1-std::exp(-kappa*term))*(v0 - theta) + /(kappa*term) + theta; + case AngledContourNoCV: + { + const Real r = rho - sigma*freq_ / (v0 + kappa*theta*term); + tanPhi_ = std::tan( + (r*freq_ < 0.0)? M_PI/12*boost::math::sign(freq_) : 0.0 + ); + } break; default: QL_FAIL("unknown control variate"); @@ -411,34 +512,52 @@ namespace QuantLib { == std::complex(0.0), "only Heston model is supported"); - const std::complex z(u, -0.5); - - std::complex phiBS; - - switch (cpxLog_) { - case AndersenPiterbarg: - case AndersenPiterbargOptCV: - phiBS = std::exp( - -0.5*vAvg_*term_*(z*z + std::complex(-z.imag(), z.real()))); - break; - case AsymptoticChF: - phiBS = std::exp(u*phi_ + psi_); - break; - default: - QL_FAIL("unknown control variate"); + constexpr std::complex i(0, 1); + + if (cpxLog_ == AngledContour || cpxLog_ == AngledContourNoCV || cpxLog_ == AsymptoticChF) { + const std::complex h_u(u, u*tanPhi_ - alpha_); + const std::complex hPrime(h_u-i); + + std::complex phiBS(0.0); + if (cpxLog_ == AngledContour) + phiBS = std::exp( + -0.5*vAvg_*term_*(hPrime*hPrime + + std::complex(-hPrime.imag(), hPrime.real()))); + else if (cpxLog_ == AsymptoticChF) + phiBS = std::exp(u*std::complex(1, tanPhi_)*phi_ + psi_); + + return std::exp(-u*tanPhi_*freq_) + *(std::exp(std::complex(0.0, u*freq_)) + *std::complex(1, tanPhi_) + *(phiBS - enginePtr_->chF(hPrime, term_))/(h_u*hPrime) + ).real()*s_alpha_; } + else if (cpxLog_ == AndersenPiterbarg || cpxLog_ == AndersenPiterbargOptCV) { + const std::complex z(u, -alpha_); + const std::complex zPrime(u, -alpha_-1); + const std::complex phiBS = std::exp( + -0.5*vAvg_*term_*(zPrime*zPrime + + std::complex(-zPrime.imag(), zPrime.real())) + ); - return (std::exp(std::complex(0.0, u*freq_)) - * (phiBS - enginePtr_->chF(z, term_)) / (u*u + 0.25)).real(); + return (std::exp(std::complex (0.0, u*freq_)) + * (phiBS - enginePtr_->chF(zPrime, term_)) / (z*zPrime) + ).real()*s_alpha_; + } + else + QL_FAIL("unknown control variate"); } Real AnalyticHestonEngine::AP_Helper::controlVariateValue() const { - if (cpxLog_ == AndersenPiterbarg || cpxLog_ == AndersenPiterbargOptCV) { + if ( cpxLog_ == AngledContour + || cpxLog_ == AndersenPiterbarg || cpxLog_ == AndersenPiterbargOptCV) { return BlackCalculator( Option::Call, strike_, fwd_, std::sqrt(vAvg_*term_)) .value(); } else if (cpxLog_ == AsymptoticChF) { + QL_REQUIRE(alpha_ == -0.5, "alpha must be equal to -0.5"); + const std::complex phiFreq(phi_.real(), phi_.imag() + freq_); using namespace ExponentialIntegral; @@ -447,35 +566,30 @@ namespace QuantLib { -2.0*Ci(-0.5*phiFreq)*std::sin(0.5*phiFreq) +std::cos(0.5*phiFreq)*(M_PI+2.0*Si(0.5*phiFreq)))).real(); } + else if (cpxLog_ == AngledContourNoCV) { + return ((alpha_ <= 0.0)? fwd_ : 0.0) + - ((alpha_ <= -1.0)? strike_ : 0.0) + -0.5*((alpha_ == 0.0)? fwd_ :0.0) + +0.5*((alpha_ == -1.0)? strike_: 0.0); + } else QL_FAIL("unknown control variate"); } std::complex AnalyticHestonEngine::chF( const std::complex& z, Time t) const { - - const Real kappa = model_->kappa(); - const Real sigma = model_->sigma(); - const Real theta = model_->theta(); - const Real rho = model_->rho(); - const Real v0 = model_->v0(); - - const Real sigma2 = sigma*sigma; - - if (sigma > 1e-4) { - const std::complex g - = kappa + rho*sigma*std::complex(z.imag(), -z.real()); - - const std::complex D = std::sqrt( - g*g + (z*z + std::complex(-z.imag(), z.real()))*sigma2); - - const std::complex G = (g-D)/(g+D); - - return std::exp(v0/sigma2*(1.0-std::exp(-D*t))/(1.0-G*std::exp(-D*t)) - *(g-D) + kappa*theta/sigma2*((g-D)*t - -2.0*std::log((1.0-G*std::exp(-D*t))/(1.0-G)))); + if (model_->sigma() > 1e-6 || model_->kappa() < 1e-8) { + return std::exp(lnChF(z, t)); } else { + const Real kappa = model_->kappa(); + const Real sigma = model_->sigma(); + const Real theta = model_->theta(); + const Real rho = model_->rho(); + const Real v0 = model_->v0(); + + const Real sigma2 = sigma*sigma; + const Real kt = kappa*t; const Real ekt = std::exp(kt); const Real e2kt = std::exp(2*kt); @@ -484,12 +598,14 @@ namespace QuantLib { return std::exp(-(((theta - v0 + ekt*((-1 + kt)*theta + v0)) *z*zpi)/ekt)/(2.*kappa)) + + (std::exp(-(kt) - ((theta - v0 + ekt *((-1 + kt)*theta + v0))*z*zpi) /(2.*ekt*kappa))*rho*(2*theta + kt*theta - v0 - kt*v0 + ekt*((-2 + kt)*theta + v0)) *(1.0 - std::complex(-z.imag(),z.real()))*z*z) /(2.*kappa*kappa)*sigma + + (std::exp(-2*kt - ((theta - v0 + ekt *((-1 + kt)*theta + v0))*z*zpi)/(2.*ekt*kappa))*z*z*zpi *(-2*rho2*squared(2*theta + kt*theta - v0 - @@ -504,8 +620,41 @@ namespace QuantLib { } std::complex AnalyticHestonEngine::lnChF( - const std::complex& z, Time T) const { - return std::log(chF(z, T)); + const std::complex& z, Time t) const { + + const Real kappa = model_->kappa(); + const Real sigma = model_->sigma(); + const Real theta = model_->theta(); + const Real rho = model_->rho(); + const Real v0 = model_->v0(); + + const Real sigma2 = sigma*sigma; + + const std::complex g + = kappa + rho*sigma*std::complex(z.imag(), -z.real()); + + const std::complex D = std::sqrt( + g*g + (z*z + std::complex(-z.imag(), z.real()))*sigma2); + + // reduce cancelation errors, see. L. Andersen and M. Lake + std::complex r(g-D); + if (g.real()*D.real() + g.imag()*D.imag() > 0.0) { + r = -sigma2*z*std::complex(z.real(), z.imag()+1)/(g+D); + } + + std::complex y; + if (D.real() != 0.0 || D.imag() != 0.0) { + y = expm1(-D*t)/(2.0*D); + } + else + y = -0.5*t; + + const std::complex A + = kappa*theta/sigma2*(r*t - 2.0*log1p(-r*y)); + const std::complex B + = z*std::complex(z.real(), z.imag()+1)*y/(1.0-r*y); + + return A+v0*B; } AnalyticHestonEngine::AnalyticHestonEngine( @@ -515,10 +664,11 @@ namespace QuantLib { VanillaOption::arguments, VanillaOption::results>(model), evaluations_(0), - cpxLog_ (Gatheral), + cpxLog_ (OptimalCV), integration_(new Integration( Integration::gaussLaguerre(integrationOrder))), - andersenPiterbargEpsilon_(Null()) { + andersenPiterbargEpsilon_(Null()), + alpha_(-0.5) { } AnalyticHestonEngine::AnalyticHestonEngine( @@ -528,24 +678,27 @@ namespace QuantLib { VanillaOption::arguments, VanillaOption::results>(model), evaluations_(0), - cpxLog_(Gatheral), + cpxLog_(OptimalCV), integration_(new Integration(Integration::gaussLobatto( relTolerance, Null(), maxEvaluations))), - andersenPiterbargEpsilon_(Null()) { + andersenPiterbargEpsilon_(1e-40), + alpha_(-0.5) { } AnalyticHestonEngine::AnalyticHestonEngine( const ext::shared_ptr& model, ComplexLogFormula cpxLog, const Integration& integration, - const Real andersenPiterbargEpsilon) + Real andersenPiterbargEpsilon, + Real alpha) : GenericModelEngine(model), evaluations_(0), cpxLog_(cpxLog), integration_(new Integration(integration)), - andersenPiterbargEpsilon_(andersenPiterbargEpsilon) { + andersenPiterbargEpsilon_(andersenPiterbargEpsilon), + alpha_(alpha) { QL_REQUIRE( cpxLog_ != BranchCorrection || !integration.isAdaptiveIntegration(), "Branch correction does not work in conjunction " @@ -556,61 +709,92 @@ namespace QuantLib { AnalyticHestonEngine::optimalControlVariate( Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho) { - if (t > 0.1 && (v0+t*kappa*theta)/sigma*std::sqrt(1-rho*rho) < 0.055) { + if (t > 0.15 && (v0+t*kappa*theta)/sigma*std::sqrt(1-rho*rho) < 0.15 + && ((kappa- 0.5*rho*sigma)*(v0 + t*kappa*theta) + + kappa*theta*std::log(4*(1-rho*rho)))/(sigma*sigma) < 0.1) { return AsymptoticChF; } else { - return AndersenPiterbargOptCV; + return AngledContour; } } - Size AnalyticHestonEngine::numberOfEvaluations() const { return evaluations_; } - void AnalyticHestonEngine::doCalculation(Real riskFreeDiscount, - Real dividendDiscount, - Real spotPrice, - Real strikePrice, - Real term, - Real kappa, Real theta, Real sigma, Real v0, Real rho, - const TypePayoff& type, - const Integration& integration, - const ComplexLogFormula cpxLog, - const AnalyticHestonEngine* const enginePtr, - Real& value, - Size& evaluations) { + Real AnalyticHestonEngine::priceVanillaPayoff( + const ext::shared_ptr& payoff, + const Date& maturity) const { - const Real ratio = riskFreeDiscount/dividendDiscount; + const ext::shared_ptr& process = model_->process(); + const Real fwd = process->s0()->value() + * process->dividendYield()->discount(maturity) + / process->riskFreeRate()->discount(maturity); - evaluations = 0; + return priceVanillaPayoff(payoff, process->time(maturity), fwd); + } - switch(cpxLog) { + Real AnalyticHestonEngine::priceVanillaPayoff( + const ext::shared_ptr& payoff, Time maturity) const { + + const ext::shared_ptr& process = model_->process(); + const Real fwd = process->s0()->value() + * process->dividendYield()->discount(maturity) + / process->riskFreeRate()->discount(maturity); + + return priceVanillaPayoff(payoff, maturity, fwd); + } + + Real AnalyticHestonEngine::priceVanillaPayoff( + const ext::shared_ptr& payoff, + Time maturity, Real fwd) const { + + Real value; + + const ext::shared_ptr& process = model_->process(); + const DiscountFactor dr = process->riskFreeRate()->discount(maturity); + + const Real strike = payoff->strike(); + const Real spot = process->s0()->value(); + QL_REQUIRE(spot > 0.0, "negative or null underlying given"); + + const DiscountFactor df = spot/fwd; + const DiscountFactor dd = dr/df; + + const Real kappa = model_->kappa(); + const Real sigma = model_->sigma(); + const Real theta = model_->theta(); + const Real rho = model_->rho(); + const Real v0 = model_->v0(); + + evaluations_ = 0; + + switch(cpxLog_) { case Gatheral: case BranchCorrection: { const Real c_inf = std::min(0.2, std::max(0.0001, - std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*term); + std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*maturity); - const Real p1 = integration.calculate(c_inf, - Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, - cpxLog, term, strikePrice, ratio, 1))/M_PI; - evaluations += integration.numberOfEvaluations(); + const Real p1 = integration_->calculate(c_inf, + Fj_Helper(kappa, theta, sigma, v0, spot, rho, this, + cpxLog_, maturity, strike, df, 1))/M_PI; + evaluations_ += integration_->numberOfEvaluations(); - const Real p2 = integration.calculate(c_inf, - Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, - cpxLog, term, strikePrice, ratio, 2))/M_PI; - evaluations += integration.numberOfEvaluations(); + const Real p2 = integration_->calculate(c_inf, + Fj_Helper(kappa, theta, sigma, v0, spot, rho, this, + cpxLog_, maturity, strike, df, 2))/M_PI; + evaluations_ += integration_->numberOfEvaluations(); - switch (type.optionType()) + switch (payoff->optionType()) { case Option::Call: - value = spotPrice*dividendDiscount*(p1+0.5) - - strikePrice*riskFreeDiscount*(p2+0.5); + value = spot*dd*(p1+0.5) + - strike*dr*(p2+0.5); break; case Option::Put: - value = spotPrice*dividendDiscount*(p1-0.5) - - strikePrice*riskFreeDiscount*(p2-0.5); + value = spot*dd*(p1-0.5) + - strike*dr*(p2-0.5); break; default: QL_FAIL("unknown option type"); @@ -620,49 +804,59 @@ namespace QuantLib { case AndersenPiterbarg: case AndersenPiterbargOptCV: case AsymptoticChF: + case AngledContour: + case AngledContourNoCV: case OptimalCV: { const Real c_inf = - std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*term)/sigma; - - const Real fwdPrice = spotPrice / ratio; + std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*maturity)/sigma; - const Real epsilon = enginePtr->andersenPiterbargEpsilon_ - *M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount); + const Real epsilon = andersenPiterbargEpsilon_ + *M_PI/(std::sqrt(strike*fwd)*dr); const ext::function uM = [&](){ - return Integration::andersenPiterbargIntegrationLimit(c_inf, epsilon, v0, term); + return Integration::andersenPiterbargIntegrationLimit( + c_inf, epsilon, v0, maturity); }; - AP_Helper cvHelper(term, fwdPrice, strikePrice, - (cpxLog == OptimalCV) - ? optimalControlVariate(term, v0, kappa, theta, sigma, rho) - : cpxLog, - enginePtr + const ComplexLogFormula finalLog = (cpxLog_ == OptimalCV) + ? optimalControlVariate(maturity, v0, kappa, theta, sigma, rho) + : cpxLog_; + + const AP_Helper cvHelper( + maturity, fwd, strike, finalLog, this, alpha_ ); const Real cvValue = cvHelper.controlVariateValue(); - const Real h_cv = integration.calculate(c_inf, cvHelper, uM) - * std::sqrt(strikePrice * fwdPrice)/M_PI; - evaluations += integration.numberOfEvaluations(); + const Real vAvg = (1-std::exp(-kappa*maturity))*(v0-theta)/(kappa*maturity) + theta; - switch (type.optionType()) + const Real scalingFactor = (cpxLog_ != OptimalCV && cpxLog_ != AsymptoticChF) + ? std::max(0.25, std::min(1000.0, 0.25/std::sqrt(0.5*vAvg*maturity))) + : Real(1.0); + + const Real h_cv = + fwd/M_PI*integration_->calculate(c_inf, cvHelper, uM, scalingFactor); + + evaluations_ += integration_->numberOfEvaluations(); + + switch (payoff->optionType()) { case Option::Call: - value = (cvValue + h_cv)*riskFreeDiscount; + value = (cvValue + h_cv)*dr; break; case Option::Put: - value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount; + value = (cvValue + h_cv - (fwd - strike))*dr; break; default: QL_FAIL("unknown option type"); } } break; - default: QL_FAIL("unknown complex log formula"); } + + return value; } void AnalyticHestonEngine::calculate() const @@ -676,35 +870,9 @@ namespace QuantLib { ext::dynamic_pointer_cast(arguments_.payoff); QL_REQUIRE(payoff, "non plain vanilla payoff given"); - const ext::shared_ptr& process = model_->process(); - - const Real riskFreeDiscount = process->riskFreeRate()->discount( - arguments_.exercise->lastDate()); - const Real dividendDiscount = process->dividendYield()->discount( - arguments_.exercise->lastDate()); - - const Real spotPrice = process->s0()->value(); - QL_REQUIRE(spotPrice > 0.0, "negative or null underlying given"); - - const Real strikePrice = payoff->strike(); - const Real term = process->time(arguments_.exercise->lastDate()); + const Date exerciseDate = arguments_.exercise->lastDate(); - doCalculation(riskFreeDiscount, - dividendDiscount, - spotPrice, - strikePrice, - term, - model_->kappa(), - model_->theta(), - model_->sigma(), - model_->v0(), - model_->rho(), - *payoff, - *integration_, - cpxLog_, - this, - results_.value, - evaluations_); + results_.value = priceVanillaPayoff(payoff, exerciseDate); } @@ -717,15 +885,14 @@ namespace QuantLib { : intAlgo_(intAlgo), gaussianQuadrature_(std::move(gaussianQuadrature)) {} AnalyticHestonEngine::Integration - AnalyticHestonEngine::Integration::gaussLobatto(Real relTolerance, - Real absTolerance, - Size maxEvaluations) { - return Integration(GaussLobatto, + AnalyticHestonEngine::Integration::gaussLobatto( + Real relTolerance, Real absTolerance, Size maxEvaluations, bool useConvergenceEstimate) { + return Integration(GaussLobatto, ext::shared_ptr( new GaussLobattoIntegral(maxEvaluations, absTolerance, relTolerance, - false))); + useConvergenceEstimate))); } AnalyticHestonEngine::Integration @@ -798,6 +965,13 @@ namespace QuantLib { new DiscreteTrapezoidIntegrator(evaluations))); } + AnalyticHestonEngine::Integration + AnalyticHestonEngine::Integration::expSinh(Real relTolerance) { + return Integration( + ExpSinh, ext::shared_ptr( + new ExpSinhIntegral(relTolerance))); + } + Size AnalyticHestonEngine::Integration::numberOfEvaluations() const { if (integrator_ != nullptr) { return integrator_->numberOfEvaluations(); @@ -812,13 +986,15 @@ namespace QuantLib { return intAlgo_ == GaussLobatto || intAlgo_ == GaussKronrod || intAlgo_ == Simpson - || intAlgo_ == Trapezoid; + || intAlgo_ == Trapezoid + || intAlgo_ == ExpSinh; } Real AnalyticHestonEngine::Integration::calculate( Real c_inf, const ext::function& f, - const ext::function& maxBound) const { + const ext::function& maxBound, + const Real scaling) const { Real retVal; @@ -831,6 +1007,11 @@ namespace QuantLib { case GaussChebyshev2nd: retVal = (*gaussianQuadrature_)(integrand1(c_inf, f)); break; + case ExpSinh: + retVal = scaling*(*integrator_)( + [scaling, f](Real x) -> Real { return f(scaling*x);}, + 0.0, std::numeric_limits::max()); + break; case Simpson: case Trapezoid: case GaussLobatto: @@ -872,10 +1053,112 @@ namespace QuantLib { const Real uMax = Brent().solve(u_Max(c_inf, epsilon), QL_EPSILON*uMaxGuess, uMaxGuess, uMaxStep); - const Real uHatMax = Brent().solve(uHat_Max(0.5*v0*t, epsilon), - QL_EPSILON*std::sqrt(uMaxGuess), - std::sqrt(uMaxGuess), 0.1*std::sqrt(uMaxGuess)); + try { + const Real uHatMaxGuess = std::sqrt(-std::log(epsilon)/(0.5*v0*t)); + const Real uHatMax = Brent().solve(uHat_Max(0.5*v0*t, epsilon), + QL_EPSILON*uHatMaxGuess, uHatMaxGuess, 0.001*uHatMaxGuess); + + return std::max(uMax, uHatMax); + } + catch (const Error&) { + return uMax; + } + } + + + void AnalyticHestonEngine::doCalculation(Real riskFreeDiscount, + Real dividendDiscount, + Real spotPrice, + Real strikePrice, + Real term, + Real kappa, Real theta, Real sigma, Real v0, Real rho, + const TypePayoff& type, + const Integration& integration, + const ComplexLogFormula cpxLog, + const AnalyticHestonEngine* const enginePtr, + Real& value, + Size& evaluations) { + + const Real ratio = riskFreeDiscount/dividendDiscount; + + evaluations = 0; + + switch(cpxLog) { + case Gatheral: + case BranchCorrection: { + const Real c_inf = std::min(0.2, std::max(0.0001, + std::sqrt(1.0-rho*rho)/sigma))*(v0 + kappa*theta*term); + + const Real p1 = integration.calculate(c_inf, + Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, + cpxLog, term, strikePrice, ratio, 1))/M_PI; + evaluations += integration.numberOfEvaluations(); + + const Real p2 = integration.calculate(c_inf, + Fj_Helper(kappa, theta, sigma, v0, spotPrice, rho, enginePtr, + cpxLog, term, strikePrice, ratio, 2))/M_PI; + evaluations += integration.numberOfEvaluations(); + + switch (type.optionType()) + { + case Option::Call: + value = spotPrice*dividendDiscount*(p1+0.5) + - strikePrice*riskFreeDiscount*(p2+0.5); + break; + case Option::Put: + value = spotPrice*dividendDiscount*(p1-0.5) + - strikePrice*riskFreeDiscount*(p2-0.5); + break; + default: + QL_FAIL("unknown option type"); + } + } + break; + case AndersenPiterbarg: + case AndersenPiterbargOptCV: + case AsymptoticChF: + case OptimalCV: { + const Real c_inf = + std::sqrt(1.0-rho*rho)*(v0 + kappa*theta*term)/sigma; + + const Real fwdPrice = spotPrice / ratio; - return std::max(uMax, uHatMax); + const Real epsilon = enginePtr->andersenPiterbargEpsilon_ + *M_PI/(std::sqrt(strikePrice*fwdPrice)*riskFreeDiscount); + + const ext::function uM = [&](){ + return Integration::andersenPiterbargIntegrationLimit(c_inf, epsilon, v0, term); + }; + + AP_Helper cvHelper(term, fwdPrice, strikePrice, + (cpxLog == OptimalCV) + ? optimalControlVariate(term, v0, kappa, theta, sigma, rho) + : cpxLog, + enginePtr + ); + + const Real cvValue = cvHelper.controlVariateValue(); + + const Real h_cv = integration.calculate(c_inf, cvHelper, uM) + * std::sqrt(strikePrice * fwdPrice)/M_PI; + evaluations += integration.numberOfEvaluations(); + + switch (type.optionType()) + { + case Option::Call: + value = (cvValue + h_cv)*riskFreeDiscount; + break; + case Option::Put: + value = (cvValue + h_cv - (fwdPrice - strikePrice))*riskFreeDiscount; + break; + default: + QL_FAIL("unknown option type"); + } + } + break; + + default: + QL_FAIL("unknown complex log formula"); + } } } diff --git a/ql/pricingengines/vanilla/analytichestonengine.hpp b/ql/pricingengines/vanilla/analytichestonengine.hpp index 47e5f2f8d83..79d1a0a9784 100644 --- a/ql/pricingengines/vanilla/analytichestonengine.hpp +++ b/ql/pricingengines/vanilla/analytichestonengine.hpp @@ -77,6 +77,10 @@ namespace QuantLib { Interest Rate Modeling, Volume I: Foundations and Vanilla Models, Atlantic Financial Press London. + L. Andersen and M. Lake, 2018 + Robust High-Precision Option Pricing by Fourier Transforms: + Contour Deformations and Double-Exponential Quadrature, + https://papers.ssrn.com/sol3/papers.cfm?abstract_id=3231626 \ingroup vanillaengines @@ -90,6 +94,9 @@ namespace QuantLib { VanillaOption::results> { public: class Integration; + class OptimalAlpha; + class AP_Helper; + enum ComplexLogFormula { // Gatheral form of characteristic function w/o control variate Gatheral, @@ -102,6 +109,10 @@ namespace QuantLib { // Gatheral form with asymptotic expansion of the characteristic function as control variate // https://hpcquantlib.wordpress.com/2020/08/30/a-novel-control-variate-for-the-heston-model AsymptoticChF, + // angled contour shift integral with control variate + AngledContour, + // angled contour shift integral w/o control variate + AngledContourNoCV, // auto selection of best control variate algorithm from above OptimalCV }; @@ -123,15 +134,18 @@ namespace QuantLib { // over the Fourier integration algorithm AnalyticHestonEngine(const ext::shared_ptr& model, ComplexLogFormula cpxLog, const Integration& itg, - Real andersenPiterbargEpsilon = 1e-8); + Real andersenPiterbargEpsilon = 1e-25, + Real alpha = -0.5); + + void calculate() const override; // normalized characteristic function std::complex chF(const std::complex& z, Time t) const; std::complex lnChF(const std::complex& z, Time t) const; - void calculate() const override; Size numberOfEvaluations() const; + [[deprecated("Use AnalyticHestonEngine::priceVanillaPayoff instead.")]] static void doCalculation(Real riskFreeDiscount, Real dividendDiscount, Real spotPrice, @@ -149,29 +163,16 @@ namespace QuantLib { Real& value, Size& evaluations); + Real priceVanillaPayoff( + const ext::shared_ptr& payoff, + const Date& maturity) const; + + Real priceVanillaPayoff( + const ext::shared_ptr& payoff, Time maturity) const; + static ComplexLogFormula optimalControlVariate( Time t, Real v0, Real kappa, Real theta, Real sigma, Real rho); - class AP_Helper { - public: - AP_Helper(Time term, - Real fwd, - Real strike, - ComplexLogFormula cpxLog, - const AnalyticHestonEngine* enginePtr); - - Real operator()(Real u) const; - Real controlVariateValue() const; - - private: - const Time term_; - const Real fwd_, strike_, freq_; - const ComplexLogFormula cpxLog_; - const AnalyticHestonEngine* const enginePtr_; - Real vAvg_; - std::complex phi_, psi_; - }; - protected: // call back for extended stochastic volatility // plus jump diffusion engines like bates model @@ -182,10 +183,15 @@ namespace QuantLib { private: class Fj_Helper; + Real priceVanillaPayoff( + const ext::shared_ptr& payoff, + const Time maturity, Real fwd) const; + + mutable Size evaluations_; const ComplexLogFormula cpxLog_; const ext::shared_ptr integration_; - const Real andersenPiterbargEpsilon_; + const Real andersenPiterbargEpsilon_, alpha_; }; @@ -197,12 +203,13 @@ namespace QuantLib { static Integration gaussChebyshev (Size integrationOrder = 128); static Integration gaussChebyshev2nd(Size integrationOrder = 128); - // for an adaptive integration algorithm Gatheral's version has to - // be used.Be aware: using a too large number for maxEvaluations might + // Gatheral's version has to be used for an adaptive integration + // algorithm .Be aware: using a too large number for maxEvaluations might // result in a stack overflow as the these integrations are based on // recursive algorithms. static Integration gaussLobatto(Real relTolerance, Real absTolerance, - Size maxEvaluations = 1000); + Size maxEvaluations = 1000, + bool useConvergenceEstimate = false); // usually these routines have a poor convergence behavior. static Integration gaussKronrod(Real absTolerance, @@ -213,13 +220,15 @@ namespace QuantLib { Size maxEvaluations = 1000); static Integration discreteSimpson(Size evaluation = 1000); static Integration discreteTrapezoid(Size evaluation = 1000); + static Integration expSinh(Real relTolerance = 1e-8); static Real andersenPiterbargIntegrationLimit( Real c_inf, Real epsilon, Real v0, Real t); Real calculate(Real c_inf, const ext::function& f, - const ext::function& maxBound = {}) const; + const ext::function& maxBound = {}, + Real scaling = 1.0) const; Real calculate(Real c_inf, const ext::function& f, @@ -233,7 +242,8 @@ namespace QuantLib { { GaussLobatto, GaussKronrod, Simpson, Trapezoid, DiscreteTrapezoid, DiscreteSimpson, GaussLaguerre, GaussLegendre, - GaussChebyshev, GaussChebyshev2nd }; + GaussChebyshev, GaussChebyshev2nd, + ExpSinh}; Integration(Algorithm intAlgo, ext::shared_ptr quadrature); @@ -244,12 +254,58 @@ namespace QuantLib { const ext::shared_ptr gaussianQuadrature_; }; - // inline + class AnalyticHestonEngine::AP_Helper { + public: + AP_Helper(Time term, Real fwd, Real strike, + ComplexLogFormula cpxLog, + const AnalyticHestonEngine* enginePtr, + const Real alpha = -0.5); + + Real operator()(Real u) const; + Real controlVariateValue() const; + + private: + const Time term_; + const Real fwd_, strike_, freq_; + const ComplexLogFormula cpxLog_; + const AnalyticHestonEngine* const enginePtr_; + const Real alpha_, s_alpha_; + Real vAvg_, tanPhi_; + std::complex phi_, psi_; + }; + + + class AnalyticHestonEngine::OptimalAlpha { + public: + OptimalAlpha( + const Time t, + const AnalyticHestonEngine* const enginePtr); + + Real operator()(Real strike) const; + std::pair alphaGreaterZero(Real strike) const; + std::pair alphaSmallerMinusOne(Real strike) const; + + Size numberOfEvaluations() const; + Real M(Real k) const; + Real k(Real x, Integer sgn) const; + Real alphaMin(Real strike) const; + Real alphaMax(Real strike) const; + + private: + std::pair findMinima(Real lower, Real upper, Real strike) const; + + const Real t_, fwd_, kappa_, theta_, sigma_, rho_; + + const Real eps_; + + const AnalyticHestonEngine* const enginePtr_; + Real km_, kp_; + mutable Size evaluations_; + }; + - inline - std::complex AnalyticHestonEngine::addOnTerm(Real, - Time, - Size) const { + inline std::complex AnalyticHestonEngine::addOnTerm( + Real, Time, Size) const { return std::complex(0,0); } } diff --git a/ql/pricingengines/vanilla/analytichestonhullwhiteengine.cpp b/ql/pricingengines/vanilla/analytichestonhullwhiteengine.cpp index 86195faef15..80481c8bae8 100644 --- a/ql/pricingengines/vanilla/analytichestonhullwhiteengine.cpp +++ b/ql/pricingengines/vanilla/analytichestonhullwhiteengine.cpp @@ -27,7 +27,9 @@ namespace QuantLib { const ext::shared_ptr& hestonModel, ext::shared_ptr hullWhiteModel, Size integrationOrder) - : AnalyticHestonEngine(hestonModel, integrationOrder), + : AnalyticHestonEngine( + hestonModel, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder)), hullWhiteModel_(std::move(hullWhiteModel)) { setParameters(); registerWith(hullWhiteModel_); @@ -38,7 +40,10 @@ namespace QuantLib { ext::shared_ptr hullWhiteModel, Real relTolerance, Size maxEvaluations) - : AnalyticHestonEngine(hestonModel, relTolerance, maxEvaluations), + : AnalyticHestonEngine( + hestonModel, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLobatto( + relTolerance, Null(), maxEvaluations)), hullWhiteModel_(std::move(hullWhiteModel)) { setParameters(); registerWith(hullWhiteModel_); diff --git a/ql/pricingengines/vanilla/batesengine.cpp b/ql/pricingengines/vanilla/batesengine.cpp index eb23662fd84..54068dee713 100644 --- a/ql/pricingengines/vanilla/batesengine.cpp +++ b/ql/pricingengines/vanilla/batesengine.cpp @@ -25,11 +25,16 @@ namespace QuantLib { BatesEngine::BatesEngine(const ext::shared_ptr & model, Size integrationOrder) - : AnalyticHestonEngine(model, integrationOrder) { } + : AnalyticHestonEngine( + model, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder)) { } BatesEngine::BatesEngine(const ext::shared_ptr& model, Real relTolerance, Size maxEvaluations) - : AnalyticHestonEngine(model, relTolerance, maxEvaluations) { } + : AnalyticHestonEngine( + model, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLobatto( + relTolerance, Null(), maxEvaluations)) { } std::complex BatesEngine::addOnTerm( Real phi, Time t, Size j) const { @@ -81,12 +86,17 @@ namespace QuantLib { BatesDoubleExpEngine::BatesDoubleExpEngine( const ext::shared_ptr & model, Size integrationOrder) - : AnalyticHestonEngine(model, integrationOrder) { } + : AnalyticHestonEngine( + model, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLaguerre(integrationOrder)) { } BatesDoubleExpEngine::BatesDoubleExpEngine( const ext::shared_ptr& model, Real relTolerance, Size maxEvaluations) - : AnalyticHestonEngine(model, relTolerance, maxEvaluations) { } + : AnalyticHestonEngine( + model, AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLobatto( + relTolerance, Null(), maxEvaluations)) { } std::complex BatesDoubleExpEngine::addOnTerm( Real phi, Time t, Size j) const { diff --git a/ql/pricingengines/vanilla/exponentialfittinghestonengine.cpp b/ql/pricingengines/vanilla/exponentialfittinghestonengine.cpp index f80d148e8c7..d9a69af46fc 100644 --- a/ql/pricingengines/vanilla/exponentialfittinghestonengine.cpp +++ b/ql/pricingengines/vanilla/exponentialfittinghestonengine.cpp @@ -23,7 +23,8 @@ #include #include #include -#include +#include +#include #include namespace QuantLib { @@ -185,12 +186,13 @@ namespace QuantLib { ExponentialFittingHestonEngine::ExponentialFittingHestonEngine( const ext::shared_ptr& model, ControlVariate cv, - Real scaling) + Real scaling, Real alpha) : GenericModelEngine(model), cv_(cv), scaling_(scaling), + alpha_(alpha), analyticEngine_(ext::make_shared(model, 1)) { if (moneyness_.empty()) { @@ -236,22 +238,22 @@ namespace QuantLib { const Real sigma = model_->sigma(); const Real rho = model_->rho(); + QL_REQUIRE(cv_ != ControlVariate::Gatheral && cv_ != ControlVariate::BranchCorrection, + "Gatheral and Branch-Correction are not supported as control-variate"); + const AnalyticHestonEngine::ComplexLogFormula analyticCV = - (cv_ == AndersenPiterbarg)? AnalyticHestonEngine::AndersenPiterbarg - : (cv_ == AndersenPiterbargOptCV)? - AnalyticHestonEngine::AndersenPiterbargOptCV - : (cv_ == AsymptoticChF)? - AnalyticHestonEngine::AsymptoticChF - : AnalyticHestonEngine::optimalControlVariate(t, v0, kappa, theta, sigma, rho); + (cv_ == ControlVariate::OptimalCV) + ? AnalyticHestonEngine::optimalControlVariate(t, v0, kappa, theta, sigma, rho) + : cv_; const AnalyticHestonEngine::AP_Helper helper( - t, fwd, strike, analyticCV, analyticEngine_.get()); + t, fwd, strike, analyticCV, analyticEngine_.get(), alpha_); const Real vAvg = (1-std::exp(-kappa*t))*(v0-theta)/(kappa*t) + theta; const Real scalingFactor = (scaling_ == Null()) ? (analyticCV != AnalyticHestonEngine::AsymptoticChF) - ? std::max(0.01, std::min(10.0, 0.25/std::sqrt(0.5*vAvg*t))) + ? std::max(0.25, std::min(1000.0, 0.25/std::sqrt(0.5*vAvg*t))) : Real(1.0) : scaling_; @@ -287,7 +289,7 @@ namespace QuantLib { s += w_i*u*helper(u*x_i); } - const Real h_cv = s * std::sqrt(strike * fwd)/M_PI; + const Real h_cv = s * fwd/M_PI; const Real cvValue = helper.controlVariateValue(); switch (payoff->optionType()) diff --git a/ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp b/ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp index 7be2d457a10..b462ad89d56 100644 --- a/ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp +++ b/ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp @@ -24,9 +24,8 @@ #ifndef quantlib_exponential_fitting_heston_engine_hpp #define quantlib_exponential_fitting_heston_engine_hpp -#include -#include #include +#include #include @@ -51,27 +50,19 @@ namespace QuantLib { VanillaOption::arguments, VanillaOption::results> { public: - enum ControlVariate { - // Gatheral form with Andersen-Piterbarg control variate - AndersenPiterbarg, - // same as AndersenPiterbarg, but a slightly better control variate - AndersenPiterbargOptCV, - // Gatheral form with asymptotic expansion of the characteristic function as control variate - AsymptoticChF, - // auto selection of best control variate algorithm from above - OptimalCV - }; + typedef AnalyticHestonEngine::ComplexLogFormula ControlVariate; explicit ExponentialFittingHestonEngine( const ext::shared_ptr& model, - ControlVariate cv = OptimalCV, - Real scaling = Null()); + ControlVariate cv = ControlVariate::OptimalCV, + Real scaling = Null(), + Real alpha = -0.5); void calculate() const override; private: const ControlVariate cv_; - const Real scaling_; + const Real scaling_, alpha_; const ext::shared_ptr analyticEngine_; static std::vector moneyness_; diff --git a/ql/termstructures/volatility/equityfx/hestonblackvolsurface.cpp b/ql/termstructures/volatility/equityfx/hestonblackvolsurface.cpp index ea8e0f4c3d1..ff259e63f3e 100644 --- a/ql/termstructures/volatility/equityfx/hestonblackvolsurface.cpp +++ b/ql/termstructures/volatility/equityfx/hestonblackvolsurface.cpp @@ -74,38 +74,23 @@ namespace QuantLib { } Volatility HestonBlackVolSurface::blackVolImpl(Time t, Real strike) const { - const ext::shared_ptr process = hestonModel_->process(); + AnalyticHestonEngine hestonEngine( + hestonModel_.currentLink(), cpxLogFormula_, integration_); + + const ext::shared_ptr& process = hestonModel_->process(); const DiscountFactor df = process->riskFreeRate()->discount(t, true); - const DiscountFactor div = process->dividendYield()->discount(t, true); - const Real spotPrice = process->s0()->value(); - const Real fwd = spotPrice - * process->dividendYield()->discount(t, true) - / process->riskFreeRate()->discount(t, true); + const Real fwd = process->s0()->value() + * process->dividendYield()->discount(t, true) / df; + const ext::shared_ptr payoff = + ext::make_shared( + fwd > strike ? Option::Put : Option::Call, strike); - const PlainVanillaPayoff payoff( - fwd > strike ? Option::Put : Option::Call, strike); + const Real npv = hestonEngine.priceVanillaPayoff(payoff, t); - const Real kappa = hestonModel_->kappa(); const Real theta = hestonModel_->theta(); - const Real rho = hestonModel_->rho(); - const Real sigma = hestonModel_->sigma(); - const Real v0 = hestonModel_->v0(); - - AnalyticHestonEngine hestonEngine( - hestonModel_.currentLink(), cpxLogFormula_, integration_); - - Real npv; - Size evaluations; - - AnalyticHestonEngine::doCalculation( - df, div, spotPrice, strike, t, - kappa, theta, sigma, v0, rho, - payoff, integration_, cpxLogFormula_, - &hestonEngine, npv, evaluations); - if (npv <= 0.0) return std::sqrt(theta); Brent solver; @@ -113,8 +98,12 @@ namespace QuantLib { const Volatility guess = std::sqrt(theta); constexpr double accuracy = std::numeric_limits::epsilon(); - return solver.solve([&](Volatility _v) { return blackValue(payoff.optionType(), strike, fwd, - t, _v, df, npv); }, - accuracy, guess, 0.01); + return solver.solve( + [&](Volatility _v) { + return blackValue( + payoff->optionType(), strike, fwd, t, _v, df, npv); + }, + accuracy, guess, 0.01 + ); } } diff --git a/ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp b/ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp index 9b435b5c917..b0821cc0079 100644 --- a/ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp +++ b/ql/termstructures/volatility/equityfx/hestonblackvolsurface.hpp @@ -35,9 +35,10 @@ namespace QuantLib { public: explicit HestonBlackVolSurface( const Handle& hestonModel, - AnalyticHestonEngine::ComplexLogFormula cpxLogFormula = AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::ComplexLogFormula cpxLogFormula = + AnalyticHestonEngine::AngledContour, AnalyticHestonEngine::Integration integration = - AnalyticHestonEngine::Integration::gaussLaguerre(164)); + AnalyticHestonEngine::Integration::gaussLaguerre(160)); DayCounter dayCounter() const override; Date maxDate() const override; diff --git a/test-suite/functions.cpp b/test-suite/functions.cpp index 5b2b3512501..10fca7e9887 100644 --- a/test-suite/functions.cpp +++ b/test-suite/functions.cpp @@ -25,6 +25,7 @@ #include #include #include +#include using namespace QuantLib; using namespace boost::unit_test_framework; @@ -293,6 +294,35 @@ BOOST_AUTO_TEST_CASE(testWeightedModifiedBesselFunctions) { } } } + +BOOST_AUTO_TEST_CASE(testExpm1) { + BOOST_TEST_MESSAGE("Testing complex valued expm1..."); + + const std::complex z = std::complex(1.2, 0.5); + BOOST_CHECK_SMALL(std::abs(std::exp(z) - 1.0 - expm1(z)), 10*QL_EPSILON); + + const std::complex calculated = expm1(std::complex(5e-6, 5e-5)); + //scipy reference value + const std::complex expected(4.998762493771078e-06,5.000024997979157e-05); + const Real tol = std::max(2.2e-14, 100*QL_EPSILON); + BOOST_CHECK_CLOSE_FRACTION(calculated.real(), expected.real(), tol); + BOOST_CHECK_CLOSE_FRACTION(calculated.imag(), expected.imag(), tol); +} + +BOOST_AUTO_TEST_CASE(testLog1p) { + BOOST_TEST_MESSAGE("Testing complex valued log1p..."); + + const std::complex z = std::complex(1.2, 0.57); + BOOST_CHECK_SMALL(std::abs(std::log(1.0+z) - log1p(z)), 10*QL_EPSILON); + + const std::complex calculated = log1p(std::complex(5e-6, 5e-5)); + //scipy reference value + const std::complex expected(5.0012374875401984e-06, 4.999974995958395e-05); + const Real tol = std::max(2.2e-14, 100*QL_EPSILON); + BOOST_CHECK_CLOSE_FRACTION(calculated.real(), expected.real(), tol); + BOOST_CHECK_CLOSE_FRACTION(calculated.imag(), expected.imag(), tol); +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END() diff --git a/test-suite/hestonmodel.cpp b/test-suite/hestonmodel.cpp index 09b87d50251..1cbd8fed317 100644 --- a/test-suite/hestonmodel.cpp +++ b/test-suite/hestonmodel.cpp @@ -1211,10 +1211,10 @@ BOOST_AUTO_TEST_CASE(testAnalyticPiecewiseTimeDependent) { const Real expected = option.NPV(); option.setPricingEngine(ext::shared_ptr( - new AnalyticPTDHestonEngine(model))); + new AnalyticPTDHestonEngine(model, 192))); const Real calculatedGatheral = option.NPV(); - if (std::fabs(calculatedGatheral-expected) > 1e-12) { + if (std::fabs(calculatedGatheral-expected) > 1e-7) { BOOST_ERROR("failed to reproduce Heston prices with Gatheral ChF" << "\n calculated: " << calculatedGatheral << "\n expected: " << expected); @@ -1224,10 +1224,10 @@ BOOST_AUTO_TEST_CASE(testAnalyticPiecewiseTimeDependent) { new AnalyticPTDHestonEngine( model, AnalyticPTDHestonEngine::AndersenPiterbarg, - AnalyticPTDHestonEngine::Integration::gaussLaguerre(164)))); + AnalyticPTDHestonEngine::Integration::gaussLobatto(1e-12, Null(), 100000)))); const Real calculatedAndersenPiterbarg = option.NPV(); - if (std::fabs(calculatedAndersenPiterbarg-expected) > 1e-8) { + if (std::fabs(calculatedAndersenPiterbarg-expected) > 1e-9) { BOOST_ERROR("failed to reproduce Heston prices Andersen-Piterbarg" << "\n calculated: " << calculatedAndersenPiterbarg << "\n expected: " << expected); @@ -1354,14 +1354,35 @@ BOOST_AUTO_TEST_CASE(testAlanLewisReferencePrices) { new AnalyticHestonEngine( model, AnalyticHestonEngine::AndersenPiterbarg, - AnalyticHestonEngine::Integration::discreteTrapezoid(92), - QL_EPSILON)); + AnalyticHestonEngine::Integration::gaussLobatto( + Null(), 1e-14, 1000000), + QL_EPSILON) + ); + + const ext::shared_ptr angledContourEngine( + new AnalyticHestonEngine( + model, + AnalyticHestonEngine::AngledContour, + AnalyticHestonEngine::Integration::gaussLobatto( + Null(), 1e-14, 1000000), + QL_EPSILON) + ); + + const ext::shared_ptr optimalCvEngine( + new AnalyticHestonEngine( + model, + AnalyticHestonEngine::OptimalCV, + AnalyticHestonEngine::Integration::gaussLobatto( + Null(), 1e-14, 1000000), + QL_EPSILON) + ); const Real strikes[] = { 80, 90, 100, 110, 120 }; const Option::Type types[] = { Option::Put, Option::Call }; const ext::shared_ptr engines[] = { laguerreEngine, gaussLobattoEngine, - cosEngine, andersenPiterbargEngine, exponentialFittingEngine }; + cosEngine, andersenPiterbargEngine, exponentialFittingEngine, + angledContourEngine, optimalCvEngine}; const Real expectedResults[][2] = { { 7.958878113256768285213263077598987193482161301733, @@ -1761,6 +1782,30 @@ BOOST_AUTO_TEST_CASE(testAllIntegrationMethods) { AnalyticHestonEngine::AndersenPiterbarg, true, expected, 1e-6, Null(), "Trapezoid with Andersen Piterbarg control variate"); + + // Angled contour shift integral + reportOnIntegrationMethodTest(option, model, + AnalyticHestonEngine::Integration::gaussLaguerre(), + AnalyticHestonEngine::AngledContour, + false, expected, tol, 128, + "Angled contour shift integral"); + + // Angled contour shift integral w/o control variate + reportOnIntegrationMethodTest(option, model, + AnalyticHestonEngine::Integration::gaussLaguerre(192), + AnalyticHestonEngine::AngledContourNoCV, + false, expected, tol, 192, + "Angled contour shift integral without control variate"); + +#ifdef QL_BOOST_HAS_EXP_SINH + // Angled contour shift integral with expSinh + reportOnIntegrationMethodTest(option, model, + AnalyticHestonEngine::Integration::expSinh(), + AnalyticHestonEngine::AngledContour, + true, expected, 1e-8, Null(), + "exp-sinh integration with angled contour shift integral"); +#endif + } BOOST_AUTO_TEST_CASE(testCosHestonCumulants) { @@ -2092,7 +2137,7 @@ BOOST_AUTO_TEST_CASE(testAndersenPiterbargPricing) { }; const ext::shared_ptr analyticEngine( - ext::make_shared(model, 178)); + ext::make_shared(model, 192)); const Date maturityDates[] = { settlementDate + Period(1, Days), @@ -2104,7 +2149,7 @@ BOOST_AUTO_TEST_CASE(testAndersenPiterbargPricing) { const Option::Type optionTypes[] = { Option::Call, Option::Put }; const Real strikes[] = { 50, 75, 90, 100, 110, 130, 150, 200}; - const Real tol = 1e-7; + const Real tol = 1e-8; for (auto maturityDate : maturityDates) { const ext::shared_ptr exercise = ext::make_shared(maturityDate); @@ -2960,7 +3005,8 @@ BOOST_AUTO_TEST_CASE(testExponentialFitting4StrikesAndMaturities) { expected = referenceValues[idx] - (fwd - strike)*df; const Real diff = std::fabs(calculated - expected); - if (diff > 1e-12) { + const Real tol = 1e-8; + if (diff > tol) { BOOST_ERROR("failed to reproduce cached extreme " "Heston model prices with exponential fitted " "Gauss-Laguerre quadrature rule" @@ -2969,7 +3015,7 @@ BOOST_AUTO_TEST_CASE(testExponentialFitting4StrikesAndMaturities) { << "\n expected : " << expected << "\n calculated: " << calculated << "\n diff : " << diff - << "\n tolerance : " << 1e-12); + << "\n tolerance : " << tol); } } } @@ -3015,17 +3061,16 @@ BOOST_AUTO_TEST_CASE(testOptimalControlVariateChoice) { } calculated = AnalyticHestonEngine::optimalControlVariate( - t, v0, kappa, theta, 0.2, rho); - if (calculated != AnalyticHestonEngine::AndersenPiterbargOptCV) { + t, v0, kappa, theta, 0.05, rho); + if (calculated != AnalyticHestonEngine::AngledContour) { BOOST_ERROR("failed to reproduce optimal control variate choice"); } calculated = AnalyticHestonEngine::optimalControlVariate( - t, 0.2, kappa, theta, sigma, rho); - if (calculated != AnalyticHestonEngine::AndersenPiterbargOptCV) { + t, 0.5, kappa, theta, sigma, rho); + if (calculated != AnalyticHestonEngine::AngledContour) { BOOST_ERROR("failed to reproduce optimal control variate choice"); } - } BOOST_AUTO_TEST_CASE(testAsymptoticControlVariate) { @@ -3210,6 +3255,126 @@ BOOST_AUTO_TEST_CASE(testLocalVolFromHestonModel) { } } +BOOST_AUTO_TEST_CASE(testOptimalAlphaKmin) { + BOOST_TEST_MESSAGE("Testing optimal Alpha k_min for characteristic function..."); + + SavedSettings backup; + + const Date todaysDate(1, January, 2023); + Settings::instance().evaluationDate() = todaysDate; + + // Example taken from figure 3 in Andersen and M. Lake, 2018 + // Robust High-Precision Option Pricing by Fourier Transforms: + const DayCounter dc = Actual365Fixed(); + const auto model = + ext::make_shared( + ext::make_shared( + Handle(flatRate(0.0, dc)), + Handle(flatRate(0.0, dc)), + Handle(ext::make_shared(150)), + 0.01, 0.1, 0.01, 2.0, 0.8 + ) + ); + + const auto engine = ext::make_shared( + model, + AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLobatto( + Null(), 1e-12, 100000) + ); + + const Real strike = 100; + + const Real alphaStar = AnalyticHestonEngine::OptimalAlpha(1.0, engine.get()) + .alphaSmallerMinusOne(strike).first; + + BOOST_CHECK_SMALL(alphaStar+3.71, 0.0051); + + const Date maturityDate = todaysDate + Period(15, Months); + + VanillaOption option( + ext::make_shared(Option::Call, strike), + ext::make_shared(maturityDate) + ); + option.setPricingEngine(engine); + const Real expected = option.NPV(); + + option.setPricingEngine( + ext::make_shared( + model, + AnalyticHestonEngine::AngledContour, + AnalyticHestonEngine::Integration::gaussLobatto( + Null(), 1e-12, 5000) + ) + ); + + if (std::abs(option.NPV() - expected)/expected > 1e-10) + BOOST_ERROR("failed to reproduce reference values with Angled-Contour"); + + option.setPricingEngine( + ext::make_shared( + model, ExponentialFittingHestonEngine::ControlVariate::OptimalCV + ) + ); + + if (std::abs(option.NPV() - expected)/expected > 1e-8) + BOOST_ERROR("failed to reproduce reference values with Angled-Contour" + "and exponential fitting"); +} + +BOOST_AUTO_TEST_CASE(testOptimalAlphaKmax) { + BOOST_TEST_MESSAGE("Testing optimal Alpha k_min for characteristic function..."); + + SavedSettings backup; + + const Date todaysDate(1, January, 2022); + Settings::instance().evaluationDate() = todaysDate; + + const DayCounter dc = Actual365Fixed(); + const Handle yTS + = Handle(flatRate(0.0, dc)); + const Handle spot + = Handle(ext::make_shared(75)); + + const Time T = 2; + const Real strike = 100; + + // case 1: kappa - sigma*rho > 0 + auto model = ext::make_shared( + ext::make_shared(yTS, yTS, spot, .1, 1.2, 0.2, 0.2, -0.8)); + auto engine = ext::make_shared(model); + Real alphaStar = AnalyticHestonEngine::OptimalAlpha(T, engine.get()) + .alphaGreaterZero(strike).first; + BOOST_CHECK_SMALL(alphaStar - 3.22615, 1e-4); + + // case 2: kappa - sigma*rho < 0, T < t_cut + model = ext::make_shared( + ext::make_shared(yTS, yTS, spot, 0.1, 1.2, 0.2, 1.5, 0.9)); + engine = ext::make_shared(model); + alphaStar = AnalyticHestonEngine::OptimalAlpha(T, engine.get()) + .alphaGreaterZero(strike).first; + BOOST_CHECK_SMALL(alphaStar - 0.31137, 1e-4); + + // case 3: kappa - sigma*rho < 0, T >= t_cut + model = ext::make_shared( + ext::make_shared(yTS, yTS, spot, 0.1, 1.2, 0.2, 2.25, 0.9)); + engine = ext::make_shared(model); + alphaStar = AnalyticHestonEngine::OptimalAlpha(T, engine.get()) + .alphaGreaterZero(strike).first; + BOOST_CHECK_SMALL(alphaStar - 0.11940, 1e-4); + + // case 4: kappa - sigma*rho == 0 + model = ext::make_shared( + ext::make_shared(yTS, yTS, spot, 0.1, 1.0, 0.2, 2.0, 0.5)); + engine = ext::make_shared(model); + alphaStar = AnalyticHestonEngine::OptimalAlpha(T, engine.get()) + .alphaGreaterZero(strike).first; + BOOST_CHECK_SMALL(alphaStar - 0.28006, 1e-4); +} +BOOST_AUTO_TEST_SUITE_END() + +BOOST_AUTO_TEST_SUITE(HestonModelExperimentalTest) + BOOST_AUTO_TEST_CASE(testAnalyticPDFHestonEngine) { BOOST_TEST_MESSAGE("Testing analytic PDF Heston engine..."); diff --git a/test-suite/integrals.cpp b/test-suite/integrals.cpp index 93e3c5a7607..fe5d5a4daa0 100644 --- a/test-suite/integrals.cpp +++ b/test-suite/integrals.cpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include #include @@ -37,6 +38,8 @@ #include #include +#include + using namespace QuantLib; using namespace boost::unit_test; @@ -187,6 +190,21 @@ BOOST_AUTO_TEST_CASE(testTanhSinh) { } #endif +#ifdef QL_BOOST_HAS_EXP_SINH +BOOST_AUTO_TEST_CASE(testExpSinh) { + BOOST_TEST_MESSAGE("Testing exp-sinh integration..."); + + const ExpSinhIntegral integrator; + testSingle(integrator, + "f(x) = Gaussian(x)", NormalDistribution(), + 0.0, std::numeric_limits::max(), 0.5); + + testSingle(integrator, + "f(x) = x*e^(-x)", [](Real x) { return x*std::exp(-x); }, + 0.0, std::numeric_limits::max(), 1.0); +} +#endif + BOOST_AUTO_TEST_CASE(testGaussLegendreIntegrator) { BOOST_TEST_MESSAGE("Testing Gauss-Legendre integrator..."); @@ -360,116 +378,117 @@ BOOST_AUTO_TEST_CASE(testExponentialIntegral) { // reference values are calculated with Mathematica or Python/mpmath const Real data[][10] = { - {0.0010000000000000000208, 0.0, 0.00099999994444444613193, 0.0, -6.330539864080593754, 0.0, -6.3295393640250381967, 0.0, 6.3315393641361493112, 0.0}, - {0.00070710678118654756077, 0.00070710678118654751747, 0.00070710682047025644818, 0.00070710674190283627304, -6.3305396140806145873, 0.785397913397448279, -6.329832507338711751, 0.7861055202179185354, 6.3312467208225174236, -0.78469130657697802259}, - {6.1232339957367660136e-20, 0.0010000000000000000208, 6.1232350162201617204e-20, 0.001000000055555557243, -6.330539364080593754, 1.570796326794896558, -6.3305398640805937539, 1.5717963267393410041, 6.330539864080593754, -1.5697963268504521119}, - {-0.00070710678118654747417, 0.00070710678118654760407, -0.00070710682047025636158, 0.00070710674190283635964, -6.3305396140806145873, 2.356194740192344837, -6.3312467208225174235, 2.3569013470128150935, 6.3298325073387117511, -2.3554871333718745805}, - {-0.0010000000000000000208, 1.2246467991473532027e-19, -0.00099999994444444613193, 1.2246465950321797194e-19, -6.330539864080593754, 3.141592653589793116, -6.3315393641361493112, 3.1415926535897931161, 6.3295393640250381967, -3.1415926535897931159}, - {-0.00070710678118654764736, -0.00070710678118654743088, -0.00070710682047025653477, -0.00070710674190283618645, -6.3305396140806145873, -2.3561947401923450819, -6.3312467208225174237, -2.3569013470128153382, 6.3298325073387117509, 2.3554871333718748256}, - {-1.8369701987210298041e-19, -0.0010000000000000000208, -1.8369705049015472569e-19, -0.001000000055555557243, -6.330539364080593754, -1.5707963267948968029, -6.3305398640805937541, -1.5717963267393412491, 6.3305398640805937538, 1.5697963268504523568}, - {0.00070710678118654738758, -0.00070710678118654769066, 0.00070710682047025627499, -0.00070710674190283644623, -6.3305396140806145873, -0.78539791339744852393, -6.3298325073387117512, -0.78610552021791878051, 6.3312467208225174234, 0.78469130657697826735}, - {0.10000000000000000555, 0.0, 0.099944461108276955702, 0.0, -1.7278683866572965838, 0.0, -1.6228128139692766136, 0.0, 1.8229239584193906159, 0.0}, - {0.07071067811865475853, 0.0707106781186547542, 0.070749950041603603669, 0.070671382625480302524, -1.7253704697591484327, 0.78289816362892975745, -1.6546990871336681256, 0.8586481132075703999, 1.7960418523846287393, -0.7171482131243632012}, - {6.123233995736766226e-18, 0.10000000000000000555, 6.1334444895968712791e-18, 0.10005557222505700112, -1.7228683861943336147, 1.5707963267948965577, -1.7278683866572965777, 1.670740787903173514, 1.7278683866572965899, -1.4708518656866196026}, - {-0.070710678118654749871, 0.07071067811865476286, -0.070749950041603595024, 0.070671382625480311198, -1.7253704697591484321, 2.3586944899608633586, -1.7960418523846287312, 2.4244444404654299234, 1.6546990871336681349, -2.2829445403822227075}, - {-0.10000000000000000555, 1.2246467991473532452e-17, -0.099944461108276955702, 1.2226067414369268591e-17, -1.7278683866572965838, 3.1415926535897931166, -1.8229239584193906159, 3.1415926535897931277, 1.6228128139692766136, -3.1415926535897931031}, - {-0.07071067811865476719, -0.070710678118654745541, -0.070749950041603612314, -0.07067138262548029385, -1.7253704697591484333, -2.3586944899608636035, -1.7960418523846287474, -2.4244444404654301511, 1.6546990871336681163, 2.2829445403822229697}, - {-1.8369701987210298678e-17, -0.10000000000000000555, -1.8400333469093536425e-17, -0.10005557222505700112, -1.7228683861943336147, -1.5707963267948968038, -1.7278683866572966021, -1.6707407879031737577, 1.7278683866572965654, 1.4708518656866198463}, - {0.070710678118654741211, -0.07071067811865477152, 0.070749950041603586379, -0.070671382625480319872, -1.7253704697591484315, -0.78289816362893000238, -1.6546990871336681441, -0.85864811320757066212, 1.7960418523846287232, 0.71714821312436342884}, - {0.25, 0.0, 0.2491335703197571641, 0.0, -0.82466306258094565309, 0.0, -0.54254326466191372953, 0.0, 1.0442826344437381945, 0.0}, - {0.17677669529663688651, 0.17677669529663687569, 0.17738935115398995617, 0.17616173766104897661, -0.80911938627521916259, 0.76977321991145556311, -0.63295764861417017313, 0.97841245803743094036, 0.98528112393626814822, -0.62363375572945104943}, - {1.5308084989341914715e-17, 0.25, 1.5468043259937507815e-17, 0.2508696848909122325, -0.79341294955282596203, 1.5707963267948965561, -0.82466306258094563794, 1.819929897114653724, 0.82466306258094566823, -1.3216627564751393958}, - {-0.17677669529663686486, 0.17677669529663689734, -0.17738935115398993475, 0.17616173766104899848, -0.80911938627521915876, 2.3718194336783375529, -0.98528112393626813018, 2.517958897860342088, 0.63295764861417019883, -2.1631801955523621542}, - {-0.25, 3.0616169978683829431e-17, -0.2491335703197571641, 3.0298246679137807606e-17, -0.82466306258094565309, 3.1415926535897931198, -1.0442826344437381945, 3.1415926535897931431, 0.54254326466191372953, -3.1415926535897930812}, - {-0.17677669529663690816, -0.17677669529663685404, -0.1773893511539899776, -0.17616173766104895473, -0.80911938627521916642, -2.3718194336783377978, -0.98528112393626816627, -2.5179588978603422901, 0.63295764861417014743, 2.163180195552362442}, - {-4.5924254968025744146e-17, -0.25, -4.6404129781529084775e-17, -0.2508696848909122325, -0.79341294955282596203, -1.5707963267948968087, -0.82466306258094569853, -1.8199298971146539613, 0.82466306258094560764, 1.3216627564751396331}, - {0.17677669529663684321, -0.17677669529663691899, 0.17738935115398991333, -0.17616173766104902036, -0.80911938627521915494, -0.769773219911455808, -0.63295764861417022453, -0.97841245803743122809, 0.98528112393626811213, 0.62363375572945125148}, - {0.5, 0.0, 0.49310741804306668916, 0.0, -0.17778407880661290134, 0.0, 0.45421990486317357992, 0.0, 0.55977359477616081175, 0.0}, - {0.35355339059327377302, 0.35355339059327375138, 0.3584268697133189464, 0.34860625536259795411, -0.11658254521497154353, 0.72290178026868503268, 0.23202371014762644078, 1.2063214162395304511, 0.46518880057756951253, -0.48946767681289259981}, - {3.0616169978683829431e-17, 0.5, 3.190788489546069124e-17, 0.50699674981966719583, -0.052776844956493615913, 1.5707963267948965502, -0.17778407880661287198, 2.0639037448379632547, 0.17778407880661293069, -1.0776889087518298763}, - {-0.35355339059327372973, 0.35355339059327379467, -0.35842686971331890493, 0.34860625536259799919, -0.11658254521497152822, 2.4186908733211080836, -0.46518880057756948275, 2.652124976776900558, -0.2320237101476263804, -1.9352712373502626237}, - {-0.5, 6.1232339957367658861e-17, -0.49310741804306668916, 5.8712695126944314728e-17, -0.17778407880661290134, 3.141592653589793131, -0.55977359477616081175, 3.1415926535897931642, -0.45421990486317357992, -3.1415926535897930366}, - {-0.35355339059327381632, -0.35355339059327370808, -0.35842686971331898787, -0.34860625536259790903, -0.11658254521497155883, -2.4186908733211083279, -0.4651888005775695423, -2.6521249767769007193, -0.23202371014762650116, 1.9352712373502629509}, - {-9.1848509936051488292e-17, -0.5, -9.5723654690421041556e-17, -0.50699674981966719583, -0.052776844956493615913, -1.5707963267948968264, -0.1777840788066129894, -2.0639037448379634696, 0.17778407880661281327, 1.0776889087518300913}, - {0.35355339059327368643, -0.35355339059327383797, 0.35842686971331886346, -0.34860625536259804427, -0.11658254521497151291, -0.72290178026868527697, 0.23202371014762632001, -1.2063214162395307784, 0.46518880057756945298, 0.48946767681289276116}, - {1.0, 0.0, 0.94608307036718301494, 0.0, 0.33740392290096813466, 0.0, 1.8951178163559367555, 0.0, 0.21938393439552027368, 0.0}, - {0.70710678118654754605, 0.70710678118654750275, 0.74519215535365930209, 0.66666481741950631398, 0.56680209825930888934, 0.53562961732242986806, 1.233466915678815284, 1.7803588648261259588, 0.099862719160197444251, -0.28997455411880742612}, - {6.1232339957367658861e-17, 1.0, 7.1960319005373041642e-17, 1.0572508753757285146, 0.83786694098020824089, 1.5707963267948965247, 0.33740392290096818619, 2.5168793971620796011, -0.33740392290096808314, -0.62471325642771357121}, - {-0.70710678118654745945, 0.70710678118654758935, -0.74519215535365923063, 0.66666481741950641427, 0.5668020982593089504, 2.605963036267363253, -0.099862719160197405024, 2.8516180994709857664, -1.2334669156788151226, -1.3612337887636670908}, - {-1.0, 1.2246467991473531772e-16, -0.94608307036718301494, 1.0305047480945999774e-16, 0.33740392290096813466, 3.1415926535897931723, -0.21938393439552027368, 3.1415926535897931934, -1.8951178163559367555, -3.1415926535897929056}, - {-0.70710678118654763265, -0.70710678118654741616, -0.74519215535365937355, -0.66666481741950621369, 0.56680209825930882828, -2.6059630362673634878, -0.099862719160197483478, -2.8516180994709858582, -1.2334669156788154453, 1.3612337887636674684}, - {-1.8369701987210297658e-16, -1.0, -2.158809570266204413e-16, -1.0572508753757285146, 0.83786694098020824089, -1.5707963267948969027, 0.33740392290096798009, -2.5168793971620797334, -0.33740392290096828924, 0.62471325642771370354}, - {0.70710678118654737286, -0.70710678118654767594, 0.74519215535365915917, -0.66666481741950651456, 0.56680209825930901147, -0.53562961732243010279, 1.2334669156788149613, -1.7803588648261263365, 0.099862719160197365796, 0.28997455411880751794}, - {1.5, 0.0, 1.3246835311721196804, 0.0, 0.47035631719539988668, 0.0, 3.301285449129797838, 0.0, 0.1000195824066326519, 0.0}, - {1.0606601717798213191, 1.0606601717798212541, 1.1839593855572855609, 0.9194789610047786165, 0.93002583013377829907, 0.22553329329273497242, 1.8495047911385570699, 2.5292224190594471214, -0.010546869128999664075, -0.16130364794487607552}, - {9.1848509936051488292e-17, 1.5, 1.3038076345658281264e-16, 1.7006525157682152449, 1.600632933361582593, 1.5707963267948964752, 0.47035631719539994775, 2.8954798579670162953, -0.4703563171953998256, -0.24611279562277693453}, - {-1.0606601717798211892, 1.060660171779821384, -1.1839593855572854849, 0.91947896100477878934, 0.93002583013377843491, 2.9160593602970581693, 0.010546869128999701077, 2.9802890056449171422, -1.8495047911385567612, -0.61237023453034594437}, - {-1.5, 1.8369701987210297658e-16, -1.3246835311721196804, 1.2215790424776667532e-16, 0.47035631719539988668, 3.1415926535897932298, -0.1000195824066326519, 3.1415926535897932111, -3.301285449129797838, -3.1415926535897926896}, - {-1.060660171779821449, -1.0606601717798211242, -1.1839593855572856369, -0.91947896100477844366, 0.93002583013377816324, -2.9160593602970583627, 0.010546869128999627072, -2.9802890056449171836, -1.8495047911385573786, 0.6123702345303462898}, - {-2.7554552980815446488e-16, -1.5, -3.9114229038065365107e-16, -1.7006525157682152449, 1.600632933361582593, -1.5707963267948970514, 0.47035631719539970344, -2.8954798579670163126, -0.47035631719540006991, 0.24611279562277695186}, - {1.0606601717798210593, -1.0606601717798215139, 1.1839593855572854089, -0.91947896100477896218, 0.93002583013377857075, -0.22553329329273516584, 1.8495047911385564526, -2.5292224190594474668, -0.010546869128999738079, 0.16130364794487611693}, - {2.0, 0.0, 1.6054129768026948486, 0.0, 0.4229808287748649957, 0.0, 4.9542343560018901634, 0.0, 0.048900510708061119567, 0.0}, - {1.4142135623730950921, 1.4142135623730950055, 1.6883194933582990243, 1.0649044719839209475, 1.1044891171909028954, -0.19981522706084708457, 2.1693935891748240916, 3.4589310472140426889, -0.039584645206981933184, -0.08229206049744467714}, - {1.2246467991473531772e-16, 2.0, 2.2208114946641807464e-16, 2.5015674333549756415, 2.4526669226469145219, 1.5707963267948963889, 0.42298082877486505138, 3.1762093035975914933, -0.42298082877486494002, 0.034616650007798203864}, - {-1.4142135623730949189, 1.4142135623730951787, -1.6883194933582989874, 1.0649044719839212109, 1.1044891171909031294, 3.3414078806506402814, 0.039584645206981962593, 3.0593005930923485567, -2.169393589174823594, 0.31733839362424952895}, - {-2.0, 2.4492935982947063545e-16, -1.6054129768026948486, 1.1135681831696612616e-16, 0.4229808287748649957, 3.1415926535897932894, -0.048900510708061119567, 3.1415926535897932219, -4.9542343560018901634, -3.1415926535897923336}, - {-1.4142135623730952653, -1.4142135623730948323, -1.6883194933582990613, -1.064904471983920684, 1.1044891171909026613, -3.3414078806506403646, 0.039584645206981903775, -3.059300593092348566, -2.1693935891748245892, -0.31733839362424937184}, - {-3.6739403974420595317e-16, -2.0, -6.6624344842268023737e-16, -2.5015674333549756415, 2.4526669226469145219, -1.5707963267948973103, 0.42298082877486482866, -3.1762093035975913914, -0.42298082877486516273, -0.03461665000779830579}, - {1.4142135623730947457, -1.4142135623730953519, 1.6883194933582989504, -1.0649044719839214744, 1.1044891171909033635, 0.19981522706084700137, 2.1693935891748230965, -3.458931047214042846, -0.039584645206981992002, 0.082292060497444686426}, - {5.0, 0.0, 1.5499312449446741373, 0.0, -0.19002974965664387862, 0.0, 40.185275355803177455, 0.0, 0.0011482955912753257973, 0.0}, - {3.5355339059327377302, 3.5355339059327375138, 3.6871508611543429709, -3.157181373909001357, -3.1547681046716106125, -2.1118502909279661892, -6.311949478580612776, 7.3697974788772077195, -0.0024132692373907452624, 0.0045042434314801607547}, - {3.0616169978683829431e-16, 5.0, 4.5436362172750887551e-15, 20.09321182569722639, 20.092063530105951065, 1.5707963267948920752, -0.19002974965664393734, 3.1207275717395707391, 0.1900297496566438199, -0.020865081850222464588}, - {-3.5355339059327372973, 3.5355339059327379467, -3.6871508611543449094, -3.1571813739090021642, -3.1547681046716114182, 5.2534429445177613695, 0.0024132692373907438925, 3.1460968970212734025, 6.3119494785806111631, 4.2282048252874106008}, - {-5.0, 6.1232339957367658861e-16, -1.5499312449446741373, -1.174343542809398959e-16, -0.19002974965664387862, 3.1415926535897932037, -0.0011482955912753257973, 3.1415926535897932376, -40.185275355803177455, -3.1415926535897750631}, - {-3.5355339059327381632, -3.5355339059327370808, -3.6871508611543410324, 3.1571813739090005499, -3.1547681046716098067, -5.2534429445177574859, 0.0024132692373907466323, -3.1460968970212733959, 6.3119494785806143889, -4.2282048252874183613}, - {-9.1848509936051488292e-16, -5.0, -1.3630908648516543815e-14, -20.09321182569722639, 20.092063530105951065, -1.5707963267949102514, -0.19002974965664370247, -3.1207275717395708086, 0.19002974965664405477, 0.020865081850222534065}, - {3.5355339059327368643, -3.5355339059327383797, 3.6871508611543468479, 3.1571813739090029713, -3.154768104671612224, 2.1118502909279700728, -6.3119494785806095501, -7.3697974788771999589, -0.0024132692373907425226, -0.004504243431480167346}, - {10.0, 0.0, 1.6583475942188740493, 0.0, -0.045456433004455372635, 0.0, 2492.2289762418777591, 0.0, 4.1569689296853242774e-6, 0.0}, - {7.0710678118654754605, 7.0710678118654750275, -3.7745175303432929615, 62.642575559232345941, 62.642571122904060317, 5.3452347019798827768, 125.28514668213645736, -7.5489559055283299711, 4.4363282856236060522e-6, -0.000079155158306803868409}, - {6.1232339957367658861e-16, 10.0, 6.7436601941373126554e-13, 1246.1144901994233444, 1246.1144860424544147, 1.5707963267942222532, -0.045456433004455405946, 3.2291439210137707199, 0.045456433004455339323, 0.087551267423977378721}, - {-7.0710678118654745945, 7.0710678118654758935, 3.7745175303433438135, 62.642575559232397046, 62.642571122904111422, -2.2036420483901403904, -4.4363282856235323221e-6, 3.1415134984314864345, -125.28514668213635515, -10.690548559118021506}, - {-10.0, 1.2246467991473531772e-15, -1.6583475942188740493, -6.6623371218859982286e-17, -0.045456433004455372635, 3.1415926535897933412, -4.1569689296853242774e-6, 3.1415926535897932385, -2492.2289762418777591, -3.1415926535870957744}, - {-7.0710678118654763265, -7.0710678118654741616, 3.7745175303432421091, -62.642575559232294835, 62.642571122904009212, 2.2036420483900386859, -4.4363282856236797822e-6, -3.1415134984314864347, -125.28514668213655957, 10.690548559118224914}, - {-1.8369701987210297658e-15, -10.0, -2.0230980582411937966e-12, -1246.1144901994233444, 1246.1144860424544147, -1.5707963267969197173, -0.045456433004455272699, -3.2291439210137705144, 0.04545643300445547257, -0.087551267423977584235}, - {7.0710678118654737286, -7.0710678118654767594, -3.7745175303433946658, -62.642575559232448152, 62.642571122904162528, -5.3452347019799844813, 125.28514668213625294, 7.5489559055281265623, 4.4363282856234585915e-6, 0.000079155158306804015139}, - {15.0, 0.0, 1.6181944437083687391, 0.0, 0.046278677674360439604, 0.0, 234955.85249076830358, 0.0, 1.9186278921478669771e-8, 0.0}, - {10.606601717798213191, 10.606601717798212541, -471.05515346873571713, -1328.2594907541181068, -1328.2594913017652105, 472.62595127241629118, -2656.5189820558856064, -942.11030841435617353, 5.4764710367805069957e-7, 1.4768856774331466472e-6}, - {9.1848509936051488292e-16, 15.0, 1.0008479153886857898e-10, 117477.92624539374493, 117477.92624537455865, 1.5707963266948118277, 0.046278677674360479423, 3.1889907705032654049, -0.046278677674360399786, 0.047398116913472073375}, - {-10.606601717798211892, 10.60660171779821384, 471.05515346873477896, -1328.2594907541203959, -1328.2594913017674995, -469.48435861882555978, -5.4764710367805350438e-7, 3.1415941304754706716, 2656.5189820558810283, -945.2519010679478431}, - {-15.0, 1.8369701987210297658e-15, -1.6181944437083687391, 7.9637292204322254463e-17, 0.046278677674360439604, 3.1415926535897933315, -1.9186278921478669771e-8, 3.1415926535897932385, -234955.85249076830358, -3.1415926531894540723}, - {-10.60660171779821449, -10.606601717798211242, 471.05515346873665531, 1328.2594907541158178, -1328.2594913017629215, 469.48435861882743613, -5.4764710367804789476e-7, -3.1415941304754706716, 2656.5189820558901845, 945.25190106794409045}, - {-2.7554552980815446488e-15, -15.0, -3.0025437461660077385e-10, -117477.92624539374493, 117477.92624537455865, -1.5707963270951509938, 0.046278677674360320148, -3.1889907705032652188, -0.04627867767436055906, -0.047398116913472259445}, - {10.606601717798210593, -10.606601717798215139, -471.05515346873384079, 1328.2594907541226849, -1328.2594913017697886, -472.62595127241441484, -2656.5189820558764503, 942.11030841435992621, 5.4764710367805630921e-7, -1.4768856774331489463e-6}, - {20.0, 0.0, 1.5482417010434398402, 0.0, 0.04441982084535331654, 0.0, 25615652.66405658882, 0.0, 9.8355252906498816904e-11, 0.0}, - {14.142135623730950921, 14.142135623730950055, 24486.68965358318626, 26235.092183360235647, 26235.092183384186671, -24485.118857281672279, 52470.184366744507202, 48973.379307191653857, -2.3951024244493564362e-8, -2.5280915398646225218e-8}, - {1.2246467991473531772e-15, 20.0, 1.4853900090407493933e-8, 12807826.332028294459, 12807826.332028294361, 1.5707963119409965288, 0.044419820845353372442, 3.1190380278383364344, -0.044419820845353260638, -0.02255462575145675408}, - {-14.142135623730949189, 14.142135623730951787, -24486.689653583186682, 26235.092183360320531, 26235.092183384271555, 24488.260449935262494, 2.3951024244493652701e-8, 3.1415926283088778398, -52470.184366744337435, 48970.23771453806322}, - {-20.0, 2.4492935982947063545e-15, -1.5482417010434398402, 1.1180354790034662799e-16, 0.04441982084535331654, 3.1415926535897931885, -9.8355252906498816904e-11, 3.1415926535897932385, -25615652.66405658882, -3.1415925941741928768}, - {-14.142135623730952653, -14.142135623730948323, -24486.689653583185838, -26235.092183360150763, 26235.092183384101787, -24488.26044993526165, 2.3951024244493476023e-8, -3.1415926283088778398, -52470.184366744676971, -48970.237714538064908}, - {-3.6739403974420595317e-15, -20.0, -4.4561700271222483452e-8, -12807826.332028294459, 12807826.332028294361, -1.5707963713565968905, 0.044419820845353148834, -3.1190380278383365344, -0.044419820845353484245, 0.022554625751456854031}, - {14.142135623730947457, -14.142135623730953519, 24486.689653583187103, -26235.092183360405415, 26235.09218338435644, 24485.118857281673122, 52470.184366744167666, -48973.37930719165217, -2.3951024244493741041e-8, 2.528091539864622434e-8}, - {24.989999999999998437, 0.0, 1.5315374843262580141, 0.0, -0.0072448862399538447391, 0.0, 2977286754.0962403117, 0.0, 5.4047414055481939457e-13, 0.0}, - {17.67059846185182207, 17.670598461851820988, -885673.16331345207799, -400248.00540652140574, -400248.00540652215797, 885674.7341097792085, -800496.01081304623667, -1771346.3266269055961, 7.5223231588352466724e-10, 3.3560562518531656391e-10}, - {1.5301961755346176992e-15, 24.989999999999998437, 2.1825789542468514196e-6, 1488643377.0481201559, 1488643377.0481201559, 1.5707941442159423724, -0.0072448862399538534498, 3.1023338111211545727, 0.0072448862399538360284, -0.039258842468638544566}, - {-17.670598461851819906, 17.670598461851823152, 885673.16331345318248, -400248.0054065240787, -400248.00540652483093, -885671.59251712672318, -7.5223231588352706351e-10, 3.1415926539253988636, 800496.01081304089075, -1771349.4682195569769}, - {-24.989999999999998437, 3.0603923510692353985e-15, -1.5315374843262580141, -1.7421457417638512695e-17, -0.0072448862399538447391, 3.1415926535897931172, -5.4047414055481939457e-13, 3.1415926535897932385, -2977286754.0962403117, -3.1415839232739762511}, - {-17.670598461851824234, -17.670598461851818824, 885673.16331345097349, 400248.00540651873277, -400248.005406519485, 885671.5925171245142, -7.5223231588352227096e-10, -3.1415926539253988636, 800496.01081305158263, 1771349.4682195613948}, - {-4.5905885266038530977e-15, -24.989999999999998437, -6.5477368627405542588e-6, -1488643377.0481201559, 1488643377.0481201559, -1.5708028745317593598, -0.0072448862399538186069, -3.1023338111211548151, 0.0072448862399538708713, 0.039258842468638787005}, - {17.670598461851817742, -17.670598461851825316, -885673.16331345428695, 400248.00540652675168, -400248.00540652750391, -885674.73410978141745, -800496.01081303554482, 1771346.3266269011781, 7.5223231588352945979e-10, -3.3560562518531458362e-10} + {1e-10, 0.0, 1.0e-10, 0.0, -22.4486352650389, 0.0, -22.4486352649389, 0.0, 22.4486352651389, 0.0}, + {7.0710678118655e-11, 7.0710678118655e-11, 7.0710678118655e-11, 7.0710678118655e-11, -22.4486352650389, 0.785398163397448, -22.4486352649682, 0.785398163468159, 22.4486352651096, -0.785398163326738}, + {3.0901699437495e-11, 9.5105651629515e-11, 3.0901699437495e-11, 9.5105651629515e-11, -22.4486352650389, 1.25663706143591, -22.448635265008, 1.25663706153102, 22.4486352650698, -1.25663706134081}, + {0.0, 1e-10, 0.0, 1.0e-10, -22.4486352650389, 1.5707963267949, -22.4486352650389, 1.5707963268949, 22.4486352650389, -1.5707963266949}, + {0.0, 1e-10, 0.0, 1.0e-10, -22.4486352650389, 1.5707963267949, -22.4486352650389, 1.5707963268949, 22.4486352650389, -1.5707963266949}, + {-8.0901699437495e-11, 5.8778525229247e-11, -8.0901699437495e-11, 5.8778525229247e-11, -22.4486352650389, 2.51327412287184, -22.4486352651198, 2.51327412293062, 22.448635264958, -2.51327412281306}, + {-1e-10, 0.0, -1.0e-10, 0.0, -22.4486352650389, 3.14159265358979, -22.4486352651389, 0.0, 22.4486352649389, -3.14159265358979}, + {-8.0901699437495e-11, -5.8778525229247e-11, -8.0901699437495e-11, -5.8778525229247e-11, -22.4486352650389, -2.51327412287184, -22.4486352651198, -2.51327412293062, 22.448635264958, 2.51327412281306}, + {0.0, -1e-10, 0.0, -1.0e-10, -22.4486352650389, -1.5707963267949, -22.4486352650389, -1.5707963268949, 22.4486352650389, 1.5707963266949}, + {3.0901699437495e-11, -9.5105651629515e-11, 3.0901699437495e-11, -9.5105651629515e-11, -22.4486352650389, -1.25663706143591, -22.448635265008, -1.25663706153102, 22.4486352650698, 1.25663706134081}, + {9.8768834059514e-11, -1.5643446504023002e-11, 9.8768834059514e-11, -1.5643446504023e-11, -22.4486352650389, -0.157079632679488, -22.4486352649402, -0.157079632695132, 22.4486352651377, 0.157079632663845}, + {0.15, 0.0, 0.149812626514082, 0.0, -1.32552404918277, 0.0, -1.16408641729839, 0.0, 1.46446167052028, 0.0}, + {0.1060660171779825, 0.1060660171779825, 0.106198510172016, 0.105933345197561, -1.31990959342105, 0.779773166034167, -1.21397624822349, 0.897221670932746, 1.42584293861861, -0.684824650588713}, + {0.0463525491562425, 0.14265847744427249, 0.0465043664443717, 0.1427686871506, -1.31535197062462, 1.25332575154654, -1.27825242518864, 1.40248660838809, 1.37065439517488, -1.11739007291224}, + {0.0, 0.15, 0.0, 0.150187626610941, -1.31427404390933, 1.5707963267949, -1.32552404918277, 1.72060895330898, 1.32552404918277, -1.42098370028081}, + {0.0, 0.15, 0.0, 0.150187626610941, -1.31427404390933, 1.5707963267949, -1.32552404918277, 1.72060895330898, 1.32552404918277, -1.42098370028081}, + {-0.1213525491562425, 0.0881677878438705, -0.121410363295163, 0.0879894647931175, -1.32164680474487, 2.51862071457814, -1.43946484971679, 2.59626744276408, 1.19687588593211, -2.41957522097486}, + {-0.15, 0.0, -0.149812626514082, 0.0, -1.32552404918277, 3.14159265358979, -1.46446167052028, 0.0, 1.16408641729839, -3.14159265358979}, + {-0.1213525491562425, -0.0881677878438705, -0.121410363295163, -0.0879894647931175, -1.32164680474487, -2.51862071457814, -1.43946484971679, -2.59626744276408, 1.19687588593211, 2.41957522097486}, + {0.0, -0.15, 0.0, -0.150187626610941, -1.31427404390933, -1.5707963267949, -1.32552404918277, -1.72060895330898, 1.32552404918277, 1.42098370028081}, + {0.0463525491562425, -0.14265847744427249, 0.0465043664443717, -0.1427686871506, -1.31535197062462, -1.25332575154654, -1.27825242518864, -1.40248660838809, 1.37065439517488, 1.11739007291224}, + {0.148153251089271, -0.0234651697560345, 0.147986276837203, -0.0233801359873959, -1.32524974813753, -0.155344509602526, -1.16622995490181, -0.182371337566645, 1.46287076355731, 0.135270572544445}, + {0.25, 0.0, 0.249133570319757, 0.0, -0.824663062580946, 0.0, -0.542543264661914, 0.0, 1.04428263444374, 0.0}, + {0.1767766952966375, 0.1767766952966375, 0.177389351153991, 0.17616173766105, -0.809119386275216, 0.769773219911456, -0.632957648614166, 0.978412458037432, 0.985281123936265, -0.623633755729451}, + {0.0772542485937375, 0.2377641290737875, 0.0779581492943877, 0.238274358309521, -0.796425249249655, 1.24741416450428, -0.745153392294084, 1.50303646097033, 0.898260598498369, -1.02852866129867}, + {0.0, 0.25, 0.0, 0.250869684890912, -0.793412949552826, 1.5707963267949, -0.824663062580946, 1.81992989711465, 0.824663062580946, -1.32166275647514}, + {0.0, 0.25, 0.0, 0.250869684890912, -0.793412949552826, 1.5707963267949, -0.824663062580946, 1.81992989711465, 0.824663062580946, -1.32166275647514}, + {-0.2022542485937375, 0.1469463130731175, -0.20252086544385, 0.146120744825161, -0.813939960005834, 2.52811043072268, -1.00626764691037, 2.64616186234439, 0.60229889383601, -2.35061809970499}, + {-0.25, 0.0, -0.249133570319757, 0.0, -0.824663062580946, 3.14159265358979, -1.04428263444374, 0.0, 0.542543264661914, -3.14159265358979}, + {-0.2022542485937375, -0.1469463130731175, -0.20252086544385, -0.146120744825161, -0.813939960005834, -2.52811043072268, -1.00626764691037, -2.64616186234439, 0.60229889383601, 2.35061809970499}, + {0.0, -0.25, 0.0, -0.250869684890912, -0.793412949552826, -1.5707963267949, -0.824663062580946, -1.81992989711465, 0.824663062580946, 1.32166275647514}, + {0.0772542485937375, -0.2377641290737875, 0.0779581492943877, -0.238274358309521, -0.796425249249655, -1.24741416450428, -0.745153392294084, -1.50303646097033, 0.898260598498369, 1.02852866129867}, + {0.246922085148785, -0.0391086162600575, 0.24614979209014, -0.0387156766342252, -0.823906068503191, -0.152275113509673, -0.546488805945054, -0.201435843693654, 1.04188216592042, 0.122428128357486}, + {1.0, 0.0, 0.946083070367183, 0.0, 0.337403922900968, 0.0, 1.89511781635594, 0.0, 0.21938393439552, 0.0}, + {0.70710678118655, 0.70710678118655, 0.745192155353662, 0.666664817419508, 0.566802098259312, 0.535629617322428, 1.23346691567882, 1.78035886482613, 0.0998627191601961, -0.289974554118806}, + {0.30901699437495, 0.95105651629515, 0.355652074843551, 0.983694298574337, 0.782614772996823, 1.09956193553216, 0.643964830804846, 2.31231301720838, -0.112533957890793, -0.475476714030747}, + {0.0, 1.0, 0.0, 1.05725087537573, 0.837866940980208, 1.5707963267949, 0.337403922900968, 2.51687939716208, -0.337403922900968, -0.624713256427714}, + {0.0, 1.0, 0.0, 1.05725087537573, 0.837866940980208, 1.5707963267949, 0.337403922900968, 2.51687939716208, -0.337403922900968, -0.624713256427714}, + {-0.80901699437495, 0.58778525229247, -0.824526943360603, 0.5349755552469, 0.491722358913221, 2.74478237579885, -0.14431784116889, 2.91012082986304, -1.43603057378731, -1.62893165104155}, + {-1.0, 0.0, -0.946083070367183, 0.0, 0.337403922900968, 3.14159265358979, -0.21938393439552, 0.0, -1.89511781635594, -3.14159265358979}, + {-0.80901699437495, -0.58778525229247, -0.824526943360603, -0.5349755552469, 0.491722358913221, -2.74478237579885, -0.14431784116889, -2.91012082986304, -1.43603057378731, 1.62893165104155}, + {0.0, -1.0, 0.0, -1.05725087537573, 0.837866940980208, -1.5707963267949, 0.337403922900968, -2.51687939716208, -0.337403922900968, 0.624713256427714}, + {0.30901699437495, -0.95105651629515, 0.355652074843551, -0.983694298574337, 0.782614772996823, -1.09956193553216, 0.643964830804846, -2.31231301720838, -0.112533957890793, 0.475476714030747}, + {0.98768834059514, -0.15643446504023, 0.939353669480516, -0.132366326809511, 0.347743692745538, -0.0857637957494435, 1.86192420379474, -0.4235071237, 0.214836056406461, 0.0577866622153682}, + {5.0, 0.0, 1.54993124494467, 0.0, -0.190029749656644, 0.0, 40.1852753558032, 0.0, 0.00114829559127533, 0.0}, + {3.53553390593275, 3.53553390593275, 3.68715086115432, -3.15718137390906, -3.15476810467167, -2.11185029092794, -6.31194947858072, 7.36979747887716, -0.00241326923739065, 0.00450424343148012}, + {1.5450849718747501, 4.75528258147575, 14.299679516973, 6.85221185491562, 6.85257226323722, -12.7303117750282, -0.931350039879264, 2.99045284011251, 0.0356665739529384, 0.0160488285537158}, + {0.0, 5.0, 0.0, 20.0932118256972, 20.0920635301059, 1.5707963267949, -0.190029749656644, 3.12072757173957, 0.190029749656644, -0.0208650818502225}, + {0.0, 5.0, 0.0, 20.0932118256972, 20.0920635301059, 1.5707963267949, -0.190029749656644, 3.12072757173957, 0.190029749656644, -0.0208650818502225}, + {-4.04508497187475, 2.93892626146235, -2.0577013528011, -1.96223940975232, -1.9637046590567, 3.61921566552724, 0.00286020292932927, 3.14261835694337, 6.84905720502975, 11.1883945116728}, + {-5.0, 0.0, -1.54993124494467, 0.0, -0.190029749656644, 3.14159265358979, -0.00114829559127533, 0.0, -40.1852753558032, -3.14159265358979}, + {-4.04508497187475, -2.93892626146235, -2.0577013528011, 1.96223940975232, -1.9637046590567, -3.61921566552724, 0.00286020292932927, -3.14261835694337, 6.84905720502975, -11.1883945116728}, + {0.0, -5.0, 0.0, -20.0932118256972, 20.0920635301059, -1.5707963267949, -0.190029749656644, -3.12072757173957, 0.190029749656644, 0.0208650818502225}, + {1.5450849718747501, -4.75528258147575, 14.299679516973, -6.85221185491562, 6.85257226323722, 12.7303117750282, -0.931350039879264, -2.99045284011251, 0.0356665739529384, -0.0160488285537158}, + {4.9384417029757, -0.7821723252011501, 1.53351371140353, 0.167535111630988, -0.252671967618136, -0.0455545136665558, 31.7637646606649, -20.6127722347705, 0.000742118122850436, 0.000971589948194675}, + {10.0, 0.0, 1.65834759421887, 0.0, -0.0454564330044554, 0.0, 2492.22897624188, 0.0, 4.15696892968532e-6, 0.0}, + {7.0710678118655, 7.0710678118655, -3.77451753034182, 62.6425755592338, 62.6425711229056, 5.34523470197841, 125.285146682139, -7.54895590552534, 4.43632828562146e-6, -7.91551583068017e-5}, + {3.0901699437495003, 9.5105651629515, 303.07292777526, -690.037761260879, -690.037754650298, -301.502129842997, -0.659900725018632, 5.27667742385125, -0.00134856502993308, 0.00415958644984393}, + {0.0, 10.0, 0.0, 1246.11449019942, 1246.11448604245, 1.5707963267949, -0.0454564330044554, 3.22914392101377, 0.0454564330044554, 0.0875512674239774}, + {0.0, 10.0, 0.0, 1246.11449019942, 1246.11448604245, 1.5707963267949, -0.0454564330044554, 3.22914392101377, 0.0454564330044554, 0.0875512674239774}, + {-8.0901699437495, 5.8778525229247, -14.6236949578037, 13.4643508624518, 13.4645870261785, 16.1946084513107, -2.79815608075126e-5, 3.14158769865141, -157.085481478947, -317.2439811058}, + {-10.0, 0.0, -1.65834759421887, 0.0, -0.0454564330044554, 3.14159265358979, -4.15696892968532e-6, 0.0, -2492.22897624188, -3.14159265358979}, + {-8.0901699437495, -5.8778525229247, -14.6236949578037, -13.4643508624518, 13.4645870261785, -16.1946084513107, -2.79815608075126e-5, -3.14158769865141, -157.085481478947, 317.2439811058}, + {0.0, -10.0, 0.0, -1246.11449019942, 1246.11448604245, -1.5707963267949, -0.0454564330044554, -3.22914392101377, 0.0454564330044554, -0.0875512674239774}, + {3.0901699437495003, -9.5105651629515, 303.07292777526, 690.037761260879, -690.037754650298, 301.502129842997, -0.659900725018632, -5.27667742385125, -0.00134856502993308, -0.00415958644984393}, + {9.8768834059514, -1.5643446504023002, 1.78956084261706, 0.114701769782499, -0.118816490702582, 0.198823504802007, 411.904076239608, -2157.22483235914, -6.48699583272709e-7, 4.66032253043785e-6}, + {25.0, 0.0, 1.53148255099996, 0.0, -0.00684859717970259, 0.0, 3005950906.52555, 0.0, 5.34889975534022e-13, 0.0}, + {17.67766952966375, 17.67766952966375, -894423.548678786, -396595.979622699, -396595.9796227, 894425.119475113, -793191.959245399, -1788847.09735757, 7.48981460647877e-10, 3.27816276287981e-10}, + {7.72542485937375, 23.77641290737875, 395787595.545024, 194501516.12134, 194501516.12134, -395787593.974227, -80.7948153607822, -39.8888851700048, 1.72503667797818e-5, 2.36415887840135e-6}, + {0.0, 25.0, 0.0, 1502975453.26277, 1502975453.26277, 1.5707963267949, -0.00684859717970259, 3.10227887779486, 0.00684859717970259, -0.0393137757949353}, + {0.0, 25.0, 0.0, 1502975453.26277, 1502975453.26277, 1.5707963267949, -0.00684859717970259, 3.10227887779486, 0.00684859717970259, -0.0393137757949353}, + {-20.22542485937375, 14.69463130731175, -19129.3494470458, 45406.0213041107, 45406.0213041213, 19130.9202433848, 5.85665949258649e-11, 3.14159265356458, -2432061.38760638, 25010638.0968068}, + {-25.0, 0.0, -1.53148255099996, 0.0, -0.00684859717970259, 3.14159265358979, -5.34889975534022e-13, 0.0, -3005950906.52555, -3.14159265358979}, + {-20.22542485937375, -14.69463130731175, -19129.3494470458, -45406.0213041107, 45406.0213041213, -19130.9202433848, 5.85665949258649e-11, -3.14159265356458, -2432061.38760638, -25010638.0968068}, + {0.0, -25.0, 0.0, -1502975453.26277, 1502975453.26277, -1.5707963267949, -0.00684859717970259, -3.10227887779486, 0.00684859717970259, 0.0393137757949353}, + {7.72542485937375, -23.77641290737875, 395787595.545024, -194501516.12134, 194501516.12134, 395787593.974227, -80.7948153607822, 39.8888851700048, 1.72503667797818e-5, -2.36415887840135e-6}, + {24.6922085148785, -3.91086162600575, 0.61973692887531, 0.318459426938049, -0.318931296543192, -0.950420524151913, -1816162045.63054, 1255955799.5082, -4.40593065675657e-13, -5.79490191675286e-13}, + {50.0, 0.0, 1.55161707248594, 0.0, -0.00562838632411631, 0.0, 1.05856368971317e+20, 0.0, 3.78326402955046e-24, 0.0}, + {35.3553390593275, 35.3553390593275, 53807668130.5995, -22948660925283.2, -22948660925283.2, -53807668129.0287, -45897321850566.4, 107615336261.199, -9.9766761181828e-21, 8.71502630154959e-18}, + {15.4508497187475, 47.5528258147575, 2.49903843573354e+18, -3.83240358282137e+18, -3.83240358282137e+18, -2.49903843573354e+18, -68343.3715391731, 77339.6040605891, 3.71621275609622e-10, 3.85406628982992e-9}, + {0.0, 50.0, 0.0, 5.29281844856585e+19, 5.29281844856585e+19, 1.5707963267949, -0.00562838632411631, 3.12241339928083, 0.00562838632411631, -0.0191792543089607}, + {0.0, 50.0, 0.0, 5.29281844856585e+19, 5.29281844856585e+19, 1.5707963267949, -0.00562838632411631, 3.12241339928083, 0.00562838632411631, -0.0191792543089607}, + {-40.4508497187475, 29.3892626146235, -57258797567.9644, -12906669326.6389, -12906669326.6389, 57258797569.5352, -8.55226617604501e-21, 3.14159265358979, 6.68228261723918e+15, -3.43017053184612e+15}, + {-50.0, 0.0, -1.55161707248594, 0.0, -0.00562838632411631, 3.14159265358979, -3.78326402955046e-24, 0.0, -1.05856368971317e+20, -3.14159265358979}, + {-40.4508497187475, -29.3892626146235, -57258797567.9644, 12906669326.6389, -12906669326.6389, -57258797569.5352, -8.55226617604501e-21, -3.14159265358979, 6.68228261723918e+15, 3.43017053184612e+15}, + {0.0, -50.0, 0.0, -5.29281844856585e+19, 5.29281844856585e+19, -1.5707963267949, -0.00562838632411631, -3.12241339928083, 0.00562838632411631, 0.0191792543089607}, + {15.4508497187475, -47.5528258147575, 2.49903843573354e+18, 3.83240358282137e+18, -3.83240358282137e+18, 2.49903843573354e+18, -68343.3715391731, -77339.6040605891, 3.71621275609622e-10, -3.85406628982992e-9}, + {49.384417029757, -7.8217232520115, -16.8292457944994, 16.9326906903424, -16.9326976506474, -18.4000381995002, 1.09489979806082e+19, -5.61228684199658e+19, -8.51344869310291e-25, 6.95142343223447e-24}, + {700.0, 0.0, 1.57199393223749, 0.0, 0.000778810012739756, 0.0, 1.45097873605256e+301, 0.0, 1.40651876623403e-307, 0.0}, + {494.974746830585, 494.974746830585, -5.39480977313549e+211, -3.7907051625115e+211, -3.7907051625115e+211, 5.39480977313549e+211, -7.58141032502299e+211, -1.0789619546271e+212, 1.26627531288803e-218, 8.89746644202181e-219}, + {216.311896062465, 665.7395614066049, 6.68861022474796e+285, -6.86204916856497e+285, -6.86204916856497e+285, -6.68861022474796e+285, 4.35129688126332e+89, -1.25283433405018e+91, 9.10599247691995e-98, -1.3494793845188e-97}, + {0.0, 700.0, 0.0, 7.2548936802628e+300, 7.2548936802628e+300, 1.5707963267949, 0.000778810012739756, 3.14279025903239, -0.000778810012739756, 0.00119760544259495}, + {0.0, 700.0, 0.0, 7.2548936802628e+300, 7.2548936802628e+300, 1.5707963267949, 0.000778810012739756, 3.14279025903239, -0.000778810012739756, 0.00119760544259495}, + {-566.311896062465, 411.449676604729, 4.13964135191794e+174, 3.47943069430311e+175, 3.47943069430311e+175, -4.13964135191794e+174, 1.39494929258574e-249, 3.14159265358979, 9.43022777090499e+242, 8.40743888884655e+242}, + {-700.0, 0.0, -1.57199393223749, 0.0, 0.000778810012739756, 3.14159265358979, -1.40651876623403e-307, 0.0, -1.45097873605256e+301, -3.14159265358979}, + {-566.311896062465, -411.449676604729, 4.13964135191794e+174, -3.47943069430311e+175, 3.47943069430311e+175, 4.13964135191794e+174, 1.39494929258574e-249, -3.14159265358979, 9.43022777090499e+242, -8.40743888884655e+242}, + {0.0, -700.0, 0.0, -7.2548936802628e+300, 7.2548936802628e+300, -1.5707963267949, 0.000778810012739756, -3.14279025903239, -0.000778810012739756, -0.00119760544259495}, + {216.311896062465, -665.7395614066049, 6.68861022474796e+285, 6.86204916856497e+285, -6.86204916856497e+285, 6.68861022474796e+285, 4.35129688126332e+89, 1.25283433405018e+91, 9.10599247691995e-98, 1.3494793845188e-97}, + {691.381838416598, -109.50412552816101, -2.38570018769502e+44, -9.72638025849046e+43, 9.72638025849046e+43, -2.38570018769502e+44, -2.15172979114587e+297, -1.50043260461905e+297, -7.44435180959991e-304, 2.26013762375079e-304} }; - - const Real tol = 1e-10; + const Real tol = 100*QL_EPSILON; for (const auto& i : data) { const Real x = i[0]; const Real y = (std::abs(i[1]) < 1e-12) ? 0.0 : i[1]; const std::complex z(x, y); - const std::complex si = Si(z); std::complex ref(i[2], i[3]); Real diff = std::abs(si-ref)/std::abs(ref); - if (diff > tol + if (diff > tol || std::isnan(diff) || (std::abs(ref.real()) < tol && std::abs(si.real()) > tol) || (std::abs(ref.imag()) < tol && std::abs(si.imag()) > tol)) { reportSiCiFail("Si", z, si, ref, diff, tol); @@ -478,7 +497,7 @@ BOOST_AUTO_TEST_CASE(testExponentialIntegral) { const std::complex ci = Ci(z); ref = std::complex(i[4], i[5]); diff = std::min(std::abs(ci-ref), std::abs(ci-ref)/std::abs(ref)); - if (diff > tol + if (diff > tol || std::isnan(diff) || (std::abs(ref.real()) < tol && std::abs(ci.real()) > tol) || (std::abs(ref.imag()) < tol && std::abs(ci.imag()) > tol)) { reportSiCiFail("Ci", z, ci, ref, diff, tol); @@ -487,7 +506,7 @@ BOOST_AUTO_TEST_CASE(testExponentialIntegral) { const std::complex ei = Ei(z); ref = std::complex(i[6], i[7]); diff = std::abs(ei-ref)/std::abs(ref); - if (diff > tol + if (diff > tol || std::isnan(diff) || (std::abs(ref.real()) < tol && std::abs(ei.real()) > tol) || (std::abs(ref.imag()) < tol && std::abs(ei.imag()) > tol)) { reportSiCiFail("Ei", z, ei, ref, diff, tol); @@ -495,9 +514,8 @@ BOOST_AUTO_TEST_CASE(testExponentialIntegral) { const std::complex e1 = E1(z); ref = std::complex(i[8], i[9]); - diff = std::min(std::abs(e1-ref), std::abs(e1-ref)/std::abs(ref)); - if (std::abs(z) < 10.0) - if (diff > tol + diff = std::abs(e1-ref)/std::abs(ref); + if (diff > 10*tol || std::isnan(diff) || (std::abs(ref.real()) < tol && std::abs(e1.real()) > tol) || (std::abs(ref.imag()) < tol && std::abs(e1.imag()) > tol)) { reportSiCiFail("E1", z, e1, ref, diff, tol); @@ -558,6 +576,67 @@ BOOST_AUTO_TEST_CASE(testRealSiCiIntegrals) { } } +BOOST_AUTO_TEST_CASE(testExponentialIntegralLimits) { + BOOST_TEST_MESSAGE("Testing limits for Ei..."); + + using namespace ExponentialIntegral; + + const Real largeValue = 0.75*std::log(0.1*QL_MAX_REAL); + + const std::complex largeValuePosImag = + Ei(std::complex(largeValue, std::numeric_limits::min())); + + const Real tol = 1000*QL_EPSILON; + + BOOST_CHECK_CLOSE(largeValuePosImag.imag(), M_PI, tol); + + BOOST_CHECK_CLOSE( + largeValuePosImag.real(), std::exp(largeValue)/largeValue, 1e3/largeValue); + + const std::complex largeValueNegImag = + Ei(std::complex(largeValue, -std::numeric_limits::min())); + + BOOST_CHECK_CLOSE(largeValueNegImag.imag(), -M_PI, tol); + BOOST_CHECK_CLOSE( + largeValueNegImag.real(), std::exp(largeValue)/largeValue, 1e3/largeValue); + + const std::complex largeValueZeroImag = + Ei(std::complex(largeValue)); + BOOST_CHECK(largeValueZeroImag.imag() == Real(0.0)); + + if (std::numeric_limits::has_infinity) { + const std::complex ei_0 = Ei(std::complex(0.0)); + BOOST_CHECK( + ei_0 == std::complex(-std::numeric_limits::infinity())); + } + + const Real smallR = QL_EPSILON*QL_EPSILON; + for (Integer x = -100; x < 100; ++x) { + const Real phi = x/100.0 * M_PI; + const std::complex z = std::polar(smallR, phi); + const std::complex ei = Ei(z); + + // principal branch + const std::complex limit_ei = M_EULER_MASCHERONI + std::log(z); + + BOOST_CHECK_CLOSE(ei.real(), limit_ei.real(), tol); + BOOST_CHECK_CLOSE(ei.imag(), limit_ei.imag(), tol); + } + + const Real largeR = largeValue; + for (Integer x = -10; x < 10; ++x) { + const Real phi = x/10.0 * M_PI; + if (std::abs(phi) > 0.5*M_PI) { + const std::complex z = std::polar(largeR, phi); + const std::complex ei = Ei(z); + + const Real limit_ei_imag = boost::math::sign(z.imag())*M_PI; + BOOST_CHECK(close_enough(ei.real(), 0.0)); + BOOST_CHECK_CLOSE(ei.imag(), limit_ei_imag, tol); + } + } +} + BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE_END() diff --git a/test-suite/marketmodel.cpp b/test-suite/marketmodel.cpp index 861e6e1404b..1611579bbfe 100644 --- a/test-suite/marketmodel.cpp +++ b/test-suite/marketmodel.cpp @@ -4457,28 +4457,24 @@ void MarketModelTest::testStochVolForwardsAndOptionlets() { Real rho = 0.0; Real v1 = v0*volCoeff*volCoeff; - - - ext::shared_ptr payoff( new PlainVanillaPayoff(Option::Call, forwardStrikes[i])); - - - Real trueValue =0.0; - Size evaluations =0; - - AnalyticHestonEngine::doCalculation( - 1.0, // no discounting - 1.0, // no discounting - todaysForwards[i] + displacement, - todaysForwards[i] + displacement, rateTimes[i], kappa, theta, - sigma, v1, rho, *payoff, - AnalyticHestonEngine::Integration::gaussLaguerre(), - // AnalyticHestonEngine::Integration::gaussLobatto(1e-8, - // 1e-8), - AnalyticHestonEngine::Gatheral, nullptr, trueValue, evaluations); - + Real trueValue = AnalyticHestonEngine( + ext::make_shared( + ext::make_shared( + Handle(flatRate(0.0, dayCounter)), + Handle(flatRate(0.0, dayCounter)), + Handle(ext::make_shared(todaysForwards[i] + displacement)), + v1, kappa, theta, sigma, rho + ) + ), + AnalyticHestonEngine::Gatheral, + AnalyticHestonEngine::Integration::gaussLaguerre(128) + ).priceVanillaPayoff( + ext::make_shared(payoff->optionType(), todaysForwards[i] + displacement), + rateTimes[i] + ); trueValue *= accruals[i] * todaysDiscounts[i + 1]; diff --git a/test-suite/normalclvmodel.cpp b/test-suite/normalclvmodel.cpp index ee9a39a412e..e37ef137bdd 100644 --- a/test-suite/normalclvmodel.cpp +++ b/test-suite/normalclvmodel.cpp @@ -515,10 +515,11 @@ BOOST_AUTO_TEST_CASE(testMoustacheGraph, *precondition(if_speed(Slow))) { } const Real expected[] = { - 0.00931214, 0.0901481, 0.138982, 0.112059, 0.0595901, - 0.0167549, -0.00906787, -0.0206768, -0.0225628, -0.0203593, - -0.016036, -0.0116629, -0.00728792, -0.00328821, - -0.000158562, 0.00502041, 0.00347706, 0.00238216, }; + 0.0093023407, 0.090099067, 0.13903052, 0.11206888, 0.059511722, + 0.016627484, -0.0091168814, -0.020725808, -0.022513829, + -0.020280851, -0.015967361, -0.011476705, -0.0071604902, + -0.003219595, -0.00014875974, 0.0049910036, 0.0034770584, + 0.0023821554}; const Real tol = 1e-5; for (Size u=0; u < n; ++u) { diff --git a/test-suite/quantlibbenchmark.cpp b/test-suite/quantlibbenchmark.cpp index 79be4b47697..472ce10d398 100644 --- a/test-suite/quantlibbenchmark.cpp +++ b/test-suite/quantlibbenchmark.cpp @@ -35,70 +35,11 @@ ./quantlib-benchmark --mp The number of floating point operations of a given test case was measured - using the perfex library, http://user.it.uu.se/~mikpe/linux/perfctr - and PAPI, http://icl.cs.utk.edu/papi - - Example results: 1. i7 7820X@3.6GHz :24192.2 mflops - 2. i7 4702HQ@2.2GHz : 6524.9 mflops - 3. i7 870@2.93GHz : 4759.2 mflops - 4. Core2 Q9300@2.5Ghz : 2272.6 mflops - 5. Core2 Q6600@2.4Ghz : 1984.0 mflops - 6. i3 540@3.1Ghz : 1755.3 mflops - 7. Raspberry Pi4@1.5GHz : 1704.2 mflops - 8. Core2 Dual@2.0Ghz : 835.9 mflops - 9. Athlon 64 X2 4400+ : 824.2 mflops - 10. Cortex-A57@2.0GHz : 821.7 mflops - 11. Core2 Dual@2.0Ghz : 754.1 mflops - 12. Pentium4 Dual@2.8Ghz : 423.8 mflops - 13. Raspberry Pi3@1.2GHz : 309.2 mflops - 14. Pentium4@3.0Ghz : 266.3 mflops - 15. PentiumIII@1.1Ghz : 146.2 mflops - 16. Alpha 2xEV68@833Mhz : 184.6 mflops - 17. Wii PowerPC 750@729MHz : 46.1 mflops - 18. Raspberry Pi ARM@700Mhz: 28.3 mflops - 19. MIPS R5000@150MHz : 12.6 mflops - 20. RISC-V on FPGA@25Mhz : 2.4 mflops - 21. Strong ARM@206Mhz : 1.4 mflops - 22. SPARC v7@25MHz : 0.78mflops - - Remarks: OS: Linux, static libs - 1. g++-6.3.0 -O3 -ffast-math -march=core-avx2 - Remark: 16 processes - 2. g++-4.8.1 -O3 -ffast-math -march=core-avx2 - Remark: eight processes - 3. gcc-4.6.3, -O3 -ffast-math -mfpmath=sse,387 -march=corei7 - Remark: eight processes - 4. icc-11.0, -gcc-version=420 -fast -fp-model fast=2 -ipo-jobs2 - Remark: four processes - 5. icc-11.0, -gcc-version=420 -fast -fp-model fast=2 -ipo-jobs2 - Remark: four processes - 6. gcc-4.4.5, -O3 -ffast-math -mfpmath=sse,387 -msse4.2 -march=core2 - Remark: four processes - 7. gcc-8.3.0, -O3 -ffast-math -mcpu=cortx-a8 -mfpu=neon-fp-armv8 - Remark: four processes - 8. icc-11.0, -gcc-version=420 -fast -fp-model fast=2 -ipo-jobs2 - Remark: two processes - 9. icc-11.0, -gcc-version=420 -xSSSE3 -O3 -ipo -no-prec-div -static - -fp-model fast=2 -ipo-jobs2, Remark: two processes - 10. clang++-6.0.1 -O2, Remark: four processes - 11. gcc-4.2.1, -O3 -ffast-math -mfpmath=sse,387 -msse3 -funroll-all-loops - Remark: two processes - 12. gcc-4.0.1, -O3 -march=pentium4 -ffast-math - -mfpmath=sse,387 -msse2 -funroll-all-loops, Remark: two processes - 13. gcc-4.9.2 -O2, Remark: four processes - 14. gcc-4.0.1, -O3 -march=pentium4 -ffast-math - -mfpmath=sse,387 -msse2 -funroll-all-loops - 15. gcc-4.1.1, -O3 -march=pentium3 -ffast-math - -mfpmath=sse,387 -msse -funroll-all-loops - 16. gcc-3.3.5, -O3 -mcpu=e67 -funroll-all-loops, Remark: two processes - 17. gcc-4.9.2, -O2 -g on a Nintendo Wii - 18. gcc-4.6.3, -O3 - 19. gcc-4-7-4, -O2 on a SGI Indy - 20. gcc-9.2, -O2 on RISC-V softcore on an Artix7 100T FPGA - 21. gcc-3.4.3, -O2 -g on a Zaurus PDA - 22. gcc-7.5.0, -O2 on a Sun SPARCstation IPC, FPU: Weitek 3170 - - This benchmark is derived from quantlibtestsuite.cpp. Please see the + using PAPI, http://icl.cs.utk.edu/papi + + Example results can be found at https://openbenchmarking.org/test/pts/quantlib + + This benchmark is derived from quantlibtestsuite.cpp. Please see the copyrights therein. */ @@ -124,10 +65,15 @@ #include #include -/* PAPI code -#include + +/* initialize PAPI on Linux + sudo sysctl -w kernel.perf_event_paranoid=0 + export PAPI_EVENTS="PAPI_TOT_INS,PAPI_FP_OPS,PAPI_FP_INS" + export PAPI_REPORT=1 */ +//#include + + /* Use BOOST_MSVC instead of _MSC_VER since some other vendors (Metrowerks, for example) also #define _MSC_VER @@ -296,7 +242,7 @@ namespace { Benchmark("EuropeanOption::FdEngines", [] { QuantLibTests::EuropeanOptionTests::testFdEngines().test_method(); }, 148.43), Benchmark("FdHestonTest::testFdmHestonAmerican", [] { QuantLibTests::FdHestonTests::testFdmHestonAmerican().test_method(); }, 234.21), Benchmark("HestonModel::DAXCalibration", [] { QuantLibTests::HestonModelTests::testDAXCalibration().test_method(); }, 555.19), - Benchmark("InterpolationTest::testSabrInterpolation", [] { QuantLibTests::InterpolationTests::testSabrInterpolation().test_method(); }, 2266.06), + Benchmark("InterpolationTest::testSabrInterpolation", [] { QuantLibTests::InterpolationTests::testSabrInterpolation().test_method(); }, 295.63), Benchmark("JumpDiffusion::Greeks", [] { QuantLibTests::JumpDiffusionTests::testGreeks().test_method(); }, 433.77), Benchmark("MarketModelCmsTest::testCmSwapsSwaptions", [] { QuantLibTests::MarketModelCmsTests::testMultiStepCmSwapsAndSwaptions().test_method(); }, 11497.73), Benchmark("MarketModelSmmTest::testMultiSmmSwaptions", [] { QuantLibTests::MarketModelSmmTests::testMultiStepCoterminalSwapsAndSwaptions().test_method(); }, 11244.95), @@ -306,28 +252,19 @@ namespace { Benchmark("ShortRateModel::Swaps", [] { QuantLibTests::ShortRateModelTests::testSwaps().test_method(); }, 454.73) }; - /* PAPI code - float real_time, proc_time, mflops; - long_long lflop, flop=0; - */ - class TimedBenchmark { public: - explicit TimedBenchmark(std::function f) : f_(std::move(f)) {} + TimedBenchmark(std::function f, const std::string name) + : f_(std::move(f)), name_(std::move(name)) {} void startMeasurement() const { - /* PAPI code - lflop = flop; - PAPI_flops(&real_time, &proc_time, &flop, &mflops); - */ + //QL_REQUIRE(PAPI_hl_region_begin(name_.c_str()) == PAPI_OK, + // "could not initialize PAPI"); } void stopMeasurement() const { - /* PAPI code - PAPI_flops(&real_time, &proc_time, &flop, &mflops); - printf("Real_time: %f Proc_time: %f Total mflop: %f\n", - real_time, proc_time, (flop-lflop)/1e6); - */ + //QL_REQUIRE(PAPI_hl_region_end(name_.c_str()) == PAPI_OK, + // "could not stop PAPI"); } double operator()() const { @@ -342,6 +279,7 @@ namespace { } private: std::function f_; + const std::string name_; }; void printResults( @@ -449,7 +387,7 @@ int main(int argc, char* argv[] ) { std::for_each(bm.begin(), bm.end(), [&runTimes](const Benchmark& iter) { runTimes.emplace_back( - iter, TimedBenchmark(iter.getTestCase())()); + iter, TimedBenchmark(iter.getTestCase(), iter.getName())()); }); printResults(nProc, runTimes); } @@ -514,7 +452,7 @@ int main(int argc, char* argv[] ) { mq.receive(&id, sizeof(unsigned), recvd_size, priority); while (id != terminateId) { - result_type a(id, TimedBenchmark(bm[id].getTestCase())()); + result_type a(id, TimedBenchmark(bm[id].getTestCase(), bm[id].getName())()); rq.send(&a, sizeof(result_type), 0); mq.receive(&id, sizeof(unsigned), recvd_size, priority); diff --git a/test-suite/riskneutraldensitycalculator.cpp b/test-suite/riskneutraldensitycalculator.cpp index c1d86d6f869..9d38c838f1a 100644 --- a/test-suite/riskneutraldensitycalculator.cpp +++ b/test-suite/riskneutraldensitycalculator.cpp @@ -588,7 +588,7 @@ BOOST_AUTO_TEST_CASE(testBlackScholesWithSkew, *precondition(if_speed(Fast))) { ext::make_shared( Handle(ext::make_shared(hestonProcess)), AnalyticHestonEngine::AndersenPiterbarg, - AnalyticHestonEngine::Integration::discreteTrapezoid(64))); + AnalyticHestonEngine::Integration::discreteTrapezoid(128))); const ext::shared_ptr timeGrid(new TimeGrid(maturity, 51)); @@ -692,7 +692,7 @@ BOOST_AUTO_TEST_CASE(testBlackScholesWithSkew, *precondition(if_speed(Fast))) { const Real calculatedLocalVol = std::exp(localVolCalc.invcdf(quantile, maturity)); - const Real localVolTol = 0.1; + const Real localVolTol = 0.2; if (std::fabs(expected - calculatedLocalVol) > localVolTol) { BOOST_FAIL("failed to match Heston and local Volatility invcdf" << "\n t: " << maturity diff --git a/test-suite/testsuite.vcxproj b/test-suite/testsuite.vcxproj index 64c6dbd624f..30da8dd6c93 100644 --- a/test-suite/testsuite.vcxproj +++ b/test-suite/testsuite.vcxproj @@ -829,4 +829,4 @@ - + \ No newline at end of file