Higher order spatial operators for Finite Difference Methods

QuantLib.vcxproj

@@ -589,6 +589,8 @@
     <ClInclude Include="ql\methods\finitedifferences\operators\fdmlinearopcomposite.hpp" />
     <ClInclude Include="ql\methods\finitedifferences\operators\fdmlinearopiterator.hpp" />
     <ClInclude Include="ql\methods\finitedifferences\operators\fdmlinearoplayout.hpp" />
+    <ClInclude Include="ql\methods\finitedifferences\operators\numericaldifferentiation.hpp" />
+    <ClInclude Include="ql\methods\finitedifferences\operators\nthorderderivativeop.hpp" />
     <ClInclude Include="ql\methods\finitedifferences\operators\firstderivativeop.hpp" />
     <ClInclude Include="ql\methods\finitedifferences\operators\ninepointlinearop.hpp" />
     <ClInclude Include="ql\methods\finitedifferences\operators\secondderivativeop.hpp" />
@@ -1752,7 +1754,6 @@
     <ClInclude Include="ql\experimental\math\moorepenroseinverse.hpp" />
     <ClInclude Include="ql\experimental\math\multidimintegrator.hpp" />
     <ClInclude Include="ql\experimental\math\multidimquadrature.hpp" />
-    <ClInclude Include="ql\experimental\math\numericaldifferentiation.hpp" />
     <ClInclude Include="ql\experimental\math\piecewisefunction.hpp" />
     <ClInclude Include="ql\experimental\math\piecewiseintegral.hpp" />
     <ClInclude Include="ql\experimental\math\polarstudenttrng.hpp" />
@@ -1881,6 +1882,8 @@
     <ClCompile Include="ql\methods\finitedifferences\operators\fdmlinearoplayout.cpp" />
     <ClCompile Include="ql\methods\finitedifferences\operators\firstderivativeop.cpp" />
     <ClCompile Include="ql\methods\finitedifferences\operators\ninepointlinearop.cpp" />
+    <ClCompile Include="ql\methods\finitedifferences\operators\numericaldifferentiation.cpp" />
+    <ClCompile Include="ql\methods\finitedifferences\operators\nthorderderivativeop.cpp" />
     <ClCompile Include="ql\methods\finitedifferences\operators\secondderivativeop.cpp" />
     <ClCompile Include="ql\methods\finitedifferences\operators\secondordermixedderivativeop.cpp" />
     <ClCompile Include="ql\methods\finitedifferences\operators\triplebandlinearop.cpp" />
@@ -2654,7 +2657,6 @@
     <ClCompile Include="ql\experimental\math\gaussiancopulapolicy.cpp" />
     <ClCompile Include="ql\experimental\math\multidimintegrator.cpp" />
     <ClCompile Include="ql\experimental\math\multidimquadrature.cpp" />
-    <ClCompile Include="ql\experimental\math\numericaldifferentiation.cpp" />
     <ClCompile Include="ql\experimental\math\piecewiseintegral.cpp" />
     <ClCompile Include="ql\experimental\math\tcopulapolicy.cpp" />
     <ClCompile Include="ql\experimental\math\zigguratrng.cpp" />

QuantLib.vcxproj.filters

@@ -3606,8 +3606,11 @@
     <ClInclude Include="ql\experimental\math\multidimquadrature.hpp">
       <Filter>experimental\math</Filter>
     </ClInclude>
-    <ClInclude Include="ql\experimental\math\numericaldifferentiation.hpp">
-      <Filter>experimental\math</Filter>
+    <ClInclude Include="ql\methods\finitedifferences\operators\numericaldifferentiation.hpp">
+      <Filter>methods\finitedifferences\operators</Filter>
+    </ClInclude>
+    <ClInclude Include="ql\methods\finitedifferences\operators\nthorderderivativeop.hpp">
+      <Filter>methods\finitedifferences\operators</Filter>
     </ClInclude>
     <ClInclude Include="ql\experimental\math\piecewisefunction.hpp">
       <Filter>experimental\math</Filter>
@@ -6341,8 +6344,11 @@
     <ClCompile Include="ql\experimental\math\multidimquadrature.cpp">
       <Filter>experimental\math</Filter>
     </ClCompile>
