Merge pull request #812.

added Heston engine based on exponentially-fitted Laguerre quadrature rule

QuantLib.vcxproj

@@ -1549,6 +1549,7 @@
     <ClInclude Include="ql\pricingengines\vanilla\bjerksundstenslandengine.hpp" />
     <ClInclude Include="ql\pricingengines\vanilla\coshestonengine.hpp" />
     <ClInclude Include="ql\pricingengines\vanilla\discretizedvanillaoption.hpp" />
+    <ClInclude Include="ql\pricingengines\vanilla\exponentialfittinghestonengine.hpp" />
     <ClInclude Include="ql\pricingengines\vanilla\fdamericanengine.hpp" />
     <ClInclude Include="ql\pricingengines\vanilla\fdbatesvanillaengine.hpp" />
     <ClInclude Include="ql\pricingengines\vanilla\fdbermudanengine.hpp" />
@@ -2555,6 +2556,7 @@
     <ClCompile Include="ql\pricingengines\vanilla\bjerksundstenslandengine.cpp" />
     <ClCompile Include="ql\pricingengines\vanilla\coshestonengine.cpp" />
     <ClCompile Include="ql\pricingengines\vanilla\discretizedvanillaoption.cpp" />
+    <ClCompile Include="ql\pricingengines\vanilla\exponentialfittinghestonengine.cpp" />
     <ClCompile Include="ql\pricingengines\vanilla\fdbatesvanillaengine.cpp" />
     <ClCompile Include="ql\pricingengines\vanilla\fdblackscholesvanillaengine.cpp" />
     <ClCompile Include="ql\pricingengines\vanilla\fdcevvanillaengine.cpp" />

QuantLib.vcxproj.filters

@@ -4309,6 +4309,9 @@
     <ClInclude Include="ql\math\integrals\momentbasedgaussianpolynomial.hpp">
       <Filter>math\integrals</Filter>
     </ClInclude>
+    <ClInclude Include="ql\pricingengines\vanilla\exponentialfittinghestonengine.hpp">
+      <Filter>pricingengines\vanilla</Filter>
+    </ClInclude>
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="ql\methods\montecarlo\brownianbridge.cpp">
@@ -6954,5 +6957,8 @@
     <ClCompile Include="ql\methods\finitedifferences\schemes\cranknicolsonscheme.cpp">
       <Filter>methods\finitedifferences\schemes</Filter>
     </ClCompile>
+    <ClCompile Include="ql\pricingengines\vanilla\exponentialfittinghestonengine.cpp">
+      <Filter>pricingengines\vanilla</Filter>
+    </ClCompile>
   </ItemGroup>
 </Project>

ql/CMakeLists.txt

