Initial layers

Near the start the fast process has not yet mixed. A short correction on the fast clock repairs the expansion there.

The expansion built from solvability describes the answer once the fast process has mixed. At the very start it has not mixed, and the expansion does not satisfy the starting condition. A short correction on the fast time scale repairs this.

The outer expansion is a quasi-steady state

Each term of the expansion is found by solving $Q w = f$. The solution describes the fast process after it has settled into equilibrium with the slowly changing forcing $f(t)$.

Borrowing a word from fluid mechanics, this part is called the outer expansion. It is a function of the slow time $t$ only.

Nothing forces the outer expansion to take the right value at $t = 0$. On the regime-switching page, for example, the ratio $\omega = (a_1 - a_2)/(a_1 + a_2)$ must start at zero, because both regimes start with the same payoff. The outer series gives, at second order,

The outer expansion is not valid for the first moments, before the switch has had time to reach its equilibrium with the forcing.

The layer lives on the fast clock

The mismatch is removed on the fast time scale. Define the fast time $\tau = t/\varepsilon$, which counts switches rather than years, and write

Substituting into the equation for $w$ and keeping the fast dynamics gives, to leading order,

where $Q_0 = \varepsilon Q$ is the generator on the fast clock. Every direction orthogonal to the stationary distribution decays under $Q_0$, so the correction $\eta$ dies out exponentially in $\tau$.

For the symmetric two-state switch it decays like $e^{-2\tau} = e^{-2\lambda t}$, so it lives for a time of about $1/(2\lambda)$ in the original units. This short interval is the initial layer.

A uniform two-state theorem

The layer can be made quantitative. Put $\varepsilon=1/\lambda$, write $a_{1,2}=m(1\pm\omega)$, and let $b=(g_1+g_2)/2$ and $q=(g_1-g_2)/2$. Then

If $q$ is real- or complex-valued and continuously differentiable on $[0,T]$, then, for sufficiently small $\varepsilon$, uniformly for every $0\leq t\leq T$,

This is a uniform first-order composite approximation. It both satisfies the starting condition and retains the small permanent shift left by the transient. The constants in the remainder may depend on $T$ and bounds for $q$ and $q'$, but not on $t$ or sufficiently small $\varepsilon$.

Unequal switching rates

Symmetry is not needed. Let the fast-clock generator and its stationary law be

with $\alpha,\beta>0$, and let the physical generator be $Q_0/\varepsilon$. Put $m=\pi_1a_1+\pi_2a_2$, $r=(a_1-a_2)/m$, $\delta=g_1-g_2$, $\bar g=\pi_1g_1+\pi_2g_2$, and $v=\pi_1\pi_2$. The two-dimensional linear system is exactly equivalent to

For $\delta\in C^1[0,T]$, real or complex valued, the uniform first-order composite is

The components follow from $a_1=m(1+\pi_2r)$ and $a_2=m(1-\pi_1r)$. Substitution leaves an $O(\varepsilon)$ residual in the ratio equation; the stable kernel has mass $\varepsilon/\rho$, so the ratio and both components have uniform $O(\varepsilon^2)$ errors.

For $\alpha=\beta=1$, one has $\rho=2$, $\delta=2q$ and $r=2\omega$, recovering the preceding formula exactly. Thus the stationary weights determine the amplitude correction, while the nonzero eigenvalue $-\rho$ determines the layer width.

There is a new feature at second order. Set $s=\pi_2-\pi_1$, $D(t)=\int_0^t\delta(u)du$, $E(t)=e^{-\rho t/\varepsilon}$ and

If $\delta\in C^2[0,T]$, then the uniform second-order ratio is

The last term is a resonant layer: the unequal stationary weights couple the leading transient back into the linear term $s\delta r$. It vanishes for a symmetric chain. The accompanying amplitude is

Indeed, applying $d/dt+\rho/\varepsilon$ to the resonant term gives exactly $-\varepsilon s\delta(0)\delta(t)E(t)/\rho$, cancelling the uncancelled layer contribution from $s\delta r$. Substitution then leaves an $O(\varepsilon^2)$ ratio residual, and the stable kernel again contributes its mass $\varepsilon/\rho$, proving a uniform $O(\varepsilon^3)$ component error. When $\alpha=\beta=1$, $s=0$, $\rho=2$, $\delta=2q$ and $r=2\omega$; the formula reduces to the symmetric second-order composite below, up to terms already inside the stated remainder.

