Regime-switching survival

A hazard rate whose mean level and volatility switch between two regimes, expanded to all orders in the inverse switching rate.

The problem comes from the second half of my 2001 Stanford thesis. The hazard rate, or short rate, is an Ornstein–Uhlenbeck process whose mean level and volatility switch between two regimes at rate $\lambda$.

The quantity wanted is the survival probability, which is also a zero-coupon bond price. When the regimes switch infinitely fast the answer is Vasicek's formula with the regime averages plugged in. The corrector adds terms in $\varepsilon = 1/\lambda$.

The expansion

The hazard follows

where $y_t \in \{1,2\}$ switches at rate $\lambda$ in each direction. The survival probability is

Write

With $\varepsilon = 1/\lambda$, the upper sign when the process starts in regime 1, and $t > 0$ fixed as $\lambda \to \infty$,

The factor $e^{-B(t)x} e^{\int_0^t \bar g}$ is Vasicek's bond price with the averaged parameters. The $\varepsilon$ term in $\log M$ is the same for both regimes and grows with maturity. The factor $1 \pm R$ depends on the starting regime, and it vanishes as $\lambda \to \infty$.

The expansion can be checked against a numerical solution. Because $\kappa$ does not switch, $u_y(t,x) = a_y(t)\,e^{-B(t)x}$, where $a$ solves a two-dimensional linear ODE driven by $\bar g \pm \tilde g$ and the switching generator. This reduction is known: Elliott and Mamon (2002) solve the case of a switching mean level, and Elliott and Siu (2009) give Markov-modulated exponential-affine bond prices in general.

Against the numerical solution, the error after zero, one and two orders falls like $\varepsilon$, $\varepsilon^2$ and $\varepsilon^3$. The convergence demo shows the three slopes, and the survival demo shows the curves. The published fast-switching expansions treat option prices and transition densities, and a search of the literature did not find the corrector for this bond price or survival probability.

The remainder is $O(\varepsilon^3)$ at each fixed $t > 0$, not uniformly down to $t = 0$. The outer series does not satisfy $u_y(0,x) = 1$: at $t = 0$ the bracket is $1 \mp \tfrac14\varepsilon^2\tilde g'(0)$, with $\tilde g'(0) = -\kappa\tilde\theta$, because the regime has not yet had time to switch.

An initial layer repairs the start. To leading order it adds to $R(t)$ the transient

which is of order $\varepsilon^2$ for $t$ up to a few multiples of $1/(2\lambda)$ and falls below $\varepsilon^3$ only once $t$ exceeds about $\log(1/\varepsilon)/(2\lambda)$. With the transient included, the remainder is $O(\varepsilon^3)$ uniformly on bounded intervals of $t$.

The survival probability written out

Every integral in $\log M$ is a combination of $\int_0^T B^k$, and $B$ is a polynomial in $E = e^{-\kappa T}$. The integrals are therefore elementary, and the survival probability to second order is

This is Vasicek's survival probability with averaged parameters, multiplied by an explicit correction. The $\varepsilon$ term in the exponent is the Green–Kubo correction, and the bracket is the memory of the starting regime. Its last term is the matched initial layer; it makes the stated remainder uniform for $T\geq0$ on every fixed bounded maturity interval.

At the parameters of the survival demo, $\kappa = 2$, $\theta = (0.15,\ 0.02)$, $\sigma^2 = (0.0555,\ 0.0055)$, $x = 0.12$, with $T = 4$, $\lambda = 20$ and a start in regime 1:

Derivation of the first-order formula

The idea page explains the same two corrections without calculus: one is half the variance of the fluctuating hazard, the other is the memory of the starting regime.

The coupled equations

By the Feynman–Kac formula the two survival probabilities, one per starting regime, solve

both starting from $u_i(0,x) = 1$. The last terms are the switching: over a short interval $dt$ the regime jumps with probability $\lambda\,dt$, and the expectation then continues from the other regime.

