/src/quantlib/ql/experimental/shortrate/generalizedhullwhite.hpp

Source (jump to first uncovered line)
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2010 SunTrust Bank
 Copyright (C) 2010, 2014 Cavit Hafizoglu

 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 generalizedhullwhite.hpp
    \brief generalized Hull-White model
*/

#ifndef quantlib_generalized_hull_white_hpp
#define quantlib_generalized_hull_white_hpp

#include <ql/experimental/shortrate/generalizedornsteinuhlenbeckprocess.hpp>
#include <ql/math/interpolation.hpp>
#include <ql/models/shortrate/onefactormodel.hpp>
#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
#include <utility>

namespace QuantLib {

    //! Parameter that holds an interpolation object
    class InterpolationParameter : public Parameter {
    private:
        class Impl final : public Parameter::Impl {
        public:
          Real value(const Array&, Time t) const override { return interpolator_(t); }
          void reset(const Interpolation& interp) { interpolator_ = interp; }
        private:
            Interpolation interpolator_;
        };
    public:
        explicit InterpolationParameter(
            Size count,
            const Constraint& constraint = NoConstraint())
        : Parameter(count,
            ext::shared_ptr<Parameter::Impl>(
                new InterpolationParameter::Impl()),
                constraint)
        { }
        void reset(const Interpolation &interp) {
            ext::shared_ptr<InterpolationParameter::Impl> impl =
                ext::dynamic_pointer_cast<InterpolationParameter::Impl>(impl_);
            if (impl != nullptr)
                impl->reset(interp);
        }
    };

