Does my Python code correctly simulate realizations of a CIR process?

Question

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?

Answer 1

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!