Simulating Mean-Reverting Processes in Python The CIR Model for Interest Rates
While Geometric Brownian Motion (GBM) is standard for stock prices, variables like interest rates or volatility often exhibit mean reversion – a tendency to drift back towards a long-term average. They also typically cannot be negative. The Cox-Ingersoll-Ross (CIR) model, a type of Square-Root Diffusion (SRD) process, captures these features. This article demonstrates how to simulate the CIR process in Python using both Euler discretization and an exact method, highlighting its use in modeling interest rates.
The Square-Root Diffusion / CIR Model
The SRD process xt is defined by the Stochastic Differential Equation (SDE): \[d x_t=\kappa\left(\theta-x_t\right) d t+\sigma \sqrt{x_t} d W_t\]
- \(κ\): Speed of mean reversion (how quickly it returns to the average).
- \(θ\): Long-term mean level.
- \(σ\): Volatility parameter.
- \(x_t\): Ensures volatility decreases as \(x_t\) approaches zero, helping maintain non-negativity (under the Feller condition \(2κθ≥σ^2\)).
- \(dW_t\): Wiener process representing randomness.
When used for the short-term interest rate, this is known as the CIR model.
Setup
We’ll use standard libraries. Assume imports
like numpy as np, numpy.random as npr, math, time,
and matplotlib.pyplot as plt are done. We’ll reuse
the print_statistics function from the previous article if
needed for comparisons.
import numpy as np
import numpy.random as npr
import math
import time
from pylab import plt, mpl
import scipy.stats as scs # For print_statistics
# (Include print_statistics function definition here if running standalone)
# def print_statistics(a1, a2): ... (see previous article)
npr.seed(200) # Use a different seed for variety
np.set_printoptions(precision=4)
print("=== SRD/CIR Simulation Section ===")
# Parameters for CIR interest rate model
x0 = 0.05 # Initial interest rate (5%)
kappa = 3.0 # Mean reversion speed
theta = 0.02 # Long-term mean interest rate (2%)
sigma_srd = 0.1 # Volatility parameter
T = 2.0 # Time horizon (years)
M = 50 # Number of time steps per path
I = 10000 # Number of paths
dt = T / M # Time interval1. Euler-Maruyama Discretization
A common simulation approach is the Euler scheme, approximating the SDE over small steps \(Δt\): \(x_{t+\Delta t} \approx x_t+\kappa\left(\theta-x_t\right) \Delta t+\sigma \sqrt{\max \left(x_t, 0\right)} \sqrt{\Delta t} Z_t\) Using \(\max(xt,0)\) is crucial to avoid errors if \(x_t\) slightly undershoots zero due to discretization.
# Euler Discretization for SRD
def srd_euler(x0, kappa, theta, sigma_srd, T, M, I):
dt = T / M
xh = np.zeros((M + 1, I)) # Helper array (can go negative)
x = np.zeros_like(xh) # Actual process (non-negative)
xh[0] = x0
x[0] = x0
for t in range(1, M + 1):
# Euler step (max(x, 0) prevents issues with negative rates in sqrt)
xh[t] = (xh[t - 1] +
kappa * (theta - np.maximum(xh[t - 1], 0)) * dt +
sigma_srd * np.sqrt(np.maximum(xh[t - 1], 0)) *
math.sqrt(dt) * npr.standard_normal(I))
x = np.maximum(xh, 0) # Ensure non-negativity
return x
x1 = srd_euler(x0, kappa, theta, sigma_srd, T, M, I)
# Plot histogram of final values (Euler)
plt.figure(figsize=(10, 6))
plt.hist(x1[-1], bins=50, label='Final Value (Euler)')
plt.xlabel('value')
plt.ylabel('frequency')
plt.title('Final SRD Values (Euler Scheme)')
plt.legend()
plt.show()
# Plot first 10 simulated paths (Euler)
plt.figure(figsize=(10, 6))
plt.plot(np.linspace(0, T, M+1), x1[:, :10], lw=1.5)
plt.xlabel('time (years)')
plt.ylabel('value (e.g., interest rate)')
plt.title('Sample SRD Paths (Euler Scheme)')
plt.show()
-
Explanation: The function iterates through time steps,
calculating the next value based on the previous one, the drift
term \(κ(θ−xt)Δt\), and the diffusion
term \(\sigma \sqrt{\max \left(x_t, 0\right)}
\sqrt{\Delta t} Z_t\). A
final np.maximum(xh, 0) ensures the stored path stays
non-negative.
2. Exact Simulation Scheme (CIR)
The CIR process benefits from an exact simulation method based on the non-central chi-square (\(χ^2\)) distribution. This avoids the discretization error of the Euler scheme. The value \(xt+Δt\) conditional on \(x_t\) is given by \(c×Y\), where \(Y\) follows a non-central \(χ^2\) distribution with specific parameters:
# Exact Simulation for SRD (CIR model)
def srd_exact(x0, kappa, theta, sigma_srd, T, M, I):
dt = T / M
x = np.zeros((M + 1, I))
x[0] = x0
for t in range(1, M + 1):
df = 4 * theta * kappa / sigma_srd ** 2 # Degrees of freedom
c = (sigma_srd ** 2 * (1 - np.exp(-kappa * dt))) / (4 * kappa) # Scaling factor
nc = np.exp(-kappa * dt) / c * x[t - 1] # Non-centrality parameter
# Sample from non-central chi-square distribution
x[t] = c * npr.noncentral_chisquare(df, nc, size=I)
return x
x2 = srd_exact(x0, kappa, theta, sigma_srd, T, M, I)
# Plot histogram of final values (Exact)
plt.figure(figsize=(10, 6))
plt.hist(x2[-1], bins=50, label='Final Value (Exact)')
plt.xlabel('value')
plt.ylabel('frequency')
plt.title('Final SRD Values (Exact Scheme)')
plt.legend()
plt.show()
# Plot first 10 simulated paths (Exact)
plt.figure(figsize=(10, 6))
plt.plot(np.linspace(0, T, M+1), x2[:, :10], lw=1.5)
plt.xlabel('time (years)')
plt.ylabel('value (e.g., interest rate)')
plt.title('Sample SRD Paths (Exact Scheme)')
plt.show()
-
Explanation: This function calculates the parameters
(df, c, nc) at each step based on
the model parameters and the previous value \(x_t−1\). It then samples
from npr.noncentral_chisquare and scales the result to get
the exact value \(x_t\).
3. Comparison: Euler vs. Exact
Let’s compare the results and performance.
# Compare statistics of Euler vs Exact final values
print("\nComparing statistics of final SRD values (Euler vs Exact):")
print_statistics(x1[-1], x2[-1])
# Timing comparison
print("\nTiming SRD methods:")
I_large = 250000 # Larger number of paths for timing
start_time = time.time()
x1_timed = srd_euler(x0, kappa, theta, sigma_srd, T, M, I_large)
print(f"Euler Time: {time.time() - start_time:.4f} seconds")
start_time = time.time()
x2_timed = srd_exact(x0, kappa, theta, sigma_srd, T, M, I_large)
print(f"Exact Time: {time.time() - start_time:.4f} seconds")
print("\nComparing statistics for larger simulation:")
print_statistics(x1_timed[-1], x2_timed[-1])
x1 = None; x2 = None; x1_timed = None; x2_timed = None # Clear memoryComparing statistics of final SRD values (Euler vs Exact):
statistic data set 1 data set 2
---------------------------------------------
size 10000.000 10000.000
min 0.004 0.005
max 0.049 0.052
mean 0.020 0.020
std 0.006 0.006
skew 0.505 0.555
kurtosis 0.303 0.402
Timing SRD methods:
Euler Time: 0.7942 seconds
Exact Time: 1.4020 seconds
Comparing statistics for larger simulation:
statistic data set 1 data set 2
---------------------------------------------
size 250000.000 250000.000
min 0.003 0.003
max 0.069 0.058
mean 0.020 0.020
std 0.006 0.006
skew 0.564 0.574
kurtosis 0.517 0.485- Discussion: The statistics are generally close, but
the Euler scheme might show slight biases, especially with
larger
dt. The Euler scheme is often faster computationally, but the exact method eliminates discretization error. The choice depends on the required accuracy and computational constraints. For high accuracy or large time steps, the exact method is preferred if available.
Use Cases in Finance
- Interest Rate Derivatives: The CIR model is widely used to price bonds, swaps, caps, floors, swaptions, and other interest rate-sensitive instruments.
- Volatility Modeling: The variance process in the Heston stochastic volatility model follows a CIR process.
- Credit Risk: Modeling the intensity (hazard rate) of default events in some credit derivative models.
Conclusion
The CIR/SRD process provides a more realistic framework than GBM for modeling mean-reverting, non-negative variables like interest rates. Python allows straightforward simulation using either the approximate Euler scheme or the more accurate exact method. Understanding these simulation techniques is crucial for pricing and risk-managing interest rate derivatives and for building more advanced models like Heston.