Removing the rate

Try $u_i = e^{-B(t)x}\,a_i(t)$. Each derivative in $x$ brings down a factor of $-B$, and the terms proportional to $x$ cancel when $B' = 1 - \kappa B$. Because $\kappa$ is the same in both regimes, one $B$ works for both, and what is left is a pair of ordinary differential equations:

with $g_i = -\kappa\theta_i B + \tfrac12\sigma_i^2B^2$. Solving this pair numerically gives the reference used throughout the examples.

Mean and ratio

Write the two solutions through their mean and their normalized difference:

Adding the two equations gives the mean, and the switching terms cancel:

Subtracting them, and dividing by $m$, gives the ratio:

The quasi-steady balance

In the equation for $\omega$ the term $-2\lambda\,\omega$ is large. Switching pulls the two regimes back together at rate $2\lambda$, while $\tilde g$ pushes them apart. When $\lambda$ is large, $\omega$ settles quickly near the value that balances the two:

Putting this into the equation for $m$ and integrating gives

With $a_{1,2} = m\,(1 \pm \omega)$ and the factor $e^{-B(t)x}$ restored, this is the first-order formula above. The same balance, carried to higher orders, gives every term.

All orders

The balance can be carried to every order, because the two equations for $m$ and $\omega$ are exact.

Expanding $\omega = \sum_n \varepsilon^n \omega_n$ gives $\omega_1 = \tilde g/2$ and

Every $\omega_n$ is a polynomial in $E = 1 - e^{-\kappa t}$, because $dE/dt = \kappa(1-E)$, so every term of $\log m$ integrates in closed form.

The outer series does not start at $\omega(0) = 0$. The fast transient that corrects it decays within about $1/(2\lambda)$, but it leaves a permanent shift in $\log m$:

The layer is solved to all orders in the span of $\tau^k e^{-2j\tau}$, with $\tau = \lambda t$. With both series the error after $n$ orders falls like $\varepsilon^{n+1}$.

The certificate verify_all_orders.py checks orders 1 to 6 against a 30-digit numerical solution. The survival demo shows the terms adding up.

Unequal switching rates

Let regime 1 switch to regime 2 at rate $\lambda_{12}$ and back at rate $\lambda_{21}$, with total rate $\gamma = \lambda_{12} + \lambda_{21}$ and stationary probabilities $p = \lambda_{21}/\gamma$ and $q = \lambda_{12}/\gamma$. Write $G_i = -\kappa\theta_i B + \tfrac12\sigma_i^2 B^2$ and $D = G_1 - G_2$. In the coordinates $s = p\,a_1 + q\,a_2$ and $r = (a_1 - a_2)/s$ the generator separates exactly, and

With $\varepsilon = 1/\gamma$ the outer series has $r_1 = D$ and

again polynomials in $E$. The initial layer decays like $e^{-\gamma t}$ and is solved to all orders in the same way. With $p = q = \tfrac12$ and $\gamma = 2\lambda$ this is the symmetric case above.

The certificate verify_general_orders.py checks orders 1 to 6 at $p = 0.3$ against a 30-digit numerical solution. The homogenize package implements both cases.

Thesis erratum: the dynamic-mean kernel

Chapter 3 of the linked 2001 thesis contains two distinct errors in its ``Dynamic Mean Vasicek'' approximation. They are historical source errata; the equations on this page already use the unambiguous convention that $\lambda$ is the jump rate out of each state.

First, the thesis uses two incompatible meanings of $\lambda$. The generator printed on p. 85 has off-diagonal rate $\lambda/2$, for which the centered state decays as $e^{-\lambda s}$ and the mean holding time is $2/\lambda$.

The numerical examples on pp. 125–128 instead call $\lambda$ the number of switches per year. That is the convention

with mean holding time $1/\lambda$. The factor of two changes prices, not just notation.

The corrected remaining-maturity formula

