homogenization

The fast-switching engine

One recursion for every model whose state factors out of an exponential-affine formula.

Each example has a slow process $x_t$, the rate, intensity, variance or count of interest, and a fast process $y_t$ that it does not observe: a Markov chain with a finite number of regimes, or a fast mean-reverting diffusion. The quantity wanted is a Feynman–Kac expectation,

$$\mathbb{E}\big[e^{-\int_0^t c\cdot x_r\,dr}\,\phi(x_t, y_t)\big].$$

When the model is exponential-affine in $x$ and the coefficient of $x$ does not depend on the regime, the state factors out and what remains is a linear system in time only:

$$u_i(t,x) = e^{-B(t)\cdot x}\,a_i(t), \qquad a' = \big(Q + \operatorname{diag} g(t)\big)\,a, \qquad a(0) = \text{the payoff in each regime}.$$

The model enters only through the functions $g_i(t)$, one per regime, and the terminal vector.

ModelSwitching parametersTime function
Vasicek, any number of regimesmean level, volatility$-\kappa\theta_i B + \tfrac12\sigma_i^2 B^2$
Gaussian factors (certificate)means, volatilities, correlations$-\sum_j \kappa_j\theta_{ji} B_j + \tfrac12\sum_{j,l}\rho_{jl,i}\sigma_{ji}\sigma_{li}B_jB_l$
CIRmean level$-\kappa\theta_i B_{\mathrm{CIR}}$
Two-name credit, CIR namesmean levels$-\kappa\,(c_1\theta_{1,i} + c_2\theta_{2,i})\,B_{\mathrm{CIR}}$
Vasicek with jumpsmean, volatility, jump intensity$-\kappa\theta_i B + \tfrac12\sigma_i^2B^2 + \ell_i\big(\tfrac{1}{1+mB} - 1\big)$
Poisson countsarrival rate$(z-1)\,\ell_i$
Hestonlong-run variance$\kappa\theta_i D(t)$, complex
Black–Scholesvolatility$iu(r - \tfrac12\sigma_i^2) - \tfrac12u^2\sigma_i^2$
Bond optionsmean level, volatility$-\kappa\theta_i \tilde B + \tfrac12\sigma_i^2\tilde B^2$, indicator terminal vector
Mertonvolatility, jump intensity$iu(r - \tfrac12\sigma_i^2) - \tfrac12u^2\sigma_i^2 + \ell_i(\phi_J(u) - 1 - iu\bar k)$
Variance gamma$\sigma$, $\nu$, $\theta$$iu(r + \omega_i) + \psi_i(u)$, the regime's Lévy exponent
Bateslong-run variance, jump intensity$\kappa\theta_i D(t) + \ell_i(\phi_J(u) - 1 - iu\bar k)$
Hull–Whitevolatility$\tfrac12\sigma_i^2 B^2$
G2++volatilities, correlation$\tfrac12\sigma_i^2B_a^2 + \tfrac12\eta_i^2B_b^2 + \rho_i\sigma_i\eta_iB_aB_b$
Equity with ratesboth parts' parametersthe rates forcing at weight $1 - iw$ plus the equity forcing at $w$
CEVvolatility scale$-z\,\sigma_i^2\,h(\tau)$, for the transform of the total variance
Fast mean-reverting factormean level and volatility, continuouslyan operator in the Hermite basis, with correlation

Deriving the recursion

The derivation uses two ideas from the Background pages: the solvability condition and the initial layer. It is written here for any number of regimes.

Step 1: separate the average from the shape

Let $\pi$ be the stationary distribution of the chain, so that $\pi Q = 0$. Split $a$ into its stationary average and its shape:

$$s = \pi\cdot a, \qquad w = \frac{a}{s} - \mathbf 1, \qquad \pi\cdot w = 0 .$$

Multiply the system $a' = (Q + \operatorname{diag} g)a$ on the left by $\pi$. The switching term drops out because $\pi Q = 0$, which leaves an equation for the average:

