Heston Model: Support for Angled Contour Shift Integrals (lballabio#1826)

Author: lballabio

QuantLib.vcxproj

@@ -1011,6 +1011,7 @@
     <ClInclude Include="ql\math\distributions\poissondistribution.hpp" />
     <ClInclude Include="ql\math\distributions\studenttdistribution.hpp" />
     <ClInclude Include="ql\math\errorfunction.hpp" />
+    <ClInclude Include="ql\math\expm1.hpp" />
     <ClInclude Include="ql\math\factorial.hpp" />
     <ClInclude Include="ql\math\fastfouriertransform.hpp" />
     <ClInclude Include="ql\math\functional.hpp" />
@@ -1019,6 +1020,7 @@
     <ClInclude Include="ql\math\integrals\all.hpp" />
     <ClInclude Include="ql\math\integrals\discreteintegrals.hpp" />
     <ClInclude Include="ql\math\integrals\exponentialintegrals.hpp" />
+    <ClInclude Include="ql\math\integrals\expsinhintegral.hpp" />
     <ClInclude Include="ql\math\integrals\filonintegral.hpp" />
     <ClInclude Include="ql\math\integrals\gaussianorthogonalpolynomial.hpp" />
     <ClInclude Include="ql\math\integrals\gaussianquadratures.hpp" />
@@ -2260,6 +2262,7 @@
     <ClCompile Include="ql\math\distributions\normaldistribution.cpp" />
     <ClCompile Include="ql\math\distributions\studenttdistribution.cpp" />
     <ClCompile Include="ql\math\errorfunction.cpp" />
+    <ClCompile Include="ql\math\expm1.cpp" />
     <ClCompile Include="ql\math\factorial.cpp" />
     <ClCompile Include="ql\math\incompletegamma.cpp" />
     <ClCompile Include="ql\math\integrals\discreteintegrals.cpp" />

QuantLib.vcxproj.filters

@@ -987,6 +987,9 @@
     <ClInclude Include="ql\math\errorfunction.hpp">
       <Filter>math</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\expm1.hpp">
+      <Filter>math</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\factorial.hpp">
       <Filter>math</Filter>
     </ClInclude>
@@ -1164,6 +1167,9 @@
     <ClInclude Include="ql\math\integrals\discreteintegrals.hpp">
       <Filter>math\integrals</Filter>
     </ClInclude>
+    <ClInclude Include="ql\math\integrals\expsinhintegral.hpp">
+      <Filter>math\integrals</Filter>
+    </ClInclude>
     <ClInclude Include="ql\math\integrals\filonintegral.hpp">
       <Filter>math\integrals</Filter>
     </ClInclude>
@@ -4850,6 +4856,9 @@
     <ClCompile Include="ql\math\errorfunction.cpp">
       <Filter>math</Filter>
     </ClCompile>
+    <ClCompile Include="ql\math\expm1.cpp">
+      <Filter>math</Filter>
+    </ClCompile>
     <ClCompile Include="ql\math\factorial.cpp">
       <Filter>math</Filter>
     </ClCompile>

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,13 +1427,15 @@ 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
     math/generallinearleastsquares.hpp
     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

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 \

ql/math/all.hpp

@@ -9,6 +9,7 @@
 #include <ql/math/bspline.hpp>
 #include <ql/math/comparison.hpp>
 #include <ql/math/errorfunction.hpp>
+#include <ql/math/expm1.hpp>
 #include <ql/math/factorial.hpp>
 #include <ql/math/fastfouriertransform.hpp>
 #include <ql/math/functional.hpp>

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
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ 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 <ql/math/expm1.hpp>
+#include <ql/math/functional.hpp>
+
+#include <cmath>
+
+namespace QuantLib {
+    std::complex<Real> expm1(const std::complex<Real>& 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<Real>(
+                exp_1*cos_1 + exp_1 + cos_1,
+                std::sin(b)*std::exp(a)
+            );
+        }
+        else {
+            return std::exp(z)-1.0;
+        }
+    }
+
+    std::complex<Real> log1p(const std::complex<Real>& z) {
+        const Real a = z.real(), b = z.imag();
+        if (std::abs(a) < 0.5 && std::abs(b) < 0.5) {
+            return std::complex<Real>(
+                0.5*std::log1p(a*a + 2*a + b*b),
+                std::arg(1.0 + z)
+            );
+        }
+        else {
+            return std::log(1.0+z);
+        }
+    }
+}

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
+ <quantlib-dev@lists.sf.net>. The license is also available online at
+ <http://quantlib.org/license.shtml>.
+
+ 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 <ql/types.hpp>
+#include <complex>
+
+namespace QuantLib {
+    std::complex<Real> expm1(const std::complex<Real>& z);
+    std::complex<Real> log1p(const std::complex<Real>& z);
+}
+#endif
+

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 \

