Logo Search packages:      
Sourcecode: quantlib version File versions  Download package

cubicspline.hpp

Go to the documentation of this file.
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
 Copyright (C) 2001, 2002, 2003 Nicolas Di Césaré
 Copyright (C) 2004 Ferdinando Ametrano

 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 cubicspline.hpp
    \brief cubic spline interpolation between discrete points
*/

#ifndef quantlib_cubic_spline_hpp
#define quantlib_cubic_spline_hpp

#include <ql/math/interpolation.hpp>
#include <ql/methods/finitedifferences/tridiagonaloperator.hpp>
#include <ql/utilities/null.hpp>
#include <vector>

namespace QuantLib {

    namespace detail {

        class CoefficientHolder {
          public:
            CoefficientHolder(Size n)
            : n_(n), primitiveConst_(n-1), a_(n-1), b_(n-1), c_(n-1),
              monotonicityAdjustments_(n) {}
            virtual ~CoefficientHolder() {}
            Size n_;
            // P[i](x) = y[i] +
            //           a[i]*(x-x[i]) +
            //           b[i]*(x-x[i])^2 +
            //           c[i]*(x-x[i])^3
            std::vector<Real> primitiveConst_, a_, b_, c_;
            std::vector<bool> monotonicityAdjustments_;
        };

        template <class I1, class I2> class CubicSplineInterpolationImpl;

    }

    //! %Cubic spline interpolation between discrete points.
    /*! It implements different type of end conditions: not-a-knot,
        first derivative value, second derivative value.

        It also implements Hyman's monotonicity constraint filter
        which ensures that in the regions of monotoniticity of the input
        (so, three successive increasing or decreasing values)
        the interpolating spline remains monotonic at the expense of the
        second derivative of the curve which will no longer be continuous
        where the filter has been applied. If the interpolating spline
        is already monotonic, the Hyman filter leaves it unchanged.

        See R. L. Dougherty, A. Edelman, and J. M. Hyman,
        "Nonnegativity-, Monotonicity-, or Convexity-Preserving CubicSpline and
        Quintic Hermite Interpolation"
        Mathematics Of Computation, v. 52, n. 186, April 1989, pp. 471-494.

        \test the correctness of the returned values is tested by
              reproducing results available in literature.
    */
00077     class CubicSplineInterpolation : public Interpolation {
      public:
00079         enum BoundaryCondition {
            //! Make second(-last) point an inactive knot
00081             NotAKnot,
            //! Match value of end-slope
00083             FirstDerivative,
            //! Match value of second derivative at end
00085             SecondDerivative,
            //! Match first and second derivative at either end
00087             Periodic,
            /*! Match end-slope to the slope of the cubic that matches
                the first four data at the respective end
            */
00091             Lagrange
        };
        /*! \pre the \f$ x \f$ values must be sorted. */
        template <class I1, class I2>
00095         CubicSplineInterpolation(const I1& xBegin, const I1& xEnd, const I2& yBegin,
                    CubicSplineInterpolation::BoundaryCondition leftCondition,
                    Real leftConditionValue,
                    CubicSplineInterpolation::BoundaryCondition rightCondition,
                    Real rightConditionValue,
                    bool monotonicityConstraint) {
            impl_ = boost::shared_ptr<Interpolation::Impl>(new
                detail::CubicSplineInterpolationImpl<I1,I2>(xBegin, xEnd, yBegin,
                                               leftCondition,
                                               leftConditionValue,
                                               rightCondition,
                                               rightConditionValue,
                                               monotonicityConstraint));
            impl_->update();
            coeffs_ =
                boost::dynamic_pointer_cast<detail::CoefficientHolder>(impl_);
        }
        const std::vector<Real>& primitiveConstants() const {
            return coeffs_->primitiveConst_;
        }
        const std::vector<Real>& aCoefficients() const { return coeffs_->a_; }
        const std::vector<Real>& bCoefficients() const { return coeffs_->b_; }
        const std::vector<Real>& cCoefficients() const { return coeffs_->c_; }
        const std::vector<bool>& monotonicityAdjustments() const {
            return coeffs_->monotonicityAdjustments_;
        }
      private:
        boost::shared_ptr<detail::CoefficientHolder> coeffs_;
    };