The matched outer formula

Simply deleting every exponential term does not give a second-order outer approximation when $\delta(0)\ne0$. Indeed, Laplace expansion of the boundary integral gives, at fixed positive maturity,

so the initial layer leaves the permanent log-amplitude shift $-\varepsilon^2v\delta(0)^2/\rho^2$. The correctly matched second-order outer formula is therefore

For every fixed $t_0>0$, this has $O(\varepsilon^3)$ error uniformly on $[t_0,T]$. More sharply, if $\delta(0)\ne0$, the same bound holds on the expanding interval

The factor two is sharp for generic components: the leading omitted layer is $-\varepsilon\delta(0)e^{-\rho t/\varepsilon}/\rho$, which becomes $O(\varepsilon^3)$ exactly at that scale. If $\delta(0)=0$, the permanent shift and leading layer vanish; the second-order layer then needs only $t\ge(\varepsilon/\rho)\log(1/\varepsilon)$.

These are outer-only statements. The composite formula above remains uniformly $O(\varepsilon^3)$ all the way to $t=0$.

The expansion continues uniformly without assuming $q(0)=0$. If $q$ is twice continuously differentiable, put $I_\varepsilon(t)=\int_0^tq(s)e^{-2s/\varepsilon}ds$. Then

Both orders of the exponential layer are needed to match $\omega(0)=0$ through second order. The integral terms retain the permanent amplitude shift accumulated while the layer decays.

The short-rate and survival examples have the simplifying cancellation $q(0)=0$. There $I_\varepsilon(t)=O(\varepsilon^2)$ uniformly, so the formula reduces, through order $\varepsilon^2$, to

Why the bound is uniform

Insert the first composite ratio into its Riccati equation. Its residual is

The error equation has the stable leading operator $d/dt+2/\varepsilon$. Convolution with its kernel has mass $\varepsilon/2$, turning this residual into an $O(\varepsilon^2)$ error uniformly, including at $t=0$. The general second composite has residual

because the $q(0)e^{-2t/\varepsilon}$ and $q'(0)e^{-2t/\varepsilon}$ derivative terms cancel their stiff contributions exactly. The same stable-kernel argument gives $O(\varepsilon^3)$ for arbitrary $q(0)$. This is the finite-dimensional boundary-layer mechanism in the primary rapid-switching expansions of Il'in, Khasminskii and Yin and Khasminskii and Yin.

The logarithmic crossover

A fixed number of switching times is enough to make the layer numerically small, but not enough to recover the next asymptotic order. When $q(0)\neq0$, the outer-only first-order error contains $\varepsilon q(0)e^{-2t/\varepsilon}/2$. To make this $O(\varepsilon^2)$ requires

When $q(0)=0$, exactly the same crossover makes the omitted second-order layer $\varepsilon^2q'(0)e^{-2t/\varepsilon}/4$ as small as the nominal $O(\varepsilon^3)$ remainder. Thus “away from a layer of width $O(\varepsilon)$” is correct for fixed positive maturity, but a uniform next-order claim needs either the exponential correction or the logarithmic factor.

Numerical certificate

Direct integration of the original two-component system over the whole interval, with the fast layer resolved, gives the following observed orders as $\varepsilon$ is halved:

For the unequal-rate case, naively deleting the second-order layer gives fixed-maturity error $\varepsilon^{2.023}$. Including the permanent amplitude shift improves this to $\varepsilon^{3.014}$ at fixed maturity and $\varepsilon^{3.011}$ uniformly beyond the sharp crossover. At the smallest tested $\varepsilon$, the naive error divided by its predicted leading $\varepsilon^2v\delta(0)^2e^{\int\bar g}/\rho^2$ coefficient is $1.006306$.

The code also evaluates the layer at $t=(\varepsilon/2)\log(1/\varepsilon)$ and obtains exactly a constant times $\varepsilon^3$ in the second case. Certificate: verify_uniform_layer.py.

Higher orders of the layer