ql/math/integrals/all.hpp

@@ -3,6 +3,7 @@
 
 #include <ql/math/integrals/discreteintegrals.hpp>
 #include <ql/math/integrals/exponentialintegrals.hpp>
+#include <ql/math/integrals/expsinhintegral.hpp>
 #include <ql/math/integrals/filonintegral.hpp>
 #include <ql/math/integrals/gausslaguerrecosinepolynomial.hpp>
 #include <ql/math/integrals/gausslobattointegral.hpp>

ql/math/integrals/exponentialintegrals.cpp

@@ -23,14 +23,12 @@
 
 #include <ql/errors.hpp>
 #include <ql/mathconstants.hpp>
+#include <ql/math/comparison.hpp>
 #include <ql/math/integrals/exponentialintegrals.hpp>
 
+#include <boost/math/special_functions/sign.hpp>
 #include <cmath>
 
-#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<Real> _Ei(
+            const std::complex<Real>& z, const std::complex<Real>& acc) {
 
-        std::complex<Real> E1(std::complex<Real> 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<Real>::has_infinity) {
+                return std::complex<Real>(-std::numeric_limits<Real>::infinity());
+            }
 
-            std::complex<Real> 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<Real> Ei(std::complex<Real> 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<Real> s(0.0), sn=z;
+            const auto match = [=](
+                const std::complex<Real>& z1, const std::complex<Real>& z2)
+                -> bool {
+                    const std::complex<Real> 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<Real>(std::exp(z)/z) + acc;
 
+            if (abs_z > 1.1*z_asym) {
+                std::complex<Real> ei = acc + std::complex<Real>(Real(0.0), sign(z.imag())*M_PI);
+                std::complex<Real> 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<Real> 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<Real>(0.0, sign(z.imag())*M_PI))
+                        - std::exp(z)/ (1.0 - z + ei);
+
+                QL_FAIL("series conversion issue for Ei(" << z << ")");
+            }
+
+            std::complex<Real> 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<Real> 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<Real>(r.real(), acc.imag());
+        }
+
+        std::complex<Real> Ei(const std::complex<Real>&z) {
+            return _Ei(z, std::complex<Real>(0.0));
+        }
+
+        std::complex<Real> E1(const std::complex<Real>& z) {
+            if (z.imag() < 0.0) {
+                return -_Ei(-z, std::complex<Real>(0.0, -M_PI));
+            }
+            else if (z.imag() > 0.0 || z.real() < 0.0) {
+                return -_Ei(-z, std::complex<Real>(0.0, M_PI));
+            }
+            else {
+                return -Ei(-z);
+            }
         }
 
         // Reference:
         // https://functions.wolfram.com/GammaBetaErf/ExpIntegralEi/introductions/ExpIntegrals/ShowAll.html
-        std::complex<Real> Si(std::complex<Real> z) {
-            const std::complex<Real> 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<Real> Si(const std::complex<Real>& z) {
+            if (std::abs(z) <= 0.2) {
+                std::complex<Real> 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<Real> i(0.0, 1.0);
+                return 0.5*i*(E1(-i*z) - E1(i*z)
+                        - std::complex<Real>(0.0, ((z.real() >= 0 && z.imag() >= 0)
+                                || (z.real() > 0 && z.imag() < 0) )? M_PI : -M_PI));
+            }
         }
 
-        std::complex<Real> Ci(std::complex<Real> z) {
+        std::complex<Real> Ci(const std::complex<Real>& z) {
             const std::complex<Real> 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<Real> 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;
         }
     }
 }