C++：实现量化相关的Interpolation插值测试实例

#include "interpolations.hpp"
#include "utilities.hpp"
#include <ql/experimental/volatility/noarbsabrinterpolation.hpp>
#include <ql/math/bspline.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/integrals/simpsonintegral.hpp>
#include <ql/math/interpolations/backwardflatinterpolation.hpp>
#include <ql/math/interpolations/bicubicsplineinterpolation.hpp>
#include <ql/math/interpolations/chebyshevinterpolation.hpp>
#include <ql/math/interpolations/cubicinterpolation.hpp>
#include <ql/math/interpolations/forwardflatinterpolation.hpp>
#include <ql/math/interpolations/kernelinterpolation.hpp>
#include <ql/math/interpolations/kernelinterpolation2d.hpp>
#include <ql/math/interpolations/lagrangeinterpolation.hpp>
#include <ql/math/interpolations/linearinterpolation.hpp>
#include <ql/math/interpolations/multicubicspline.hpp>
#include <ql/math/interpolations/sabrinterpolation.hpp>
#include <ql/math/kernelfunctions.hpp>
#include <ql/math/optimization/levenbergmarquardt.hpp>
#include <ql/math/randomnumbers/sobolrsg.hpp>
#include <ql/math/richardsonextrapolation.hpp>
#include <ql/tuple.hpp>
#include <ql/utilities/dataformatters.hpp>
#include <ql/utilities/null.hpp>
#include <cmath>
#include <utility>

using namespace QuantLib;
using namespace boost::unit_test_framework;

namespace {
   

    std::vector<Real> xRange(Real start, Real finish, Size points) {
   
        std::vector<Real> x(points);
        Real dx = (finish-start)/(points-1);
        for (Size i=0; i<points-1; i++)
            x[i] = start+i*dx;
        x[points-1] = finish;
        return x;
    }

    std::vector<Real> gaussian(const std::vector<Real>& x) {
   
        std::vector<Real> y(x.size());
        for (Size i=0; i<x.size(); i++)
            y[i] = std::exp(-x[i]*x[i]);
        return y;
    }

    std::vector<Real> parabolic(const std::vector<Real>& x) {
   
        std::vector<Real> y(x.size());
        for (Size i=0; i<x.size(); i++)
            y[i] = -x[i]*x[i];
        return y;
    }

    template <class I, class J>
    void checkValues(const char* type,
                     const CubicInterpolation& cubic,
                     I xBegin, I xEnd, J yBegin) {
   
        Real tolerance = 2.0e-15;
        while (xBegin != xEnd) {
   
            Real interpolated = cubic(*xBegin);
            if (std::fabs(interpolated-*yBegin) > tolerance) {
   
                BOOST_ERROR(type << " interpolation failed at x = " << *xBegin
                            << std::scientific
                            << "\n    interpolated value: " << interpolated
                            << "\n    expected value:     " << *yBegin
                            << "\n    error:              "
                            << std::fabs(interpolated-*yBegin));
            }
            ++xBegin; ++yBegin;
        }
    }

    void check1stDerivativeValue(const char* type,
                                 const CubicInterpolation& cubic,
                                 Real x,
                                 Real value) {
   
        Real tolerance = 1.0e-14;
        Real interpolated = cubic.derivative(x);
        Real error = std::fabs(interpolated-value);
        if (error > tolerance) {
   
            BOOST_ERROR(type << " interpolation first derivative failure\n"
                        << "at x = " << x
                        << "\n    interpolated value: " << interpolated
                        << "\n    expected value:     " << value
                        << std::scientific
                        << "\n    error:              " << error);
        }
    }

    void check2ndDerivativeValue(const char* type,
                                 const CubicInterpolation& cubic,
                                 Real x,
                                 Real value) {
   
        Real tolerance = 1.0e-13;
        Real interpolated = cubic.secondDerivative(x);
        Real error = std::fabs(interpolated-value);
        if (error > tolerance) {
   
            BOOST_ERROR(type << " interpolation second derivative failure\n"
                        << "at x = " << x
                        << "\n    interpolated value: " << interpolated
                        << "\n    expected value:     " << value
                        << std::scientific
                        << "\n    error:              " << error);
        }
    }

