Improve CIR Discretization

ql/models/shortrate/onefactormodels/coxingersollross.hpp

@@ -25,8 +25,7 @@
 #define quantlib_cox_ingersoll_ross_hpp
 
 #include <ql/models/shortrate/onefactormodel.hpp>
-#include <ql/stochasticprocess.hpp>
-#include <ql/processes/eulerdiscretization.hpp>
+#include <ql/processes/coxingersollrossprocess.hpp>
 
 namespace QuantLib {
 
@@ -58,6 +57,7 @@ namespace QuantLib {
 
         ext::shared_ptr<Lattice> tree(const TimeGrid& grid) const override;
 
+
         class Dynamics;
       protected:
         Real A(Time t, Time T) const override;
@@ -70,41 +70,17 @@ namespace QuantLib {
 
       private:
         class VolatilityConstraint;
-        class HelperProcess;
 
         Parameter& theta_;
         Parameter& k_;
         Parameter& sigma_;
         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<discretization>(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<StochasticProcess1D>(
-                        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
-

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 {
     }
 
 }
-

ql/models/shortrate/onefactormodels/extendedcoxingersollross.hpp

@@ -75,20 +75,19 @@ 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 {
       public:
         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
-

ql/processes/coxingersollrossprocess.hpp

@@ -25,6 +25,7 @@
 #define quantlib_coxingersollross_process_hpp
 
 #include <ql/stochasticprocess.hpp>
+#include <ql/math/distributions/normaldistribution.hpp>
 
 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