Numerical SDE Methods for Interest Rate Dynamics

This post shows how to simulate short-rate models for interest rates and use them in Monte Carlo pricing/valuation. We cover discretization choices (Euler–Maruyama, Milstein, exact), correlation handling, and practical tips for stability and accuracy.

Scenarios:

Path‑dependent discounting for cashflows: $DF = \exp\big(-\int r_t dt\big)$.
Pricing/valuation under short‑rate models: Vasicek/OU, CIR, and a note on Hull–White.
Monthly or finer time steps with correlated factors.

Models

Vasicek (Ornstein–Uhlenbeck on $r$): $\mathrm{d}r = \kappa(\theta - r)\mathrm{d}t + \sigma \mathrm{d}W_t$.
- Mean‑reverting, Gaussian; can be negative; has exact discretization in discrete time.
CIR (Cox–Ingersoll–Ross): $\mathrm{d}r = \kappa(\theta - r) \mathrm{d}t + \sigma\sqrt{r} \mathrm{d}W_t$.
- Mean‑reverting, strictly positive if Feller: $2\kappa\theta \ge \sigma^2$; has exact transition via noncentral $\chi^2$.
Hull–White (time‑dependent Vasicek): $\mathrm{d}r = a(t)\big(b(t) - r\big) \mathrm{d}t + \sigma(t) \mathrm{d}W_t$.
- Matches an initial yield curve; can simulate with Euler or semi‑exact integrator.

Discretization methods

Euler–Maruyama (EM): strong order 0.5; simple and widely used.
Milstein: strong order 1.0 for scalar SDEs; adds a diffusion derivative term (useful for √r in CIR).
Exact schemes (when available):
- Vasicek/OU in discrete time has a closed form.
- CIR has an exact transition: $r_{t+\Delta} = c \cdot \chi^2_{\nu}(\lambda)$ for appropriate $(\nu, \lambda, c)$.

Guidance:

Prefer exact discretizations when available (Vasicek, CIR) for accuracy/stability.
For CIR with EM/Milstein, ensure positivity (e.g., full truncation, Alfonsi scheme) if not using exact sampling.
Use small $\Delta t$ for path‑dependent discounting; check convergence of results to $\Delta t \to 0$.

Correlation

For multiple factors (e.g., multi‑curve, FX joint modeling), correlate Brownian shocks at each time step:

Build correlation matrix $C$.
Compute its Cholesky $L$.
At each step, draw $u \sim \mathcal{N}(0, I)$ and set $Z = L u$. Use $Z$ components in each factor’s update.

Python examples

Below are compact NumPy examples for Vasicek (exact) and CIR (exact sampling). Replace parameters with calibrated values.

Exact Vasicek simulation and bond pricing

import numpy as np

# Horizon and grid
T_y = 5.0            # years
steps_per_year = 12
dt = 1.0/steps_per_year
T = int(T_y / dt)
N = 50_000           # paths

# Vasicek params
kappa = 0.8
theta = 0.02
sigma = 0.01
r0 = 0.03

# Precompute exact discrete coefficients
exp_kdt = np.exp(-kappa*dt)
mean_coef = theta * (1 - exp_kdt)
var_coef = (sigma**2) * (1 - exp_kdt**2) / (2*kappa)

r = np.full(N, r0)
accum_int = np.zeros(N)  # integral of r_t dt

for _ in range(T):
    z = np.random.normal(size=N)
    # Exact update: r_{t+dt} = theta + (r_t - theta)e^{-k dt} + sqrt(var) * z
    r = theta + (r - theta) * exp_kdt + np.sqrt(var_coef) * z
    accum_int += r * dt

# Discount factor to horizon (pathwise)
DF = np.exp(-accum_int)
price_mc = DF.mean()
print({'vasicek_mc_ZCB_price': float(price_mc)})

Note: Vasicek has a closed‑form zero‑coupon bond price. Use it to sanity‑check Monte Carlo.

Exact CIR simulation via noncentral chi‑square

