/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ /* Copyright (C) 2004 Ferdinando Ametrano Copyright (C) 2001, 2002, 2003 Nicolas Di Césaré Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl 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/reference/license.html>. 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/FiniteDifferences/tridiagonaloperator.hpp> #include <ql/Utilities/null.hpp> #include <vector> namespace QuantLib { namespace detail { template <class I1, class I2> class CubicSplineImpl; class CoefficientHolder { public: CoefficientHolder(Size n) : n_(n), primitiveConst_(n-1), a_(n-1), b_(n-1), c_(n-1) {} 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_; }; } //! %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 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 Cubic 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. */ 00074 class CubicSpline : public Interpolation { public: 00076 enum BoundaryCondition { //! Make second(-last) point an inactive knot 00078 NotAKnot, //! Match value of end-slope 00080 FirstDerivative, //! Match value of second derivative at end 00082 SecondDerivative, //! Match first and second derivative at either end 00084 Periodic, /*! Match end-slope to the slope of the cubic that matches the first four data at the respective end */ 00088 Lagrange }; protected: private: boost::shared_ptr<detail::CoefficientHolder> coeffs_; public: /*! \pre the \f$ x \f$ values must be sorted. */ template <class I1, class I2> 00096 CubicSpline(const I1& xBegin, const I1& xEnd, const I2& yBegin, CubicSpline::BoundaryCondition leftCondition, Real leftConditionValue, CubicSpline::BoundaryCondition rightCondition, Real rightConditionValue, bool monotonicityConstraint) { impl_ = boost::shared_ptr<Interpolation::Impl>( new detail::CubicSplineImpl<I1,I2>( xBegin, xEnd, yBegin, leftCondition, leftConditionValue, rightCondition, rightConditionValue, monotonicityConstraint)); coeffs_ = boost::dynamic_pointer_cast<detail::CoefficientHolder>(impl_); } 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_; } }; // convenience classes //! %Cubic spline with monotonicity constraint 00126 class MonotonicCubicSpline : public CubicSpline { public: /*! \pre the \f$ x \f$ values must be sorted. */ template <class I1, class I2> 00130 MonotonicCubicSpline(const I1& xBegin, const I1& xEnd, const I2& yBegin, CubicSpline::BoundaryCondition leftCondition, Real leftConditionValue, CubicSpline::BoundaryCondition rightCondition, Real rightConditionValue) : CubicSpline(xBegin,xEnd,yBegin, leftCondition,leftConditionValue, rightCondition,rightConditionValue, true) {} }; //! %Cubic spline with null second derivative at end points 00143 class NaturalCubicSpline : public CubicSpline { public: /*! \pre the \f$ x \f$ values must be sorted. */ template <class I1, class I2> 00147 NaturalCubicSpline(const I1& xBegin, const I1& xEnd, const I2& yBegin) : CubicSpline(xBegin,xEnd,yBegin, SecondDerivative, 0.0, SecondDerivative, 0.0, false) {} }; //! Natural cubic spline with monotonicity constraint 00156 class NaturalMonotonicCubicSpline : public CubicSpline { public: /*! \pre the \f$ x \f$ values must be sorted. */ template <class I1, class I2> 00160 NaturalMonotonicCubicSpline(const I1& xBegin, const I1& xEnd, const I2& yBegin) : CubicSpline(xBegin,xEnd,yBegin, SecondDerivative, 0.0, SecondDerivative, 0.0, true) {} }; //! %cubic-spline interpolation factory and traits 00169 class Cubic { public: Cubic(CubicSpline::BoundaryCondition leftCondition = CubicSpline::SecondDerivative, Real leftConditionValue = 0.0, CubicSpline::BoundaryCondition rightCondition = CubicSpline::SecondDerivative, Real rightConditionValue = 0.0, bool monotonicityConstraint = false) : leftCondition_(leftCondition), leftValue_(leftConditionValue), rightCondition_(rightCondition), rightValue_(rightConditionValue), monotone_(monotonicityConstraint) {} template <class I1, class I2> Interpolation interpolate(const I1& xBegin, const I1& xEnd, const I2& yBegin) const { return CubicSpline(xBegin,xEnd,yBegin, leftCondition_,leftValue_, rightCondition_,rightValue_, monotone_); } enum { global = 1 }; private: CubicSpline::BoundaryCondition leftCondition_; Real leftValue_; CubicSpline::BoundaryCondition rightCondition_; Real rightValue_; bool monotone_; }; namespace detail { template <class I1, class I2> class CubicSplineImpl : public Interpolation::templateImpl<I1,I2>, public CoefficientHolder { public: CubicSplineImpl(const I1& xBegin, const I1& xEnd, const I2& yBegin, CubicSpline::BoundaryCondition leftCondition, Real leftConditionValue, CubicSpline::BoundaryCondition rightCondition, Real rightConditionValue, bool monotonicityConstraint) : Interpolation::templateImpl<I1,I2>(xBegin,xEnd,yBegin), CoefficientHolder(xEnd-xBegin), monotone_(false), constrained_(monotonicityConstraint), leftType_(leftCondition), rightType_(rightCondition), leftValue_(leftConditionValue), rightValue_(rightConditionValue) { calculate(); } void calculate() { 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 CubicSpline::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 CubicSpline::FirstDerivative: L.setFirstRow(1.0, 0.0); tmp[0] = leftValue_; break; case CubicSpline::SecondDerivative: L.setFirstRow(2.0, 1.0); tmp[0] = 3.0*S[0] - leftValue_*dx[0]/2.0; break; case CubicSpline::Periodic: case CubicSpline::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 CubicSpline::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 CubicSpline::FirstDerivative: L.setLastRow(0.0, 1.0); tmp[n_-1] = rightValue_; break; case CubicSpline::SecondDerivative: L.setLastRow(1.0, 2.0); tmp[n_-1] = 3.0*S[n_-2] + rightValue_*dx[n_-2]/2.0; break; case CubicSpline::Periodic: case CubicSpline::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); 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; monotone_ = 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; monotone_ = 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; monotone_ = 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 monotone_, constrained_; CubicSpline::BoundaryCondition leftType_, rightType_; Real leftValue_, rightValue_; }; } } #endif