    void checkNotAKnotCondition(const char* type,
                                const CubicInterpolation& cubic) {
   
        Real tolerance = 1.0e-14;
        const std::vector<Real>& c = cubic.cCoefficients();
        if (std::fabs(c[0]-c[1]) > tolerance) {
   
            BOOST_ERROR(type << " interpolation failure"
                        << "\n    cubic coefficient of the first"
                        << " polinomial is " << c[0]
                        << "\n    cubic coefficient of the second"
                        << " polinomial is " << c[1]);
        }
        Size n = c.size();
        if (std::fabs(c[n-2]-c[n-1]) > tolerance) {
   
            BOOST_ERROR(type << " interpolation failure"
                        << "\n    cubic coefficient of the 2nd to last"
                        << " polinomial is " << c[n-2]
                        << "\n    cubic coefficient of the last"
                        << " polinomial is " << c[n-1]);
        }
    }

    void checkSymmetry(const char* type,
                       const CubicInterpolation& cubic,
                       Real xMin) {
   
        Real tolerance = 1.0e-15;
        for (Real x = xMin; x < 0.0; x += 0.1) {
   
            Real y1 = cubic(x), y2 = cubic(-x);
            if (std::fabs(y1-y2) > tolerance) {
   
                BOOST_ERROR(type << " interpolation not symmetric"
                            << "\n    x = " << x
                            << "\n    g(x)  = " << y1
                            << "\n    g(-x) = " << y2
                            << "\n    error:  " << std::fabs(y1-y2));
            }
        }
    }

    template <class F>
    class errorFunction {
   
      public:
        errorFunction(F f) : f_(std::move(f)) {
   }
        Real operator()(Real x) const {
   
            Real temp = f_(x)-std::exp(-x*x);
            return temp*temp;
        }
      private:
        F f_;
    };

    template <class F>
    errorFunction<F> make_error_function(const F& f) {
   
        return errorFunction<F>(f);
    }

    Real multif(Real s, Real t, Real u, Real v, Real w) {
   
        return std::sqrt(s * std::sinh(std::log(t)) +
                         std::exp(std::sin(u) * std::sin(3 * v)) +
                         std::sinh(std::log(v * w)));
    }

    Real epanechnikovKernel(Real u){
   

        if(std::fabs(u)<=1){
   
            return (3.0/4.0)*(1-u*u);
        }else{
   
            return 0.0;
        }
    }

}


/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
   SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
   http://math.lanl.gov/~mac/papers/numerics/H83.pdf
*/
void InterpolationTest::testSplineErrorOnGaussianValues() {
   

    BOOST_TEST_MESSAGE("Testing spline approximation on Gaussian data sets...");

    Size points[]                = {
         5,      9,     17,     33 };

    // complete spline data from the original 1983 Hyman paper
    Real tabulatedErrors[]     = {
    3.5e-2, 2.0e-3, 4.0e-5, 1.8e-6 };
    Real toleranceOnTabErr[]   = {
    0.1e-2, 0.1e-3, 0.1e-5, 0.1e-6 };

    // (complete) MC spline data from the original 1983 Hyman paper
    // NB: with the improved Hyman filter from the Dougherty, Edelman, and
    //     Hyman 1989 paper the n=17 nonmonotonicity is not filtered anymore
    //     so the error agrees with the non MC method.
    Real tabulatedMCErrors[]   = {
    1.7e-2, 2.0e-3, 4.0e-5, 1.8e-6 };
    Real toleranceOnTabMCErr[] = {
    0.1e-2, 0.1e-3, 0.1e-5, 0.1e-6 };

    SimpsonIntegral integral(1e-12, 10000);
    std::vector<Real> x, y;

    // still unexplained scale factor needed to obtain the numerical
    // results from the paper
    Real scaleFactor = 1.9;

    for (Size i=0; i<LENGTH(points); i++) {
   
        Size n = points[i];
        std::vector<Real> x = xRange(-1.7, 1.9, n);
        std::vector<Real> y = gaussian(x);

        // Not-a-knot
        CubicInterpolation f(x.begin(), x.end(), y.begin(),
                             CubicInterpolation::Spline, false,
                             CubicInterpolation::NotAKnot, Null<Real>(),
                             CubicInterpolation::NotAKnot, Null<Real>());
        f.update();
        Real result = std::sqrt(integral(make_error_function(f), -1.7, 1.9));
        result /= scaleFactor;
        if (std::fabs(result-tabulatedErrors[i]) > toleranceOnTabErr[i])
            BOOST_ERROR("Not-a-knot spline interpolation "
                        << "\n    sample points:      " << n
                        << "\n    norm of difference: " << result
                        << "\n    it should be:       " << tabulatedErrors[i]);

        // MC not-a-knot
        f = CubicInterpolation(x.begin(), x.end(), y.begin(),
                               CubicInterpolation::Spline, true,
                               CubicInterpolation::NotAKnot, Null<Real>(),
                               CubicInterpolation::NotAKnot, Null<Real>());
        f.update();
        result = std::sqrt(integral(make_error_function(f), -1.7, 1.9));
        result /= scaleFactor;
        if (std::fabs(result-tabulatedMCErrors[i]) > toleranceOnTabMCErr[i])
            BOOST_ERROR("MC Not-a-knot spline interpolation "
                        << "\n    sample points:      " << n
                        << "\n    norm of difference: " << result
                        << "\n    it should be:       "
                        << tabulatedMCErrors[i]);
    }

}

