Short-Rate Interest Rate Models › Simulation of Short-Rate Paths
Quantitative FinanceAll chapters Short-Rate Interest Rate Models › Simulation of Short-Rate PathsOne section · Sections in this chapter6 / 13 ↓
Reading tools Definitions & inline links
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Let's implement both models and visualize their behavior. We'll start by simulating multiple paths to understand how mean reversion shapes rate dynamics.
import numpy as np
# Set random seed for reproducibility
np.random.seed(42)
# Model parameters
kappa = 0.5 # Speed of mean reversion
theta = 0.05 # Long-run mean (5%)
sigma = 0.02 # Volatility
r0 = 0.02 # Initial rate (2%)
# Simulation parameters
T = 10 # Time horizon (years)
n_steps = 2500 # Number of time steps
n_paths = 1000 # Number of simulation paths
dt = T / n_steps
# Time grid
t = np.linspace(0, T, n_steps + 1)Vasicek Simulation
The Vasicek model admits two approaches to simulation. The simpler approach uses Euler-Maruyama discretization which approximates the continuous-time SDE with discrete time steps. The discretization follows the general pattern of replacing with , with , and with where is standard normal:
The components of this discrete recursion are:
- rate at the next time step
- : time step size
- : standard normal random variable,
- : current rate
- : Vasicek parameters
Since the Vasicek model has an exact solution, we can also simulate using the exact transition density. This approach is superior because it produces samples from exactly the correct distribution regardless of time step size, eliminating discretization error. The implementation below uses this exact simulation method.
def simulate_vasicek(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None):
"""Simulate Vasicek model paths using exact transition density."""
if seed is not None:
np.random.seed(seed)
dt = T / n_steps
rates = np.zeros((n_paths, n_steps + 1))
rates[:, 0] = r0
for i in range(n_steps):
# Exact mean and variance of r(t+dt) given r(t)
exp_decay = np.exp(-kappa * dt)
mean = rates[:, i] * exp_decay + theta * (1 - exp_decay)
variance = (sigma**2 / (2 * kappa)) * (1 - np.exp(-2 * kappa * dt))
std = np.sqrt(variance)
# Sample from normal distribution
rates[:, i + 1] = mean + std * np.random.randn(n_paths)
return rates
# Simulate Vasicek paths
vasicek_paths = simulate_vasicek(
r0, kappa, theta, sigma, T, n_steps, n_paths, seed=42
)CIR Simulation
For the CIR model, exact simulation is more complex due to the non-central chi-squared transition distribution. While exact methods exist (involving transformations to chi-squared random variables), practitioners often use the Euler-Maruyama scheme with a reflection or truncation to handle potential negative values:
Each term in this discretization corresponds to its continuous-time counterpart:
- standard normal random variable
- : diffusion term dependent on current rate level
- : rate at next step
- : current rate
- : CIR parameters
- : time step size
The challenge with this discretization is that the right-hand side can become negative even though the continuous-time process stays positive. This occurs because the discrete approximation does not perfectly capture the behavior of the square-root diffusion near zero. We apply a full truncation scheme where we take before computing the square root, and again after adding the shock. This is one of several approaches in the literature; alternatives include reflection schemes and implicit methods.
def simulate_cir(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=None):
"""Simulate CIR model paths using Euler-Maruyama with truncation."""
if seed is not None:
np.random.seed(seed)
dt = T / n_steps
rates = np.zeros((n_paths, n_steps + 1))
rates[:, 0] = r0
for i in range(n_steps):
r_current = np.maximum(
rates[:, i], 0
) # Truncation for numerical stability
drift = kappa * (theta - r_current) * dt
diffusion = sigma * np.sqrt(r_current * dt) * np.random.randn(n_paths)
rates[:, i + 1] = r_current + drift + diffusion
rates[:, i + 1] = np.maximum(rates[:, i + 1], 0) # Ensure non-negative
return rates
# Check Feller condition parameters
feller_lhs = 2 * kappa * theta
feller_rhs = sigma**2
condition_met = feller_lhs >= feller_rhs
# Simulate CIR paths
cir_paths = simulate_cir(r0, kappa, theta, sigma, T, n_steps, n_paths, seed=42)Feller condition: 2κθ = 0.0500 >= σ² = 0.0004 Condition satisfied: True
The Feller condition is satisfied with these parameters, so the CIR model should keep rates strictly positive.
Visualizing Path Behavior
Let's plot sample paths from both models to see mean reversion in action.


The models show clear mean reversion. Paths starting at 2% tend to drift upward toward the 5% long-run mean. The key visual difference is that some Vasicek paths dip below zero, while CIR paths remain positive throughout.
Distribution Analysis
Let's examine the distribution of rates at the 10-year horizon.

# Calculate statistics
vasicek_terminal = vasicek_paths[:, -1]
cir_terminal = cir_paths[:, -1]
# Theoretical mean for both models
exp_decay = np.exp(-kappa * T)
theoretical_mean = r0 * exp_decay + theta * (1 - exp_decay)
# Vasicek statistics
v_mean = vasicek_terminal.mean()
v_std = vasicek_terminal.std()
v_min = vasicek_terminal.min()
v_neg_count = (vasicek_terminal < 0).sum()
v_neg_pct = (vasicek_terminal < 0).mean() * 100
# CIR statistics
c_mean = cir_terminal.mean()
c_std = cir_terminal.std()
c_min = cir_terminal.min()Terminal Rate Statistics (T = 10 years): -------------------------------------------------- Theoretical mean: 4.9798% Vasicek: Simulated mean: 5.0281% Std deviation: 2.0835% Minimum rate: -0.5242% Negative paths: 5 (0.5%) CIR: Simulated mean: 4.9922% Std deviation: 0.4670% Minimum rate: 3.8295%
Both models converge toward the theoretical mean of approximately 5%, as expected from the mean-reversion dynamics. The Vasicek distribution is symmetric and normal, while the CIR distribution shows slight right-skewness and is bounded below by zero.