fittedbonddiscountcurve.hpp

Go to the documentation of this file.

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2007 Allen Kuo
 Copyright (C) 2009 Ferdinando Ametrano
 Copyright (C) 2015 Andres Hernandez
 
 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 fittedbonddiscountcurve.hpp
    \brief discount curve fitted to a set of bonds
*/
 
#ifndef quantlib_fitted_bond_discount_curve_hpp
#define quantlib_fitted_bond_discount_curve_hpp
 
#include <ql/termstructures/yield/bondhelpers.hpp>
#include <ql/math/optimization/method.hpp>
#include <ql/patterns/lazyobject.hpp>
#include <ql/math/array.hpp>
#include <ql/utilities/clone.hpp>
#include <ql/math/optimization/constraint.hpp>
 
namespace QuantLib {
 
    //! Discount curve fitted to a set of fixed-coupon bonds
    /*! This class fits a discount function \f$ d(t) \f$ over a set of
        bonds, using a user defined fitting method. The discount
        function is fit in such a way so that all cashflows of all
        input bonds, when discounted using \f$ d(t) \f$, will
        reproduce the set of input bond prices in an optimized
        sense. Minimized price errors are weighted by the inverse of
        their respective bond duration.
 
        The FittedBondDiscountCurve class acts as a generic wrapper,
        while its inner class FittingMethod provides the
        implementation details. Developers thus need only derive new
        fitting methods from the latter.
 
        \warning The method can be slow if there are many bonds to
                 fit. Speed also depends on the particular choice of
                 fitting method chosen and its convergence properties
                 under optimization.  See also todo list for
                 BondDiscountCurveFittingMethod.
 
        \todo refactor the bond helper class so that it is pure
              virtual and returns a generic bond or its cash
              flows. Derived classes would include helpers for
              fixed-rate and zero-coupon bonds. In this way, both
              bonds and bills can be used to fit a discount curve
              using the exact same machinery. At present, only
              fixed-coupon bonds are supported. An even better way to
              move forward might be to get rate helpers to return
              cashflows, in which case this class could be used to fit
              any set of cash flows, not just bonds.
 
        \todo add more fitting diagnostics: smoothness, standard
              deviation, student-t test, etc. Generic smoothness
              method may be useful for smoothing splines fitting. See
              Fisher, M., D. Nychka and D. Zervos: "Fitting the term
              structure of interest rates with smoothing splines."
              Board of Governors of the Federal Reserve System,
              Federal Resere Board Working Paper, 95-1.
 
        \todo add extrapolation routines
 
        \ingroup yieldtermstructures
    */
    class FittedBondDiscountCurve : public YieldTermStructure,
                                    public LazyObject {
      public:
        class FittingMethod;
 
        //! \name Constructors
        //@{
        //! reference date based on current evaluation date
        FittedBondDiscountCurve(Natural settlementDays,
                                const Calendar& calendar,
                                std::vector<ext::shared_ptr<BondHelper> > bonds,
                                const DayCounter& dayCounter,
                                const FittingMethod& fittingMethod,
                                Real accuracy = 1.0e-10,
                                Size maxEvaluations = 10000,
                                Array guess = Array(),
                                Real simplexLambda = 1.0,
                                Size maxStationaryStateIterations = 100);
 
        //! curve reference date fixed for life of curve
        FittedBondDiscountCurve(const Date& referenceDate,
                                std::vector<ext::shared_ptr<BondHelper> > bonds,
                                const DayCounter& dayCounter,
                                const FittingMethod& fittingMethod,
                                Real accuracy = 1.0e-10,
                                Size maxEvaluations = 10000,
                                Array guess = Array(),
                                Real simplexLambda = 1.0,
                                Size maxStationaryStateIterations = 100);
 
        //! don't fit, use precalculated parameters
        FittedBondDiscountCurve(Natural settlementDays,
                                const Calendar& calendar,
                                const FittingMethod& fittingMethod,
                                Array parameters,
                                Date maxDate,
                                const DayCounter& dayCounter);
 
        //! don't fit, use precalculated parameters
        FittedBondDiscountCurve(const Date& referenceDate,
                                const FittingMethod& fittingMethod,
                                Array parameters,
                                Date maxDate,
                                const DayCounter& dayCounter);
        //@}
 
        //! \name Inspectors
        //@{
        //! total number of bonds used to fit the yield curve
        Size numberOfBonds() const;
        //! the latest date for which the curve can return values
        Date maxDate() const override;
        //! class holding the results of the fit
        const FittingMethod& fitResults() const;
        //@}
 
        //! \name Other utilities
        //@{
        /*! This allows to try out multiple guesses and avoid local minima */
        void resetGuess(const Array& guess);
        //@}
 
        //! \name Observer interface
        //@{
        void update() override;
        //@}
 
      private:
        void setup();
        void performCalculations() const override;
        DiscountFactor discountImpl(Time) const override;
        // target accuracy level to be used in the optimization routine
        Real accuracy_;
        // max number of evaluations to be used in the optimization routine
        Size maxEvaluations_;
        // sets the scale in the (Simplex) optimization routine
        Real simplexLambda_;
        // max number of evaluations where no improvement to solution is made
        Size maxStationaryStateIterations_;
        // a guess solution may be passed into the constructor to speed calcs
        Array guessSolution_;
        mutable Date maxDate_;
        std::vector<ext::shared_ptr<BondHelper> > bondHelpers_;
        Clone<FittingMethod> fittingMethod_;
    };
 
 
    //! Base fitting method used to construct a fitted bond discount curve
    /*! This base class provides the specific methodology/strategy
        used to construct a FittedBondDiscountCurve.  Derived classes
        need only define the virtual function discountFunction() based
        on the particular fitting method to be implemented, as well as
        size(), the number of variables to be solved for/optimized. The
        generic fitting methodology implemented here can be termed
        nonlinear, in contrast to (typically faster, computationally)
        linear fitting method.
 
        Optional parameters for FittingMethod include an Array of
        weights, which will be used as weights to each bond. If not given
        or empty, then the bonds will be weighted by inverse duration
 
        An optional Array may be provided as an L2 regularizor in this case
        a L2 (gaussian) penalty is applied to each parameter starting from the 
        initial guess. This is the same as giving a Gaussian prior on the parameters
 
        \todo derive the special-case class LinearFittingMethods from
              FittingMethod. A linear fitting to a set of basis
              functions \f$ b_i(t) \f$ is any fitting of the form
              \f[
              d(t) = \sum_{i=0} c_i b_i(t)
              \f]
              i.e., linear in the unknown coefficients \f$ c_i
              \f$. Such a fitting can be reduced to a linear algebra
              problem \f$ Ax = b \f$, and for large numbers of bonds,
              would typically be much faster computationally than the
              generic non-linear fitting method.
 
        \warning some parameters to the Simplex optimization method
                 may need to be tweaked internally to the class,
                 depending on the fitting method used, in order to get
                 proper/reasonable/faster convergence.
    */
    class FittedBondDiscountCurve::FittingMethod {
        friend class FittedBondDiscountCurve;
        // internal class
        class FittingCost;
      public:
        virtual ~FittingMethod() = default;
        //! total number of coefficients to fit/solve for
        virtual Size size() const = 0;
        //! output array of results of optimization problem
        Array solution() const;
        //! final number of iterations used in the optimization problem
        Integer numberOfIterations() const;
        //! final value of cost function after optimization
        Real minimumCostValue() const;
        //! error code of the optimization
        EndCriteria::Type errorCode() const;
        //! clone of the current object
        virtual std::unique_ptr<FittingMethod> clone() const = 0;
        //! return whether there is a constraint at zero
        bool constrainAtZero() const;
        //! return weights being used
        Array weights() const;
        //! return l2 penalties being used
        Array l2() const;
        //! return optimization method being used
        ext::shared_ptr<OptimizationMethod> optimizationMethod() const;
        //! return optimization contraint
        const Constraint& constraint() const;
        //! open discountFunction to public
        DiscountFactor discount(const Array& x, Time t) const;
      protected:
        //! constructors
        FittingMethod(bool constrainAtZero = true,
                      const Array& weights = Array(),
                      ext::shared_ptr<OptimizationMethod> optimizationMethod =
                          ext::shared_ptr<OptimizationMethod>(),
                      Array l2 = Array(),
                      Real minCutoffTime = 0.0,
                      Real maxCutoffTime = QL_MAX_REAL,
                      Constraint constraint = NoConstraint());
        //! rerun every time instruments/referenceDate changes
        virtual void init();
        //! discount function called by FittedBondDiscountCurve
        virtual DiscountFactor discountFunction(const Array& x,
                                                Time t) const = 0;
 