/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
   SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
   http://math.lanl.gov/~mac/papers/numerics/H83.pdf
*/
void InterpolationTest::testSplineOnGaussianValues() {
   

    BOOST_TEST_MESSAGE("Testing spline interpolation on a Gaussian data set...");

    Real interpolated, interpolated2;
    Size n = 5;

    std::vector<Real> x(n), y(n);
    Real x1_bad=-1.7, x2_bad=1.7;

    for (Real start = -1.9, j=0; j<2; start+=0.2, j++) {
   
        x = xRange(start, start+3.6, n);
        y = gaussian(x);

        // Not-a-knot spline
        CubicInterpolation f(x.begin(), x.end(), y.begin(),
                             CubicInterpolation::Spline, false,
                             CubicInterpolation::NotAKnot, Null<Real>(),
                             CubicInterpolation::NotAKnot, Null<Real>());
        f.update();
        checkValues("Not-a-knot spline", f,
                    x.begin(), x.end(), y.begin());
        checkNotAKnotCondition("Not-a-knot spline", f);
        // bad performance
        interpolated = f(x1_bad);
        interpolated2= f(x2_bad);
        if (interpolated>0.0 && interpolated2>0.0 ) {
   
            BOOST_ERROR("Not-a-knot spline interpolation "
                        << "bad performance unverified"
                        << "\nat x = " << x1_bad
                        << " interpolated value: " << interpolated
                        << "\nat x = " << x2_bad
                        << " interpolated value: " << interpolated
                        << "\n at least one of them was expected to be < 0.0");
        }

        // MC not-a-knot spline
        f = CubicInterpolation(x.begin(), x.end(), y.begin(),
                               CubicInterpolation::Spline, true,
                               CubicInterpolation::NotAKnot, Null<Real>(),
                               CubicInterpolation::NotAKnot, Null<Real>());
        f.update();
        checkValues("MC not-a-knot spline", f,
                    x.begin(), x.end(), y.begin());
        // good performance
        interpolated = f(x1_bad);
        if (interpolated<0.0) {
   
            BOOST_ERROR("MC not-a-knot spline interpolation "
                        << "good performance unverified\n"
                        << "at x = " << x1_bad
                        << "\ninterpolated value: " << interpolated
                        << "\nexpected value > 0.0");
        }
        interpolated = f(x2_bad);
        if (interpolated<0.0) {
   
            BOOST_ERROR("MC not-a-knot spline interpolation "
                        << "good performance unverified\n"
                        << "at x = " << x2_bad
                        << "\ninterpolated value: " << interpolated
                        << "\nexpected value > 0.0");
        }
    }
}


