Three regimes

A short rate whose mean level and volatility move among three regimes that switch quickly.

The model

The short rate $x_t$ is an Ornstein–Uhlenbeck process whose parameters are set by a hidden regime $y_t \in \{1, 2, 3\}$:

$$dx_t = \kappa\,(\theta_{y_t} - x_t)\,dt + \sigma_{y_t}\,dW_t .$$

The regime is a continuous-time Markov chain with generator $s\,Q_1$,

$$Q_1 = \begin{pmatrix} -3 & 2 & 1 \\ 1 & -2 & 1 \\ 0.5 & 1.5 & -2 \end{pmatrix},$$

so that from regime 1 it jumps to regime 2 at rate $2s$ and to regime 3 at rate $s$, and so on. The chain is not reversible.

There are two time scales: the rate reverts towards its current level over $1/\kappa$, and the regime changes every $1/s$ or so. Switching is fast when $s \gg \kappa$.

Parameters: $\kappa = 2$, $\theta = (0.20, 0.08, 0.02)$ and $\sigma^2 = (0.04, 0.01, 0.0025)$.

The quantity

The zero-coupon bond price

$$u_i(t,x) = \mathbb{E}\big[e^{-\int_0^t x_r\,dr} \mid x_0 = x,\ y_0 = i\big]$$

for maturity $t$, which is also the survival probability when $x$ is a default intensity.

Averaging

When the regime switches infinitely fast the rate sees only the stationary mixture. The chain's stationary distribution is $\pi = (0.208, 0.458, 0.333)$, and the limit is Vasicek's model with mean level $\bar\theta = \pi\cdot\theta$ and variance $\bar s = \pi\cdot\sigma^2$.

Reduction to a linear system

By the Feynman–Kac formula the prices $u_i(t,x)$, one for each starting regime, solve the coupled equations

$$\begin{aligned} \partial_t u_i = {}& \kappa(\theta_i - x)\,\partial_x u_i + \tfrac12\sigma_i^2\,\partial_{xx}u_i - x\,u_i \\ & + \sum_j Q_{ij}\,u_j , \end{aligned}$$

starting from $u_i(0,x) = 1$. The last term is the regime switching: it couples the equation for regime $i$ to the others at the switching rates. Try the Vasicek form $u_i = e^{-B(t)x}a_i(t)$. Then $\partial_x u_i = -B\,u_i$ and $\partial_{xx}u_i = B^2u_i$, and after dividing by $e^{-Bx}$ the terms proportional to $x$ are

$$-B'\,x = \kappa B\,x - x .$$

They cancel for every regime when $B' = 1 - \kappa B$, that is $B(t) = (1 - e^{-\kappa t})/\kappa$. Because $\kappa$ is the same in every regime, one $B$ works for all of them. What remains involves $t$ only:

$$a' = \big(s\,Q_1 + \operatorname{diag} g(t)\big)\,a, \quad a(0) = \mathbf 1, \qquad g_i(t) = -\kappa\theta_i B(t) + \tfrac12\sigma_i^2 B(t)^2 .$$

The expansion

