The likelihood of a switching short rate

A short rate observed at discrete times, with a level and volatility switched by a hidden regime much faster than the sampling: the first-order stationary-start transition density, what it identifies, and the cone-corrected test for switching it implies.

The model

The short rate follows a Vasicek model whose mean level and variance depend on a hidden Markov chain $y_t$:

$$dx_t = \kappa\,\big(\theta_{y_t} - x_t\big)\,dt + \sigma_{y_t}\,dW_t .$$

The rate is observed every $\Delta$, giving $x_0, x_1, \dots, x_n$. The chain is not observed. It switches many times between observations, so $\lambda\Delta \gg 1$ for switching rates $\lambda$.

There are three time scales: the regimes last about $1/\lambda$, observations come every $\Delta$, and the rate reverts over $1/\kappa$.

Parameters for the checks: $\kappa = 0.5$, $\theta = (0.08,\ 0.02)$, $\sigma^2 = (0.0009,\ 0.0001)$, monthly observations $\Delta = 1/12$.

The quantity

The calculation below gives the transition density with the regime drawn from its stationary law $\pi$ at the start of each interval and summed at the end. Its sum

$$\ell_C = \sum_{k} \log p_\pi\big(x_{k+1} \mid x_k\big)$$

is a stationary-start transition composite likelihood. The distinction matters: after the regime is hidden, $x_t$ alone is generally not Markov, because the past observations filter the current regime. Fast mixing makes endpoint regime memory exponentially small, but that fact by itself does not make the displayed product the exact observed-data likelihood.

A full likelihood propagates the hidden-state filter. The fixed-panel theorem below shows when the two likelihoods nevertheless agree through first order; it does not make the observations exactly Markov or give a uniform long-panel equivalence.

When the composite likelihood is first-order exact