/* See J. M. Hyman, "Accurate monotonicity preserving cubic interpolation"
   SIAM J. of Scientific and Statistical Computing, v. 4, 1983, pp. 645-654.
   http://math.lanl.gov/~mac/papers/numerics/H83.pdf
*/
void InterpolationTest::testSplineOnRPN15AValues() {
   

    BOOST_TEST_MESSAGE("Testing spline interpolation on RPN15A data set...");

    const Real RPN15A_x[] = {
   
        7.99,       8.09,       8.19,      8.7,
        9.2,     10.0,     12.0,     15.0,     20.0
    };
    const Real RPN15A_y[] = {
   
        0.0, 2.76429e-5, 4.37498e-5, 0.169183,
        0.469428, 0.943740, 0.998636, 0.999919, 0.999994
    };

    Real interpolated;

    // Natural spline
    CubicInterpolation f = CubicInterpolation(
                                    std::begin(RPN15A_x), std::end(RPN15A_x),
                                    std::begin(RPN15A_y),
                                    CubicInterpolation::Spline, false,
                                    CubicInterpolation::SecondDerivative, 0.0,
                                    CubicInterpolation::SecondDerivative, 0.0);
    f.update();
    checkValues("Natural spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    check2ndDerivativeValue("Natural spline", f,
                            *std::begin(RPN15A_x), 0.0);
    check2ndDerivativeValue("Natural spline", f,
                            *(std::end(RPN15A_x)-1), 0.0);
    // poor performance
    Real x_bad = 11.0;
    interpolated = f(x_bad);
    if (interpolated<1.0) {
   
        BOOST_ERROR("Natural spline interpolation "
                    << "poor performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value > 1.0");
    }


    // Clamped spline
    f = CubicInterpolation(std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y),
                           CubicInterpolation::Spline, false,
                           CubicInterpolation::FirstDerivative, 0.0,
                           CubicInterpolation::FirstDerivative, 0.0);
    f.update();
    checkValues("Clamped spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    check1stDerivativeValue("Clamped spline", f,
                            *std::begin(RPN15A_x), 0.0);
    check1stDerivativeValue("Clamped spline", f,
                            *(std::end(RPN15A_x)-1), 0.0);
    // poor performance
    interpolated = f(x_bad);
    if (interpolated<1.0) {
   
        BOOST_ERROR("Clamped spline interpolation "
                    << "poor performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value > 1.0");
    }


    // Not-a-knot spline
    f = CubicInterpolation(std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y),
                           CubicInterpolation::Spline, false,
                           CubicInterpolation::NotAKnot, Null<Real>(),
                           CubicInterpolation::NotAKnot, Null<Real>());
    f.update();
    checkValues("Not-a-knot spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    checkNotAKnotCondition("Not-a-knot spline", f);
    // poor performance
    interpolated = f(x_bad);
    if (interpolated<1.0) {
   
        BOOST_ERROR("Not-a-knot spline interpolation "
                    << "poor performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value > 1.0");
    }


    // MC natural spline values
    f = CubicInterpolation(std::begin(RPN15A_x), std::end(RPN15A_x),
                           std::begin(RPN15A_y),
                           CubicInterpolation::Spline, true,
                           CubicInterpolation::SecondDerivative, 0.0,
                           CubicInterpolation::SecondDerivative, 0.0);
    f.update();
    checkValues("MC natural spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    // good performance
    interpolated = f(x_bad);
    if (interpolated>1.0) {
   
        BOOST_ERROR("MC natural spline interpolation "
                    << "good performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value < 1.0");
    }


    // MC clamped spline values
    f = CubicInterpolation(std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y),
                           CubicInterpolation::Spline, true,
                           CubicInterpolation::FirstDerivative, 0.0,
                           CubicInterpolation::FirstDerivative, 0.0);
    f.update();
    checkValues("MC clamped spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    check1stDerivativeValue("MC clamped spline", f,
                            *std::begin(RPN15A_x), 0.0);
    check1stDerivativeValue("MC clamped spline", f,
                            *(std::end(RPN15A_x)-1), 0.0);
    // good performance
    interpolated = f(x_bad);
    if (interpolated>1.0) {
   
        BOOST_ERROR("MC clamped spline interpolation "
                    << "good performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value < 1.0");
    }


    // MC not-a-knot spline values
    f = CubicInterpolation(std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y),
                           CubicInterpolation::Spline, true,
                           CubicInterpolation::NotAKnot, Null<Real>(),
                           CubicInterpolation::NotAKnot, Null<Real>());
    f.update();
    checkValues("MC not-a-knot spline", f,
                std::begin(RPN15A_x), std::end(RPN15A_x), std::begin(RPN15A_y));
    // good performance
    interpolated = f(x_bad);
    if (interpolated>1.0) {
   
        BOOST_ERROR("MC clamped spline interpolation "
                    << "good performance unverified\n"
                    << "at x = " << x_bad
                    << "\ninterpolated value: " << interpolated
                    << "\nexpected value < 1.0");
    }
}

