Improve CIR Discretization

[mmencke]

Using the Quadratic Exponential scheme of Leif Andersen to avoid "explosions" in the short-rate when the square root process gets close to zero.

I have made a small program to showcase the issue. I have attached a plot of the terminal distribution of the discretized processes compared to the true distribution in the Cox-Ingersoll-Ross model. The plot on the left is of the approach that is in QuantLib today and the plot on the right is of my suggested code.

Run-time is given below:

Simulate Square Root: Run time = 592.542 ms
Quadratic Exponential: Run time = 713.863 ms

Code

//
//  main.cpp
//  ImproveCir
//
//  Created by Magnus Mencke on 10/06/2021.
//  Copyright © 2021 Magnus Mencke. All rights reserved.
//

#include <iostream>
#include <fstream>

#include <ql/quantlib.hpp>

using namespace QuantLib;

void writeCsv(const Matrix data, std::string file_name) {
    std::ofstream myfile;
    myfile.open (file_name);
    
    Size n =data.rows();
    Size m = data.columns();
    
    for(int i=0;i<n;i++) {
        for(int j=0;j<m;j++) {
            myfile << data[i][j];
            if(j+1<m) {
                myfile << ", ";
            }
        }
        myfile << "\n";
    }
    
    myfile.close();
}

class CirHelperProcess : public StochasticProcess1D {
public:
    CirHelperProcess( Real k, Real sigma, Real y0, Real theta)
    : 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_;
};

int main() {
    BigNatural seed = 1;
    Size nTimeSteps = 252;
    Size nSamples = 10000;
    Statistics statistics;
    
    bool brownianBridge = false;
    
    // Parameters from [Brigo, Morini and Pallavicini (2013), p. 125]
    Real kappa=0.4;
    Real mu = 0.026;
    Real v = 0.14;
    Real y0 = 0.0165;
    
    TimeGrid timeGrid(1,nTimeSteps);
    
    PseudoRandom::rsg_type rsg = PseudoRandom::make_sequence_generator(nTimeSteps, seed);
    
    Matrix cirDistribution(nSamples,2);
    
    auto start = std::chrono::steady_clock::now();
    auto end = std::chrono::steady_clock::now();
    auto diff = end - start;
    
    /*******************************************************************/
    /* SIMULATE SQUARE ROOT                                            */
    /*******************************************************************/
    
    start = std::chrono::steady_clock::now();
    
    ext::shared_ptr<StochasticProcess1D> cirProcessSqrt(new CirHelperProcess(kappa, v, y0, mu));
    
    ext::shared_ptr<SingleVariate<PseudoRandom>::path_generator_type> cirRsgSqrt(new SingleVariate<PseudoRandom>::path_generator_type(cirProcessSqrt,timeGrid, rsg, brownianBridge));
    
    for(int i=0;i<nSamples;i++) {
        Path path = (cirRsgSqrt->next()).value;
        cirDistribution[i][0] =std::pow(path[nTimeSteps],2);
    }
    
    end = std::chrono::steady_clock::now();
    
    diff = end - start;
    
    std::cout << "Simulate Square Root: Run time = " << std::chrono::duration <double, std::milli> (diff).count() << " ms" << std::endl;
    
    /*******************************************************************/
    /* QUADRATIC EXPONENTIAL SCHEME                                    */
    /*******************************************************************/
    
    start = std::chrono::steady_clock::now();
    
    ext::shared_ptr<StochasticProcess1D> cirProcessQe(new CoxIngersollRossProcess(kappa, v, y0, mu));
    
    ext::shared_ptr<SingleVariate<PseudoRandom>::path_generator_type> cirRsgQe(new SingleVariate<PseudoRandom>::path_generator_type(cirProcessQe,timeGrid, rsg, brownianBridge));
    
    for(int i=0;i<nSamples;i++) {
        Path path = (cirRsgQe->next()).value;
        cirDistribution[i][1] =path[nTimeSteps];
    }
    
    end = std::chrono::steady_clock::now();
    
    diff = end - start;
    
    std::cout << "Quadratic Exponential: Run time = " << std::chrono::duration <double, std::milli> (diff).count() << " ms" << std::endl;
    
    
    writeCsv(cirDistribution,"cirDistribution.csv");
    
    return 0;
}

[boring-cyborg]

Thanks for opening this pull request! It might take a while before we look at it, so don't worry if there seems to be no feedback. We'll get to it.

[CLAassistant]

All committers have signed the CLA.

[coveralls]

Coverage decreased (-0.003%) to 71.101% when pulling 383f260 on mmencke:master into 956d572 on lballabio:master.

[boring-cyborg]

Congratulations on your first merged pull request!

[lballabio]

Thank you!