Monte Carlo Simulation
Exploring the stochastic nature of short-term interest rates.
Model Parameters
Adjust variables to recalculate the model.
Zero-Coupon Bond Yield Curve
The theoretical term structure for zero-coupon bonds based on the affine property of the Vasicek model.
Statistical Probability
Likelihood of the rate outcome using the Ornstein-Uhlenbeck transition density
Expected Rate
-
Uncertainty
-
95% Confidence
-
The Mathematics
The core logic driving the simulation.
1. The Vasicek Model SDE
The Vasicek model describes the evolution of interest rates as a Stochastic Differential Equation (SDE). It captures the tension between random market shocks and a long-term economic equilibrium.
$$dr_t = a(b - r_t)dt + \sigma dW_t$$
- $r_t$: The instantaneous short rate at time $t$.
- $a$: Speed of Reversion. The gravitational pull back to equilibrium.
- $b$: Long-Term Mean. The equilibrium level.
- $\sigma$: Volatility. The magnitude of the random shocks.
- $dW_t$: Wiener Process (Brownian Motion). $dW_t \sim N(0, dt)$.
2. The Ornstein-Uhlenbeck Solution
This SDE is a specific case of the Ornstein-Uhlenbeck process and can be solved explicitly using Itō's Lemma:
$$r_t = r_0 e^{-at} + b(1 - e^{-at}) + \sigma \int_0^t e^{-a(t-s)} dW_s$$
- Deterministic Part: $r_0 e^{-at} + b(1 - e^{-at})$ (This creates the "Expected Path" line).
- Stochastic Part: $\sigma \int_0^t e^{-a(t-s)} dW_s$ (This creates the "random noise" in the simulation).
3. Bond Pricing (Affine Term)
The model produces a closed-form solution for bond prices. The price of a Zero-Coupon Bond $P(t,T)$ maturing at time $T$ is:
$$P(t, T) = A(t, T) e^{-B(t, T) r_t}$$
This relationship implies that the Yield Curve is determined entirely by the short rate $r_t$ and the model parameters.
The Implementation
The Python code constructing the models.
1. Euler-Maruyama Discretisation
Since continuous time ($dt$) is impossible to model on a discrete computer, we approximate the Vasicek SDE using the Euler-Maruyama method.
def simulate_vasicek(r0, a, b, sigma, T=1.0, dt=0.01, num_sims=1):
.
.
.
# Scale random shocks by sqrt(dt) as variance grows linearly with time.
shocks = sigma * np.sqrt(dt) * np.random.normal(size=(N, num_sims))
# Iterate through time steps to apply the Euler-Maruyama update rule.
for t in range(1, N):
drift = a * (b - rates[t-1, :]) * dt
rates[t, :] = rates[t-1, :] + drift + shocks[t, :]
return time, rates2. Probability Distribution
To generate the probability distributions, we do not need to run simulations. We use the statistical properties of the normal distribution derived directly from the model:
def calculate_future_distribution(r0, a, b, sigma, t):
expected_mean = r0 * np.exp(-a * t) + b * (1 - np.exp(-a * t))
# Variance (derived from Ito Isometry)
variance = (sigma**2 / (2 * a)) * (1 - np.exp(-2 * a * t))
std_dev = np.sqrt(variance)
return expected_mean, std_dev