        //! constrains discount function to unity at \f$ T=0 \f$, if true
        bool constrainAtZero_;
        //! internal reference to the FittedBondDiscountCurve instance
        FittedBondDiscountCurve* curve_;
        //! solution array found from optimization, set in calculate()
        Array solution_;
        //! optional guess solution to be passed into constructor.
        /*! The idea is to use a previous solution as a guess solution to
            the discount curve, in an attempt to speed up calculations.
        */
        Array guessSolution_;
        //! base class sets this cost function used in the optimization routine
        ext::shared_ptr<FittingCost> costFunction_;
      private:
        // curve optimization called here- adjust optimization parameters here
        void calculate();
        // array of normalized (duration) weights, one for each bond helper
        Array weights_;
        // array of l2 penalties one for each parameter
        Array l2_;
        // whether or not the weights should be calculated internally
        bool calculateWeights_;
        // total number of iterations used in the optimization routine
        // (possibly including gradient evaluations)
        Integer numberOfIterations_;
        // final value for the minimized cost function
        Real costValue_;
        // error code returned by OptimizationMethod::minimize()
        EndCriteria::Type errorCode_ = EndCriteria::None;
        // optimization method to be used, if none provided use Simplex
        ext::shared_ptr<OptimizationMethod> optimizationMethod_;
        // optimization constraint, if none provided use NoConstraint
        Constraint constraint_;
        // flat extrapolation of instantaneous forward before / after cutoff
        Real minCutoffTime_, maxCutoffTime_;
    };
 
