Generalized Hull White model with non-constant parameters
Locked
 
|
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----------------------- |
Re: Generalized Hull White model with non-constant parameters
|
|
On Tue, 2009-11-10 at 10:31 -0800, 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. Thanks, Javit. I'll add the files to the library as soon as I get some time. I've a few questions, namely: 1) are the files attached here the same as the ones in the Sourceforge patch tracker, or is one version more recent than the other? 2) it would be nice to have a test case for the new classes to add to the test suite. Do you think you can find some time to code it? 3) and finally, there's the usual matter of copyright attribution. Should the copyright be attributed to you or your company (if any?) And if you're not self-employed, are you sure that your company is ok with you contributing code? Later, Luigi -- Hanlon's Razor: Never attribute to malice that which is adequately explained by stupidity. ------------------------------------------------------------------------------ Let Crystal Reports handle the reporting - Free Crystal Reports 2008 30-Day trial. Simplify your report design, integration and deployment - and focus on what you do best, core application coding. Discover what's new with Crystal Reports now. http://p.sf.net/sfu/bobj-july _______________________________________________ QuantLib-dev mailing list [hidden email] https://lists.sourceforge.net/lists/listinfo/quantlib-dev |
|
|
1) The code is the same as in Sourceforge. Actually, in my sourceforge post, I provided a link (not the source code) to my original message in the quantlib-dev group.
For 3), I'm working on the copyright thing with our group manager. It's going to take sometime. Meanwhile, I would like to tidy up the code and write a good example case. For 2), there is a concern. Sadruddin Rejeb and/or StatPro Italia srl didn't refer to any material about the trinomial tree generation method used in the code. Comparing to the documents I have used in coding, the trinomial tree structures are different by observation. I need some information about the trinomial tree generation method in the trinomialtree class. Or instead, I may create another trinomialtree class. In either caser, it would be better to know more about the current trinomialtree class. Thank you, Javit
|
|
|
In reply to this post by javit
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 GeneralizedHW.zipGeneralizedHW.zipGeneralizedHW.zip |
«
Return to quantlib-dev
|
1 view|%1 views
| Free forum by Nabble | Edit this page |