Why We Model Rates as Random

The problem with deterministic rates

Under the annuity measure $\mathbb{Q}^A$ of Chapter 03 the forward swap rate is a martingale, so $\mathbb{E}^{A}\bigl[S_{\alpha,\beta}(T_\alpha)\bigr] = S_{\alpha,\beta}(0)$, which at an at-the-money quote is the strike: plugging the mean in returns $\max(K - K,\,0) = 0$ for every cell on that surface. Convexity of $\max(\cdot,0)$ forces $\mathbb{E}^{A}\bigl[\max(S_{\alpha,\beta}(T_\alpha) - K,\,0)\bigr] > 0$, and no bootstrap produces the difference.

The Wiener process

The simplest object that captures unpredictable, continuous movement is the Wiener process $W_t$. Its increments are independent across non-overlapping intervals and normally distributed:

The variance equals $dt$, not $dt^2$: shocks grow as the square root of time. This is why daily vol scales to annual vol by $\sqrt{252}$, and why longer-dated options are more expensive.

From noise to a short rate model

A Wiener process is just noise. A rate model adds a drift, where the rate is pulled, and a diffusion, how much noise surrounds that pull:

The choice of $\mu$ and $\sigma$ is the model. Mean reversion toward a long-run level $b$ at speed $a$, with constant volatility $\sigma$, gives the Vasicek model ([1], Section 3.2) under the risk-neutral measure $\mathbb{Q}$:

Above $b$ the drift pulls down, below $b$ it pulls up: Ornstein-Uhlenbeck applied to rates. Both constants have a unit you can feel. At $a = 0.15$ the pull closes half the gap to $b$ in $\ln 2 / a = 4.62$ years, and the fan stops widening once that pull balances the noise, at $\sigma/\sqrt{2a} = 1.46\%$; five years reaches $\sigma\sqrt{(1 - e^{-2aT})/2a} = 1.29\%$.

The Vasicek bond price formula

Vasicek's advantage is a closed form for the bond price. Risk-neutral pricing gives

Under Vasicek $r$ is Markov under $\mathbb{Q}$ with constant $(a, b, \sigma)$ and $r_T \mid r_t$ is Gaussian, so $\int_t^T r_s\,ds$ is a limit of sums of Gaussians and Gaussian too. The expectation of a lognormal is exponential-affine:

$B(t,T) = \int_0^{T-t} e^{-au}\,du$ and $-\partial \ln P(t,T)/\partial r_t = B(t,T)$: the bond's sensitivity to today's short rate, accumulated over the remaining life and damped by the mean reversion (derivation in [1], Section 3.2.1). It also lets the martingale test below reprice analytically at any simulated $r_s$, with no second Monte Carlo nested inside the first.

At $a = 0.15$, $b = 4\%$, $\sigma = 0.80\%$, $r_0 = 4.33\%$:

$B(0,5) = \frac{1 - e^{-0.15 \times 5}}{0.15} = 3.5176$, and $\ln A(0,5) = \bigl(b - \frac{\sigma^2}{2a^2}\bigr)(B - T) - \frac{\sigma^2}{4a}B^2 = -0.0585$, giving $P(0,5) = e^{-0.0585} \times e^{-3.5176 \times 0.0433} = 0.8099$.

Monte Carlo pricing

Start at $r_0$, draw increments, accumulate $\exp(-\int_0^T r_s\,ds)$ along the path, and average the discounted payoff:

The exponential factor is the path-wise stochastic discount factor, kept in integral form throughout to avoid confusion with the deterministic curve $P(0,T)$. $h(r_T)$ is the payoff at maturity.

The martingale test

Test 1: discount factor alone

Can the simulation reproduce today's $P(0,T)$ from the stochastic discount factor alone? Under $\mathbb{Q}$ the path-wise factor must average to the deterministic price:

The identity does not care which short-rate model you feed it, so this runs on the engine that ships: xvafoundations.models.HullWhite1F on the Chapter 01 SOFR curve, every path starting at that curve's own $r(0) = 4.3282\%$.

Then plot the difference, not the two curves. On a bond-price axis running 97.9% down to 66.7% a 10bp error is a tenth of a pixel: an overlay cannot fail visibly, so it is not a test.

Why the accumulator is not a trapezoid on the short rate

Exact simulation buys an unbiased short rate. It does not buy an unbiased money-market account on a finite grid. Write $r_t = x_t + \alpha(t)$: the deterministic $\alpha$ carries $f(0,t)$, which log-linear interpolation of $\ln P(0,T)$ leaves piecewise constant with a jump at every curve pillar. Trapezoid $r$ itself and each panel straddling a jump misses $\tfrac{1}{2}\Delta t$ times it, so the rule drops to first order, and the straddling-panel count is the pillar count, which does not shrink with $\Delta t$.

The two pillars inside the first year drop $f$ by 13.7bp and 56.2bp, and $0.025 \times 69.8 = 1.75$bp predicts the 1.74bp of accumulation lost, or the 1.67bp in the first column once multiplied by $P(0,1) = 0.9596$. Halving the step halves each column rather than quartering it: first order, and no path count touches it.