There is a separate reversal of time in Eq. (3.7). Let $\theta_{1,2}=\bar\theta\pm\delta\tilde\theta$ and $\sigma_{1,2}^2=\bar s\pm\delta\tilde s$, with fixed $\lambda,\kappa,T>0$. Define

The argument of $B$ is the maturity remaining after calendar time $s$. If $P_i^\delta(T,x)$ is the exact switched bond price, differentiation of the two-state affine ODE at $\delta=0$ gives, as $\delta\to0$ with all displayed parameters fixed,

Indeed, writing the time-to-maturity regime factors as a column $a^\delta(t)$ and $c=(1,-1)^\top$, their exact affine reduction is

At $\delta=0$, $a^0(t)=a_0(t)\mathbf1$. Therefore $z(t)=a_0(t)^{-1}\partial_\delta a^\delta(t)|_{\delta=0}$ satisfies $\dot z=Qz+\tilde g(t)c$, $z(0)=0$.

Since $Qc=-2\lambda c$, variation of constants and $s=T-t$ give $z(T)=cH_\lambda(T)$. The common short-rate factor $e^{-B(T)x}$ cancels from the price ratio.

The remainder is cubic because changing $\delta$ to $-\delta$ swaps the two states, so the log ratio is odd. Equivalently, conditional parameter averaging gives

and the deterministic time-inhomogeneous Vasicek formula is $P_i^{\rm DM}=P_{\rm av}\exp(\pm\delta H_\lambda)$. Thus the corrected dynamic-mean approximation has the exact first derivative in the regime contrast at every fixed switching rate and maturity. It is not claimed to be exact at finite contrast.

The required integrals are closed form. For $p=1,2$ and $c>0$,

At a removable resonance $c=j\kappa$, the quotient is interpreted as $Te^{-cT}$. Eq. (3.7) instead uses $e^{-\lambda s}B(s)$: it combines the p. 85 half-rate convention with elapsed rather than remaining maturity.

Fixed contrast under fast switching

The small-$\delta$ statement above holds at every fixed switching rate. A complementary result keeps the regime contrast finite and lets $\lambda\to\infty$.

Put $D_\lambda(T)=H_\lambda(T)$ with the unscaled contrasts $\tilde\theta,\tilde s$, and write the exact affine amplitudes as $a_1=m(1+\omega)$ and $a_2=m(1-\omega)$. Their equations are exactly

The linear convolution $D_\lambda$ solves the first equation with the nonlinear term $-\tilde g\omega^2$ removed. For fixed $T>0$, integration by parts gives

Because $\omega=O(\lambda^{-1})$, variation of constants in the difference equation gives $\omega(T)-D_\lambda(T)=O(\lambda^{-3})$. Expanding $\log((1+\omega)/(1-\omega))$ and integrating the equation for $m$ therefore yield

For the analytic Vasicek forcing, these remainder orders are sharp and their first omitted coefficients are explicit. Since $\tilde g(0)=0$, endpoint Laplace expansion and variation of constants give

Using $2\operatorname{arctanh}\omega=2\omega+2\omega^3/3+O(\omega^5)$ therefore sharpens the relative-price result to

The common level needs one extra endpoint term: exactly

Integration by parts in the first term contributes $-\tilde g(T)^2/(8\lambda^2)$, and the final logarithm contributes the same amount. Hence

Thus the corrected dynamic mean captures the exact regime log-price contrast through second order in $1/\lambda$, even at finite contrast. It does not give the full price to first order: it omits the common positive Green–Kubo convexity term in the second display. This separates a relative-price result from a common-level approximation.

The endpoint expansions and the common-level expansion above have fixed $T,\kappa,\theta_i,\sigma_i$ as $\lambda\to\infty$. In particular, no maturity-uniform claim is being made for the common Green–Kubo remainder.

A nonasymptotic bound, uniform in maturity