Write $Q=mQ_0$ and let $K_m(x,x')$ be the matrix transition density: its $(i,j)$ entry is the joint density of $x'$ and ending regime $j$, given $x$ and starting regime $i$. For fixed $\Delta>0$, suppose uniformly on a compact observation set that

$$K_m=p_0\,\mathbf1\pi^\top+m^{-1}K_1+O(m^{-2}), \qquad K_m\mathbf1=p_0\mathbf1+m^{-1}h+O(m^{-2}),$$

with $p_0$ bounded away from zero. This is the regime-resolved version of the fast-switching expansion; the exponentially small endpoint transient is absorbed by the remainder.

If the incoming filter is $\nu_k=\pi+m^{-1}r_k+O(m^{-2})$, where $r_k\mathbf1=0$, its predictive density differs from the stationary-start density only at second order:

$$\nu_kK_m\mathbf1-\pi K_m\mathbf1 =m^{-1}r_k(p_0\mathbf1)+O(m^{-2})=O(m^{-2}).$$

Moreover the Bayes update $\nu_{k+1}=\nu_kK_m/(\nu_kK_m\mathbf1)$ is again $\pi+O(m^{-1})$. Induction therefore gives, for a fixed panel of $N$ transitions started from $\pi$,

$$\ell_F-\ell_C=O(Nm^{-2}),$$

where $\ell_F$ is the exact filtered log likelihood. An arbitrary initial regime prior changes the first predictive density by $O(m^{-1})$; after that observation the filter is $\pi+O(m^{-1})$, so the remaining gap is $O(Nm^{-2})$.

For a growing panel the displayed bound tends to zero only if $N/m^2\to0$, and it is $o(m^{-1})$ only if $N/m\to0$. No long-panel equivalence is asserted.

The second-order likelihood coefficient is explicit. For transition $k$, write

$$K_{m,k}=p_k\mathbf1\pi^\top+m^{-1}K_{1,k}+O(m^{-2}),\qquad h_k=K_{1,k}\mathbf1.$$

For a stationary start set $r_0=0$ and propagate

$$r_{k+1}=\frac{\pi K_{1,k}-(\pi h_k)\pi}{p_k}.$$

The incoming $r_k$ drops out of this recursion because $r_k\mathbf1=0$: the leading kernel maps every starting regime to the same stationary ending law. Consequently, for fixed $N$ and a uniform kernel expansion,

$$\ell_F-\ell_C=\frac1{m^2}\sum_{k=0}^{N-1}\frac{r_kh_k}{p_k}+O(Nm^{-3}).$$

For an arbitrary initial prior $\nu$, the first coefficient is $(\nu-\pi)h_0/p_0$, while the filter entering the second interval has $r_1=\{\nu K_{1,0}-(\nu h_0)\pi\}/p_0$. After the next update, its first-order correction no longer depends on $\nu$. Thus the theorem identifies not just the rates but the pathwise coefficients and the exact point at which initial-prior memory disappears at each order.

Averaging

When the chain switches infinitely fast the rate is Vasicek with averaged parameters, observed as an AR(1):

$$x_{k+1} = \bar\theta + (x_k - \bar\theta)\,e^{-\kappa\Delta} + \sqrt v\,Z_k, \qquad v = \bar s\,\frac{1 - e^{-2\kappa\Delta}}{2\kappa}, \qquad Z_k \sim N(0, 1),$$

with $\bar\theta = \pi\cdot\theta$ and $\bar s = \pi\cdot\sigma^2$. The transition density is the Gaussian $f(x_{k+1})$ with mean $m_k = \bar\theta + (x_k - \bar\theta)e^{-\kappa\Delta}$ and variance $v$.

Reduction to a linear system

The Laplace transform of the next observation is exponential-affine. Try $\mathbb{E}\big[e^{-c\,x_{k+1}} \mid x_k,\ y\big] = e^{-B(\Delta)\,x_k}\,a_y(\Delta)$. With no discounting the coefficient solves $B' = -\kappa B$ from $B(0) = c$, so

$$B(t) = c\,e^{-\kappa t},$$

the same in every regime. What remains is

$$a' = \big(Q + \operatorname{diag} g\big)\,a, \qquad a(0) = \mathbf 1, \qquad g_i(t) = -\kappa\,\theta_i\,B(t) + \tfrac12\,\sigma_i^2\,B(t)^2 .$$

The first-order transition composite likelihood

Let $u = \theta - \bar\theta$ and $w = \sigma^2 - \bar s$ be the fluctuations of level and variance across regimes. The first-order correction to $\log(\pi\cdot a)$ is the Green–Kubo integral of the fluctuating forcing, as on the yield-curve page:

$$\int_0^\Delta \mathcal K\big(\tilde g, \tilde g\big)\,dt, \qquad \tilde g = -\kappa\,u\,B + \tfrac12\,w\,B^2, \qquad \mathcal K(f, h) = -\pi\cdot\big(f\,Q^{\#}h\big).$$

Since $B = c\,e^{-\kappa t}$, this is a polynomial in $c$:

$$\int_0^\Delta \mathcal K(\tilde g, \tilde g)\,dt = A_2\,c^2 + A_3\,c^3 + A_4\,c^4,$$ $$\begin{aligned} A_2 &= \kappa^2\,\mathcal K(u, u)\,\frac{1 - e^{-2\kappa\Delta}}{2\kappa}, \\ A_3 &= -\kappa\,\mathcal K_{\mathrm{sym}}(u, w)\,\frac{1 - e^{-3\kappa\Delta}}{3\kappa}, \\ A_4 &= \tfrac14\,\mathcal K(w, w)\,\frac{1 - e^{-4\kappa\Delta}}{4\kappa}. \end{aligned}$$

A factor $c^n$ on the transform is the $n$-th derivative of the Gaussian density, and that derivative is a Hermite polynomial times the density:

$$f^{(n)}(x) = (-1)^n\,v^{-n/2}\,\mathrm{He}_n(z)\,f(x), \qquad z = \frac{x - m_k}{\sqrt v}.$$

So the transition density and the stationary-start transition composite likelihood are, to first order,

$$\begin{aligned} p(x_{k+1} \mid x_k) &= f(x_{k+1})\Big(1 + \frac{A_2}{v}\,\mathrm{He}_2(z_k) - \frac{A_3}{v^{3/2}}\,\mathrm{He}_3(z_k) + \frac{A_4}{v^2}\,\mathrm{He}_4(z_k)\Big), \\ \ell_C &= \sum_k \Big(\log f(x_{k+1}) + \frac{A_2}{v}\,\mathrm{He}_2(z_k) - \frac{A_3}{v^{3/2}}\,\mathrm{He}_3(z_k) + \frac{A_4}{v^2}\,\mathrm{He}_4(z_k)\Big), \end{aligned}$$ $$\mathrm{He}_2(z) = z^2 - 1, \qquad \mathrm{He}_3(z) = z^3 - 3z, \qquad \mathrm{He}_4(z) = z^4 - 6z^2 + 3 .$$

The terms have plain meanings. $A_2$ widens the density, as extra variance would.

$A_3$ skews it, through the covariance between a regime's level and its volatility. $A_4$ fattens its tails, through the variance of the volatility.

What the stationary-start transition density identifies

The cone-corrected test for switching

With no switching the matrix is zero. The scores for $A_3$ and $A_4$ at zero are sums of Hermite polynomials of the AR(1) residuals, whose sample means are the residual skewness and excess kurtosis:

$$\frac{\partial\ell_C}{\partial A_3} = -\frac{1}{v^{3/2}}\sum_k \mathrm{He}_3(z_k), \qquad \frac{\partial\ell_C}{\partial A_4} = \frac{1}{v^2}\sum_k \mathrm{He}_4(z_k).$$

Since $\mathbb{E}[\mathrm{He}_3^2] = 6$ and $\mathbb{E}[\mathrm{He}_4^2] = 24$ under the null, put

$$Z_3=\sqrt{n}\,\hat S/\sqrt6,\qquad Z_4=\sqrt{n}\,\hat K/\sqrt{24}.$$

The ordinary Jarque–Bera statistic is $Z_3^2+Z_4^2$; it tests both positive and negative excess kurtosis. That is not the score statistic for this model.

Here $A_4\geq0$. After the variance correction $A_2$ is profiled into the Gaussian variance, the closure of the local alternative cone in the skewness–kurtosis coordinates is $\mathbb R\times\mathbb R_+$.

Indeed, for every $A_4>0$ and either sign of $A_3$, the positive-semidefinite constraint $A_3^2\leq4A_2A_4$ can be met by choosing $A_2$ large enough, while $A_2$ is absorbed by the interior variance nuisance parameter. Projecting the Gaussian score onto that cone gives

$$\mathrm{LM}_+=Z_3^2+(Z_4^+)^2 =n\left\{\frac{\hat S^2}{6}+\frac{\max(\hat K,0)^2}{24}\right\}.$$

Thus negative sample excess kurtosis contributes zero, not its square. Under the Gaussian null, with regular estimation of the Gaussian AR parameters, Hermite orthogonality makes the third- and fourth-order scores asymptotically independent standard normals and orthogonal to the nuisance scores. Therefore

$$\mathrm{LM}_+\ \Longrightarrow\ \tfrac12\chi_1^2+\tfrac12\chi_2^2.$$

This is a chi-bar-square law, with 95% critical value $5.138381$. It is the projection geometry of boundary tests described by Self and Liang (1987) and Shapiro (1988). Ordinary Jarque–Bera remains the right omnibus score test for unrestricted nonnormality; $\mathrm{LM}_+$ is the directional test implied by the fast-switching variance-mixture model.

General half-space theorem

The independent Hermite coordinates are convenient, but neither independence nor one-dimensionality of the unrestricted score is essential. Let $Z=(X,Y)\in\mathbb R^p\times\mathbb R$ be the unconstrained Gaussian local estimator after profiling the regular nuisance parameters, with positive-definite covariance

$$\Sigma=\begin{pmatrix}A&b\\b^\top&s\end{pmatrix},\qquad s>0,$$

and local cone $\mathbb R^p\times\mathbb R_+$. Equivalently, $\Sigma^{-1}$ is the Fisher metric.

Its Schur complement $C=A-bb^\top/s$ is positive definite. Put

$$V=Y/\sqrt{s},\qquad U=C^{-1/2}(X-bY/s).$$

Under the null, $U\sim N_p(0,I_p)$ and $V\sim N(0,1)$ are independent. To see the projection explicitly, the block inverse of $\Sigma$ has upper-left block $C^{-1}$ and upper-right block $-C^{-1}b/s$.

If $Y<0$, the Karush–Kuhn–Tucker equations put the projected point on the face $Y=0$ and set its unrestricted coordinate to $X-bY/s$; if $Y\geq0$, the point is already feasible. Its squared Fisher norm is therefore

$$\mathrm{LM}_{\Sigma}=\lVert U\rVert^2+\max(V,0)^2 =\begin{cases} Z^\top\Sigma^{-1}Z,& Y\geq0,\\[2pt] (X-bY/s)^\top C^{-1}(X-bY/s),&Y<0. \end{cases}$$

Consequently, for every $p\geq1$ and every positive-definite $\Sigma$,

$$\mathrm{LM}_{\Sigma}\ \Longrightarrow\ \tfrac12\chi_p^2+\tfrac12\chi_{p+1}^2.$$

The proof is just the block-Gaussian regression identity: $X-bY/s$ is independent of $Y$ with covariance $C$. What is invariant is the Fisher-metric projection, not marginal standardization followed by coordinatewise clipping.

The latter generally has the wrong null law when the coordinates are correlated. The skewness–kurtosis statistic above is the $p=1$, $\Sigma=I$ specialization supplied by Hermite orthogonality.

Several one-sided directions: the orthogonal case

There is a similarly explicit result for a product orthant, but it needs an additional assumption. After profiling nuisances and residualizing the unrestricted scores, suppose the canonical Gaussian score is $(U,V)\in\mathbb R^p\times\mathbb R^q$, where $U\sim N_p(0,I_p)$, $V\sim N_q(0,I_q)$, and the two blocks are independent.

Thus the constrained canonical coordinates are Fisher-orthogonal. For the cone $\mathbb R^p\times\mathbb R_+^q$, metric projection is coordinatewise in this canonical system and

$$\mathrm{LM}_{p,q}=\lVert U\rVert^2+\sum_{r=1}^q(V_r^+)^2 \ \Longrightarrow\ 2^{-q}\sum_{j=0}^q {q\choose j}\chi_{p+j}^2.$$

Indeed, the $q$ signs are independent fair coins. Conditional on exactly $j$ positive coordinates, the sum contains $p+j$ independent squared standard normals; averaging over the binomial number of positive coordinates proves the formula.

Correlation between the unrestricted and constrained blocks does not invalidate it: regress the unrestricted block on the constrained block first. The constrained covariance itself must remain diagonal, so that its positive rescaling preserves the orthant.

The binomial weights are not universal once $q\geq2$. Correlation within the constrained block changes the cone angles and therefore the chi-bar-square weights, as in the multivariate one-sided tests of Kudô (1963) and the cone formulation of Shapiro (1988).

The two-constraint case remains explicit even with unrestricted scores. Partition the profiled Gaussian covariance as

$$\Sigma=\begin{pmatrix}A&B\\B^\top&S\end{pmatrix},\qquad C=A-BS^{-1}B^\top,$$

where $X\in\mathbb R^p$ is unrestricted and $Y\in\mathbb R^2$ is constrained to the positive orthant. Then $U=C^{-1/2}(X-BS^{-1}Y)$ is independent of $Y$, and the projected Fisher norm splits as

$$\mathrm{LM}_{p,2}=\lVert U\rVert^2+ \left\lVert\Pi_{\mathbb R_+^2}^{S^{-1}}(Y)\right\rVert_{S^{-1}}^2.$$

Let $\rho=S_{12}/\sqrt{S_{11}S_{22}}$ be the marginal correlation of the constrained block. Positive diagonal rescaling reduces $S$ to a correlation matrix without changing the orthant. Direct two-dimensional projection gives

$$w_0=\frac14-\frac{\arcsin\rho}{2\pi},\qquad w_1=\frac12,\qquad w_2=\frac14+\frac{\arcsin\rho}{2\pi},$$

and hence the exact null law

$$\mathrm{LM}_{p,2}\ \Longrightarrow\ w_0\chi_p^2+w_1\chi_{p+1}^2+w_2\chi_{p+2}^2.$$

The cross block $B$ changes the residualized unrestricted score but not these weights. In particular, replacing $\rho$ by the correlation of $Y$ conditional on $X$ is wrong: minimizing over the unrestricted cone coordinates leaves the marginal metric $S^{-1}$ on $Y$, not the conditional covariance. This theorem includes the $(1/4,1/2,1/4)$ law when $\rho=0$ and separates the universal one-constraint result from its genuinely covariance-dependent multi-constraint extension.

Three correlated one-sided directions: exact weights

The first higher-dimensional correlated case is still closed form. Let $Y\in\mathbb R^3$ be the constrained block in the same partition, let $R$ be the correlation matrix of its marginal covariance $S$, and let $R^\#$ be the correlation matrix obtained by standardizing $S^{-1}$. Define

$$\Omega(M)=\frac18+\frac{\arcsin M_{12}+\arcsin M_{13}+\arcsin M_{23}}{4\pi}.$$

Then the four chi-bar-square weights, ordered by face dimension, are

$$w_0=\Omega(R^\#),\qquad w_1=\frac12-\Omega(R),\qquad w_2=\frac12-\Omega(R^\#),\qquad w_3=\Omega(R),$$ $$\mathrm{LM}_{p,3}\ \Longrightarrow\ \sum_{j=0}^3w_j\chi_{p+j}^2.$$

To prove the endpoint formulas, positively rescale $Y$ so its covariance is $R$. Projection lands in the interior exactly when $Y>0$. Expanding $\mathbf1_{\{Y>0\}}=\prod_{i=1}^3(1+\operatorname{sgn}Y_i)/8$, central symmetry kills the odd sign moments, and the bivariate quadrant identity of Sheppard (1899) gives $\mathbb E[\operatorname{sgn}Y_i\operatorname{sgn}Y_j]=2\arcsin(R_{ij})/\pi$.

Hence $w_3=\Pr(Y>0)=\Omega(R)$. Projection lands at the vertex exactly when $S^{-1}Y\leq0$; this vector has covariance $S^{-1}$, so symmetry gives $w_0=\Omega(R^\#)$. The conic Gauss–Bonnet identities $w_0+w_2=w_1+w_3=1/2$ give the two remaining weights.

As before, $B$ affects only the independent residualized unrestricted block. The precision correlation enters only the vertex weight: neither $R^\#$ nor $R$ may be replaced by the correlation of $Y$ conditional on the unrestricted scores.

When $R=I_3$, both endpoint probabilities are $1/8$ and the formula reduces to the orthogonal binomial weights $(1,3,3,1)/8$. For four or more correlated constraints, the face-probability characterization remains valid but this elementary arcsine closure no longer follows from the sign expansion.

Arbitrarily many correlated directions: an exact face formula

The individual face probabilities factor in every dimension, giving an active-set form of the multivariate one-sided weights studied by Kudô (1963) and the cone probabilities of Shapiro (1988). Continue to use the marginal constrained covariance $S$, let $A\subseteq\{1,\ldots,q\}$ be the coordinates active in the projection, and let $I=A^c$. Write

$$S_{A\mid I}=S_{AA}-S_{AI}S_{II}^{-1}S_{IA},\qquad \omega(C)=\Pr\{N(0,C)>0\},$$

with the empty-dimensional orthant probability equal to one. Then the probability of landing in the relative interior of face $A$ is exactly

$$w_A=\omega(S_{A\mid I})\,\omega(S_{II}^{-1}),\qquad w_j=\sum_{|A|=j}w_A,$$ $$\mathrm{LM}_{p,q}\ \Longrightarrow\ \sum_{j=0}^q w_j\chi_{p+j}^2.$$

Here $S_{A\mid I}=S_{AA}$ when $I$ is empty, and the second factor is one when $I$ is empty. To prove the factorization, put $P=S^{-1}$ and solve the KKT equations on face $A$. The active projected coordinates are

$$Z_A=Y_A+P_{AA}^{-1}P_{AI}Y_I =Y_A-S_{AI}S_{II}^{-1}Y_I,$$

while the inactive multiplier condition is $S_{II}^{-1}Y_I<0$. The first vector is the Gaussian regression residual with covariance $S_{A\mid I}$ and is independent of $Y_I$; the second has covariance $S_{II}^{-1}$. Central symmetry turns its negative-orthant probability into $\omega(S_{II}^{-1})$, proving the product.

The same formula exposes an exact covariance-inversion polarity. Put $P=S^{-1}$ and apply the face formula under $P$ to the complementary active set $I$. The block-inverse identities

$$P_{I\mid A}=S_{II}^{-1},\qquad P_{AA}^{-1}=S_{A\mid I}$$

exchange the two orthant factors, so face by face

$$w_A(S)=w_{A^c}(S^{-1}),\qquad w_j(S)=w_{q-j}(S^{-1}).$$ $$H_S(t):=\sum_{j=0}^q w_j(S)t^j =t^q H_{S^{-1}}(t^{-1}).$$

Thus covariance inversion reverses, rather than preserves, the chi-bar dimensions. In particular the mean face dimensions satisfy $\sum_j j\,w_j(S)+\sum_j j\,w_j(S^{-1})=q$.

This is the covariance-coordinate version of conic intrinsic-volume polarity, $v_j(C)=v_{q-j}(C^\circ)$; see Amelunxen and Lotz (2017). It is an exact algebraic identity, not an asymptotic or numerical symmetry.

This use of a conditional covariance inside a specified constrained face should not be confused with replacing the entire marginal constrained covariance $S$ by $\operatorname{Cov}(Y\mid X)$ after looking at the unrestricted scores. The latter remains wrong. The cross block $B$ disappears only because $X-BS^{-1}Y$ is first residualized and contributes an independent $\chi_p^2$ term.

The formula is exact for every finite $q$; its computational cost is the evaluation of Gaussian orthant probabilities and a sum over $2^q$ faces. It reduces to the preceding $q=2$ and $q=3$ formulas and to binomial weights when $S$ is diagonal.

Local alternatives and power

The same projection gives an exact local-asymptotic power calculation. Suppose a triangular sequence of alternatives, after the same nuisance-parameter profiling, satisfies

$$\left(Z_3,Z_4\right)\ \Longrightarrow\ \left(G_3+\delta_3,G_4+\delta_4\right),\qquad G_3,G_4\stackrel{\mathrm{iid}}{\sim}N(0,1),\qquad \delta_4\geq0.$$

Then, by continuous mapping,

$$\mathrm{LM}_+\ \Longrightarrow\ (G_3+\delta_3)^2+\max(G_4+\delta_4,0)^2.$$

The general theorem also gives an exact one-dimensional power formula. If the raw local mean is $(\delta_X,\delta_Y)$, put

$$\eta=C^{-1/2}(\delta_X-b\delta_Y/s),\qquad d=\delta_Y/\sqrt{s}, \qquad \Lambda=\lVert\eta\rVert^2.$$

Then $\lVert U\rVert^2$ is noncentral chi-square with $p$ degrees of freedom and noncentrality $\Lambda$, independently of $V\sim N(d,1)$. Writing $\overline F_{p,\Lambda}$ for its survival function, the limiting power at critical value $c$ is

$$\beta_p(\eta,d;c)=\Phi(-d)\overline F_{p,\Lambda}(c) +\int_0^{\sqrt c}\overline F_{p,\Lambda}(c-y^2)\phi(y-d)\,dy +1-\Phi(\sqrt c-d).$$

For $p=1$, $\overline F_{1,\delta_3^2}(a^2)=\Pr\{|G+\delta_3|>a\}$, so this is exactly the skewness–kurtosis formula. It is generally not a noncentral chi-bar-square with fixed one-half weights: when $d>0$, the boundary-face probability is $\Phi(-d)$ rather than $1/2$.

At $(\eta,d)=(0,0)$ the formula reduces exactly to the chi-bar-square null law. The result is confined to the local Gaussian score experiment; it does not assert finite-sample optimality, a rate relating $n$ to switching speed, or equivalence to the full filtered likelihood test.

This is an information-matrix-type construction, related to the full-likelihood test that Carrasco, Hu and Ploberger (2014) show to be optimal for Markov switching, building on White (1982) and Chesher (1984). Because the density used here is a stationary-start composite likelihood, no corresponding full-data optimality claim is made. What the continuous-time model adds is the structure of the coefficients: they are Green–Kubo integrals with explicit time factors, bound by the inequality above.

Tests for switching are usually non-standard because the switching parameters are not identified under the null (Hansen 1992; Garcia 1998; Cho and White 2007). Here the first-order working parameters are Green–Kubo combinations: $A_3$ and $A_4$ are identified by the transition shape, while $A_2$ is profiled with the Gaussian variance. The null remains the vertex of a positive-semidefinite cone, so the kurtosis direction is one-sided.

Unconstrained moment diagnostics follow from the same expansion:

$$\hat A_3 = -\frac{\hat S\,v^{3/2}}{6}, \qquad \hat A_4 = \frac{\hat K\,v^2}{24} .$$

A negative $\hat A_4$ diagnoses a direction outside the model; it is not a feasible switching estimate. A constrained estimator should optimize the full positive-semidefinite Green–Kubo matrix and the Gaussian variance jointly, rather than clipping these two formulas separately.

Checks

The transition density against its exact value by Fourier inversion, as a share of the peak density:

switching rate $\lambda$regime changes per observation $\lambda\Delta$Gaussianfirst order
201.714.37%5.40%
403.37.64%2.05%
806.73.67%0.55%

The first-order error falls by about four for each doubling of $\lambda$ once $\lambda\Delta$ is several units, as a second-order remainder should. With fewer than two regime changes per observation the expansion is poor.

The two three-state chains of the yield-curve page, built to share averages and Green–Kubo matrix, observed quarterly with $\kappa = 1$:

speed factorchain A vs Gaussianchain A vs chain B
201.751%0.1293%
400.876%0.0320%
800.437%0.0079%

Their densities differ at second order, while each differs from the Gaussian at first order.

A separate cone-test certificate simulates 20,000 Gaussian samples of size 4,000. At the chi-bar-square 95% cutoff, $\mathrm{LM}_+$ rejects 5.46% while the unprojected square rejects 7.45%; sample excess kurtosis is negative 52.29% of the time. For a standardized uniform sample, ordinary Jarque–Bera is 11964.19 because excess kurtosis is $-1.19819$, whereas the cone statistic is only 0.41.

An independent simulation of two million limiting Gaussian score pairs agrees with the power integral within $9.6\times10^{-4}$. At local shifts $(\delta_3,\delta_4)=(1,0),(0,1),(1,1),(0,2)$, the exact directional powers are respectively $0.141681$, $0.177730$, $0.283486$, and $0.490980$; the size-matched ordinary Jarque–Bera powers are $0.132710$, $0.132710$, $0.225545$, and $0.415427$.

The same certificate tests correlated Gaussian coordinates with correlation $0.75$. The canonical formula and the direct piecewise Fisher-metric projection agree within $1.42\times10^{-14}$.

At the same chi-bar cutoff, the metric statistic rejects $5.0056\%$ under the null, whereas marginal standardization followed by naive clipping rejects $5.9701\%$. For a raw mean shift corresponding to canonical shift $(1,1)$, simulated power is $0.283067$ against the exact value $0.283486$.

For the general theorem, the certificate also uses three unrestricted coordinates and one one-sided coordinate with a dense $4\times4$ covariance. The 95% point of $\tfrac12\chi_3^2+\tfrac12\chi_4^2$ is $8.761053$.

Across 1.5 million null draws, rejection is $4.9966\%$ and the canonical statistic agrees with the direct piecewise Fisher projection within $4.98\times10^{-14}$. At canonical mean $(0.6,-0.8,0.4,0.9)$, the noncentral-$\chi^2$ integral gives power $0.202280$ and simulation gives $0.203155$.

For the orthant extension, one million draws with $p=2$, $q=3$, dense cross-block covariance, and a diagonal constrained covariance block give face frequencies $(0.125228,0.374093,0.374995,0.125684)$ against exact binomial weights $(1,3,3,1)/8$. The 95% critical value is $8.946572$ and null rejection is $5.0178\%$. A second million-draw experiment deliberately correlates two constrained scores at $\rho=0.70$.

Its exact weights are $(0.126592,0.5,0.373408)$; the simulated face frequencies are $(0.126491,0.499893,0.373616)$. With $p=2$, its correctly calibrated 95% critical value $8.410794$ gives $5.0320\%$ rejection. Incorrectly imposing the independent-binomial cutoff $8.022108$ raises rejection to $5.9396\%$.

The conditional constrained correlation is $-0.011608$ and its cutoff gives $5.9540\%$ rejection, directly checking that the theorem uses marginal $S$.

A third million-draw experiment uses $p=2$, three correlated constraints, and a dense cross block. The exact weights $(0.059323,0.320310,0.440677,0.179690)$ agree with simulated face frequencies $(0.059344,0.320617,0.440516,0.179523)$ to $3.1\times10^{-4}$.

The exact critical value $9.305539$ gives $4.9824\%$ rejection, while the independent-binomial cutoff $8.946572$ gives $5.7718\%$. Thus the new formula is numerically calibrated and the binomial shortcut is again rejected.

Finally, a dense four-constraint experiment checks all 16 active sets in the general formula. Its exact weights are $(0.018400,0.156882,0.363697,0.343118,0.117902)$ and the simulated frequencies are $(0.018262,0.157558,0.362940,0.343282,0.117958)$, with maximum dimension-weight error $7.57\times10^{-4}$; the maximum error over the 16 individual face probabilities is $5.33\times10^{-4}$. The exact 95% cutoff $10.398896$ rejects at $5.0215\%$; imposing the independent-binomial weights rejects at $6.2549\%$.

Repeating the only numerical four-dimensional orthant integrations with four randomizations changes any weight by at most $7.1\times10^{-8}$, and the even/odd Gauss–Bonnet sums equal $1/2$ within $3.5\times10^{-8}$. Under covariance inversion, the maximum complementary-face and reversed-dimension discrepancies are $2.72\times10^{-8}$, at the scale of that integration, while the two block identities agree to $6.66\times10^{-16}$.

A separate matrix-kernel certificate propagates the exact two-state filter by Fourier inversion for a fixed six-transition path. Doubling $m$ through $40,80,160,320$ gives observed orders $2.045$ for the stationary-start log-likelihood gap, $2.043$ for the specified-start tail after the first observation, and $1.019$ for the filter's distance from $\pi$. The specified-start first-observation gap has order $0.906$, confirming the stated $O(m^{-1})$ initial-prior caveat.

Refining both Fourier and time grids changes the fastest-rate log gap by less than $4\times10^{-15}$. The kernel formula predicts coefficients $-2.529716$, $-0.183148$, and $-2.530169$ for the stationary gap, specified-start first interval, and specified-start tail. Richardson extrapolation of the exact filtered likelihood gives $-2.528356$, $-0.183174$, and $-2.528853$, respectively; the largest discrepancy is $1.4\times10^{-3}$.

Moments of the AR(1) residuals from 200,000 exactly simulated monthly observations, against the first-order predictions $-6A_3/v^{3/2}$, $24A_4/v^2$ and $1 + 2A_2/v$:

$\lambda$skewnesspredictedexcess kurtosispredictedvariance ratiopredicted
60+0.0821+0.09300.33730.38421.00991.0075
120+0.0361+0.04650.17910.19211.00481.0037
240+0.0249+0.02320.09210.09611.00031.0019

Scope

The expansion needs many regime changes between observations. When regimes last longer than the sampling interval, the first-order cancellation above is unavailable and the exact likelihood needs the filter of Hamilton (1989). Even in the fast-switching limit, the equivalence theorem is pointwise on compact observation sets with predictive densities bounded away from zero; it is not a tail-uniform likelihood theorem.

A switching variance alone gives the classic normal variance mixture of Clark (1973). For closed-form likelihood expansions of diffusions without regimes see Aït-Sahalia (2002).

Certificates: transition expansion, filtered likelihood, and cone test.