    //! Generalized Hull-White model class.
    /*! This class implements the standard Black-Karasinski model defined by
        \f[
        d f(r_t) = (\theta(t) - \alpha f(r_t))dt + \sigma dW_t,
        \f]
        where \f$ alpha \f$ and \f$ sigma \f$ are piecewise linear functions.

        \ingroup shortrate
    */
    class GeneralizedHullWhite : public OneFactorAffineModel,
                                 public TermStructureConsistentModel {
      public:

        GeneralizedHullWhite(
            const Handle<YieldTermStructure>& yieldtermStructure,
            const std::vector<Date>& speedstructure,
            const std::vector<Date>& volstructure,
            const std::vector<Real>& speed,
            const std::vector<Real>& vol,
            const std::function<Real(Real)>& f = {},
            const std::function<Real(Real)>& fInverse = {});

        template <class SpeedInterpolationTraits,class VolInterpolationTraits>
        GeneralizedHullWhite(
            const Handle<YieldTermStructure>& yieldtermStructure,
            const std::vector<Date>& speedstructure,
            const std::vector<Date>& volstructure,
            const std::vector<Real>& speed,
            const std::vector<Real>& vol,
            const SpeedInterpolationTraits &speedtraits,
            const VolInterpolationTraits &voltraits,
            const std::function<Real(Real)>& f = {},
            const std::function<Real(Real)>& fInverse = {}) :
            OneFactorAffineModel(2), TermStructureConsistentModel(yieldtermStructure),
            speedstructure_(speedstructure), volstructure_(volstructure),
            a_(arguments_[0]), sigma_(arguments_[1]),
            f_(f), fInverse_(fInverse)
        {
            initialize(yieldtermStructure,speedstructure,volstructure,
                speed,vol,speedtraits,voltraits,f,fInverse);
        }

        ext::shared_ptr<ShortRateDynamics> dynamics() const override {
            QL_FAIL("no defined process for generalized Hull-White model, "
                    "use HWdynamics()");
        }

        ext::shared_ptr<Lattice> tree(const TimeGrid& grid) const override;

        //Analytical calibration of HW

        GeneralizedHullWhite(
                  const Handle<YieldTermStructure>& yieldtermStructure,
                  Real a = 0.1, Real sigma = 0.01);


        ext::shared_ptr<ShortRateDynamics> HWdynamics() const;

        //! Only valid under Hull-White model
        Real discountBondOption(Option::Type type,
                                Real strike,
                                Time maturity,
                                Time bondMaturity) const override;

        //! vector to pass to 'calibrate' to fit only volatility
        std::vector<bool> fixedReversion() const;

      protected:
        //Analytical calibration of HW
        Real a() const { return a_(0.0); }
        Real sigma() const { return sigma_(0.0); }
        void generateArguments() override;
        Real A(Time t, Time T) const override;
        Real B(Time t, Time T) const override;
        Real V(Time t, Time T) const;

      private:

        class Dynamics;
        class Helper;
        class FittingParameter;// for analytic HW fitting

        std::vector<Date> speedstructure_;
        std::vector<Date> volstructure_;
        std::vector<Time> speedperiods_;
        std::vector<Time> volperiods_;
        Interpolation speed_;
        Interpolation vol_;

        std::function<Real (Time)> speed() const;
        std::function<Real (Time)> vol() const;

        Parameter& a_;
        Parameter& sigma_;
        Parameter phi_;

        std::function<Real(Real)> f_;
        std::function<Real(Real)> fInverse_;

        static Real identity(Real x) {
            return x;
        }

        template <class SpeedInterpolationTraits,class VolInterpolationTraits>
        void initialize(const Handle<YieldTermStructure>& yieldtermStructure,
            const std::vector<Date>& speedstructure,
            const std::vector<Date>& volstructure,
            const std::vector<Real>& speed,
            const std::vector<Real>& vol,
            const SpeedInterpolationTraits &speedtraits,
            const VolInterpolationTraits &voltraits,
            const std::function<Real(Real)>& f,
            const std::function<Real(Real)>& fInverse)
        {
            QL_REQUIRE(speedstructure.size()==speed.size(),
                "mean reversion inputs inconsistent");
            QL_REQUIRE(volstructure.size()==vol.size(),
                "volatility inputs inconsistent");
            if (!f_)
                f_ = identity;
            if (!fInverse_)
                fInverse_ = identity;

            DayCounter dc = yieldtermStructure->dayCounter();
            Date ref = yieldtermStructure->referenceDate();
            for (auto i : speedstructure)
                speedperiods_.push_back(dc.yearFraction(ref, i));
            for (auto i : volstructure)
                volperiods_.push_back(dc.yearFraction(ref, i));

            // interpolator x points to *periods_ vector, y points to
            // the internal Array in the parameter
            InterpolationParameter atemp(speedperiods_.size(), NoConstraint());
            a_ = atemp;
            for (Size i=0; i<speedperiods_.size(); i++)
                a_.setParam(i, speed[i]);
            speed_ = speedtraits.interpolate(speedperiods_.begin(),
                speedperiods_.end(),a_.params().begin());
            speed_.enableExtrapolation();
            atemp.reset(speed_);

            InterpolationParameter sigmatemp(volperiods_.size(), PositiveConstraint());
            sigma_ = sigmatemp;
            for (Size i=0; i<volperiods_.size(); i++)
                sigma_.setParam(i, vol[i]);
            vol_ = voltraits.interpolate(volperiods_.begin(),
                volperiods_.end(),sigma_.params().begin());
            vol_.enableExtrapolation();
            sigmatemp.reset(vol_);

            generateArguments();
            registerWith(yieldtermStructure);
        }
    };

