Menu
▾
▴
Re: [Quantlib-dev] Generalized Hull White model with non-constant parameters
|
From: javit <ca...@vi...> - 2009-12-04 16:22:23
|
I cleared some bugs in my code and modified the BermudanSwaption.cpp test
file to include a test case. I uploaded the new source files with this
message. The generalized HW model has a better fit to market data as
expected.
The attached zip file also contains the filled and signed QuantLib copyright
template.
There are some minor things need to address:
1- LevenberMarquardt algorithm doesn't work for the generalizedHW class, but
Conjugate Gradient algorithm works. This is strange. I appreciate any help
in understanding this issue.
2- fInverse function is hardcoded to the exponential function, for now.
Therefore, the results are close to the Black Karasinski model, the only
exception is that the parameters can be made time-dependent and can be as
many as one wishes.
I'll deal with this problem in the next version.
3- I abused the YieldTermStructure class a little to define the
time-dependent parameters. Any suggestions/help are welcome.
4- Due to the structure of the YieldTermStructure class, model parameters
should be non-negative. Therefore, the additions, I mentioned below, to the
parameter class are necessary to use these classes.
5- I used the 9.7 version of Quantlib and boost version 1.39 to run my
tests. The code should work with the newer version.
Please let me know if there is anything else I could do. I look forward to
reading your comments.
Thank you,
Javit
javit wrote:
>
> 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-----------------------
>
>
http://old.nabble.com/file/p26635940/GeneralizedHW.zip GeneralizedHW.zip
http://old.nabble.com/file/p26635940/GeneralizedHW.zip GeneralizedHW.zip
http://old.nabble.com/file/p26635940/GeneralizedHW.zip GeneralizedHW.zip
-----
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-tp26287370p26635940.html
Sent from the quantlib-dev mailing list archive at Nabble.com.
|
