From 383f2603861960cfd64b9a8658f5234d7f16b8c3 Mon Sep 17 00:00:00 2001 From: mmencke <57640398+mmencke@users.noreply.github.com> Date: Thu, 10 Jun 2021 23:53:54 +0200 Subject: [PATCH] Improve CIR Discretization Using the Quadratic Exponential scheme of Leif Andersen to avoid "explosions" in the short-rate when the square root process gets close to zero. --- .../onefactormodels/coxingersollross.hpp | 41 ++++-------------- .../extendedcoxingersollross.cpp | 3 +- .../extendedcoxingersollross.hpp | 10 ++--- ql/processes/coxingersollrossprocess.hpp | 42 ++++++++++++++++++- 4 files changed, 54 insertions(+), 42 deletions(-) diff --git a/ql/models/shortrate/onefactormodels/coxingersollross.hpp b/ql/models/shortrate/onefactormodels/coxingersollross.hpp index 0c27c311a77..1843e48ef1b 100644 --- a/ql/models/shortrate/onefactormodels/coxingersollross.hpp +++ b/ql/models/shortrate/onefactormodels/coxingersollross.hpp @@ -25,8 +25,7 @@ #define quantlib_cox_ingersoll_ross_hpp #include -#include -#include +#include namespace QuantLib { @@ -58,6 +57,7 @@ namespace QuantLib { ext::shared_ptr tree(const TimeGrid& grid) const override; + class Dynamics; protected: Real A(Time t, Time T) const override; @@ -70,7 +70,6 @@ namespace QuantLib { private: class VolatilityConstraint; - class HelperProcess; Parameter& theta_; Parameter& k_; @@ -78,33 +77,10 @@ namespace QuantLib { Parameter& r0_; }; - class CoxIngersollRoss::HelperProcess : public StochasticProcess1D { - public: - HelperProcess(Real theta, Real k, Real sigma, Real y0) - : y0_(y0), theta_(theta), k_(k), sigma_(sigma) { - discretization_ = - ext::shared_ptr(new EulerDiscretization); - } - - Real x0() const override { return y0_; } - Real drift(Time, Real y) const override { - return (0.5*theta_*k_ - 0.125*sigma_*sigma_)/y - - 0.5*k_*y; - } - Real diffusion(Time, Real) const override { return 0.5 * sigma_; } - - private: - Real y0_, theta_, k_, sigma_; - }; - //! %Dynamics of the short-rate under the Cox-Ingersoll-Ross model - /*! The state variable \f$ y_t \f$ will here be the square-root of the - short-rate. It satisfies the following stochastic equation - \f[ - dy_t=\left[ - (\frac{k\theta }{2}+\frac{\sigma ^2}{8})\frac{1}{y_t}- - \frac{k}{2}y_t \right] d_t+ \frac{\sigma }{2}dW_{t} - \f]. + /*! The short-rate is simulated directly using the Quadratic Exponential + discretization scheme as described in Leif Andersen, + Efficient Simulation of the Heston Stochastic Volatility Model. */ class CoxIngersollRoss::Dynamics : public OneFactorModel::ShortRateDynamics { @@ -114,14 +90,13 @@ namespace QuantLib { Real sigma, Real x0) : ShortRateDynamics(ext::shared_ptr( - new HelperProcess(theta, k, sigma, std::sqrt(x0)))) {} + new CoxIngersollRossProcess(k, sigma, x0, theta))) {} - Real variable(Time, Rate r) const override { return std::sqrt(r); } - Real shortRate(Time, Real y) const override { return y * y; } + Real variable(Time, Rate r) const override { return r; } + Real shortRate(Time, Real y) const override { return y; } }; } #endif - diff --git a/ql/models/shortrate/onefactormodels/extendedcoxingersollross.cpp b/ql/models/shortrate/onefactormodels/extendedcoxingersollross.cpp index bf3a59086bc..fbd792b5ab2 100644 --- a/ql/models/shortrate/onefactormodels/extendedcoxingersollross.cpp +++ b/ql/models/shortrate/onefactormodels/extendedcoxingersollross.cpp @@ -28,7 +28,7 @@ namespace QuantLib { Real theta, Real k, Real sigma, Real x0, bool withFellerConstraint) : CoxIngersollRoss(x0, theta, k, sigma, withFellerConstraint), - TermStructureConsistentModel(termStructure) { + TermStructureConsistentModel(termStructure){ generateArguments(); } @@ -102,4 +102,3 @@ namespace QuantLib { } } - diff --git a/ql/models/shortrate/onefactormodels/extendedcoxingersollross.hpp b/ql/models/shortrate/onefactormodels/extendedcoxingersollross.hpp index 04b9ccb0b92..af8d3f1be8f 100644 --- a/ql/models/shortrate/onefactormodels/extendedcoxingersollross.hpp +++ b/ql/models/shortrate/onefactormodels/extendedcoxingersollross.hpp @@ -75,11 +75,10 @@ namespace QuantLib { //! Short-rate dynamics in the extended Cox-Ingersoll-Ross model /*! The short-rate is here \f[ - r_t = \varphi(t) + y_t^2 + r_t = \varphi(t) + y_t \f] where \f$ \varphi(t) \f$ is the deterministic time-dependent - parameter used for term-structure fitting and \f$ y_t \f$ is the - state variable, the square-root of a standard CIR process. + parameter used for term-structure fitting and \f$ y_t \f$ is a standard CIR process. */ class ExtendedCoxIngersollRoss::Dynamics : public CoxIngersollRoss::Dynamics { @@ -87,8 +86,8 @@ namespace QuantLib { Dynamics(Parameter phi, Real theta, Real k, Real sigma, Real x0) : CoxIngersollRoss::Dynamics(theta, k, sigma, x0), phi_(std::move(phi)) {} - Real variable(Time t, Rate r) const override { return std::sqrt(r - phi_(t)); } - Real shortRate(Time t, Real y) const override { return y * y + phi_(t); } + Real variable(Time t, Rate r) const override { return r - phi_(t); } + Real shortRate(Time t, Real y) const override { return y + phi_(t); } private: Parameter phi_; @@ -153,4 +152,3 @@ namespace QuantLib { #endif - diff --git a/ql/processes/coxingersollrossprocess.hpp b/ql/processes/coxingersollrossprocess.hpp index de08e680c65..9b19f416abb 100644 --- a/ql/processes/coxingersollrossprocess.hpp +++ b/ql/processes/coxingersollrossprocess.hpp @@ -25,6 +25,7 @@ #define quantlib_coxingersollross_process_hpp #include +#include namespace QuantLib { @@ -38,6 +39,7 @@ namespace QuantLib { */ class CoxIngersollRossProcess : public StochasticProcess1D { public: + CoxIngersollRossProcess(Real speed, Volatility vol, Real x0 = 0.0, @@ -53,7 +55,10 @@ namespace QuantLib { Real volatility() const; Real level() const; Real variance(Time t0, Real x0, Time dt) const override; - + Real evolve (Time t0, + Real x0, + Time dt, + Real dw) const override; private: Real x0_, speed_, level_; Volatility volatility_; @@ -95,6 +100,41 @@ namespace QuantLib { return std::sqrt(variance(t,x0,dt)); } + inline Real CoxIngersollRossProcess::evolve (Time t0, + Real x0, + Time dt, + Real dw) const { + Real result; + // This is the QuadraticExponential scheme + // for details see Leif Andersen, + // Efficient Simulation of the Heston Stochastic Volatility Model + const Real ex = std::exp(-speed_*dt); + + const Real m = level_+(x0-level_)*ex; + const Real s2 = x0*volatility_*volatility_*ex/speed_*(1-ex) + + level_*volatility_*volatility_/(2*speed_)*(1-ex)*(1-ex); + const Real psi = s2/(m*m); + + if (psi <= 1.5) { + const Real b2 = 2/psi-1+std::sqrt(2/psi*(2/psi-1)); + const Real b = std::sqrt(b2); + const Real a = m/(1+b2); + + result = a*(b+dw)*(b+dw); + } + else { + const Real p = (psi-1)/(psi+1); + const Real beta = (1-p)/m; + + const Real u = CumulativeNormalDistribution()(dw); + + result = ((u <= p) ? 0.0 : std::log((1-p)/(1-u))/beta); + } + + + return result; + } + } #endif