Thread: [Quantlib-dev] Generalized Hull White model with non-constant parameters

SourceForge Headquarters 1320 Columbia Street Suite 310 San Diego, CA 92101 +1 (858) 422-6466

I modified the Black-Karasinski and the Ornstein-Uhlenbcek classes to build a
generalized Hull-White model with time-dependent drift and volatility
parameters. It consists of four files pasted at the bottom of this message: 
GeneralizedOUprocess.hpp, 
GeneralizedOUprocess.cpp, 
GeneralizedHW.hpp and
GeneralizedHW.cpp. 

These files need an addition/modification to the parameter class. Change the
NoConstraint() condition in PiecewiseConstantParameter class to
PositiveConstraint(). Or define another class in the parameter.hpp and name
it PiecewiseConstantParameter2. This is how I did it. I changed the code in
PiecewiseConstantParameter class

from --------
     public:
        PiecewiseConstantParameter(const std::vector<Time>& times)
        : Parameter(times.size()+1,
                    boost::shared_ptr<Parameter::Impl>(
                                 new
PiecewiseConstantParameter2::Impl(times)),
                    NoConstraint())
to ----------
     public:
        PiecewiseConstantParameter2(const std::vector<Time>& times)
        : Parameter(times.size(),
                    boost::shared_ptr<Parameter::Impl>(
                                 new
PiecewiseConstantParameter2::Impl(times)),
                    PositiveConstraint())

You will obviously nedd to make some changes to the paths in the include
commands.

I would like these to be added to the quantlib library in the future. Please
let me know what you think and also what I should do to make it a
non-experimental contribution.

Thank you,
Javit

----------- GeneralizedOUprocess.hpp --------------

#ifndef quantlib_generalized_ou_process_hpp
#define quantlib_generalized_ou_process_hpp

#include <ql/stochasticprocess.hpp>
#include <ql/models/shortrate/onefactormodel.hpp>
#include <ql/termstructures/yieldtermstructure.hpp>

#include <ql/termstructures/yield/zerocurve.hpp>
#include <ql/math/interpolations/linearinterpolation.hpp>
#include <ql/math/interpolations/loginterpolation.hpp>

namespace QuantLib {

    //! Piecewise linear Ornstein-Uhlenbeck process class
    /*! This class describes the Ornstein-Uhlenbeck process governed by
        \f[
            dx = a (level - x_t) dt + \sigma dW_t
        \f]

        \ingroup processes

		where the coefficients a and sigma are piecewise linear.
    */
    class GeneralizedOUprocess : public StochasticProcess1D {
      public:
		  GeneralizedOUprocess(const Handle<YieldTermStructure>& speedTS, 
			  const Handle<YieldTermStructure>& vol, 
			  Real x0 = 0.0, 
			  Real level = 0.0);
         //! \name StochasticProcess interface
        //@{
        Real x0() const;

        Real drift(Time t, Real x) const;
        Real diffusion(Time t, Real x) const;

        Handle<YieldTermStructure> speed() const;
        Handle<YieldTermStructure> volatility() const;
        Real level() const;

        Real expectation(Time t0, Real x0, Time dt) const;
        Real stdDeviation(Time t0, Real x0, Time dt) const;
        Real variance(Time t0, Real x0, Time dt) const;
        //@}
      private:
        Real x0_, level_;
		Handle<YieldTermStructure> speed_;
		Handle<YieldTermStructure> volatility_;
    };

}

#endif

----------- GeneralizedOUprocess.cpp --------------

#include "C:\Documents and Settings\****\My Documents\Visual Studio
2008\Projects\GeneralizedHW\GeneralizedHW\generalizedOUprocess.hpp"

namespace QuantLib {

	GeneralizedOUprocess::GeneralizedOUprocess(const
Handle<YieldTermStructure>& speedTS,
		const Handle<YieldTermStructure>& vol,
		Real x0,
		Real level)
		: x0_(x0), speed_(speedTS), level_(level), volatility_(vol) {

					QL_REQUIRE(x0 >= 0.0, "negative initial data given");
			QL_REQUIRE(level >= 0.0, "negative level given"); 
	}

	Real GeneralizedOUprocess::x0() const {
		return x0_;
	}

	Handle<YieldTermStructure> GeneralizedOUprocess::speed() const {
		return speed_;
	}

	Handle<YieldTermStructure> GeneralizedOUprocess::volatility() const {
		return volatility_;
	}

	Real GeneralizedOUprocess::drift(Time, Real x) const {
        return 0;
    }

	Real GeneralizedOUprocess::diffusion(Time, Real) const {
        return 0;
    }

	Real GeneralizedOUprocess::level() const {
		return level_;
	}

