I've written the following function which should simulate realizations of a CIR process:

def cir_simulations(alpha: float, mu: float, sigma: float, delta_t: float, num_steps: int, num_sims: int):

    sims_shape = (num_steps, num_sims)
    sims = np.empty(sims_shape)
    std_norm_rand = np.random.standard_normal(sims_shape)

    prev_r = mu

    for r_t, rnd in zip(sims, std_norm_rand):
        delta_r = alpha * (mu - prev_r) * delta_t + sigma * np.sqrt(prev_r) * np.sqrt(delta_t) * rnd[:]
        r_t[:] = prev_r = prev_r + delta_r

    return sims

Can anyone confirm that my code will correctly simulate a CIR process?

As a secondary question, are there any improvements you'd make to this code?

  • $\begingroup$ I don't see any obvious errors. $\endgroup$ – noob2 Oct 25 '20 at 20:31
  • 2
    $\begingroup$ Your conditional distribution around zero will be very wrong if the Feller condition doesn't hold, and slightly wrong if it does. The question you need to ask yourself is what do you do if prev_r < 0. $\endgroup$ – river_rat Oct 26 '20 at 22:42

I would like to address your second question, namely:

are there any improvements you'd make to this code?

Yes, there are a lot of improvements, but it depends on what you want to compute at the end:

  1. Using a Euler-Maruyama scheme as a discretization scheme is the most basic choice. Maybe you can investigate about different discretization schemes. A good book to look into it is "Numerical Solution of Stochastic Differential Equations" by Kloeden and Platen.

  2. If you decide to implement a better algorithm for the SDE discretization, you will see that Python will have a really bad performance. In that sense, you might want to inspect Julia programming language and its DifferentialEquations.jl library (which is probably the best library for differential equations, including stochastic differential equations, in the world). This library (actually, one of its dependences, the StochasticDiffEq.jl library) has a tonne of methods for SDEs simulations. Learning Julia's syntax is really easy and you can achieve C performance using syntax similar to Python or MatLab.

Hope this helps!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.