Menu
▾
▴
[Quantlib-dev] Generalized Hull White model with non-constant parameters
|
From: javit <ca...@vi...> - 2009-11-10 18:31:14
|
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.
|