	Real GeneralizedOUprocess::expectation(const Time t, Real x0, Time dt)
const {

		Real speed;
			if (t!=0) {
				speed = - std::log(speed_->discount(t))/t;
			} else { speed = 0.000001;}

			return level_ + (x0 - level_) * std::exp(-speed*dt);
	}

	Real GeneralizedOUprocess::stdDeviation(Time t, Real x0,Time dt) const {
			return std::sqrt(variance(t,x0,dt));
	}

	Real GeneralizedOUprocess::variance(const Time t, Real x, Time dt) const {
			Real speed;
			Volatility vol;
		if (t!=0) {
				speed = - std::log(speed_->discount(t))/t;
				vol =  - std::log(volatility_->discount(t))/(t);
			} else {
				speed = 0.000001;
				vol = 0.00001;
		}

		if (speed < std::sqrt(QL_EPSILON)) {
			// algebraic limit for small speed
			return vol*vol*dt;
		} else {
			return 0.5*vol*vol/speed*
				(1.0 - std::exp(-2.0*speed*dt));
		}
	}

}

--------- GeneralizedHW.hpp -------------

#ifndef quantlib_GeneralizedHW_hpp
#define quantlib_GeneralizedHW_hpp

#include <ql/models/shortrate/onefactormodel.hpp>
#include "C:\Documents and Settings\****\My Documents\Visual Studio
2008\Projects\GeneralizedHW\GeneralizedHW\generalizedOUprocess.hpp"

namespace QuantLib {

	//! 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 GeneralizedHW : public OneFactorModel,
		public TermStructureConsistentModel {

		// TO DO: Build your fInverse class according to your model. The template
class T
		// should overload the () operator to return a function value, f inverse.

	public:

		GeneralizedHW(const Handle<YieldTermStructure>& yieldtermStructure,
			const std::vector<Date>& speedstructure,
			const std::vector<Date>& volstructure);

		boost::shared_ptr<ShortRateDynamics> dynamics() const {
			QL_FAIL("no defined process for generalized Hull-White model");
		} 

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

	private:

		class Dynamics;
		class Helper;

		std::vector<Date> speedstructure_;
		std::vector<Date> volstructure_;

		Handle<YieldTermStructure> speed() const;
		Handle<YieldTermStructure> vol() const;

		Parameter& a_;
		Parameter& sigma_;
		Parameter phi_;

	};

	//! 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 not be determined analytically)
	used for term-structure fitting and  x_t is the state
	variable following an Ornstein-Uhlenbeck process.
	*/

	class GeneralizedHW::Dynamics : public GeneralizedHW::ShortRateDynamics {
	public:
		Dynamics(const Parameter& fitting, const Handle<YieldTermStructure>&
alpha,
			const Handle<YieldTermStructure>& sigma)
			: ShortRateDynamics(boost::shared_ptr<StochasticProcess1D>(
			new GeneralizedOUprocess(alpha, sigma))), fitting_(fitting) {}

		Real variable(Time t, Rate r) const {
			return std::log(r) - fitting_(t);
		}

		Real shortRate(Time t, Real x) const {
			return std::exp(x + fitting_(t));
		}

	private:
		Parameter fitting_;	
	};

}

#endif

----------------- GeneralizedHW.cpp -------------
#include "C:\Documents and Settings\****\My Documents\Visual Studio
2008\Projects\GeneralizedHW\GeneralizedHW\GeneralizedHW.hpp"
#include <ql/methods/lattices/trinomialtree.hpp>
#include <ql/math/solvers1d/brent.hpp>
#include <ql/termstructures/yield/zerocurve.hpp>
#include <ql/math/solvers1d/bisection.hpp>

namespace QuantLib {

	/* Private function used by solver to determine time-dependent parameter
	df(r) = [theta(t) - a(t) f(r)]dt + sigma(t) dz
	dg = [theta(t) - a(t) g(t)] dt
	dx = -a(t) x dt + sigma(t) dz
	x = f(r) - g(t)

	*/
	//Change the overloaded operator to change the model by changing the
function below
	//fInverse is a user-chosen function. When fInverse = exp(), the model
becomes Black-Karasinski model.

	Real fInverse_(Real x) { 
		return std::exp(x);
	}

	class GeneralizedHW::Helper {
	public:
		Helper(const Size i, const Real xMin, const Real dx,
			const Real discountBondPrice,
			const boost::shared_ptr<ShortRateTree>& tree)
			: size_(tree->size(i)),
			dt_(tree->timeGrid().dt(i)),
			xMin_(xMin), dx_(dx),
			statePrices_(tree->statePrices(i)),
			discountBondPrice_(discountBondPrice){}

		Real operator()(const Real theta) const {
			Real value = discountBondPrice_;
			Real x = xMin_;
			for (Size j=0; j<size_; j++) {
				Real discount = std::exp(- fInverse_(theta+x)*dt_);
				//std::cout <<"fInverse  " <<fInverse_(theta+x)<<std::endl;
				value -= statePrices_[j]*discount;
				x += dx_;
			}

			return value;
		}