/* Blossey, Frigyik, Farnum "A Note On CubicSpline Splines"
   Applied Linear Algebra and Numerical Analysis AMATH 352 Lecture Notes
   http://www.amath.washington.edu/courses/352-winter-2002/spline_note.pdf
*/
void InterpolationTest::testSplineOnGenericValues() {
   

    BOOST_TEST_MESSAGE("Testing spline interpolation on generic values...");

    const Real generic_x[] = {
    0.0, 1.0, 3.0, 4.0 };
    const Real generic_y[] = {
    0.0, 0.0, 2.0, 2.0 };
    const Real generic_natural_y2[] = {
    0.0, 1.5, -1.5, 0.0 };

    Real interpolated, error;
    Size i, n = LENGTH(generic_x);
    std::vector<Real> x35(3);

    // Natural spline
    CubicInterpolation f(std::begin(generic_x), std::end(generic_x),
                         std::begin(generic_y),
                         CubicInterpolation::Spline, false,
                         CubicInterpolation::SecondDerivative,
                         generic_natural_y2[0],
                         CubicInterpolation::SecondDerivative,
                         generic_natural_y2[n-1]);
    f.update();
    checkValues("Natural spline", f,
                std::begin(generic_x), std::end(generic_x), std::begin(generic_y));
    // cached second derivative
    for (i=0; i<n; i++) {
   
        interpolated = f.secondDerivative(generic_x[i]);
        error = interpolated - generic_natural_y2[i];
        if (std::fabs(error)>3e-16) {
   
            BOOST_ERROR("Natural spline interpolation "
                        << "second derivative failed at x=" << generic_x[i]
                        << "\ninterpolated value: " << interpolated
                        << "\nexpected value:     " << generic_natural_y2[i]
                        << "\nerror:              " << error);
        }
    }
    x35[1] = f(3.5);


    // Clamped spline
    Real y1a = 0.0, y1b = 0.0;
    f = CubicInterpolation(std::begin(generic_x), std::end(generic_x), std::begin(generic_y),
                    CubicInterpolation::Spline, false,
                    CubicInterpolation::FirstDerivative, y1a,
                    CubicInterpolation::FirstDerivative, y1b);
    f.update();
    checkValues("Clamped spline", f,
                std::begin(generic_x), std::end(generic_x), std::begin(generic_y));
    check1stDerivativeValue("Clamped spline", f,
                            *std::begin(generic_x), 0.0);
    check1stDerivativeValue("Clamped spline", f,
                            *(std::end(generic_x)-1), 0.0);
    x35[0] = f(3.5);


    // Not-a-knot spline
    f = CubicInterpolation(std::begin(generic_x), std::end(generic_x), std::begin(generic_y),
                           CubicInterpolation::Spline, false,
                           CubicInterpolation::NotAKnot, Null<Real>(),
                           CubicInterpolation::NotAKnot, Null<Real>());
    f.update();
    checkValues("Not-a-knot spline", f,
                std::begin(generic_x), std::end(generic_x), std::begin(generic_y));
    checkNotAKnotCondition("Not-a-knot spline", f);

    x35[2] = f(3.5);

    if (x35[0]>x35[1] || x35[1]>x35[2]) {
   
        BOOST_ERROR("Spline interpolation failure"
                    << "\nat x = " << 3.5
                    << "\nclamped spline    " << x35[0]
                    << "\nnatural spline    " << x35[1]
                    << "\nnot-a-knot spline " << x35[2]
                    << "\nvalues should be in increasing order");
    }
}


void InterpolationTest::testSimmetricEndConditions() {
   

    BOOST_TEST_MESSAGE("Testing symmetry of spline interpolation "
                       "end-conditions...");

    Size n = 9;

    std::vector<Real> x, y;
    x = xRange(-1.8, 1.8, n);
    y = gaussian(x);

    // Not-a-knot spline
    CubicInterpolation f(x.begin(), x.end(), y.begin(),
                         CubicInterpolation::Spline, false,
                         CubicInterpolation::NotAKnot, Null<Real>(),
                         CubicInterpolation::NotAKnot, Null<Real>());
    f.update();
    checkValues("Not-a-knot spline", f,
                x.begin(), x