The relative-price statement has a stronger form that needs neither differentiability nor a fixed contrast. Let $q:[0,\infty)\to\mathbb R$ be bounded and continuous, put $M=\|q\|_\infty$ and $\eta=M/(2\lambda)$, and define

Define the invariant-radius function

For every finite $\eta\ge0$, with no small-contrast assumption, all maturities simultaneously satisfy

This is a finite-rate inequality, not asymptotic notation. Its proof is a direct scalar comparison. At $\omega=\rho(\eta)$ the largest possible outward derivative is

and the analogous inward inequality holds at $-\rho(\eta)$. Thus the Riccati flow leaves $[-\rho(\eta),\rho(\eta)]$ invariant. If $(\mathcal Vh)(T)=\int_0^T e^{-2\lambda(T-s)}h(s)\,ds$, variation of constants gives exactly

Finally, $\log(a_1/a_2)=2\operatorname{arctanh}\omega$, and $\operatorname{arctanh}x-x$ is increasing for $x\ge0$. The triangle inequality gives the displayed bound. As $\eta\downarrow0$, its right-hand side is $\tfrac83\eta^3+O(\eta^5)$, so it retains the cubic perturbative rate while remaining finite and valid for every contrast-to-switching ratio.

The bound also resolves a joint scaling that the fixed-contrast expansion does not cover. For a family $q_\lambda$ with $M_\lambda=\|q_\lambda\|_\infty$ and $M_\lambda/\lambda\to0$,

This corollary concerns only the regime log-price contrast. It neither removes nor bounds the separate common-level Green–Kubo correction.

Unequal transition rates: the quadratic term returns

The cubic rate above is a symmetry cancellation, not a generic two-state fact. Let the transition rates be $r_{12},r_{21}>0$, let $g_i(t)$ be the two affine amplitude forcings, put $q=g_1-g_2$, and write $\ell=\log(a_1/a_2)$. Directly from the two amplitude equations,

Set $s=r_{12}+r_{21}$ and let the unequal-rate dynamic-mean approximation solve

For a constant value $z$, define the positive equilibrium ratio

If $|q(t)|\le M$, the interval $[\log y_*(-M),\log y_*(M)]$ is invariant. With

every maturity simultaneously obeys the explicit finite-rate bound

Indeed, after subtracting the linear term the exact nonlinearity is

Variation of constants and $|\ell|\le L$ give the bound. Since $L=O(M/s)$, the first term is generally quadratic:

Thus fixed unequal rate proportions give a quadratic error, whereas $r_{12}=r_{21}$ removes the $\cosh\ell-1$ term exactly and restores the cubic order. Under the symmetric convention $r_{12}=r_{21}=\lambda$ and $q=2\tilde g$, this is the same cancellation used above. The statement is again only about the relative log price, not the common price level.

The quadratic order is sharp, and constant forcing gives its signed coefficient. Write

The equilibrium equation is $x-\sinh\ell_*+\alpha(\cosh\ell_*-1)=0$. Since its derivative in $\ell_*$ is $-1$ at the origin, the analytic implicit-function theorem gives, uniformly for $\alpha$ in every compact subset of $(-1,1)$,

The equilibrium is actually explicit to every order. With the square-root and logarithm branches selected continuously from $x=0$,

The generating function for the Legendre polynomials $P_n$ therefore gives the complete series

The only nearest inverse-function branch points are $x_\pm=\alpha\pm i\sqrt{1-\alpha^2}$, and both have modulus one. Hence the convergence radius is exactly one for every pair of positive transition rates: unequal rate proportions change every coefficient but not the natural scale $|q|<r_{12}+r_{21}$. Since $|P_n(\alpha)|\le1$ for real $|\alpha|\le1$, truncating after degree $N$ has the explicit uniform bound

The symmetry $\ell_{-\alpha}(-x)=-\ell_\alpha(x)$ gives $c_n(-\alpha)=(-1)^{n+1}c_n(\alpha)$. At equal rates this removes every even coefficient and reduces the formula to $\ell_0(x)=\operatorname{asinh}x$. This is an exact constant-forcing theorem; it does not replace the preceding maturity-uniform bound for time-dependent forcing.