    // convenience classes

    //! %Cubic spline with monotonicity constraint
00129     class MonotonicCubicSpline : public CubicSplineInterpolation {
      public:
        /*! \pre the \f$ x \f$ values must be sorted. */
        template <class I1, class I2>
00133         MonotonicCubicSpline(const I1& xBegin, const I1& xEnd,
                             const I2& yBegin,
                             CubicSplineInterpolation::BoundaryCondition leftCondition,
                             Real leftConditionValue,
                             CubicSplineInterpolation::BoundaryCondition rightCondition,
                             Real rightConditionValue)
        : CubicSplineInterpolation(xBegin,xEnd,yBegin,
                      leftCondition,leftConditionValue,
                      rightCondition,rightConditionValue,
                      true) {}
    };

    //! %Cubic spline with null second derivative at end points
00146     class NaturalCubicSpline : public CubicSplineInterpolation {
      public:
        /*! \pre the \f$ x \f$ values must be sorted. */
        template <class I1, class I2>
00150         NaturalCubicSpline(const I1& xBegin, const I1& xEnd,
                           const I2& yBegin)
        : CubicSplineInterpolation(xBegin,xEnd,yBegin,
                      SecondDerivative, 0.0,
                      SecondDerivative, 0.0,
                      false) {}
    };

    //! Natural cubic spline with monotonicity constraint
00159     class NaturalMonotonicCubicSpline : public CubicSplineInterpolation {
      public:
        /*! \pre the \f$ x \f$ values must be sorted. */
        template <class I1, class I2>
00163         NaturalMonotonicCubicSpline(const I1& xBegin, const I1& xEnd,
                                    const I2& yBegin)
        : CubicSplineInterpolation(xBegin,xEnd,yBegin,
                      SecondDerivative, 0.0,
                      SecondDerivative, 0.0,
                      true) {}
    };

    //! %Cubic spline interpolation factory and traits
00172     class CubicSpline {
      public:
        CubicSpline(CubicSplineInterpolation::BoundaryCondition leftCondition
                        = CubicSplineInterpolation::SecondDerivative,
                    Real leftConditionValue = 0.0,
                    CubicSplineInterpolation::BoundaryCondition rightCondition
                        = CubicSplineInterpolation::SecondDerivative,
                    Real rightConditionValue = 0.0,
                    bool monotonicityConstraint = true)
        : leftType_(leftCondition), rightType_(rightCondition),
          leftValue_(leftConditionValue), rightValue_(rightConditionValue),
          monotonic_(monotonicityConstraint) {}
        template <class I1, class I2>
        Interpolation interpolate(const I1& xBegin, const I1& xEnd,
                                  const I2& yBegin) const {
            return CubicSplineInterpolation(xBegin, xEnd, yBegin,
                                            leftType_,leftValue_,
                                            rightType_,rightValue_,
                                            monotonic_);
        }
        enum { global = 1 };
        enum { requiredPoints = 3 };
      private:
        CubicSplineInterpolation::BoundaryCondition leftType_, rightType_;
        Real leftValue_, rightValue_;
        bool monotonic_;
    };


    namespace detail {

        template <class I1, class I2>
        class CubicSplineInterpolationImpl : public Interpolation::templateImpl<I1,I2>,
                                public CoefficientHolder {
          public:
            CubicSplineInterpolationImpl(const I1& xBegin, const I1& xEnd, const I2& yBegin,
                            CubicSplineInterpolation::BoundaryCondition leftCondition,
                            Real leftConditionValue,
                            CubicSplineInterpolation::BoundaryCondition rightCondition,
                            Real rightConditionValue,
                            bool monotonicityConstraint)
            : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin),
              CoefficientHolder(xEnd-xBegin),
              constrained_(monotonicityConstraint),
              leftType_(leftCondition), rightType_(rightCondition),
              leftValue_(leftConditionValue),
              rightValue_(rightConditionValue) {}