At higher orders the layer feels the slow forcing too. The forcing is evaluated at $t = \varepsilon\tau$, which is expanded in a Taylor series around $t = 0$, and the products of decaying exponentials with these polynomials produce terms of the form

with $\mu$ an eigenvalue of $Q_0$ or a sum of eigenvalues. This family is closed under the operations the expansion needs: products, derivatives, integrals from $0$ to $\tau$, and solving $d\eta/d\tau = Q_0\eta + f$. So the layer can be computed order by order in closed form, just like the outer series.

The permanent part of the layer

Most of the layer disappears after a few multiples of $1/\lambda$. Only one part survives.

The quantity $s$, the stationary average of the answer, is an integral over time. The transient contributes a fixed amount to that integral before it dies out, which leaves a permanent shift in $\log s$.

On the regime-switching page this shift first appears at order $\varepsilon^4$, where $\omega_2(0) = -\tilde g'(0)/4$ is the second outer term of the ratio at the start:

It is small because the mismatch it repairs is already of order $\varepsilon^2$, and the forcing it meets is itself small near $t = 0$. Without it, however, the expansion stops improving at order $\varepsilon^3$: adding more outer terms does not reduce the error below the size of the missing shift. The convergence demo shows the orders that result when both are included.

The layer to second order for two states

For the two-state switch the whole layer can be written down. With $q = \tilde g$ and the ratio $\omega$ of the regime-switching page, $\omega' = q\,(1 - \omega^2) - 2\omega/\varepsilon$ from $\omega(0) = 0$, and the outer terms $\omega_1 = q/2$, $\omega_2 = -q'/4$ miss the start by $\tfrac12\varepsilon q(0) - \tfrac14\varepsilon^2 q'(0)$. The layer returns each to zero on the fast clock:

uniformly on $[0, T]$. The stiff parts of the two exponentials cancel in the equation, leaving a source of order $\varepsilon^2$ acting through a kernel of mass $\varepsilon/2$, which is where the order of the remainder comes from. Integrating $q\,\omega$ gives the mean, with $I(t) = \int_0^t q(s)\,e^{-2s/\varepsilon}\,ds$:

The two terms in $I$ are permanent shifts of the kind described above: $I$ settles at about $\tfrac12\varepsilon q(0)$, so the first is of order $\varepsilon^2$ and the second of order $\varepsilon^3$. When $q(0) = 0$, as for the affine survival examples, both vanish and the first shift is the $\varepsilon^4$ term. When $q(0) \ne 0$, as for counts, jumps and Black–Scholes, the first is the term $-\tfrac14\varepsilon^2\tilde g^{\,2}(1 - e^{-2\lambda T})$ that the exact count formula carries.

The certificate takes $q(t) = 0.4 + 0.3\sin 2t - 0.2\,t\,e^{-t}$ on $[0, 3]$ and measures the largest error as $\varepsilon$ halves from $0.2$ to $0.025$. The outer series is out by order $\varepsilon$ in $\omega$, which costs order $\varepsilon^2$ in $\log m$ because the miss lasts a time of order $\varepsilon$; the first composite is second order in both; the second composite is third order in both, with measured orders $3.00$ and $2.99$.

Certificate: verify_layer.py.

Layers at order zero

If the payoff itself depends on the final state of the fast process, the mismatch is of order one. This happens for an option on a bond whose price depends on the regime at expiry.

The unequal-rate coordinates make this case explicit. Let the initial vector be $c=(c_1,c_2)$ and suppose

Put $D(t)=\int_0^t\delta(u)du$, $E(t)=e^{-\rho t/\varepsilon}$ and $I_\varepsilon(t)=\int_0^t\delta(u)E(u)du$. For bounded $\delta\in C^1[0,T]$, the uniform first composite is

The first term $r_0E$ is the order-zero layer. The next two layer terms are not optional when the rates are unequal or the initial contrast is large: they are the first nonlinear Duhamel correction.

Applying $d/dt+\rho/\varepsilon$ to the displayed ratio leaves an $O(\varepsilon)$ residual; the stable kernel then gives the stated $O(\varepsilon^2)$ error. The components remain $a_1=m(1+\pi_2r)$ and $a_2=m(1-\pi_1r)$, so the formula reproduces $c_1,c_2$ exactly at $t=0$.

