A fast mean-reverting factor
The continuous version of regime switching, with correlation, to all orders in $\sqrt\varepsilon$.
Replace the chain by a fast Ornstein–Uhlenbeck factor $dY = -\varepsilon^{-1}Y\,dt + \sqrt{2/\varepsilon}\,dZ$, with standard normal stationary law, driving the mean level and volatility of the rate: $dx = \kappa(\theta(Y) - x)\,dt + \sigma(Y)\,dW$, with $d\langle W, Z\rangle = \rho\,dt$. This is the setting of fast mean-reverting stochastic volatility for interest rates. Because $\kappa$ does not depend on $Y$, $u = e^{-B(t)x}a(t,y)$ and, with $\delta = \sqrt\varepsilon$,
$$a_t = \Big(\frac{1}{\delta^2}\mathcal L + \frac1\delta\,G_{-1}(t) + G_0(t)\Big)a,$$ $$G_{-1} = -\sqrt2\,\rho\,B(t)\,\sigma(y)\,\partial_y, \qquad G_0 = -\kappa\theta(y)B(t) + \tfrac12\sigma(y)^2B(t)^2,$$where $\mathcal L = -y\,\partial_y + \partial_{yy}$. In the Hermite basis $\mathcal L$ is $\operatorname{diag}(0, -1, -2, \dots)$, the constant function plays the role of the vector of ones, and the expectation under the stationary law plays the role of $\pi$. Multiplication by a polynomial in $y$ couples neighbouring modes, so every order of the expansion involves finitely many of them.
The engine carries over with $G_{-1}$ and $G_0$ in place of $\operatorname{diag} g$: writing $v = a/\pi a$, the outer terms are $w_n = \mathcal L^{\#}\big([w']_{n-2} - [F_{-1}]_{n-1} - [F_0]_{n-2}\big)$ with $F_p(w) = G_p v - v\,\pi G_p v$. The first correction, of order $\sqrt\varepsilon$, needs the correlation:
$$\log\frac{\pi a}{\pi a\big|_{\rho = 0}} = -\sqrt{2\varepsilon}\,\rho\int_0^t B(r)\;\mathbb{E}\big[\sigma(Y)\,\partial_y\phi_r(Y)\big]\,dr + O(\varepsilon), \qquad \mathcal L\phi_r = \bar g(r) - g(r,\cdot).$$For a check with an exact answer, take $\theta(y) = \theta_0 + \theta_1 y$ and $\sigma(y) = \sigma_0 + \sigma_1 y$. Then $a = e^{A + C_1y + C_2y^2}$ exactly, with $A$, $C_1$ and $C_2$ solving three ODEs. The table gives the error at $t = 1$, $y = 0.3$, after each order in $\sqrt\varepsilon$, for $\kappa = 1$, $\theta = 0.05 + 0.03y$, $\sigma = 0.25 + 0.2y$ and $\rho = -0.7$.
| 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}$. Without correlation only even orders appear and the expansion is in $\varepsilon$. A Monte Carlo simulation of the two-factor model agrees with the exact solution; see verify_fast_factor.py. The engines are fastswitch_op.py and fastswitch_gen.py.