            void update() {

                TridiagonalOperator L(n_);
                Array tmp(n_);
                std::vector<Real> dx(n_-1), S(n_-1);

                Size i=0;
                dx[i] = this->xBegin_[i+1] - this->xBegin_[i];
                S[i] = (this->yBegin_[i+1] - this->yBegin_[i])/dx[i];
                for (i=1; i<n_-1; ++i) {
                    dx[i] = this->xBegin_[i+1] - this->xBegin_[i];
                    S[i] = (this->yBegin_[i+1] - this->yBegin_[i])/dx[i];

                    L.setMidRow(i, dx[i], 2.0*(dx[i]+dx[i-1]), dx[i-1]);
                    tmp[i] = 3.0*(dx[i]*S[i-1] + dx[i-1]*S[i]);
                }

                /**** BOUNDARY CONDITIONS ****/

                // left condition
                switch (leftType_) {
                  case CubicSplineInterpolation::NotAKnot:
                    // ignoring end condition value
                    L.setFirstRow(dx[1]*(dx[1]+dx[0]),
                                  (dx[0]+dx[1])*(dx[0]+dx[1]));
                    tmp[0] = S[0]*dx[1]*(2.0*dx[1]+3.0*dx[0]) +
                             S[1]*dx[0]*dx[0];
                    break;
                  case CubicSplineInterpolation::FirstDerivative:
                    L.setFirstRow(1.0, 0.0);
                    tmp[0] = leftValue_;
                    break;
                  case CubicSplineInterpolation::SecondDerivative:
                    L.setFirstRow(2.0, 1.0);
                    tmp[0] = 3.0*S[0] - leftValue_*dx[0]/2.0;
                    break;
                  case CubicSplineInterpolation::Periodic:
                  case CubicSplineInterpolation::Lagrange:
                    // ignoring end condition value
                    QL_FAIL("this end condition is not implemented yet");
                  default:
                    QL_FAIL("unknown end condition");
                }

                // right condition
                switch (rightType_) {
                  case CubicSplineInterpolation::NotAKnot:
                    // ignoring end condition value
                    L.setLastRow(-(dx[n_-2]+dx[n_-3])*(dx[n_-2]+dx[n_-3]),
                                 -dx[n_-3]*(dx[n_-3]+dx[n_-2]));
                    tmp[n_-1] = -S[n_-3]*dx[n_-2]*dx[n_-2] -
                                 S[n_-2]*dx[n_-3]*(3.0*dx[n_-2]+2.0*dx[n_-3]);
                    break;
                  case CubicSplineInterpolation::FirstDerivative:
                    L.setLastRow(0.0, 1.0);
                    tmp[n_-1] = rightValue_;
                    break;
                  case CubicSplineInterpolation::SecondDerivative:
                    L.setLastRow(1.0, 2.0);
                    tmp[n_-1] = 3.0*S[n_-2] + rightValue_*dx[n_-2]/2.0;
                    break;
                  case CubicSplineInterpolation::Periodic:
                  case CubicSplineInterpolation::Lagrange:
                    // ignoring end condition value
                    QL_FAIL("this end condition is not implemented yet");
                  default:
                    QL_FAIL("unknown end condition");
                }

                // solve the system
                tmp = L.solveFor(tmp);

                std::fill(monotonicityAdjustments_.begin(),
                          monotonicityAdjustments_.end(), false);
                if (constrained_) {
                    Real correction;
                    Real pm, pu, pd, M;
                    for (i=0; i<n_; ++i) {
                        if (i==0) {
                            if (tmp[i]*S[0]>0.0) {
                                correction = tmp[i]/std::fabs(tmp[i]) *
                                    std::min<Real>(std::fabs(tmp[i]),
                                                   std::fabs(3.0*S[0]));
                            } else {
                                correction = 0.0;
                            }
                            if (correction!=tmp[i]) {
                                tmp[i] = correction;
                                monotonicityAdjustments_[i] = true;
                            }
                        } else if (i==n_-1) {
                            if (tmp[i]*S[n_-2]>0.0) {
                                correction = tmp[i]/std::fabs(tmp[i]) *
                                    std::min<Real>(std::fabs(tmp[i]),
                                                   std::fabs(3.0*S[n_-2]));
                            } else {
                                correction = 0.0;
                            }
                            if (correction!=tmp[i]) {
                                tmp[i] = correction;
                                monotonicityAdjustments_[i] = true;
                            }
                        } else {
                            pm=(S[i-1]*dx[i]+S[i]*dx[i-1])/
                                (dx[i-1]+dx[i]);
                            M = 3.0 * std::min(std::min(std::fabs(S[i-1]),
                                                        std::fabs(S[i])),
                                               std::fabs(pm));
                            if (i>1) {
                                if ((S[i-1]-S[i-2])*(S[i]-S[i-1])>0.0) {
                                    pd=(S[i-1]*(2.0*dx[i-1]+dx[i-2])
                                        -S[i-2]*dx[i-1])/
                                        (dx[i-2]+dx[i-1]);
                                    if (pm*pd>0.0 && pm*(S[i-1]-S[i-2])>0.0) {
                                        M = std::max<Real>(M, 1.5*std::min(
                                                std::fabs(pm),std::fabs(pd)));
                                    }
                                }
                            }
                            if (i<n_-2) {
                                if ((S[i]-S[i-1])*(S[i+1]-S[i])>0.0) {
                                    pu=(S[i]*(2.0*dx[i]+dx[i+1])-S[i+1]*dx[i])/
                                        (dx[i]+dx[i+1]);
                                    if (pm*pu>0.0 && -pm*(S[i]-S[i-1])>0.0) {
                                        M = std::max<Real>(M, 1.5*std::min(
                                                std::fabs(pm),std::fabs(pu)));
                                    }
                                }
                            }
                            if (tmp[i]*pm>0.0) {
                                correction = tmp[i]/std::fabs(tmp[i]) *
                                    std::min(std::fabs(tmp[i]), M);
                            } else {
                                correction = 0.0;
                            }
                            if (correction!=tmp[i]) {
                                tmp[i] = correction;
                                monotonicityAdjustments_[i] = true;
                            }
                        }
                    }
                }

                for (i=0; i<n_-1; ++i) {
                    a_[i] = tmp[i];
                    b_[i] = (3.0*S[i] - tmp[i+1] - 2.0*tmp[i])/dx[i];
                    c_[i] = (tmp[i+1] + tmp[i] - 2.0*S[i])/(dx[i]*dx[i]);
                }

                primitiveConst_[0] = 0.0;
                for (i=1; i<n_-1; ++i) {
                    primitiveConst_[i] = primitiveConst_[i-1]
                        + dx[i-1] *
                        (this->yBegin_[i-1] + dx[i-1] *
                         (a_[i-1]/2.0 + dx[i-1] *
                          (b_[i-1]/3.0 + dx[i-1] * c_[i-1]/4.0)));
                }
            }
            Real value(Real x) const {
                Size j = this->locate(x);
                Real dx = x-this->xBegin_[j];
                return this->yBegin_[j] + dx*(a_[j] + dx*(b_[j] + dx*c_[j]));
            }
            Real primitive(Real x) const {
                Size j = this->locate(x);
                Real dx = x-this->xBegin_[j];
                return primitiveConst_[j]
                    + dx*(this->yBegin_[j] + dx*(a_[j]/2.0
                    + dx*(b_[j]/3.0 + dx*c_[j]/4.0)));
            }
            Real derivative(Real x) const {
                Size j = this->locate(x);
                Real dx = x-this->xBegin_[j];
                return a_[j] + (2.0*b_[j] + 3.0*c_[j]*dx)*dx;
            }
            Real secondDerivative(Real x) const {
                Size j = this->locate(x);
                Real dx = x-this->xBegin_[j];
                return 2.0*b_[j] + 6.0*c_[j]*dx;
            }
          private:
            bool constrained_;
            CubicSplineInterpolation::BoundaryCondition leftType_, rightType_;
            Real leftValue_, rightValue_;
        };

    }

}

#endif

Generated by  Doxygen 1.6.0   Back to index