    //! Short-rate dynamics in the generalized Hull-White model
    /*! The short-rate is here

        f(r_t) = x_t + g(t)

        where g is the deterministic time-dependent
        parameter (which can't be determined analytically)
        used for initial term-structure fitting and  x_t is the state
        variable following an Ornstein-Uhlenbeck process.

        In this version, the function f may also be defined as a piece-wise linear
        function and can be calibrated to the away-from-the-money instruments.

    */
    class GeneralizedHullWhite::Dynamics
        : public GeneralizedHullWhite::ShortRateDynamics {
      public:
        Dynamics(Parameter fitting,
                 const std::function<Real(Time)>& alpha,
                 const std::function<Real(Time)>& sigma,
                 std::function<Real(Real)> f,
                 std::function<Real(Real)> fInverse)
        : ShortRateDynamics(ext::shared_ptr<StochasticProcess1D>(
              new GeneralizedOrnsteinUhlenbeckProcess(alpha, sigma))),
          fitting_(std::move(fitting)), _f_(std::move(f)), _fInverse_(std::move(fInverse)) {}

        //classical HW dynamics
        Dynamics(Parameter fitting, Real a, Real sigma)
        : GeneralizedHullWhite::ShortRateDynamics(
              ext::shared_ptr<StochasticProcess1D>(new OrnsteinUhlenbeckProcess(a, sigma))),
          fitting_(std::move(fitting)), _f_(identity()), _fInverse_(identity()) {}

        Real variable(Time t, Rate r) const override { return _f_(r) - fitting_(t); }

        Real shortRate(Time t, Real x) const override { return _fInverse_(x + fitting_(t)); }

      private:
        Parameter fitting_;
        std::function<Real(Real)> _f_;
        std::function<Real(Real)> _fInverse_;
        struct identity {
            Real operator()(Real x) const {return x;};
        };
    };

    //! Analytical term-structure fitting parameter \f$ \varphi(t) \f$.
    /*! \f$ \varphi(t) \f$ is analytically defined by
        \f[
            \varphi(t) = f(t) + \frac{1}{2}[\frac{\sigma(1-e^{-at})}{a}]^2,
        \f]
        where \f$ f(t) \f$ is the instantaneous forward rate at \f$ t \f$.
    */
    class GeneralizedHullWhite::FittingParameter
        : public TermStructureFittingParameter {
      private:
        class Impl final : public Parameter::Impl {
          public:
            Impl(Handle<YieldTermStructure> termStructure, Real a, Real sigma)
            : termStructure_(std::move(termStructure)), a_(a), sigma_(sigma) {}

            Real value(const Array&, Time t) const override {
                Rate forwardRate =
                    termStructure_->forwardRate(t, t, Continuous, NoFrequency);
                Real temp = a_ < std::sqrt(QL_EPSILON) ?
                            Real(sigma_*t) :
                            Real(sigma_*(1.0 - std::exp(-a_*t))/a_);
                return (forwardRate + 0.5*temp*temp);
            }

          private:
            Handle<YieldTermStructure> termStructure_;
            Real a_, sigma_;
        };
      public:
        FittingParameter(const Handle<YieldTermStructure>& termStructure,
                         Real a, Real sigma)
        : TermStructureFittingParameter(ext::shared_ptr<Parameter::Impl>(
                      new FittingParameter::Impl(termStructure, a, sigma))) {}
    };

    // Analytic fitting dynamics
    inline ext::shared_ptr<OneFactorModel::ShortRateDynamics>
    GeneralizedHullWhite::HWdynamics() const {
        return ext::shared_ptr<ShortRateDynamics>(
          new Dynamics(phi_, a(), sigma()));
    }

    namespace detail {
        template <class I1, class I2>
        class LinearFlatInterpolationImpl;
    }

    //! %Linear interpolation between discrete points with flat extapolation
    /*! \ingroup interpolations */
    class LinearFlatInterpolation : public Interpolation {
      public:
        /*! \pre the \f$ x \f$ values must be sorted. */
        template <class I1, class I2>
        LinearFlatInterpolation(const I1& xBegin, const I1& xEnd,
                            const I2& yBegin) {
            impl_ = ext::shared_ptr<Interpolation::Impl>(new
                detail::LinearFlatInterpolationImpl<I1,I2>(xBegin, xEnd,
                                                       yBegin));
            impl_->update();
        }
    };