    // inline
 
    inline Size FittedBondDiscountCurve::numberOfBonds() const {
        return bondHelpers_.size();
    }
 
    inline Date FittedBondDiscountCurve::maxDate() const {
        calculate();
        return maxDate_;
    }
 
    inline const FittedBondDiscountCurve::FittingMethod&
    FittedBondDiscountCurve::fitResults() const {
        calculate();
        return *fittingMethod_;
    }
 
    inline void FittedBondDiscountCurve::update() {
        YieldTermStructure::update();
        LazyObject::update();
    }
 
    inline void FittedBondDiscountCurve::setup() {
        for (auto& bondHelper : bondHelpers_)
            registerWith(bondHelper);
    }
 
    inline DiscountFactor FittedBondDiscountCurve::discountImpl(Time t) const {
        calculate();
        return fittingMethod_->discount(fittingMethod_->solution_, t);
    }
 
    inline Integer
    FittedBondDiscountCurve::FittingMethod::numberOfIterations() const {
        return numberOfIterations_;
    }
 
    inline
    Real FittedBondDiscountCurve::FittingMethod::minimumCostValue() const {
        return costValue_;
    }
 
    inline 
    EndCriteria::Type FittedBondDiscountCurve::FittingMethod::errorCode() const {
        return errorCode_;
    }
 
    inline Array FittedBondDiscountCurve::FittingMethod::solution() const {
        return solution_;
    }
    
    inline bool FittedBondDiscountCurve::FittingMethod::constrainAtZero() const {
        return constrainAtZero_;
    }
    
    inline Array FittedBondDiscountCurve::FittingMethod::weights() const {
        return weights_;
    }
 
    inline Array FittedBondDiscountCurve::FittingMethod::l2() const {
        return l2_;
    }
 
    inline ext::shared_ptr<OptimizationMethod> 
    FittedBondDiscountCurve::FittingMethod::optimizationMethod() const {
        return optimizationMethod_;
    }
 
    inline const Constraint& FittedBondDiscountCurve::FittingMethod::constraint() const {
        return constraint_;
    }
 
    inline DiscountFactor FittedBondDiscountCurve::FittingMethod::discount(const Array& x, Time t) const {
        if (t < minCutoffTime_) {
            // flat fwd extrapolation before min cutoff time
            return std::exp(std::log(discountFunction(x, minCutoffTime_)) / minCutoffTime_ * t);
        } else if (t > maxCutoffTime_) {
            // flat fwd extrapolation after max cutoff time
            return discountFunction(x, maxCutoffTime_) *
                   std::exp((std::log(discountFunction(x, maxCutoffTime_ + 1E-4)) -
                             std::log(discountFunction(x, maxCutoffTime_))) *
                            1E4 * (t - maxCutoffTime_));
        } else {
            return discountFunction(x, t);
        }
    }
}
 
#endif