A fast mean-reverting factor
The continuous version of regime switching, with correlation, derived term by term.
In most of the examples a Markov chain switches the parameters of a rate. Here a continuous factor does the same job: it moves the mean level and the volatility smoothly and quickly. The correction terms follow one at a time from the idea of averaging and the solvability condition.
The model
The short rate $x_t$ reverts to a level, and diffuses with a volatility, that both depend on a hidden factor $Y_t$:
$$dx_t = \kappa\,\big(\theta(Y_t) - x_t\big)\,dt + \sigma(Y_t)\,dW_t .$$The factor is a fast Ornstein–Uhlenbeck process,
$$dY_t = -\frac{1}{\varepsilon}\,Y_t\,dt + \sqrt{\frac{2}{\varepsilon}}\,dZ_t, \qquad d\langle W, Z\rangle = \rho\,dt .$$Its stationary law is standard normal whatever the value of $\varepsilon$. It forgets its starting point over a time of order $\varepsilon$, while the rate reverts over $1/\kappa$. So $\varepsilon \ll 1/\kappa$ separates the two time scales.
The correlation $\rho$ ties the shocks to the rate to the shocks to the factor. This is the setting of fast mean-reverting stochastic volatility for interest rates, with the mean level allowed to move as well.
For a check with a closed-form answer, both functions are taken to be linear:
$$\theta(y) = \theta_0 + \theta_1\,y, \qquad \sigma(y) = \sigma_0 + \sigma_1\,y .$$Parameters: $\kappa = 1$, $\theta_0 = 0.05$, $\theta_1 = 0.03$, $\sigma_0 = 0.25$, $\sigma_1 = 0.2$, $\rho = -0.7$.
The linear $\sigma$ is a signed diffusion loading rather than a volatility. It is negative when $y < -\sigma_0/\sigma_1 = -1.25$, which the stationary law visits about 11% of the time, and there the instantaneous correlation between the shocks to the rate and to the factor is $-\rho$ rather than $\rho$.
The non-negative loading $|\sigma(y)|$ has the same $\sigma^2$, so the same averaged model, but a different order-$\sqrt\varepsilon$ term, because that term is linear in the loading: at $t = 1$ its coefficient is $0.000122$ with the signed loading and $0.000035$ with $|\sigma|$. The closed forms and the checks below use the signed loading.
The quantity
The bond price for maturity $t$, starting from rate $x$ and factor value $y$:
$$u(t, x, y) = \mathbb{E}\big[e^{-\int_0^t x_r\,dr} \mid x_0 = x,\ Y_0 = y\big].$$Averaging
When the factor is infinitely fast, the rate sees only its stationary law. The limit is Vasicek's model with the averaged mean level and variance,
$$\bar\theta = \mathbb{E}[\theta(Y)] = \theta_0, \qquad \bar\sigma^2 = \mathbb{E}[\sigma(Y)^2] = \sigma_0^2 + \sigma_1^2,$$with the expectations taken under the standard normal law. The corrections below are what the average leaves out.
Step 1: the rate factors out
By the Feynman–Kac formula, $u$ solves a partial differential equation in $t$, $x$ and $y$. Try the Vasicek form in $x$:
$$u(t,x,y) = e^{-B(t)\,x}\,a(t,y), \qquad B(t) = \frac{1 - e^{-\kappa t}}{\kappa}.$$Every term that involves $x$ cancels, because $B$ solves $B' = 1 - \kappa B$ and does not depend on $y$. What remains is an equation for $a$ alone. With $\delta = \sqrt\varepsilon$ it reads
$$\partial_t a = \frac{1}{\delta^2}\,\mathcal L a + \frac1\delta\,G_{-1}(t)\,a + g(t,y)\,a, \qquad a(0,y) = 1 .$$Each operator has a source.
- $\mathcal L = -y\,\partial_y + \partial_{yy}$ is the generator of the factor on its own clock. It is multiplied by $1/\delta^2 = 1/\varepsilon$ because the factor is fast.
- The correlation term comes from the cross derivative $\rho\,\sigma(y)\sqrt{2/\varepsilon}\;\partial_x\partial_y u$. Acting on $e^{-Bx}$, the $\partial_x$ produces $-B$, which leaves
- The multiplication term collects the drift and the variance of the rate:
For linear $\theta$ and $\sigma$ this is a quadratic in $y$. Writing $y^2 = \mathrm{He}_2(y) + 1$, with the Hermite polynomials $\mathrm{He}_1 = y$ and $\mathrm{He}_2 = y^2 - 1$,
$$g = (g_0 + g_2) + g_1\,\mathrm{He}_1 + g_2\,\mathrm{He}_2,$$where the three coefficient functions are
$$g_0 = -\kappa\theta_0 B + \tfrac12\sigma_0^2B^2, \qquad g_1 = -\kappa\theta_1 B + \sigma_0\sigma_1B^2, \qquad g_2 = \tfrac12\sigma_1^2B^2 .$$The stationary average of $g$ is $\bar g = g_0 + g_2$.
Step 2: split off the average
Let $s(t) = \mathbb{E}[a(t,Y)]$ be the stationary average of $a$, and $v = a/s$ its shape, so that $\mathbb{E}[v] = 1$.
Take the stationary average of the equation for $a$. The term $\mathcal L a$ averages to zero, because the standard normal law is stationary. Dividing by $s$ gives the slow equation:
$$\frac{s'}{s} = \frac1\delta\,\mathbb{E}\big[G_{-1}v\big] + \mathbb{E}\big[g\,v\big].$$Subtracting $v$ times this from the equation for $a/s$ gives the equation for the shape:
$$\partial_t v = \frac{1}{\delta^2}\,\mathcal L v + \frac1\delta\Big(G_{-1}v - v\,\mathbb{E}[G_{-1}v]\Big) + \Big(g\,v - v\,\mathbb{E}[g\,v]\Big).$$Both brackets average to zero by construction. That is what will make every step solvable.
Step 3: expand the shape
Write the shape as a power series in $\delta$, with each correction averaging to zero:
$$v = 1 + \delta\,w_1 + \delta^2\,w_2 + \delta^3\,w_3 + \cdots, \qquad \mathbb{E}[w_n] = 0 .$$Substitute and collect powers of $\delta$. Each power gives an equation $\mathcal L w_n = f_n$. By the Fredholm alternative it can be solved when $\mathbb{E}[f_n] = 0$, and the solution divides each Hermite coefficient of $f_n$ by minus its index:
$$\mathcal L\,\mathrm{He}_k = -k\,\mathrm{He}_k \quad\Longrightarrow\quad \mathcal L^{-1}\,\mathrm{He}_k = -\frac{\mathrm{He}_k}{k}.$$Order $\delta^{-1}$
The only terms are $\mathcal L w_1$ and the correlation term acting on the constant $1$:
$$\mathcal L\,w_1 = -\big(G_{-1}1 - \mathbb{E}[G_{-1}1]\big) = 0,$$since $\partial_y 1 = 0$. So $w_1 = 0$. The factor has no effect at order $\delta$ on the shape.
Order $\delta^{0}$
Now $g$ enters:
$$\mathcal L\,w_2 = -\big(g - \bar g\big) = -g_1\,\mathrm{He}_1 - g_2\,\mathrm{He}_2 .$$The right side averages to zero, so it can be solved by dividing by $-1$ and $-2$:
$$w_2 = g_1\,\mathrm{He}_1 + \tfrac12\,g_2\,\mathrm{He}_2 .$$This is the corrector. It says how the bond price depends on the current factor value: a high factor raises the mean level and the volatility, and the corrector records by how much.
Order $\delta^{1}$
The correlation term now acts on $w_2$:
$$\mathcal L\,w_3 = -\big(G_{-1}w_2 - \mathbb{E}[G_{-1}w_2]\big).$$Write $c(t) = \sqrt2\,\rho\,B(t)$, so that $G_{-1} = -c\,\sigma(y)\,\partial_y$. Since $\partial_y w_2 = g_1 + g_2\,y$,
$$G_{-1}w_2 = -c\,(\sigma_0 + \sigma_1 y)(g_1 + g_2 y) = -c\Big[(\sigma_0g_1 + \sigma_1g_2) + (\sigma_0g_2 + \sigma_1g_1)\,\mathrm{He}_1 + \sigma_1g_2\,\mathrm{He}_2\Big].$$Removing the average and dividing by $-1$ and $-2$ gives
$$w_3 = -c\,(\sigma_0g_2 + \sigma_1g_1)\,\mathrm{He}_1 - \tfrac12\,c\,\sigma_1g_2\,\mathrm{He}_2 .$$Step 4: the slow equation, order by order
Substitute the expansion of $v$ into $s'/s = \delta^{-1}\mathbb{E}[G_{-1}v] + \mathbb{E}[g\,v]$ and collect powers of $\delta$. The Hermite polynomials are orthogonal with $\mathbb{E}[\mathrm{He}_k^2] = k!$, which makes each average a short sum:
$$\mathbb{E}[g] = \bar g, \qquad \mathbb{E}[G_{-1}w_2] = -c\,(\sigma_0g_1 + \sigma_1g_2),$$ $$\mathbb{E}[g\,w_2] = g_1^2 + g_2^2, \qquad \mathbb{E}[G_{-1}w_3] = c^2\big(\sigma_0^2g_2 + \sigma_0\sigma_1g_1 + \sigma_1^2g_2\big).$$Integrating in time, the stationary average of the bond factor is
$$\begin{aligned} \log s(t) = &\int_0^t \bar g \;-\; \delta\int_0^t c\,(\sigma_0g_1 + \sigma_1g_2) \\ &+ \delta^2\int_0^t \Big[g_1^2 + g_2^2 + c^2\big(\sigma_0^2g_2 + \sigma_0\sigma_1g_1 + \sigma_1^2g_2\big)\Big] + O(\delta^3). \end{aligned}$$Every coefficient is a polynomial in $B$:
$$\bar g = -\kappa\theta_0 B + \tfrac12(\sigma_0^2 + \sigma_1^2)B^2, \qquad g_1 = -\kappa\theta_1 B + \sigma_0\sigma_1 B^2, \qquad g_2 = \tfrac12\sigma_1^2 B^2, \qquad c = \sqrt2\,\rho\,B .$$So every integral is a combination of $I_k = \int_0^t B^k$, and the exponent is
Each term has a reading.
- The first is the averaged Vasicek model.
- The second, of order $\sqrt\varepsilon$, exists only with correlation, and it is linear in $\rho\,\sigma(y)$. Where $\rho\,\sigma(y) < 0$, as at $\rho = -0.7$ and $y$ near zero, a shock that raises the rate lowers the factor, and with it the mean level and the loading; below $y = -1.25$ the loading changes sign and so does the co-movement. This is the correlation correction of fast mean-reverting stochastic volatility.
- The third, of order $\varepsilon$, contains the fluctuation term $g_1^2 + g_2^2$. As on the idea page, it is half the variance of the fluctuating integral of $g$, and it is present even without correlation.
Step 5: the bond price
Putting the pieces together, and writing the Hermite polynomials out,
with $\log s(t)$ from the card above. The dependence on the current factor value $y$ first appears at order $\varepsilon$.
At $t = 1$, $y = 0.3$, $\varepsilon = 0.0025$ and the parameters above:
| quantity | value |
|---|---|
| $B$ | 0.632121 |
| $I_1$, $I_2$, $I_3$, $I_4$ | 0.367879, 0.168091, 0.0838978, 0.0439824 |
| coefficients of $1$, $\delta$, $\delta^2$ in $\log s$ | −0.00977930, 0.000122386, 0.0000999259 |
| the $\delta$ and $\delta^2$ terms at $\delta = 0.05$ | 0.00000611930, 0.000000249815 |
| coefficients of $\varepsilon$, $\varepsilon^{3/2}$ in the bracket | −0.00333158, −0.0000418965 |
| the $\varepsilon$ and $\varepsilon^{3/2}$ terms at $\varepsilon = 0.0025$ | −0.00000832895, −0.00000000523706 |
| bond factor $u\,e^{Bx}$ | value | error |
|---|---|---|
| averaged model | 0.99026837 | 2.1e-6 |
| through $\delta$ | 0.99027443 | 8.1e-6 |
| through $\delta^2$ | 0.99026643 | 1.3e-7 |
| numerical (Riccati) | 0.99026629 |
At these parameters the $\delta$ term moves the averaged value the wrong way before the $\delta^2$ term corrects it, which is the expansion in powers of $\sqrt\varepsilon$ at work.
All orders
Every further order repeats steps 3 and 4. At order $\delta^{n-2}$ the equation is
$$\mathcal L\,w_n = w_{n-2}' - \big[F_{-1}\big]_{n-1} - \big[F_0\big]_{n-2},$$where $F_p$ collects the terms of order $\delta^k$ in the brackets of the shape equation:
$$F_{-1}(v) = G_{-1}v - v\,\mathbb{E}[G_{-1}v], \qquad F_0(v) = g\,v - v\,\mathbb{E}[g\,v].$$The right side always averages to zero, so every step is solvable. Because $g$ is a polynomial in $y$, each $w_n$ involves only finitely many Hermite polynomials, and each is a combination of powers of $e^{-\kappa t}$. The computation is exact at every order; only the truncation of the series introduces error.
The shape starts at $v(0,y) = 1$, and the outer terms respect this for a while: $g$ vanishes at $t = 0$ because $B(0) = 0$, so $w_2(0) = w_3(0) = 0$. The first mismatch is at $w_4$, where the time derivative $w_2'$ enters.
From there an initial layer in the fast time $t/\varepsilon$ restores the starting condition. It is solved in the eigenbasis of $\mathcal L$, which here is the Hermite basis itself.
Results
For linear $\theta$ and $\sigma$ the solution also has the form $a = e^{A + C_1y + C_2y^2}$, with $A$, $C_1$ and $C_2$ solving three ordinary differential equations. Their numerical solution is the reference.
The terms derived above agree with the code to $10^{-17}$. The table gives the error in the bond factor at $t = 1$, $y = 0.3$, after each order in $\sqrt\varepsilon$.
| epsilon | order 0 | order 1 | order 2 | order 3 | order 4 | order 5 | order 6 |
|---|---|---|---|---|---|---|---|
| 0.01 | 2.1e-5 | 3.3e-5 | 1.1e-6 | 1.8e-7 | 1.2e-8 | 1.4e-9 | 1.2e-9 |
| 0.0025 | 2.1e-6 | 8.1e-6 | 1.3e-7 | 1.2e-8 | 3.5e-10 | 1.9e-11 | 1.5e-11 |
| 0.000625 | 1.0e-6 | 2.0e-6 | 1.6e-8 | 7.4e-10 | 1.0e-11 | 2.8e-13 | 2.2e-13 |
Each quartering of $\varepsilon$ halves $\sqrt\varepsilon$ and divides the order-$n$ error by about $2^{n+1}$, until round-off near $10^{-13}$. A Monte Carlo simulation of the two-factor model agrees with the reference solution; see verify_fast_factor.py. The engines are fastswitch_op.py and fastswitch_gen.py.