    //! %Linear-interpolation with flat-extrapolation factory and traits
    /*! \ingroup interpolations */
    class LinearFlat {
      public:
        template <class I1, class I2>
        Interpolation interpolate(const I1& xBegin, const I1& xEnd,
                                  const I2& yBegin) const {
            return LinearFlatInterpolation(xBegin, xEnd, yBegin);
        }
        static const bool global = false;
        static const Size requiredPoints = 1;
    };

    namespace detail {
        template <class I1, class I2>
        class LinearFlatInterpolationImpl
            : public Interpolation::templateImpl<I1,I2> {
          public:
            LinearFlatInterpolationImpl(const I1& xBegin, const I1& xEnd,
                                    const I2& yBegin)
            : Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin,
                                    LinearFlat::requiredPoints),
              primitiveConst_(xEnd-xBegin), s_(xEnd-xBegin) {}
            void update() override {
                primitiveConst_[0] = 0.0;
                for (Size i=1; i<Size(this->xEnd_-this->xBegin_); ++i) {
                    Real dx = this->xBegin_[i]-this->xBegin_[i-1];
                    s_[i-1] = (this->yBegin_[i]-this->yBegin_[i-1])/dx;
                    primitiveConst_[i] = primitiveConst_[i-1]
                        + dx*(this->yBegin_[i-1] +0.5*dx*s_[i-1]);
                }
            }
            Real value(Real x) const override {
                if (x <= this->xMin())
                    return this->yBegin_[0];
                if (x >= this->xMax())
                    return *(this->yBegin_+(this->xEnd_-this->xBegin_)-1);
                Size i = this->locate(x);
                return this->yBegin_[i] + (x-this->xBegin_[i])*s_[i];
            }
            Real primitive(Real x) const override {
                Size i = this->locate(x);
                Real dx = x-this->xBegin_[i];
                return primitiveConst_[i] +
                    dx*(this->yBegin_[i] + 0.5*dx*s_[i]);
            }
            Real derivative(Real x) const override {
                if (!this->isInRange(x))
                    return 0;
                Size i = this->locate(x);
                return s_[i];
            }
            Real secondDerivative(Real) const override { return 0.0; }

          private:
            std::vector<Real> primitiveConst_, s_;
        };
    }

}


#endif

Coverage Report

Created: 2025-08-11 06:28

