Limits of Euler–Maruyama for Simulating CIR Interest Rate Paths
Summary
The document presents a Python function that generates paths for a Cox–Ingersoll–Ross (CIR) process. It updates the rate at each time step using the mean-reversion drift and a random shock scaled by the square root of the prior rate and the time increment. The answer identifies this approach as a basic Euler–Maruyama discretization and recommends considering alternative schemes for stochastic differential equations when the application requires greater numerical accuracy.
The response offers general guidance rather than checking whether the posted code is correct. In particular, it does not assess how the discretization handles the CIR process’s nonnegative boundary, compare schemes, or report simulation tests or error estimates. It also notes that more sophisticated simulation methods may be slow in Python and points to Julia’s differential-equation ecosystem as an option. The best scheme and implementation depend on the quantities being estimated, so the exchange is a starting point for method selection rather than a complete simulation recipe.
Key ideas
- The posted simulation uses an Euler–Maruyama time-step update for the CIR process.
- Euler–Maruyama is a basic discretization and may not suit every simulation objective.
- Alternative stochastic differential equation schemes can be considered for different accuracy needs.
- The answer does not verify the code or analyze how it preserves nonnegative rates.
Tags
Full text
# Does my Python code correctly simulate realizations of a CIR process?
# Does my Python code correctly simulate realizations of a CIR process?
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 by rvignolo (score 1)
https://quant.stackexchange.com/a/58925
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:
- 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.
- 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!Shown in full with attribution under the source's licence. Licence: CC BY-SA 4.0 (Stack Exchange)
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.