	private:
		Size size_;
		Time dt_;
		Real xMin_, dx_;
		const Array& statePrices_;
		Real discountBondPrice_;
	};

	GeneralizedHW::GeneralizedHW(const Handle<YieldTermStructure>&
yieldtermStructure,
		const std::vector<Date>& speedstructure,
		const std::vector<Date>& volstructure) : OneFactorModel(2),
TermStructureConsistentModel(yieldtermStructure),
		speedstructure_(speedstructure),
		volstructure_(volstructure),
		a_(arguments_[0]), sigma_(arguments_[1]){

			std::vector<Real>speedperiods;
			speedperiods.push_back(0.0);
			for (Size i=0;i<speedstructure.size()-1;i++)
				speedperiods.push_back(
				(speedstructure[i+1]-speedstructure[i])/365.0);

			a_ = PiecewiseConstantParameter2(speedperiods);

			std::vector<Real>volperiods;
			volperiods.push_back(0.0);
			for (Size i=0;i<volstructure.size()-1;i++)
				volperiods.push_back(
				(volstructure[i+1]-volstructure[i])/365.0);

			sigma_ = PiecewiseConstantParameter2(volperiods);

			a_.setParam(0,0.0001);
			sigma_.setParam(0,0.0001);
			for (Size i=1; i< a_.size();i++){

				a_.setParam(i,0.01*i);

			}

			for (Size i=1; i< sigma_.size();i++){

				sigma_.setParam(i,0.01*i);

			}

			registerWith(yieldtermStructure);

	}

	boost::shared_ptr<Lattice> GeneralizedHW::tree(const TimeGrid& grid) const{

		TermStructureFittingParameter phi(termStructure());
		boost::shared_ptr<ShortRateDynamics> numericDynamics(
			new Dynamics(phi, speed(), vol()));
		boost::shared_ptr<TrinomialTree> trinomial(
			new TrinomialTree(numericDynamics->process(), grid)); 
		boost::shared_ptr<ShortRateTree> numericTree(
			new ShortRateTree(trinomial, numericDynamics, grid));
		typedef TermStructureFittingParameter::NumericalImpl NumericalImpl;
		boost::shared_ptr<NumericalImpl> impl =
			boost::dynamic_pointer_cast<NumericalImpl>(phi.implementation());

		impl->reset();
		Real value = 1.0;
		Real vMin = -50.0;
		Real vMax = 50.0;

		extern std::vector<Real> shifts;

		for (Size i=0; i<(grid.size() - 1); i++) {
			vMin = -50.0;
			vMax = 50.0;
			Real discountBond = termStructure()->discount(grid[i+1]);
			Real xMin = trinomial->underlying(i, 0);
			Real dx = trinomial->dx(i);
			Helper finder(i, xMin, dx, discountBond, numericTree);
			value = 0.5*(vMin + vMax);
			Brent s1d;
			s1d.setMaxEvaluations(1000);
			value =s1d.solve(finder, QL_EPSILON, value, vMin, vMax);
			impl->set(grid[i], value);
			shifts.push_back(value);
		}
		return numericTree;

	}

	Handle<YieldTermStructure> GeneralizedHW::speed() const {

		std::vector<Real> speedvals;
		speedvals.push_back(0.000001);
		for (Size i=0;i<a_.size()-1;i++)
			speedvals.push_back(
			a_(
			(speedstructure_[i+1]-speedstructure_[i])/365.0
			- 0.00001));

		Handle<YieldTermStructure> speed_(boost::shared_ptr<YieldTermStructure>(
			new InterpolatedZeroCurve<LogLinear>(speedstructure_, speedvals,
Actual365Fixed())));

		return speed_; 
	}

	Handle<YieldTermStructure> GeneralizedHW::vol() const {

		std::vector<Real> volvals;
		volvals.push_back(0.000001);
		for (Size i=0;i<sigma_.size()-1;i++)
			volvals.push_back(
			sigma_(
			(speedstructure_[i+1]-speedstructure_[i])/365.0
			- 0.00001));

		Handle<YieldTermStructure> vol_(boost::shared_ptr<YieldTermStructure>(
			new InterpolatedZeroCurve<LogLinear>(volstructure_, volvals,
Actual365Fixed())));

		return vol_; 
	}

}
--------------- End of GeneralizedHW.cpp-----------------------

-----
Cavit (Javit) Hafizoglu
mailto:jav...@su... mailto:jav...@su... 
-- 
View this message in context: http://old.nabble.com/Generalized-Hull-White-model-with-non-constant-parameters-tp26287370p26287370.html
Sent from the quantlib-dev mailing list archive at Nabble.com.