Consequently, the signed dynamic-mean equilibrium error is $\ell_*-D_*=\tfrac12\alpha x^2+O(x^3)$ whenever the rates are unequal. If the rates are equal, $\alpha=0$ and the leading error is instead $-x^3/6+O(x^5)$. This expansion is a constant-forcing sharpness result; the preceding nonasymptotic bound, not this local series, is the claim for general bounded time-dependent forcing.

The all-order certificate independently reverts the implicit equilibrium equation through degree $20$ and compares those coefficients with the Legendre recurrence. Across seven signed asymmetries and $x=\pm0.35$, the maximum value error is $3.347\times10^{-12}$, the maximum equilibrium-equation residual is $2.349\times10^{-12}$, and the largest error is $0.171530$ of the displayed uniform remainder bound.

A finite-contrast certificate with $\kappa=1.3$, $T=4$, $\theta=(0.18,0.04)$ and $\sigma=(0.16,0.05)$ gives:

The last rate-doubling orders are $2.9986$ for the log-ratio residual, $0.9980$ for the omitted common shift, and $2.0037$ after subtracting its displayed leading term. At $\lambda=128$, the scaled remainders are $\lambda^3R_\lambda=4.84610\times10^{-5}$ versus the predicted $-\tilde g(T)^3/6=4.84281\times10^{-5}$, and $\lambda^2\{C_\lambda-(2\lambda)^{-1}\int\tilde g^2\}=-1.09963\times10^{-3}$ versus $-\tilde g(T)^2/4=-1.09674\times10^{-3}$. After subtracting these endpoint coefficients, the measured orders are $4.0637$ for the contrast residual (the $32$-to-$64$ doubling, before the ODE solve reaches its precision floor) and $2.9920$ for the common residual.

Checks against the printed examples

Both examples use $\bar\theta=0.25$, $\tilde\theta=0.10$ and essentially constant volatility. The table compares the magnitude of the thesis's dynamic-mean correction, the corrected value $\sinh(|H_\lambda|)$, and the exact two-state affine ODE correction $|a_2-a_1|/(2a_{\rm av})$.

The corrected kernel is close to the exact switched calculation in both cases; the printed kernel is on the wrong side by factors of about two. The p. 128 table also reverses the sign of the correction although regime 1 has the higher mean level and therefore the cheaper bond. The time-scale row on p. 127 is carried over from the preceding example: the stated parameters imply $1/\kappa=0.04$ and $T=5$, not $0.5$ and $15$.

The certificate checks the closed forms against quadrature to $4.7\times10^{-14}$, verifies the cubic contrast remainder with measured orders $2.999910$ and $2.999950$, verifies the fixed-contrast fast-switching orders and sharp endpoint coefficients above, checks the all-ratio maturity-uniform inequality on 260 constant and bounded time-dependent cases with $0.01\le\eta\le100$ (largest error/bound ratio $0.960384$), recovers order $1.198669$ when $M_\lambda\asymp\lambda^{0.6}$, and checks the unequal-rate theorem on 200 additional bounded forcings (largest error/bound ratio $0.780043$). Constant forcing gives orders $1.999923$ for unequal rates and $3.000000$ for equal rates; the scaled unequal-rate coefficient is $0.324982574$ versus the predicted $0.325000000$. An 80-digit check over seven signed rate asymmetries verifies the displayed fourth-order series with a fifth-order residual and recovers the quadratic coefficient for both rate orderings and the symmetric cubic coefficient.

The certificate also solves both exact switched ODEs. The affine bond structure is the one introduced by Vasicek (1977); the exact Markov-modulated mean-level model is treated by Elliott and Mamon (2002). Details and the source-page evidence are recorded in Issue #68.

A note setting out the derivation is in preparation. Comments and corrections are welcome through GitHub issues.