So the library does not accumulate that way. $\int_0^t f(0,u)\,du = -\ln P(0,t)$ holds by definition of the instantaneous forward, under every interpolation and every pillar set, so $\int_0^t \alpha$ is closed form, that identity plus a Hull-White convexity term, and only the smooth $x$ is quadratured. No panel straddles a pillar. On the same seven maturities at 200 steps the residual bias is at most 0.0003bp, and it falls by exactly 4.000 at each of four halvings: second order, measured. The trapezoid survives under the name money_market_account_trapezoid, where the model-risk sweep still measures it; there the switch took the worst bias from 11.39bp to 0.012bp.

Test 2: the tower property

The stronger check: at any intermediate $s$, the path-wise discount factor times the conditional bond price must average back to today's $P(0,T)$. That is the tower property of conditional expectation, applied to pricing:

Here $P(s,T)$ is the Vasicek analytical price at the simulated $r_s$. On a path where rates rose you discount harder, but the bond you still hold is cheaper by the offsetting amount, so the product does not move: the two effects are the same effect seen twice.

Python implementation

The teaching version of both tests: exact Vasicek step, no library dependency, 50,000 paths. Its Test 1 targets the Vasicek analytical $P(0,T)$, not the Chapter 01 curve, which constant-parameter Vasicek cannot reproduce; the chart above ran the same identity on HullWhite1F, which can.

import math

import torch

torch.set_default_dtype(torch.float64)


class VasicekModel:
    """Vasicek short rate under Q:  dr_t = a (b - r_t) dt + sigma dW_t."""

    def __init__(self, r0=0.0433, a=0.15, b=0.04, sigma=0.008):
        self.r0, self.a, self.b, self.sigma = r0, a, b, sigma

    def bond_price(self, r_t, t, T):
        """Analytical Vasicek ZCB price P(t, T) given r_t."""
        tau = T - t
        B = (1.0 - math.exp(-self.a * tau)) / self.a
        ln_A = ((self.b - self.sigma ** 2 / (2 * self.a ** 2)) * (B - tau)
                - self.sigma ** 2 * B ** 2 / (4 * self.a))
        return torch.exp(ln_A - B * torch.as_tensor(r_t))

    def simulate(self, T, n_steps, n_paths, seed=137):
        """Exact Vasicek step. Returns (n_paths, n_steps + 1) short rates and
        the trapezoidal integral of r accumulated along each path, matching
        xvafoundations.pricing.money_market_account_trapezoid."""
        torch.manual_seed(seed)
        dt = T / n_steps
        decay = math.exp(-self.a * dt)
        shift = self.b * (1.0 - decay)
        vol = self.sigma * math.sqrt(
            (1.0 - math.exp(-2.0 * self.a * dt)) / (2.0 * self.a))

        rates = torch.zeros(n_paths, n_steps + 1)
        rates[:, 0] = self.r0
        for i in range(n_steps):
            rates[:, i + 1] = rates[:, i] * decay + shift + vol * torch.randn(n_paths)

        integrals = torch.zeros_like(rates)
        integrals[:, 1:] = torch.cumsum(
            0.5 * (rates[:, :-1] + rates[:, 1:]) * dt, dim=1)
        return rates, integrals


model = VasicekModel()
T_max, n_steps = 10.0, 400
rates, integrals = model.simulate(T=T_max, n_steps=n_steps, n_paths=50_000)

# Test 1: E^Q[exp(-int_0^T r ds)] = P(0,T).  Constant-parameter Vasicek has a
# smooth drift, so the trapezoid stays second order here and the residual at
# 400 steps is Monte Carlo noise, not bias.  Point the same rule at a
# bootstrapped curve and it drops to first order, because f(0,t) jumps at every
# pillar: that is the table above.  The library avoids it by integrating alpha
# in closed form and quadraturing only x.
for T in [0.5, 1.0, 2.0, 3.0, 5.0, 7.0, 10.0]:
    k = round(T * n_steps / T_max)
    sim = torch.exp(-integrals[:, k]).mean()
    err_bp = (sim - model.bond_price(model.r0, 0.0, T)).item() * 1e4
    print(f"T={T:4.1f}Y  sim={sim:.6f}  err={err_bp:+.2f}bp")

# Test 2: E^Q[exp(-int_0^s r du) * P(s,5)] = P(0,5)
P05 = model.bond_price(model.r0, 0.0, 5.0)
for s in [0.5, 1.0, 1.5, 2.0, 3.0, 4.0]:
    k = round(s * n_steps / T_max)
    avg = (torch.exp(-integrals[:, k]) * model.bond_price(rates[:, k], s, 5.0)).mean()
    print(f"s={s:4.1f}Y  E[disc*P]={avg:.6f}  err={(avg - P05).item() * 1e4:+.2f}bp")

Why this matters for XVA

Every exposure number in the XVA series averages over paths this machinery produces, so understating rate volatility understates exposure and therefore CVA. XVA-01 turns the paths into an exposure profile; XVA-03 derives the CVA integral from it.