$$\frac{s'}{s} = \bar g + \pi\cdot(g\,w), \qquad \bar g = \pi\cdot g .$$

Subtracting this from the equation for $a/s$ gives the equation for the shape. Writing $Q = Q_0/\varepsilon$, with the scale $\varepsilon = n/(-\operatorname{tr} Q)$ used by the code, the reciprocal of the mean exit rate over the $n$ states,

$$\varepsilon\,w' = Q_0\,w + \varepsilon\,F(w), \qquad F(w) = g\,(\mathbf 1 + w) - (\mathbf 1 + w)\,\big(\bar g + \pi\cdot(g\,w)\big).$$

Any positive scale would serve, provided $Q_0$ is held fixed as the chain speeds up. This one equals the mean holding time only when every state leaves at the same rate; for a symmetric two-state chain it is $1/\lambda$.

The right side of the shape equation always averages to zero under $\pi$. That is what makes the next step work.

Step 2: expand the shape

Write $w = \varepsilon w_1 + \varepsilon^2 w_2 + \cdots$ and match powers of $\varepsilon$. At order $\varepsilon^n$ the shape equation reads

$$Q_0\,w_n = w_{n-1}' - F_{n-1},$$

where $F_{n-1}$ collects the terms of $F$ of order $\varepsilon^{n-1}$. By the Fredholm alternative this has a solution because the right side averages to zero, and the solution with $\pi\cdot w_n = 0$ is

$$w_n = Q_0^{\#}\big(w_{n-1}' - F_{n-1}\big),$$

with $Q_0^{\#}$ the group inverse. Each order therefore costs one application of a fixed matrix to terms already known.

Step 3: the first order

At the first order $w_0 = 0$ and $F_0 = g - \bar g\,\mathbf 1$, so, with $Q^{\#} = \varepsilon\,Q_0^{\#}$ the group inverse of $Q$ itself,

$$\varepsilon\,w_1 = -Q^{\#}\big(g - \bar g\big).$$

Putting this into the equation for $s$ gives

$$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).$$

The correction in the exponent has a direct meaning. Because $Q^{\#}f = -\int_0^\infty e^{Q\tau}f\,d\tau$ when $\pi\cdot f = 0$,

$$-\pi\cdot\big(g\,Q^{\#}(g - \bar g)\big) = \int_0^\infty \mathbb{E}_\pi\big[\big(g(y_0) - \bar g\big)\big(g(y_\tau) - \bar g\big)\big]\,d\tau .$$

This is a Green–Kubo integral of the autocorrelation of $g$ along the chain. For real $g$ it is the integrated autocovariance, half the long-run variance rate of $\int g$, as on the idea page, now for any chain.

For complex $g$, as in the Fourier examples, the product carries no conjugate. The integral is then a second cumulant rather than a variance, and it can be complex. The factor in brackets is the memory of the starting regime.

Step 4: closed form

When every $g_i$ is a sum of exponentials $e^{-\alpha t}$, products and derivatives of such sums are again such sums. So every $w_n$ is a sum of exponentials, and every term of $\log s$ integrates in closed form. Otherwise $g_i$ is held as a Chebyshev series, which is also closed under these operations.

Step 5: the initial layer

The outer series does not match the terminal vector at $t = 0$. The mismatch is removed on the fast clock $\tau = t/\varepsilon$:

$$\eta = w - w_{\mathrm{outer}}, \qquad \frac{d\eta}{d\tau} = Q_0\,\eta + \varepsilon\,\big[F(w_{\mathrm{outer}}+\eta) - F(w_{\mathrm{outer}})\big], \qquad \eta(0) = w(0) - w_{\mathrm{outer}}(0).$$

It decays at the rates of the chain and is solved order by order in a Schur basis of $Q_0$ on the mean-zero subspace, by back substitution through the triangular form, in the span of $\tau^k e^{\mu\tau}$. A generator with a Jordan block has no eigenvector basis, but it has a triangular one, and the powers of $\tau$ come out of the back substitution. The result is $a_i = s\,(1 + w_i + \eta_i)$, with error $O(\varepsilon^{N+1})$ after $N$ orders.

Exact constant forcing with unequal rates

For a two-state generator $m\left(\begin{smallmatrix}-a&a\\b&-b\end{smallmatrix}\right)$ and constant complex forcing $(g_1,g_2)$, the outer/layer split can be written without assuming $a=b$. Put $\kappa=a+b$, $\vartheta=(a-b)/\kappa$, $\varepsilon=(m\kappa)^{-1}$, $\delta=g_1-g_2$, $z=\varepsilon\delta$, and $\bar g_\pi=(bg_1+ag_2)/\kappa$. Then

$$\phi_\pm=e^{\bar g_\pi T}\left[ \frac12\left(1+\frac{1\pm z}{d}\right)e^{(d-1+\vartheta z)T/(2\varepsilon)}+ \frac12\left(1-\frac{1\pm z}{d}\right)e^{-(d+1-\vartheta z)T/(2\varepsilon)} \right],\qquad d=\sqrt{1-2\vartheta z+z^2}.$$

The two branch points $\vartheta\pm i\sqrt{1-\vartheta^2}$ have modulus one. Consequently the slow mode has the same radius-one expansion as in the symmetric chain, now in forcing contrast divided by total switching rate. More explicitly, put

$$c_1=\frac{T\delta(1-\vartheta^2)}4,\quad c_2=\frac{T\delta\vartheta(1-\vartheta^2)}4,\quad a_{1,\pm}=\frac{\vartheta\pm1}{2},\quad a_{2,\pm}=\frac{3\vartheta^2-1}{4}\pm\frac\vartheta2.$$

The slow multiplier after $e^{\bar g_\pi T}$ has quadratic Taylor polynomial

$$P_{2,\pm}(z)=1+(a_{1,\pm}+c_1)z+ \left(a_{2,\pm}+c_2+a_{1,\pm}c_1+\tfrac12c_1^2\right)z^2.$$

Cauchy's estimate on nested disks $|z|\leq r<\rho<1$ makes its remainder $O(z^3)$. Also $|d(z)|\geq1-\rho$, so $|B_\pm(z)|\leq(1-\rho)^{-1}$; continuity and absence of a branch point give $\eta_{\vartheta,\rho}=\min_{|z|\leq\rho}\operatorname{Re}d(z)>0$. Hence the layer is at most

$$\frac{e^{\operatorname{Re}\bar g_\pi T}}{1-\rho} \exp\left[-\frac{\eta_{\vartheta,\rho}+1-|\vartheta|\rho}{2\varepsilon}T\right].$$

This proves a fixed-frequency $O((m\kappa)^{-3})$ quadratic error plus an exponential layer. The claim is pointwise in a Fourier frequency: using it under an inversion integral still requires $|g_1(u)-g_2(u)|/(m\kappa)\leq\rho$ over the approximated frequency window. This is the matrix-transform form of the two-state occupation problem studied by Pedler (1971).

Fourier inversion is a separate estimate

For Lewis inversion, if the exact shifted transform and an approximation differ by at most $E(u)$ on $[0,R]$, then stopping the approximate integral at $R$ gives

$$|C-\widehat C_R|\leq\frac{\sqrt{S_0K}e^{-rT}}\pi\left[ \int_0^R\frac{E(u)}{u^2+1/4}\,du+ M_{1/2}\{\pi-2\arctan(2R)\}\right], \qquad M_{1/2}=\mathbb E[e^{X/2}].$$

The second term follows from $|\phi(u-i/2)|\leq M_{1/2}$ and is only $O(R^{-1})$. Consequently a cubic compact-frequency transform error does not by itself imply a cubic price error. An admissible window $R=O(\lambda)$ gives only an $O(\lambda^{-1})$ model-independent tail certificate; the diffusive window $R=O(\sqrt\lambda)$ gives only $O(\lambda^{-1/2})$.

Better price rates need model-specific transform decay. The switched variance-gamma certificate proves the sharper, generator-uniform tail $O(R^{-(2T/\nu_{\max}+1)})$; for its posted parameters this is $O(R^{-5})$.

Other inputs

Time-dependent parameters, such as a Hull–White mean level $\theta(t)$ fitted to a yield curve, change only the functions $g_i$, which the engine then holds as Chebyshev series. A terminal exponential payoff $e^{-\beta x_T}$ changes only the starting value of $B$. A continuous fast factor in place of the chain is treated on its own page.

Scope

The expansion assumes an irreducible fast chain, a fixed horizon $t$ as the switching speeds up, and forcing terms $g_i$ smooth in time. The initial layer is part of the expansion. Under these conditions the certificates observe errors falling by one power of $\varepsilon$ per order.

Maturities that grow with the switching rate are a different limit, governed by the principal eigenvalue on the yield-curve page. For option prices the Fourier integral is truncated where the expansion stops applying, which is a choice of the examples rather than part of the theory.

Code and checks