@@ -719,6 +719,7 @@ set(QuantLib_SRC
     pricingengines/vanilla/bjerksundstenslandengine.cpp
     pricingengines/vanilla/coshestonengine.cpp
     pricingengines/vanilla/discretizedvanillaoption.cpp
+    pricingengines/vanilla/exponentialfittinghestonengine.cpp
     pricingengines/vanilla/fdbatesvanillaengine.cpp
     pricingengines/vanilla/fdblackscholesvanillaengine.cpp
     pricingengines/vanilla/fdcevvanillaengine.cpp
@@ -1982,6 +1983,7 @@ set(QuantLib_HDR
     pricingengines/vanilla/bjerksundstenslandengine.hpp
     pricingengines/vanilla/coshestonengine.hpp
     pricingengines/vanilla/discretizedvanillaoption.hpp
+    pricingengines/vanilla/exponentialfittinghestonengine.hpp
     pricingengines/vanilla/fdamericanengine.hpp
     pricingengines/vanilla/fdbatesvanillaengine.hpp
     pricingengines/vanilla/fdbermudanengine.hpp

ql/pricingengines/vanilla/Makefile.am

@@ -21,6 +21,7 @@ this_include_HEADERS = \
     bjerksundstenslandengine.hpp \
     coshestonengine.hpp \
     discretizedvanillaoption.hpp \
+    exponentialfittinghestonengine.hpp \
     hestonexpansionengine.hpp \
     integralengine.hpp \
     jumpdiffusionengine.hpp \
@@ -69,6 +70,7 @@ cpp_files = \
     bjerksundstenslandengine.cpp \
     coshestonengine.cpp \
     discretizedvanillaoption.cpp \
+    exponentialfittinghestonengine.cpp \
     hestonexpansionengine.cpp \
     integralengine.cpp \
     jumpdiffusionengine.cpp \

ql/pricingengines/vanilla/all.hpp

@@ -18,6 +18,7 @@
 #include <ql/pricingengines/vanilla/bjerksundstenslandengine.hpp>
 #include <ql/pricingengines/vanilla/coshestonengine.hpp>
 #include <ql/pricingengines/vanilla/discretizedvanillaoption.hpp>
+#include <ql/pricingengines/vanilla/exponentialfittinghestonengine.hpp>
 #include <ql/pricingengines/vanilla/hestonexpansionengine.hpp>
 #include <ql/pricingengines/vanilla/integralengine.hpp>
 #include <ql/pricingengines/vanilla/jumpdiffusionengine.hpp>

ql/pricingengines/vanilla/analytichestonengine.cpp

@@ -404,7 +404,7 @@ namespace QuantLib {
         const AnalyticHestonEngine* const enginePtr_;
     };
 
-    std::complex<Real> AnalyticHestonEngine::lnChF(
+    std::complex<Real> AnalyticHestonEngine::chF(
         const std::complex<Real>& z, Time t) const {
 
         const Real kappa = model_->kappa();
@@ -415,22 +415,50 @@ namespace QuantLib {
 
         const Real sigma2 = sigma*sigma;
 
-        const std::complex<Real> g
-            = kappa + rho*sigma*std::complex<Real>(z.imag(), -z.real());
+        if (sigma > 1e-4) {
+            const std::complex<Real> g
+                = kappa + rho*sigma*std::complex<Real>(z.imag(), -z.real());
 
-        const std::complex<Real> D = std::sqrt(
-            g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2);
+            const std::complex<Real> D = std::sqrt(
+                g*g + (z*z + std::complex<Real>(-z.imag(), z.real()))*sigma2);
 
-        const std::complex<Real> G = (g-D)/(g+D);
+            const std::complex<Real> G = (g-D)/(g+D);
 
-        return 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)));
+            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))));
+        }
+        else {
+            const Real kt = kappa*t;
+            const Real ekt = std::exp(kt);
+            const Real e2kt = std::exp(2*kt);
+            const Real rho2 = rho*rho;
+            const std::complex<Real> zpi = z + std::complex<Real>(0.0, 1.0);
+
+            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<Real>(-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*square<Real>()(2*theta + kt*theta - v0 -
+                    kt*v0 + ekt*((-2 + kt)*theta + v0))
+                  *z*z*zpi + 2*kappa*v0*(-zpi
+                    + e2kt*(zpi + 4*rho2*z) - 2*ekt*(2*rho2*z
+                    + kt*(zpi + rho2*(2 + kt)*z))) + kappa*theta*(zpi + e2kt
+                *(-5.0*zpi - 24*rho2*z+ 2*kt*(zpi + 4*rho2*z)) +
+                4*ekt*(zpi + 6*rho2*z + kt*(zpi + rho2*(4 + kt)*z)))))
+                /(16.*square<Real>()(square<Real>()(kappa)))*sigma2;
+        }
     }
 
-    std::complex<Real> AnalyticHestonEngine::chF(
+    std::complex<Real> AnalyticHestonEngine::lnChF(
         const std::complex<Real>& z, Time T) const {
-        return std::exp(lnChF(z, T));
+        return std::log(chF(z, T));
     }
 
     AnalyticHestonEngine::AnalyticHestonEngine(