The certificate uses $(c_1,c_2)=(1.4,0.7)$ and $(\alpha,\beta)=(1.7,0.4)$. Immediate stationary projection has error order $0.000$; retaining only $r_0E$ gives order $1.009$; the full first composite gives order $2.005$ and matches the initial vector within $2.2\times10^{-16}$. Indicator payoffs are included whenever $m_0\ne0$.

For the error statement, take real bounded forcing and nonnegative terminal data, not both zero; more generally it is enough that the exact weighted mean stay bounded away from zero. The singular case $m_0=0$ requires direct vector coordinates rather than this ratio chart. The bond-options example is the motivating application.

An all-order Volterra composite

The first correction is not an isolated calculation. Define the stable convolution

and start from the order-zero layer $r^{[0]}(t)=r_0E(t)$. Successive nonlinear Duhamel iterates are

Reconstruct the mean and the two components by

for every fixed $k$ and $T$. Indeed, on any fixed bounded ratio ball the nonlinear source has Lipschitz constant $L=\|\delta\|_\infty(|s|+2vR)$, whereas $\|{\cal V}_\varepsilon\|_{\infty\to\infty}\le\varepsilon/\rho$. Thus the Volterra map contracts by at most $\varepsilon L/\rho$.

The starting error is $O(\varepsilon)$, so each iterate supplies one further power. Integrating the exact mean equation preserves the same order. This is a constructive two-state, bounded-horizon specialization of the full asymptotic-series framework of Khasminskii and Yin.

The same argument supplies an error certificate that does not use the exact solution. Put $M=\|\delta\|_\infty$, choose $R>|r_0|$, and suppose

Then the Volterra map preserves the radius-$R$ ball and is a contraction there. Banach's successive- approximation estimate gives, for $k\ge1$,

This is directly computable. If $G=\|\bar g\|_\infty$ and $w=\max(\pi_1,\pi_2)$, the component error is at most the right-hand side multiplied by

The estimate is the integral-equation contraction argument introduced in Banach's original 1922 paper, here with the fast-kernel norm making the contraction constant explicitly $O(\varepsilon)$.

The independent certificate solves the original two-component system and the scalar iterates separately. Iterates $k=0,1,2$ have measured uniform component-error orders $1.010$, $2.006$, and $2.997$; every iterate reproduces the terminal vector to numerical precision. With $R=2$, the maximum certified contraction factor is $0.038252$ and the invariant-ball margin is at least $1.071635$.

The computable bounds dominate every measured error and have orders $2.016$ and $3.011$ for $k=1,2$. The explicit first composite above is the closed asymptotic form of the $k=1$ iterate through its stated $O(\varepsilon^2)$ remainder.

Removing the singular ratio chart

The restriction $m_0\ne0$ belongs only to the ratio coordinate, not to the two-state problem. Put $b(t)=\int_0^t\bar g(u)du$, $m=\pi_1a_1+\pi_2a_2$, $d=a_1-a_2$, and remove the common growth by

These nonsingular slow and fast coordinates satisfy

For arbitrary real or complex terminal data, including $m_0=0$, eliminate $M$ and iterate the linear Volterra map

Reconstruct

If $M_\delta=\|\delta\|_\infty$ and

the radius-$R$ ball is invariant and the same Banach estimate applies with $q_\varepsilon^{\rm dir}$. The component multiplier is now simply

Thus the $k$th iterate has uniform component error $O(\varepsilon^{k+1})$ without any positivity or nonvanishing-mean assumption. The ratio recursion remains useful when available because it evolves one normalized contrast, but it is no longer needed for the theorem.

The singular-chart certificate takes $c=(1,-\pi_1/\pi_2)$, so $m_0=0$ exactly. Iterates $k=0,1,2$ have measured orders $1.008$, $2.021$, and $3.040$.

The computable bounds have orders $2.016$ and $3.028$, dominate every measured error, and use a maximum contraction factor $0.028468$. Every iterate still reproduces both terminal components to numerical precision.

The direct theorem for any finite chain

The nonsingular construction is not intrinsically two-state. Let $Q_0$ be any irreducible finite-state generator with stationary row vector $\pi^\top$, and set $\Pi={\bf1}\pi^\top$, $P=I-\Pi$. Assume that, in one fixed vector norm,