-    <ClCompile Include="ql\experimental\math\numericaldifferentiation.cpp">
-      <Filter>experimental\math</Filter>
+    <ClCompile Include="ql\methods\finitedifferences\operators\numericaldifferentiation.cpp">
+      <Filter>methods\finitedifferences\operators</Filter>
+    </ClCompile>
+    <ClCompile Include="ql\methods\finitedifferences\operators\nthorderderivativeop.cpp">
+      <Filter>methods\finitedifferences\operators</Filter>
     </ClCompile>
     <ClCompile Include="ql\experimental\math\piecewiseintegral.cpp">
       <Filter>experimental\math</Filter>

ql/experimental/callablebonds/callablebond.cpp

@@ -23,7 +23,7 @@
 #include <ql/cashflows/cashflowvectors.hpp>
 #include <ql/termstructures/yield/zerospreadedtermstructure.hpp>
 #include <ql/math/solvers1d/brent.hpp>
-#include <ql/experimental/math/numericaldifferentiation.hpp>
+
 namespace QuantLib {
 
     CallableBond::CallableBond(Natural settlementDays,

ql/experimental/math/Makefile.am

@@ -21,7 +21,6 @@ this_include_HEADERS = \
     moorepenroseinverse.hpp \
     multidimintegrator.hpp \
     multidimquadrature.hpp \
-    numericaldifferentiation.hpp \
     particleswarmoptimization.hpp \
     piecewisefunction.hpp \
     piecewiseintegral.hpp \
@@ -37,7 +36,6 @@ cpp_files = \
     gaussiannoncentralchisquaredpolynomial.cpp \
     multidimintegrator.cpp \
     multidimquadrature.cpp \
-    numericaldifferentiation.cpp \
     particleswarmoptimization.cpp \
     piecewiseintegral.cpp \
     tcopulapolicy.cpp \

ql/experimental/math/all.hpp

@@ -18,7 +18,6 @@
 #include <ql/experimental/math/moorepenroseinverse.hpp>
 #include <ql/experimental/math/multidimintegrator.hpp>
 #include <ql/experimental/math/multidimquadrature.hpp>
-#include <ql/experimental/math/numericaldifferentiation.hpp>
 #include <ql/experimental/math/particleswarmoptimization.hpp>
 #include <ql/experimental/math/piecewisefunction.hpp>
 #include <ql/experimental/math/piecewiseintegral.hpp>

ql/experimental/shortrate/generalizedhullwhite.cpp

@@ -129,7 +129,7 @@ namespace QuantLib {
         Real lnEt = integrateMeanReversion(speed_,0,t);
         Real Et = exp(lnEt);
         Real B = 0;
-        Size N = std::min<Size>((T-t)*365, 2000);
+        Size N = std::min<Size>(Size((T-t)*365), 2000);
         if (N==0) N=1;
         Real dt = 0.5*(T-t)/N;
         Real a,b,c,_t,total=0;
@@ -152,7 +152,7 @@ namespace QuantLib {
         // Gurrieri et al, equation (37)
         Real lnEt = integrateMeanReversion(speed_,0,t);
         Real V = 0,Eu;
-        Size N = std::min<Size>((T-t)*365, 2000);
+        Size N = std::min<Size>(Size((T-t)*365), 2000);
         if (N==0) N=1;
         Real dt = 0.5*(T-t)/N;
         Real a,b,c,_t,lnE=lnEt;

ql/math/matrixutilities/sparsematrix.hpp

@@ -46,7 +46,7 @@ FOR A PARTICULAR PURPOSE.  See the license for more details.
 #pragma GCC diagnostic ignored "-Wunused-local-typedefs"
 #endif
 
-#if BOOST_VERSION > 106300
+#if BOOST_VERSION == 106400
 #include <boost/serialization/array_wrapper.hpp>
 #endif
 
@@ -70,7 +70,7 @@ namespace QuantLib {
         SparseMatrixReference;
 
     inline Disposable<Array> prod(const SparseMatrix& A, const Array& x) {
-        Array b(x.size());
+        Array b(x.size(), 0.0);
 
         for (Size i=0; i < A.filled1()-1; ++i) {
             const Size begin = A.index1_data()[i];

ql/methods/finitedifferences/meshers/fdmblackscholesmesher.cpp

@@ -53,7 +53,7 @@ namespace QuantLib {
                     dividendSchedule[i]->amount()
                 ) );
 
-        const Size intermediateTimeSteps = std::max<Size>(2, 24.0*maturity);
+        const Size intermediateTimeSteps = std::max<Size>(2, Size(24.0*maturity));
         for (Size i=0; i < intermediateTimeSteps; ++i)
             intermediateSteps.push_back(
                 std::make_pair((i+1)*(maturity/intermediateTimeSteps), 0.0));

ql/methods/finitedifferences/operators/Makefile.am

@@ -18,6 +18,8 @@ this_include_HEADERS = \
 	fdmlinearoplayout.hpp \
 	firstderivativeop.hpp \
 	ninepointlinearop.hpp \
+	nthorderderivativeop.hpp \
+	numericaldifferentiation.hpp \
 	secondderivativeop.hpp \
 	secondordermixedderivativeop.hpp \
 	triplebandlinearop.hpp
@@ -34,6 +36,8 @@ cpp_files = \
 	fdmornsteinuhlenbeckop.cpp \
 	firstderivativeop.cpp \
 	ninepointlinearop.cpp \
+	nthorderderivativeop.cpp \
+	numericaldifferentiation.cpp \
 	secondderivativeop.cpp \
 	secondordermixedderivativeop.cpp \
 	triplebandlinearop.cpp

ql/methods/finitedifferences/operators/all.hpp

@@ -15,6 +15,8 @@
 #include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
 #include <ql/methods/finitedifferences/operators/firstderivativeop.hpp>
 #include <ql/methods/finitedifferences/operators/ninepointlinearop.hpp>
+#include <ql/methods/finitedifferences/operators/nthorderderivativeop.hpp>
+#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
 #include <ql/methods/finitedifferences/operators/secondderivativeop.hpp>
 #include <ql/methods/finitedifferences/operators/secondordermixedderivativeop.hpp>
 #include <ql/methods/finitedifferences/operators/triplebandlinearop.hpp>

ql/methods/finitedifferences/operators/firstderivativeop.hpp

@@ -33,8 +33,9 @@ namespace QuantLib {
     class FirstDerivativeOp : public TripleBandLinearOp {
       public:
         FirstDerivativeOp(Size direction,
-                           const ext::shared_ptr<FdmMesher>& mesher);
+                          const ext::shared_ptr<FdmMesher>& mesher);
     };
+
 }
 
 #endif

ql/methods/finitedifferences/operators/nthorderderivativeop.cpp

@@ -0,0 +1,98 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2018 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 NthOrderDerivativeOp.cpp
+    \brief n-th order derivative linear operator
+*/
+
+#include <ql/qldefines.hpp>
+#ifndef QL_NO_UBLAS_SUPPORT
+
+#include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
+#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
+#include <ql/methods/finitedifferences/operators/nthorderderivativeop.hpp>
+
+#include <set>
+
+namespace QuantLib {
+
+    NthOrderDerivativeOp::NthOrderDerivativeOp(
+        Size direction, Size order, Integer nPoints,
+        const ext::shared_ptr<FdmMesher>& mesher)
+    : m_(mesher->layout()->size(), mesher->layout()->size()) {
+
+        const Integer hPoints = nPoints/2;
+        const bool isEven = (nPoints == 2*hPoints);
+
+        const ext::shared_ptr<FdmLinearOpLayout> layout = mesher->layout();
+        const FdmLinearOpIterator endIter = layout->end();
+
+        Array xValues = mesher->locations(direction);
+        std::set<Real> tmp(xValues.begin(), xValues.end());
+        xValues = Array(tmp.begin(), tmp.end()); //unique vector
+
+        const Integer nx(layout->dim()[direction]);
+
+        QL_REQUIRE(Integer(xValues.size()) == nx,
+            "inconsistent set of grid values in direction " << direction);
+
+        QL_REQUIRE(nPoints > 1 && Integer(nPoints) <= nx,
+             "inconsistent number of points");
+
+        Array xOffsets(nPoints);
+        const boost::function<Real(Real)> emptyFct;
+
+        for (FdmLinearOpIterator iter = layout->begin(); iter!=endIter; ++iter) {
+            const Integer ix = Integer(iter.coordinates()[direction]);
+            const Integer offset = std::max(0, hPoints - ix)
+                - std::max(0, hPoints - (nx-((isEven)? 0 : 1) - ix));
+
+            const Integer ilx = ix - hPoints + offset;
+
+            for (Integer j=0; j < nPoints; ++j) {
+                const Integer idx = ilx + j;
+                xOffsets[j] = xValues[idx] - xValues[ix];
+            }
+
+            const Array weights =
+                NumericalDifferentiation(emptyFct, order, xOffsets).weights();
+
+            const Size i = iter.index();
+            for (Integer j=0; j < nPoints; ++j) {
+                const Size k = layout->neighbourhood(iter, direction, ilx - ix + j);
+
+                m_(i, k) = weights[j];
+            }
+        }
+    }
+
+    Disposable<NthOrderDerivativeOp::array_type>
+    NthOrderDerivativeOp::apply(const array_type& r) const {
+        return prod(m_, r);
+    }
+
+
+    Disposable<SparseMatrix> NthOrderDerivativeOp::toMatrix() const {
+        SparseMatrix tmp(m_);
+        return tmp;
+    }
+
+}
+
+#endif

ql/methods/finitedifferences/operators/nthorderderivativeop.hpp

@@ -0,0 +1,50 @@
+/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
+
+/*
+ Copyright (C) 2018 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 nthorderderivativeop.hpp
+    \brief n-th order derivative linear operator
+*/
+
+#ifndef quantlib_nth_order_derivative_op_hpp
+#define quantlib_nth_order_derivative_op_hpp
+
+#include <ql/qldefines.hpp>
+#ifndef QL_NO_UBLAS_SUPPORT
+
+#include <ql/methods/finitedifferences/meshers/fdmmesher.hpp>
+#include <ql/methods/finitedifferences/operators/fdmlinearop.hpp>
+
+namespace QuantLib {
+
+    class NthOrderDerivativeOp : public FdmLinearOp {
+      public:
+        NthOrderDerivativeOp(
+            Size direction, Size order, Integer nPoints,
+            const ext::shared_ptr<FdmMesher>& mesher);
+
+        Disposable<array_type> apply(const array_type& r) const;
+        Disposable<SparseMatrix> toMatrix() const;
+
+      private:
+        SparseMatrix m_;
+    };
+}
+
+#endif
+#endif

ql/experimental/math/numericaldifferentiation.cpp → ql/methods/finitedifferences/operators/numericaldifferentiation.cpp

@@ -19,7 +19,7 @@
 
 /*! \file numericaldifferentiation.cpp */
 
-#include <ql/experimental/math/numericaldifferentiation.hpp>
+#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
 
 #ifndef QL_EXTRA_SAFETY_CHECKS
 #define BOOST_DISABLE_ASSERTS 1

test-suite/Makefile.am

@@ -97,6 +97,7 @@ QL_TESTS = \
 	money.hpp money.cpp \
 	noarbsabr.hpp noarbsabr.cpp \
 	normalclvmodel.hpp normalclvmodel.cpp \
+	nthorderderivativeop.hpp nthorderderivativeop.cpp \
 	nthtodefault.hpp nthtodefault.cpp \
 	numericaldifferentiation.hpp numericaldifferentiation.cpp \
 	observable.hpp observable.cpp \

test-suite/fdmlinearop.cpp

@@ -29,7 +29,6 @@
 #include <ql/processes/hullwhiteprocess.hpp>
 #include <ql/processes/blackscholesprocess.hpp>
 #include <ql/processes/hybridhestonhullwhiteprocess.hpp>
-#include <ql/experimental/math/numericaldifferentiation.hpp>
 #include <ql/math/interpolations/bilinearinterpolation.hpp>
 #include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
 #include <ql/math/interpolations/cubicinterpolation.hpp>
@@ -53,6 +52,7 @@
 #include <ql/methods/finitedifferences/operators/fdmblackscholesop.hpp>
 #include <ql/methods/finitedifferences/utilities/fdmmesherintegral.hpp>
 #include <ql/methods/finitedifferences/utilities/fdminnervaluecalculator.hpp>
+#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearop.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearoplayout.hpp>
 #include <ql/methods/finitedifferences/operators/fdmlinearopcomposite.hpp>

test-suite/hestonmodel.cpp

@@ -53,8 +53,8 @@
 #include <ql/math/optimization/differentialevolution.hpp>
 #include <ql/time/period.hpp>
 #include <ql/quotes/simplequote.hpp>
-#include <ql/experimental/math/numericaldifferentiation.hpp>
 #include <ql/experimental/exoticoptions/analyticpdfhestonengine.hpp>
+#include <ql/methods/finitedifferences/operators/numericaldifferentiation.hpp>
 
 #include <boost/math/special_functions/fpclassify.hpp>