import numpy as np
from numpy.random import default_rng
try:
    from scipy.stats import ncx2
    have_scipy = True
except Exception:
    have_scipy = False

T_y = 5.0
steps_per_year = 12
dt = 1.0/steps_per_year
T = int(T_y / dt)
N = 50_000

kappa = 1.2
theta = 0.03
sigma = 0.12
r0 = 0.02

gamma = 4 * kappa * theta / (sigma**2)  # degrees of freedom ν
c = (sigma**2 * (1 - np.exp(-kappa*dt))) / (4*kappa)

r = np.full(N, r0)
acc = np.zeros(N)
rng = default_rng()

for _ in range(T):
    if have_scipy:
        lam = (4 * kappa * np.exp(-kappa*dt) * r) / (sigma**2 * (1 - np.exp(-kappa*dt)))
        r = c * ncx2.rvs(df=gamma, nc=lam, size=N, random_state=rng)
    else:
        # Fallback: positivity‑preserving Milstein (approximate)
        z = rng.standard_normal(N)
        r = np.maximum(0.0, r + kappa*(theta - r)*dt + sigma*np.sqrt(np.maximum(r,0))*np.sqrt(dt)*z \
                        + 0.25*sigma**2*dt*(z**2 - 1))
    acc += r * dt

DF = np.exp(-acc)
price_mc = DF.mean()
print({'cir_mc_ZCB_price': float(price_mc)})

Tips:

Ensure Feller: 2κθ ≥ σ² for strict positivity. If violated, exact sampler still works; EM can misbehave.
Use variance reduction (antithetic variates, control variates with analytical bond price) to cut Monte Carlo noise.

Correlated two‑factor Vasicek (sketch)

import numpy as np

dt = 1/12
T = 120
N = 100_000

kap = np.array([0.8, 1.1])
th  = np.array([0.02, 0.015])
sig = np.array([0.010, 0.008])

rho = 0.5
Corr = np.array([[1.0, rho],[rho,1.0]])
L = np.linalg.cholesky(Corr)

exp_kdt = np.exp(-kap*dt)
var_coef = (sig**2) * (1 - np.exp(-2*kap*dt)) / (2*kap)

r = np.full((2,N), [0.03, 0.02])
acc = np.zeros((2,N))

for _ in range(T):
    u = np.random.normal(size=(2,N))
    Z = L @ u
    for i in (0,1):
        r[i] = th[i] + (r[i] - th[i]) * np.exp(-kap[i]*dt) + np.sqrt(var_coef[i]) * Z[i]
        acc[i] += r[i]*dt

short_rate = r[0] + r[1]  # example combination

Convergence and accuracy

Check stability: reduce $\Delta t$ until results (e.g., ZCB price) stop changing materially.
Use exact transitions where possible; otherwise prefer Milstein over EM for diffusions with state‑dependent vol.
Track strong vs weak error based on your goal (path accuracy vs distribution accuracy).

Calibration notes

Vasicek/Hull–White: fit $\kappa, \theta, \sigma$ (and $a(t), b(t)$) to match historical vol/mean and today’s curve. Hull–White matches the initial term structure by construction.
CIR: calibrate to caps/floors or the term structure of vol; respect positivity constraints.

Extensions

Time‑dependent parameters (Hull–White): implement exact OU with time‑varying drift using integrating factors, or keep small Δt and EM.
HJM/LMM for forward‑rate dynamics; Sobol sequences for quasi‑MC; American exercise (Longstaff–Schwartz) with short‑rate paths.

This provides a practical toolkit to simulate interest rates with SDEs and apply them in Monte Carlo valuation and risk.

Historical reconstruction for knowledge retention & reference; not investment advice.

Numerical SDE Methods for Interest Rate Dynamics#

Models#

Discretization methods#

Correlation#

Python examples#

Exact Vasicek simulation and bond pricing#

Exact CIR simulation via noncentral chi‑square#

Correlated two‑factor Vasicek (sketch)#

Convergence and accuracy#

Calibration notes#

Extensions#