Line	Count	Source (jump to first uncovered line)
1		/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3		/*
4		Copyright (C) 2010 SunTrust Bank
5		Copyright (C) 2010, 2014 Cavit Hafizoglu
6
7		This file is part of QuantLib, a free-software/open-source library
8		for financial quantitative analysts and developers - http://quantlib.org/
9
10		QuantLib is free software: you can redistribute it and/or modify it
11		under the terms of the QuantLib license. You should have received a
12		copy of the license along with this program; if not, please email
13		<quantlib-dev@lists.sf.net>. The license is also available online at
14		<http://quantlib.org/license.shtml>.
15
16		This program is distributed in the hope that it will be useful, but WITHOUT
17		ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18		FOR A PARTICULAR PURPOSE. See the license for more details.
19		*/
20
21		/*! \file generalizedhullwhite.hpp
22		\brief generalized Hull-White model
23		*/
24
25		#ifndef quantlib_generalized_hull_white_hpp
26		#define quantlib_generalized_hull_white_hpp
27
28		#include <ql/experimental/shortrate/generalizedornsteinuhlenbeckprocess.hpp>
29		#include <ql/math/interpolation.hpp>
30		#include <ql/models/shortrate/onefactormodel.hpp>
31		#include <ql/processes/ornsteinuhlenbeckprocess.hpp>
32		#include <utility>
33
34		namespace QuantLib {
35
36		//! Parameter that holds an interpolation object
37		class InterpolationParameter : public Parameter {
38		private:
39		class Impl final : public Parameter::Impl {
40		public:
41	0	Real value(const Array&, Time t) const override { return interpolator_(t); }
42	0	void reset(const Interpolation& interp) { interpolator_ = interp; }
43		private:
44		Interpolation interpolator_;
45		};
46		public:
47		explicit InterpolationParameter(
48		Size count,
49		const Constraint& constraint = NoConstraint())
50	0	: Parameter(count,
51	0	ext::shared_ptr<Parameter::Impl>(
52	0	new InterpolationParameter::Impl()),
53	0	constraint)
54	0	{ }
55	0	void reset(const Interpolation &interp) {
56	0	ext::shared_ptr<InterpolationParameter::Impl> impl =
57	0	ext::dynamic_pointer_cast<InterpolationParameter::Impl>(impl_);
58	0	if (impl != nullptr)
59	0	impl->reset(interp);
60	0	}
61		};
62
63		//! Generalized Hull-White model class.
64		/*! This class implements the standard Black-Karasinski model defined by
65		\f[
66		d f(r_t) = (\theta(t) - \alpha f(r_t))dt + \sigma dW_t,
67		\f]
68		where \f$ alpha \f$ and \f$ sigma \f$ are piecewise linear functions.
69
70		\ingroup shortrate
71		*/
72		class GeneralizedHullWhite : public OneFactorAffineModel,
73		public TermStructureConsistentModel {
74		public:
75
76		GeneralizedHullWhite(
77		const Handle<YieldTermStructure>& yieldtermStructure,
78		const std::vector<Date>& speedstructure,
79		const std::vector<Date>& volstructure,
80		const std::vector<Real>& speed,
81		const std::vector<Real>& vol,
82		const std::function<Real(Real)>& f = {},
83		const std::function<Real(Real)>& fInverse = {});
84
85		template <class SpeedInterpolationTraits,class VolInterpolationTraits>
86		GeneralizedHullWhite(
87		const Handle<YieldTermStructure>& yieldtermStructure,
88		const std::vector<Date>& speedstructure,
89		const std::vector<Date>& volstructure,
90		const std::vector<Real>& speed,
91		const std::vector<Real>& vol,
92		const SpeedInterpolationTraits &speedtraits,
93		const VolInterpolationTraits &voltraits,
94		const std::function<Real(Real)>& f = {},
95		const std::function<Real(Real)>& fInverse = {}) :
96		OneFactorAffineModel(2), TermStructureConsistentModel(yieldtermStructure),
97		speedstructure_(speedstructure), volstructure_(volstructure),
98		a_(arguments_[0]), sigma_(arguments_[1]),
99		f_(f), fInverse_(fInverse)
100		{
101		initialize(yieldtermStructure,speedstructure,volstructure,
102		speed,vol,speedtraits,voltraits,f,fInverse);
103		}
104
105	0	ext::shared_ptr<ShortRateDynamics> dynamics() const override {
106	0	QL_FAIL("no defined process for generalized Hull-White model, "
107	0	"use HWdynamics()");
108	0	}
109
110		ext::shared_ptr<Lattice> tree(const TimeGrid& grid) const override;
111
112		//Analytical calibration of HW
113
114		GeneralizedHullWhite(
115		const Handle<YieldTermStructure>& yieldtermStructure,
116		Real a = 0.1, Real sigma = 0.01);
117
118
119		ext::shared_ptr<ShortRateDynamics> HWdynamics() const;
120
121		//! Only valid under Hull-White model
122		Real discountBondOption(Option::Type type,
123		Real strike,
124		Time maturity,
125		Time bondMaturity) const override;
126
127		//! vector to pass to 'calibrate' to fit only volatility
128		std::vector<bool> fixedReversion() const;
129
130		protected:
131		//Analytical calibration of HW
132	0	Real a() const { return a_(0.0); }
133	0	Real sigma() const { return sigma_(0.0); }
134		void generateArguments() override;
135		Real A(Time t, Time T) const override;
136		Real B(Time t, Time T) const override;
137		Real V(Time t, Time T) const;
138
139		private:
140
141		class Dynamics;
142		class Helper;
143		class FittingParameter;// for analytic HW fitting
144
145		std::vector<Date> speedstructure_;
146		std::vector<Date> volstructure_;
147		std::vector<Time> speedperiods_;
148		std::vector<Time> volperiods_;
149		Interpolation speed_;
150		Interpolation vol_;
151
152		std::function<Real (Time)> speed() const;
153		std::function<Real (Time)> vol() const;
154
155		Parameter& a_;
156		Parameter& sigma_;
157		Parameter phi_;
158
159		std::function<Real(Real)> f_;
160		std::function<Real(Real)> fInverse_;
161
162	0	static Real identity(Real x) {
163	0	return x;
164	0	}
165
166		template <class SpeedInterpolationTraits,class VolInterpolationTraits>
167		void initialize(const Handle<YieldTermStructure>& yieldtermStructure,
168		const std::vector<Date>& speedstructure,
169		const std::vector<Date>& volstructure,
170		const std::vector<Real>& speed,
171		const std::vector<Real>& vol,
172		const SpeedInterpolationTraits &speedtraits,
173		const VolInterpolationTraits &voltraits,
174		const std::function<Real(Real)>& f,
175		const std::function<Real(Real)>& fInverse)
176	0	{
177	0	QL_REQUIRE(speedstructure.size()==speed.size(),
178	0	"mean reversion inputs inconsistent");
179	0	QL_REQUIRE(volstructure.size()==vol.size(),
180	0	"volatility inputs inconsistent");
181	0	if (!f_)
182	0	f_ = identity;
183	0	if (!fInverse_)
184	0	fInverse_ = identity;
185
186	0	DayCounter dc = yieldtermStructure->dayCounter();
187	0	Date ref = yieldtermStructure->referenceDate();
188	0	for (auto i : speedstructure)
189	0	speedperiods_.push_back(dc.yearFraction(ref, i));
190	0	for (auto i : volstructure)
191	0	volperiods_.push_back(dc.yearFraction(ref, i));
192
193		// interpolator x points to *periods_ vector, y points to
194		// the internal Array in the parameter
195	0	InterpolationParameter atemp(speedperiods_.size(), NoConstraint());
196	0	a_ = atemp;
197	0	for (Size i=0; i<speedperiods_.size(); i++)
198	0	a_.setParam(i, speed[i]);
199	0	speed_ = speedtraits.interpolate(speedperiods_.begin(),
200	0	speedperiods_.end(),a_.params().begin());
201	0	speed_.enableExtrapolation();
202	0	atemp.reset(speed_);
203
204	0	InterpolationParameter sigmatemp(volperiods_.size(), PositiveConstraint());
205	0	sigma_ = sigmatemp;
206	0	for (Size i=0; i<volperiods_.size(); i++)
207	0	sigma_.setParam(i, vol[i]);
208	0	vol_ = voltraits.interpolate(volperiods_.begin(),
209	0	volperiods_.end(),sigma_.params().begin());
210	0	vol_.enableExtrapolation();
211	0	sigmatemp.reset(vol_);
212
213	0	generateArguments();
214	0	registerWith(yieldtermStructure);
215	0	} Unexecuted instantiation: void QuantLib::GeneralizedHullWhite::initialize<QuantLib::LinearFlat, QuantLib::LinearFlat>(QuantLib::Handle<QuantLib::YieldTermStructure> const&, std::__1::vector<QuantLib::Date, std::__1::allocator<QuantLib::Date> > const&, std::__1::vector<QuantLib::Date, std::__1::allocator<QuantLib::Date> > const&, std::__1::vector<double, std::__1::allocator<double> > const&, std::__1::vector<double, std::__1::allocator<double> > const&, QuantLib::LinearFlat const&, QuantLib::LinearFlat const&, std::__1::function<double (double)> const&, std::__1::function<double (double)> const&) Unexecuted instantiation: void QuantLib::GeneralizedHullWhite::initialize<QuantLib::BackwardFlat, QuantLib::BackwardFlat>(QuantLib::Handle<QuantLib::YieldTermStructure> const&, std::__1::vector<QuantLib::Date, std::__1::allocator<QuantLib::Date> > const&, std::__1::vector<QuantLib::Date, std::__1::allocator<QuantLib::Date> > const&, std::__1::vector<double, std::__1::allocator<double> > const&, std::__1::vector<double, std::__1::allocator<double> > const&, QuantLib::BackwardFlat const&, QuantLib::BackwardFlat const&, std::__1::function<double (double)> const&, std::__1::function<double (double)> const&)
216		};
217
218		//! Short-rate dynamics in the generalized Hull-White model
219		/*! The short-rate is here
220
221		f(r_t) = x_t + g(t)
222
223		where g is the deterministic time-dependent
224		parameter (which can't be determined analytically)
225		used for initial term-structure fitting and x_t is the state
226		variable following an Ornstein-Uhlenbeck process.
227
228		In this version, the function f may also be defined as a piece-wise linear
229		function and can be calibrated to the away-from-the-money instruments.
230
231		*/
232		class GeneralizedHullWhite::Dynamics
233		: public GeneralizedHullWhite::ShortRateDynamics {
234		public:
235		Dynamics(Parameter fitting,
236		const std::function<Real(Time)>& alpha,
237		const std::function<Real(Time)>& sigma,
238		std::function<Real(Real)> f,
239		std::function<Real(Real)> fInverse)
240	0	: ShortRateDynamics(ext::shared_ptr<StochasticProcess1D>(
241	0	new GeneralizedOrnsteinUhlenbeckProcess(alpha, sigma))),
242	0	fitting_(std::move(fitting)), _f_(std::move(f)), _fInverse_(std::move(fInverse)) {}
243
244		//classical HW dynamics
245		Dynamics(Parameter fitting, Real a, Real sigma)
246		: GeneralizedHullWhite::ShortRateDynamics(
247		ext::shared_ptr<StochasticProcess1D>(new OrnsteinUhlenbeckProcess(a, sigma))),
248	0	fitting_(std::move(fitting)), _f_(identity()), _fInverse_(identity()) {}
249
250	0	Real variable(Time t, Rate r) const override { return _f_(r) - fitting_(t); }
251
252	0	Real shortRate(Time t, Real x) const override { return _fInverse_(x + fitting_(t)); }
253
254		private:
255		Parameter fitting_;
256		std::function<Real(Real)> _f_;
257		std::function<Real(Real)> _fInverse_;
258		struct identity {
259	0	Real operator()(Real x) const {return x;};
260		};
261		};
262
263		//! Analytical term-structure fitting parameter \f$ \varphi(t) \f$.
264		/*! \f$ \varphi(t) \f$ is analytically defined by
265		\f[
266		\varphi(t) = f(t) + \frac{1}{2}[\frac{\sigma(1-e^{-at})}{a}]^2,
267		\f]
268		where \f$ f(t) \f$ is the instantaneous forward rate at \f$ t \f$.
269		*/
270		class GeneralizedHullWhite::FittingParameter
271		: public TermStructureFittingParameter {
272		private:
273		class Impl final : public Parameter::Impl {
274		public:
275		Impl(Handle<YieldTermStructure> termStructure, Real a, Real sigma)
276	0	: termStructure_(std::move(termStructure)), a_(a), sigma_(sigma) {}
277
278	0	Real value(const Array&, Time t) const override {
279	0	Rate forwardRate =
280	0	termStructure_->forwardRate(t, t, Continuous, NoFrequency);
281	0	Real temp = a_ < std::sqrt(QL_EPSILON) ?
282	0	Real(sigma_*t) :
283	0	Real(sigma_(1.0 - std::exp(-a_t))/a_);
284	0	return (forwardRate + 0.5temptemp);
285	0	}
286
287		private:
288		Handle<YieldTermStructure> termStructure_;
289		Real a_, sigma_;
290		};
291		public:
292		FittingParameter(const Handle<YieldTermStructure>& termStructure,
293		Real a, Real sigma)
294	0	: TermStructureFittingParameter(ext::shared_ptr<Parameter::Impl>(
295	0	new FittingParameter::Impl(termStructure, a, sigma))) {}
296		};
297
298		// Analytic fitting dynamics
299		inline ext::shared_ptr<OneFactorModel::ShortRateDynamics>
300	0	GeneralizedHullWhite::HWdynamics() const {
301	0	return ext::shared_ptr<ShortRateDynamics>(
302	0	new Dynamics(phi_, a(), sigma()));
303	0	}
304
305		namespace detail {
306		template <class I1, class I2>
307		class LinearFlatInterpolationImpl;
308		}
309
310		//! %Linear interpolation between discrete points with flat extapolation
311		/! \ingroup interpolations /
312		class LinearFlatInterpolation : public Interpolation {
313		public:
314		/! \pre the \f$ x \f$ values must be sorted. /
315		template <class I1, class I2>
316		LinearFlatInterpolation(const I1& xBegin, const I1& xEnd,
317	0	const I2& yBegin) {
318	0	impl_ = ext::shared_ptr<Interpolation::Impl>(new
319	0	detail::LinearFlatInterpolationImpl<I1,I2>(xBegin, xEnd,
320	0	yBegin));
321	0	impl_->update();
322	0	}
323		};
324
325		//! %Linear-interpolation with flat-extrapolation factory and traits
326		/! \ingroup interpolations /
327		class LinearFlat {
328		public:
329		template <class I1, class I2>
330		Interpolation interpolate(const I1& xBegin, const I1& xEnd,
331	0	const I2& yBegin) const {
332	0	return LinearFlatInterpolation(xBegin, xEnd, yBegin);
333	0	}
334		static const bool global = false;
335		static const Size requiredPoints = 1;
336		};
337
338		namespace detail {
339		template <class I1, class I2>
340		class LinearFlatInterpolationImpl
341		: public Interpolation::templateImpl<I1,I2> {
342		public:
343		LinearFlatInterpolationImpl(const I1& xBegin, const I1& xEnd,
344		const I2& yBegin)
345	0	: Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin,
346	0	LinearFlat::requiredPoints),
347	0	primitiveConst_(xEnd-xBegin), s_(xEnd-xBegin) {}
348	0	void update() override {
349	0	primitiveConst_[0] = 0.0;
350	0	for (Size i=1; i<Size(this->xEnd_-this->xBegin_); ++i) {
351	0	Real dx = this->xBegin_[i]-this->xBegin_[i-1];
352	0	s_[i-1] = (this->yBegin_[i]-this->yBegin_[i-1])/dx;
353	0	primitiveConst_[i] = primitiveConst_[i-1]
354	0	+ dx(this->yBegin_[i-1] +0.5dx*s_[i-1]);
355	0	}
356	0	}
357	0	Real value(Real x) const override {
358	0	if (x <= this->xMin())
359	0	return this->yBegin_[0];
360	0	if (x >= this->xMax())
361	0	return *(this->yBegin_+(this->xEnd_-this->xBegin_)-1);
362	0	Size i = this->locate(x);
363	0	return this->yBegin_[i] + (x-this->xBegin_[i])*s_[i];
364	0	}
365	0	Real primitive(Real x) const override {
366	0	Size i = this->locate(x);
367	0	Real dx = x-this->xBegin_[i];
368	0	return primitiveConst_[i] +
369	0	dx(this->yBegin_[i] + 0.5dx*s_[i]);
370	0	}
371	0	Real derivative(Real x) const override {
372	0	if (!this->isInRange(x))
373	0	return 0;
374	0	Size i = this->locate(x);
375	0	return s_[i];
376	0	}
377	0	Real secondDerivative(Real) const override { return 0.0; }
378
379		private:
380		std::vector<Real> primitiveConst_, s_;
381		};
382		}
383
384		}
385
386
387		#endif