Interest rates are not stocks—they mean-revert. A 15% rate pulls back toward some long-run level; a negative rate pushes up. The right model encodes this pull while keeping rates from reaching unrealistically extreme negative values. (Mildly negative nominal rates have occurred in practice, but the Vasicek model's Gaussian tails can push rates far below zero.)
Vasicek model
The Vasicek model is the Ornstein-Uhlenbeck process applied to interest rates: dr = κ(θ − r)dt + σdW. The speed of mean reversion κ controls how fast rates snap back to the long-run mean θ. The drawback: rates can go negative because the noise term is additive.
importmathdef vasicek_step(r, kappa, theta, sigma, dt, z):
return r + kappa*(theta - r)*dt + sigma*math.sqrt(dt)*z
kappa, theta, sigma, dt = 0.5, 0.05, 0.02, 0.01
shocks = [0.3, -0.5, 0.1, -0.2, 0.4, -0.8, 0.2, -0.1, 0.6, -0.3]
r = 0.08
path = []
for z in shocks:
r = vasicek_step(r, kappa, theta, sigma, dt, z)
path.append(r)
print("Starting rate: 0.08")
print("Path:", " ".join(f"{x:.5f}"for x in path))
# Analytical expected valuedef vasicek_mean(r0, t):
return theta + (r0 - theta)*math.exp(-kappa*t)
print(f"E[r(1)] from r0=0.08: {vasicek_mean(0.08, 1.0):.6f}")
print(f"E[r(5)] from r0=0.08: {vasicek_mean(0.08, 5.0):.6f}")
CIR model
Cox-Ingersoll-Ross fixes the negativity problem: dr = κ(θ − r)dt + σ√r dW. The √r term shrinks volatility as rates approach zero, preventing them from going negative (provided 2κθ ≥ σ², the Feller condition). Same mean-reversion, better boundary behavior.
Scheme
; CIR model: dr = kappa*(theta - r)*dt + sigma*sqrt(r)*dW; Key difference from Vasicek: volatility scales with sqrt(r)
(define (cir-step r kappa theta sigma dt z)
(let ((new-r (+ r (* kappa (- theta r) dt)
(* sigma (sqrt (max r 0)) (sqrt dt) z))))
(max new-r 0))) ; floor at zero for numerical safety
(define kappa 0.5)
(define theta 0.05)
(define sigma 0.1)
(define dt 0.01)
; Feller condition: 2*kappa*theta >= sigma^2
(define feller (* 2 kappa theta))
(define sigma-sq (* sigma sigma))
(display "2*kappa*theta = ") (display feller) (newline)
(display "sigma^2 = ") (display sigma-sq) (newline)
(display "Feller satisfied? ") (display (>= feller sigma-sq)) (newline)
; Simulate from r0=0.01 (near zero โ this is where CIR matters)
(define shocks (list 0.5-1.20.8-0.31.0-0.60.20.9-0.40.1))
(define (simulate-cir r shocks)
(if (null? shocks) (list)
(let ((next (cir-step r kappa theta sigma dt (car shocks))))
(cons next (simulate-cir next (cdr shocks))))))
(define path (simulate-cir 0.01 shocks))
(display "CIR path from r0=0.01: ") (newline)
(for-each (lambda (r) (display (exact->inexact r)) (display " ")) path)
(newline)
(display "All non-negative? ")
(display (not (memq #f (map (lambda (r) (>= r 0)) path))))
Python
importmathdef cir_step(r, kappa, theta, sigma, dt, z):
new_r = r + kappa*(theta-r)*dt + sigma*math.sqrt(max(r,0))*math.sqrt(dt)*z
returnmax(new_r, 0)
kappa, theta, sigma, dt = 0.5, 0.05, 0.1, 0.01
feller = 2*kappa*theta
print(f"2*kappa*theta = {feller}, sigma^2 = {sigma**2}")
print(f"Feller satisfied? {feller >= sigma**2}")
shocks = [0.5, -1.2, 0.8, -0.3, 1.0, -0.6, 0.2, 0.9, -0.4, 0.1]
r = 0.01
path = []
for z in shocks:
r = cir_step(r, kappa, theta, sigma, dt, z)
path.append(r)
print("CIR path from r0=0.01:")
print(" ".join(f"{x:.5f}"for x in path))
print(f"All non-negative? {all(x >= 0 for x in path)}")
HJM framework
Heath-Jarrow-Morton models the entire forward rate curve instead of a single short rate. The forward rate f(t,T) is the rate locked in at time t for instantaneous borrowing at T. The drift is not free: no-arbitrage forces it to equal σ(t,T) times the integral of σ from t to T. This drift condition is the heart of HJM.
Scheme
; HJM drift condition: mu(t,T) = sigma(t,T) * integral_t^T sigma(t,s) ds; For constant volatility sigma, the drift is sigma^2 * (T-t)
(define (hjm-drift sigma t cap-t)
(* sigma sigma (- cap-t t)))
; Forward rate evolution: df(t,T) = mu(t,T)*dt + sigma*dW
(define (forward-rate-step f sigma t cap-t dt z)
(+ f (* (hjm-drift sigma t cap-t) dt) (* sigma (sqrt dt) z)))
(define sigma 0.01) ; forward rate volatility
(define dt 0.25) ; quarterly steps; Evolve 1-year forward rate over 4 quarters
(define shocks (list 0.5-0.30.8-0.2))
(define f0 0.05) ; initial 1-year forward rate
(define (evolve f t cap-t shocks)
(if (null? shocks) (list)
(let ((next (forward-rate-step f sigma t cap-t dt (car shocks))))
(cons next (evolve next (+ t dt) cap-t (cdr shocks))))))
(define path (evolve f0 01.0 shocks))
(display "Forward rate f(t,1) path: ") (newline)
(for-each (lambda (r) (display (exact->inexact r)) (display " ")) path)
(newline)
; Bond price: P(t,T) = exp(-integral_t^T f(t,s) ds); Approximate: P โ exp(-f * (T-t))
(define (bond-price f tau)
(exp (* (- f) tau)))
(display "Bond price P(0,1) at f=5%: ")
(display (exact->inexact (bond-price 0.051.0)))
Python
importmathdef hjm_drift(sigma, t, T):
return sigma**2 * (T - t)
def forward_rate_step(f, sigma, t, T, dt, z):
return f + hjm_drift(sigma, t, T)*dt + sigma*math.sqrt(dt)*z
sigma, dt = 0.01, 0.25
shocks = [0.5, -0.3, 0.8, -0.2]
f, t, T = 0.05, 0, 1.0
path = []
for z in shocks:
f = forward_rate_step(f, sigma, t, T, dt, z)
t += dt
path.append(f)
print("Forward rate f(t,1) path:")
print(" ".join(f"{x:.6f}"for x in path))
# Bond price approximation
bond_price = math.exp(-0.05 * 1.0)
print(f"Bond price P(0,1) at f=5%: {bond_price:.6f}")
Pricing caps and floors
A caplet pays max(r − K, 0) at maturity: insurance against rates rising above the strike K. A floorlet pays max(K − r, 0): insurance against rates falling. A cap is a strip of caplets; a floor is a strip of floorlets. In the Vasicek model, the distribution of future rates is normal, so caplet prices have Black-like closed forms.
With Vasicek dynamics, the price of a zero-coupon bond has a closed-form solution: P(t,T) = A(τ)e^(-B(τ)r), where τ = T−t, B(τ) = (1−e^(-κτ))/κ, and A depends on θ, σ, κ. This gives the entire yield curve from a single state variable r.