Short-Rate Interest Rate Models › Simulation of Short-Rate Paths

Michael BrenndoerferDecember 6, 20257 min read
Quantitative FinanceAll chapters
Part 7 of 15Interest Rate Models and Derivatives
Chapter 38Short-Rate Interest Rate Models › Simulation of Short-Rate Paths
Reading section 6 of 13 · Simulation of Short-Rate Paths · 7 minSection 6 of 13 · 7 min

Book contents

Section 6 of 13 · 7 min read

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.

In[8]:
Code
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 drdr with Δr\Delta r, dtdt with Δt\Delta t, and dWdW with ΔtZ\sqrt{\Delta t} Z where ZZ is standard normal:

rt+dt=rt+κ(θ−rt)dt+σdt Zr_{t+dt} = r_t + \kappa(\theta - r_t)dt + \sigma\sqrt{dt}\,Z

The components of this discrete recursion are:

  • rt+dtr_{t+dt} rate at the next time step
  • dtdt: time step size
  • ZZ: standard normal random variable, Z∼N(0,1)Z \sim \mathcal{N}(0,1)
  • rtr_t: current rate
  • κ,θ,σ\kappa, \theta, \sigma: 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.

In[9]:
Code
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:

rt+dt=rt+κ(θ−rt)dt+σrtdt Zr_{t+dt} = r_t + \kappa(\theta - r_t)dt + \sigma\sqrt{r_t}\sqrt{dt}\,Z

Each term in this discretization corresponds to its continuous-time counterpart:

  • ZZ standard normal random variable
  • rt\sqrt{r_t}: diffusion term dependent on current rate level
  • rt+dtr_{t+dt}: rate at next step
  • rtr_t: current rate
  • κ,θ,σ\kappa, \theta, \sigma: CIR parameters
  • dtdt: 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 max⁡(rt,0)\max(r_t, 0) 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.

In[10]:
Code
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)
Out[11]:
Console
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.

Out[12]:
Visualization
Line plot showing 50 Vasicek interest rate paths reverting to 5% over 10 years.
Vasicek model simulation with 50 sample paths showing mean reversion toward θ = 5%. The normal distribution of shocks allows rates to become negative, as observed in several paths crossing below zero.
Line plot showing 50 CIR interest rate paths reverting to 5% over 10 years.
CIR model simulation with 50 sample paths showing mean reversion toward θ = 5%. The square-root diffusion process ensures rates remain non-negative, with volatility decreasing as rates approach zero.

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.

Out[13]:
Visualization
Histogram comparing Vasicek and CIR rate distributions at 10-year horizon.
Terminal rate distributions at T=10 years. Vasicek produces a symmetric normal distribution centered near θ, while CIR produces a right-skewed distribution bounded at zero.
In[14]:
Code
# 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()
Out[15]:
Console
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.