which always holds in finite dimension for some $M\ge1$ and $\gamma>0$. For

put $\bar g=\pi^\top g$, $H=G-\bar g I$, $b(t)=\int_0^t\bar g$, and decompose $u=e^{-b}a=m{\bf1}+z$, where $m=\pi^\top u$ and $z=Pu$. The exact equations are

Thus, for any continuous centered path $z$, define

Start with the exact order-zero layer $z^{[0]}=e^{tQ_0/\varepsilon}z_0$, iterate $z^{[k+1]}=\mathcal T_\varepsilon z^{[k]}$, and reconstruct $a^{[k]}=e^b(m_{z^{[k]}}{\bf1}+z^{[k]})$. No eigenvector basis and no division by $m$ is used.

Here is a finite-$\varepsilon$ stopping theorem. Put $h=\|H\|_\infty$, $\ell=\|\pi^\top H\|_\infty$, $c=\|{\bf1}\|$, and $B=\|\bar g\|_\infty$. If

then $\mathcal T_\varepsilon$ preserves the radius-$R$ ball and contracts there by $q_\varepsilon$. Therefore, for $k\ge1$,

The right-hand side is completely computable. Since $q_\varepsilon=O(\varepsilon)$ and the first Duhamel correction is $O(\varepsilon)$, every fixed iterate satisfies $\|a-a^{[k]}\|_\infty=O(\varepsilon^{k+1})$ on a bounded horizon.

This includes nonreversible and defective generators, arbitrary signed or complex terminal vectors, and paths on which the stationary component vanishes. It is a bounded-horizon statement; neither $T$ nor the iteration count is allowed to grow silently with $1/\varepsilon$.

The certificate uses an irreducible nonreversible three-state generator with a Jordan block, $Q_0=\Pi-I+0.1uv^\top$, where $u=(1,-1,0)^\top$ and $v=(1,1,-2)^\top$. Since the nilpotent term $N=0.1uv^\top$ has $N^2=0$, the centered semigroup is exactly $e^{-t}(P+tN)$; the proof uses the explicit bound $M=1+2\|N\|/e=1.254874$ and $\gamma=1/2$.

For signed initial data with zero stationary mean, iterates $k=0,1,2$ have measured component-error orders $0.995$, $1.983$, and $2.956$. Closed-form forcing bounds (rather than sampled suprema) give computable error-bound orders $2.203$ and $3.185$; those bounds dominate every measured error and have maximum contraction factor $0.393246$.

An explicit first composite for arbitrary terminal data

The first direct iterate can also be reduced to a closed matched formula. Keep $R=Q_0^\#$, $E_\tau=e^{\tau Q_0}$, $m_0=\pi^\top c$, $z_0=Pc$, $h(t)=H(t){\bf1}$, and

The order-zero layer interacts once with the forcing through

This is the upper-right block of $\exp\{\tau\left(\begin{smallmatrix}Q_0&PH(0)\\0&Q_0\end{smallmatrix}\right)\}z_0$. Define

For $g\in C^1[0,T]$ and arbitrary real or complex $c$,

uniformly on a fixed bounded horizon. Every displayed correction cancels appropriately at $t=0$, so the terminal vector is reproduced exactly even when $m_0=0$. The three first-order effects are distinct: the slow Poisson corrector, the layer that matches it at the terminal point, and the block-exponential interaction of the order-zero layer with $H(0)$.

The proof follows by integration by parts on the centered semigroup: $\int_0^t\pi^\top H(s)E_{s/\varepsilon}z_0\,ds =\varepsilon\{\pi^\top H(t)E_{t/\varepsilon}Rz_0-\pi^\top H(0)Rz_0\} +O(\varepsilon^2)$. Freezing $H(\varepsilon\tau)$ at $H(0)$ in the fast convolution also costs $O(\varepsilon^2)$; the remaining slow convolution gives the two matched Poisson terms. This is the arbitrary-terminal-data first-order specialization of the matched finite-chain expansions of Khasminskii, Yin and Zhang, expressed with Meyer's group inverse.