Write $Q = s\,Q_1$ and $\bar g = \pi\cdot g$, and let $Q^{\#}$ be the group inverse of $Q$: the matrix that inverts $Q$ on vectors with $\pi\cdot v = 0$ and returns such vectors. To first order in $1/s$,

$$a_i(t) \approx \exp\Big(\int_0^t \big[\bar g - \pi\cdot\big(g\,Q^{\#}(g - \bar g)\big)\big]\Big)\,\Big(1 - \big[Q^{\#}(g - \bar g)\big]_i(t)\Big).$$

Higher orders follow from one recursion. With $s_\pi = \pi\cdot a$ and $w = a/s_\pi - \mathbf 1$, each order of $w$ applies $Q^{\#}$ once to terms already known, and an initial layer, solved in the eigenbasis of $Q$, restores $a(0) = \mathbf 1$. The engine page gives the recursion in full.

The bond price written out

For a chain of any size the expansion is written with the group inverse $Q^{\#}$ of the generator and its stationary law $\pi$. Let $u = \theta - \pi\cdot\theta$ and $v = \sigma^2 - \pi\cdot\sigma^2$ be the fluctuations of level and variance across regimes, so the forcing splits as

$$g_i = \bar g + \tilde g_i, \qquad \bar g = -\kappa\,(\pi\cdot\theta)\,B + \tfrac12(\pi\cdot\sigma^2)\,B^2, \qquad \tilde g_i = -\kappa\,u_i\,B + \tfrac12\,v_i\,B^2 .$$

First order

The first correction is the Green–Kubo integral. Since $\tilde g$ is linear in $u$ and $v$, it is a combination of $I_2$, $I_3$ and $I_4$:

$$\begin{aligned} u_i(T, x) = \;&\exp\Big(-B\,x - \kappa(\pi\cdot\theta)\,I_1 + \tfrac12(\pi\cdot\sigma^2)\,I_2 + c_2 I_2 + c_3 I_3 + c_4 I_4\Big) \\ &\times\Big(1 + \kappa B\,[Q^{\#}u]_i - \tfrac12 B^2\,[Q^{\#}v]_i\Big) + O(|Q|^{-2}), \end{aligned}$$ $$\begin{aligned} c_2 &= \kappa^2\,\mathcal K(u, u), \qquad c_3 = -\kappa\,\mathcal K_{\mathrm{sym}}(u, v), \qquad c_4 = \tfrac14\,\mathcal K(v, v), \\ \mathcal K(f, h) &= -\pi\cdot\big(f\,Q^{\#}h\big), \qquad \mathcal K_{\mathrm{sym}}(u, v) = \tfrac12\big(\mathcal K(u, v) + \mathcal K(v, u)\big), \\ I_k &= \int_0^T B^k = \frac{1}{\kappa^k}\Big(T + \sum_{j=1}^{k}\binom{k}{j}(-1)^j\,\frac{1 - e^{-j\kappa T}}{j\kappa}\Big), \qquad B = \frac{1 - e^{-\kappa T}}{\kappa}. \end{aligned}$$

Second order

Because $B' = 1 - \kappa B$, every vector in the recursion is a polynomial in $B$:

$$\begin{aligned} w_1 &= -Q^{\#}\tilde g, \\ w_2 &= Q^{\#}\big(w_1' - \tilde g\circ w_1 + \pi\cdot(\tilde g\circ w_1)\,\mathbf 1\big), \qquad w_1' = (1 - \kappa B)\,\partial_B w_1 . \end{aligned}$$

For every fixed $T>0$, the second-order exponent is $\int_0^T \pi\cdot(\tilde g\circ w_2)$, a polynomial in $B$ integrated term by term, and the bracket gains $w_2(T)$:

$$u_i(T, x) = \exp\Big(\cdots + \sum_{k} d_k\,I_k\Big)\Big(1 + \kappa B\,[Q^{\#}u]_i - \tfrac12 B^2\,[Q^{\#}v]_i + [w_2(T)]_i\Big) + O(|Q|^{-3}),\qquad T>0\ \text{fixed},$$ $$\pi\cdot(\tilde g\circ w_2) = \sum_k d_k\,B^k .$$

The remainder is $O(|Q|^{-3})$ at fixed $T > 0$. At $T = 0$ the bracket is $1 + [w_2(0)]_i$ with $w_2(0) = \kappa\,(Q^{\#})^2 u$, and an initial layer of size $|Q|^{-2}$, decaying at the rates of $Q$, restores $u_i(0,x) = 1$; the engine includes it.

The numbers

For the chain above multiplied by 8, at $T = 1$:

$$\pi = (0.2083, 0.4583, 0.3333), \qquad Q^{\#} = \begin{pmatrix}-0.0265 & 0.0126 & 0.0139 \\ 0.0048 & -0.0187 & 0.0139 \\ 0.0100 & 0.0178 & -0.0278\end{pmatrix} .$$

quantity	value
$\pi\cdot\theta$, $\pi\cdot\sigma^2$	0.085000, 0.013750
$c_2$, $c_3$, $c_4$	0.000630208, -6.48438e-05, 1.72119e-06
$B$; $I_1, \dots, I_4$	0.432332; 0.283834, 0.0951891, 0.0341266, 0.0126963
averaged exponent	-0.047597
Green–Kubo exponent $c_2I_2 + c_3I_3 + c_4I_4$	5.77979e-05
$d_1, \dots, d_6$	-2.42e-05, 5.21e-05, -9.29e-06, 6.11e-07, -2.24e-08, 4.49e-10
second-order exponent $\sum d_kI_k$	-2.223e-06
memory $\kappa B[Q^{\#}u]_i - \frac12B^2[Q^{\#}v]_i$	-0.0033837, -0.000228824, 0.00242945
$w_2(T)$	4.363e-05, 3.338e-06, -3.186e-05

start	averaged	first order	second order	numerical
regime 1	0.95351767	0.95034618	0.95038566	0.95038809
regime 2	0.95351767	0.95335458	0.95335564	0.95335647
regime 3	0.95351767	0.95588943	0.95585693	0.95585386

The second-order formula agrees with the engine's outer series to rounding, and its error against the numerical solution is 3.1e-06. A maturity-uniform formula also includes the engine's initial layer; the initial-layers page proves that correction explicitly in the symmetric two-state case.

Results

The error in $a_1(1)$, starting in regime 1, after each order, against a 30-digit numerical solution:

scale s	order 0	order 1	order 2	order 3	order 4	order 5	order 6
4	6.2e-3	1.8e-4	2.1e-5	3.7e-6	5.1e-7	1.0e-7	1.6e-8
8	3.1e-3	4.2e-5	2.4e-6	2.2e-7	1.4e-8	1.4e-9	1.1e-10
16	1.6e-3	1.0e-5	2.9e-7	1.3e-8	4.2e-10	2.2e-11	7.9e-13

Each doubling of $s$ divides the order-$n$ error by about $2^{n+1}$.

A Monte Carlo check samples 400,000 paths of the three-state chain exactly, from exponential holding times. Given its path the integrated rate is Gaussian, so each path gives the bond price in closed form. At $s = 4$ and $s = 8$ the Monte Carlo mean of $a_1(1)$ agrees with the numerical solution within two standard errors. Certificate: verify_option_mc.py.