On the defective nonreversible three-state certificate with signed zero-mean terminal data, the uniform error has observed order $2.004210$ and the terminal discrepancy is below $10^{-12}$. The calculation uses no eigendecomposition.

An explicit first composite for a finite chain

For the regime-independent terminal vector $a(0)=m_0{\bf1}$, the first finite-chain composite can be written without performing a Picard iteration. Let $Q_0^\#$ be the group inverse, normalized by

and put $h(t)=H(t){\bf1}=g(t)-\bar g(t){\bf1}$ and

If $g\in C^1[0,T]$ is real or complex valued, then uniformly for $0\le t\le T$,

The two corrector terms cancel at $t=0$, so the formula returns the terminal vector exactly. The first is the quasi-steady Poisson corrector; the second is its initial layer. The scalar term $K$ is the accumulated feedback of that corrector into the stationary component.

For completeness, substitution gives a centered residual of order $\varepsilon$: it is the sum of $-\varepsilon m_0Q_0^\#h'$, the $\varepsilon K h$ term, and $-PH$ applied to the displayed centered corrector. Convolution against the centered semigroup contributes another factor $\varepsilon$. The scalar residual is $-\varepsilon m_0\pi^\top H(t)e^{tQ_0/\varepsilon}Q_0^\#h(0)$, whose integral is $O(\varepsilon^2)$.

A coupled Gronwall estimate therefore gives the stated uniform component error. This is the finite-chain Poisson-equation form of the matched expansion; the group-inverse normalization follows Meyer's primary finite-chain treatment.

Deleting the exponential term produces the first outer formula. It is only uniformly $O(\varepsilon)$ because it misses the terminal condition, but it is $O(\varepsilon^2)$ on every fixed interval $[t_0,T]$, $t_0>0$. Under the displayed semigroup bound, the same order holds on the expanding interval

The defective three-state certificate uses the same nonnormal generator as above. The unmatched outer error has measured order $1.000000$, while the explicit composite and the outer formula after the logarithmic crossover have orders $2.002876$ and $2.002876$. The composite matches all three terminal components to $2.2\times10^{-16}$.

The explicit second composite

The same Poisson construction can be carried through one more order without diagonalizing $Q_0$. Write $R=Q_0^\#$, $E_\tau=e^{\tau Q_0}$, $H_0=H(0)$, and define

The second boundary layer, with $\tau=t/\varepsilon$, is

The convolution is explicit and numerically stable: it is the upper-right block of

If $g\in C^2[0,T]$ is real or complex valued, then the second composite is

uniformly for $0\le t\le T$. This is a fixed bounded-horizon statement. The constant $K_2(0)=-\pi^\top H_0R^2h(0)$ is not an arbitrary normalization: it is the permanent scalar trace of the first boundary layer.

Indeed $B_2(0)=-z_2(0)-K_2(0){\bf1}$, so both first- and second-order corrections cancel at the terminal point. Omitting only this trace generically leaves $O(\varepsilon^2)$ error even at a fixed positive maturity.

For the error estimate, the displayed terms make the scalar outer residual vanish through second order and leave an $O(\varepsilon^2)$ centered residual; convolution against the centered semigroup gains the third power of $\varepsilon$. The layer equations are $(\partial_\tau-Q_0)B_1=0$ and $(\partial_\tau-Q_0)B_2=H_0B_1$. Replacing $H_0$ by $H(t)$ leaves only a localized residual because $H(t)-H_0=O(\varepsilon\tau)$.

A coupled variation-of-constants estimate then gives the uniform remainder. This is the concrete second-order specialization of the matched expansions for singular Markov generators developed by Khasminskii, Yin and Zhang, with the Poisson solves represented by Meyer's group inverse.

Deleting $B_1$ and $B_2$ gives a second outer formula whose full-interval error is still only $O(\varepsilon)$ because it misses the terminal condition. Under the semigroup bound above, it regains $O(\varepsilon^3)$ on the expanding interval

The defective three-state certificate measures orders $2.002876$ and $2.960961$ for the first and second composites. The second outer formula has order $1.001913$ on the full interval and $2.833557$ after the doubled crossover; omitting the permanent trace has fixed-maturity order $2.001789$. The second composite matches all three terminal components to $2.2\times10^{-16}$.