← summary

The direction of a regime cycle in the Heston smile

A chain or continuous periodic factor that cycles in one direction contributes a commutator correction invisible to averages and symmetric covariances.

The smile from fast volatility

When volatility fluctuates fast and independently of the stock, the option price is Black–Scholes averaged over integrated variance (Hull and White 1987). Expanding that average around the mean variance gives a smile at first order in the fluctuation time scale.

Sircar and Papanicolaou (1999), section 3.3, derive it for a stationary ergodic volatility driver whose integrated fluctuations obey a central limit theorem, which a finite-state Markov chain does. Ergodicity alone does not fix the rate: long memory changes the power of the small parameter. The implied volatility is quadratic in log-moneyness, with a curvature set by the variance of integrated variance. Fouque, Lorig and Sircar (2016) carry the fast-factor expansion to second order with correlation and a terminal layer.

In that work the fast driver is a one-dimensional diffusion or a generic ergodic process whose correlation function enters through a variance. A stationary scalar diffusion on an interval or the line is reversible under the usual zero-current boundary conditions. Dimension alone is not enough: a diffusion on a periodic domain can carry stationary current, break detailed balance, and have an antisymmetric Green–Kubo matrix.

A continuous cycle in closed form

Let the fast factor be drifted Brownian motion on the circle,

$$dY_t=c\,dt+\sqrt{2D}\,dW_t\pmod {2\pi},\qquad D>0.$$

Its stationary law is uniform, but its stationary current is $c/(2\pi)$ and time reversal changes $c$ to $-c$. For every Fourier mode $n\ge1$, use the centered feature pair $\phi_{n,c}(y)=\cos(ny)$ and $\phi_{n,s}(y)=\sin(ny)$. Fourier orthogonality makes every cross-mode entry vanish, while the block for mode $n$ is exactly

$$K_n=\frac1{2\{D^2n^4+c^2n^2\}} \begin{pmatrix}Dn^2&cn\\-cn&Dn^2\end{pmatrix}.$$

Indeed, the off-diagonal correlation is $\mathbb E[\cos(nY_0)\sin(nY_t)]=\tfrac12e^{-Dn^2t}\sin(cnt)$; the reverse ordering has the opposite sign. Thus $c\mapsto-c$ preserves $K_n^{\rm sym}$ and flips $K_n^{\rm anti}$ for every mode.

If $\cos Y$ and $\sin Y$ multiply two fixed operators $A_c,A_s$, their first correction is

$$\frac{D}{2(D^2+c^2)}(A_c^2+A_s^2) +\frac{c}{2(D^2+c^2)}[A_c,A_s].$$

For $D=1,c=2$, the block is $K_1=\left(\begin{smallmatrix}0.1&0.2\\-0.2&0.1\end{smallmatrix}\right)$. Direct correlation quadrature agrees below $9\times10^{-16}$.

A periodic CTMC discretization converges to this block at rate 2.001, and an independent coupled-operator evolution gives a 2.002 rate for the residual after the full first-order rule. Keeping only the symmetric block leaves a first-order error; the forward-minus-reverse residual after the commutator correction converges at rate 2.000.

Continuous-factor certificate: verify_circle_factor.py. The physical interpretation of current and broken detailed balance is standard in periodic Brownian-motor models (Reimann 2002).

The theorem for every stationary Markov driver

Let $P_t$ be a stationary Markov semigroup on the real space $L^2(\pi)$, and let $P_t^{\leftarrow}=P_t^\dagger$ be its stationary time reversal. For centered observables $f,g$ whose correlations are integrable, define

$$K^{\rightarrow}(f,g)=\int_0^\infty\langle f,P_tg\rangle_\pi\,dt.$$

Semigroup adjointness gives, without a diffusion or finite-state assumption,

$$K^{\leftarrow}(f,g)=K^{\rightarrow}(g,f).$$

Thus stationary time reversal preserves the symmetric Green–Kubo form and flips the antisymmetric form. This is the integrated-correlation consequence of the general Markov time-reversal framework of Nagasawa (1964).

Exact occupation-law blindness at every switching rate

The preceding result has an exact, non-asymptotic counterpart. Under a stationary start, a process and its stationary time reversal have the same law for the occupation measure $O_T(B)=\int_0^T\mathbf 1_{\{Y_s\in B\}}\,ds$, because reversing a path preserves the total time it spends in every set. Hence, for every measurable state observable $f$ for which the transform exists,

$$\int_0^T f(Y_s^{\leftarrow})\,ds\ \overset d=\ \int_0^T f(Y_s)\,ds.$$

For a finite-state chain in row-generator convention, this follows directly from $Q^{\leftarrow}=\Pi^{-1}Q^\top\Pi$. If $D_z$ is any diagonal real or complex state potential, then

$$\pi e^{T(Q^{\leftarrow}+D_z)}\mathbf 1 =\mathbf 1^\top e^{T(Q^\top+D_z)}\pi^\top =\pi e^{T(Q+D_z)}\mathbf 1.$$

This is the full joint Laplace/characteristic-transform identity for the state occupation times. Thus no stationary payoff depending only on total state occupations can identify cycle direction at any switching rate, not just at first Green–Kubo order. The tilted generators are similar to transposes, so their principal eigenvalues—and therefore scalar long-time Feynman–Kac growth rates—are also identical.

Temporal ordering is the missing information. For a time-dependent loading $g$, path reversal gives

$$\mathbb E_{\pi,Q^{\leftarrow}}\exp\!\left\{\int_0^Tg(s,Y_s)\,ds\right\} =\mathbb E_{\pi,Q}\exp\!\left\{\int_0^Tg(T-s,Y_s)\,ds\right\}.$$

Constant or time-symmetric loadings remain blind, whereas an asymmetric loading schedule can distinguish the two directions. Endpoint conditioning, ordered events such as first-default races, and noncommuting switched operators can do so as well. The finite-state certificate checks the constant-potential transform to $7.51\times10^{-16}$, the principal eigenvalue to $8.88\times10^{-16}$, and the reversed-schedule identity to $4.44\times10^{-16}$; applying the same asymmetric two-stage schedule leaves a directional gap of $0.01976$.

There is also a precise converse. Suppose $-L$ is invertible on the centered subspace with bounded inverse, as for an irreducible finite chain and for the usual spectral-gap setting.

Then $K(f,g)=\langle f,(-L)^{-1}g\rangle_\pi$ is symmetric for every centered pair if and only if $(-L)^{-1}$, hence $L$, is self-adjoint—equivalently, the process is reversible. This full-space statement must not be inferred from a selected feature matrix: a one-feature Green–Kubo matrix is always symmetric, even for an irreversible process.

Two generic features detect every finite-state irreversible chain

The preceding caveat is sharp for one feature but not generic for two. For an irreducible finite-state chain, let $\mathcal H_0$ be the real centered $L^2(\pi)$ space, put $R=(-L)^{-1}|_{\mathcal H_0}$, and define its skew-adjoint part

$$J=\operatorname{Anti}R=\tfrac12(R-R^\dagger).$$

Choose an orthonormal basis of $\mathcal H_0$ and draw two observables $F,G$ with independent standard Gaussian coordinates in that basis. Their antisymmetric Green–Kubo signal is

$$\Xi=K^{\rm anti}(F,G)=\langle F,JG\rangle_\pi.$$

Then the chain is irreversible if and only if $J\ne0$, and in that case

$$\Pr(\Xi=0)=0,\qquad \mathbb E\Xi^2=\|J\|_{\rm HS}^2,$$ $$\mathbb E e^{it\Xi} =\det(I+t^2J^\dagger J)^{-1/2}.$$

Indeed, symmetry of $R$ is equivalent to symmetry of its inverse $-L$, hence to reversibility. If $J\ne0$, its kernel is a proper subspace, so $JG\ne0$ almost surely; conditional on such a $G$, the scalar $\langle F,JG\rangle$ is a nondegenerate Gaussian and is therefore nonzero almost surely. The moment and characteristic-function formulas follow by conditioning first on $G$ and evaluating a Gaussian quadratic form.

Thus two generic continuous probes identify finite-state irreversibility almost surely, although a deliberately chosen pair can still miss it. This is a population operator statement, not a finite-sample test or a lower bound on the signal for prescribed features.

The squared-signal identity is a bilinear analogue of stochastic trace probing; compare Hutchinson (1989). In the four-state certificate, $\|J\|_{\rm HS}=0.216796492$ and the exact second moment is $0.047000719$; 400,000 independent feature pairs give $0.047160377$. At $t=0.75/\|J\|_{\rm op}$, the exact characteristic function is $0.64$ and the simulation gives $0.640300132$ with imaginary part below $6.3\times10^{-4}$.

A finite random-probe certificate

The almost-sure statement has a quantitative finite-probe version. Let $s_i$ be the singular values of $J$, put $S_2=\sum_i s_i^2=\|J\|_{\rm HS}^2$ and $S_4=\sum_i s_i^4$, and retain independent standard Gaussian probes. Direct expansion of the bilinear Gaussian form gives

$$\mathbb E\Xi^4=3S_2^2+6S_4.$$

Because $J$ is real and skew-adjoint, its nonzero singular values occur in equal pairs. Hence $S_4/S_2^2\le1/2$, so $\mathbb E\Xi^4\le6S_2^2$ and $\operatorname{Var}(\Xi^2)\le5S_2^2$. Applying the inequality of Paley and Zygmund (1932) to $\Xi^2$ yields, for $0\le\theta<1$,

$$\Pr\!\left(|\Xi|\ge\theta\|J\|_{\rm HS}\right) \ge \frac{(1-\theta^2)^2}{6}.$$

At $\theta=1/2$, each pair detects a signal of at least half the root-mean-square scale with probability at least $3/32$. Thus $m$ independent pairs miss that threshold with probability at most $(29/32)^m$; in particular, 47 pairs give the explicit bound $(29/32)^{47}=0.0097872<1\%$. The empirical energy $m^{-1}\sum_{k=1}^m\Xi_k^2$ is unbiased for $\|J\|_{\rm HS}^2$, with relative variance at most $5/m$ and consequently the Chebyshev bound

$$\Pr\!\left(\left|m^{-1}\sum_{k=1}^m\Xi_k^2/S_2-1\right|\ge\varepsilon\right) \le \frac{5}{m\varepsilon^2}.$$

Action probes need five times fewer samples

The bilinear design is appropriate when an oracle returns only scalar entries $\langle x,Jy\rangle$. If the full vector action $Jx$ is available, the ordinary quadratic trace probe is sharper.

For $x\sim N(0,I)$ put $Z_x=\lVert Jx\rVert^2$. Then

$$\mathbb E Z_x=S_2,\qquad \operatorname{Var}(Z_x)=2S_4\le S_2^2.$$

The final inequality again uses the paired nonzero singular values of a real skew-adjoint matrix. Thus the empirical action-probe energy is unbiased with relative variance at most $1/m$, compared with $5/m$ for squared bilinear Gaussian probes. Paley–Zygmund also gives

$$\Pr\!\left(\lVert Jx\rVert\ge\theta\lVert J\rVert_{\rm HS}\right) \ge { (1-\theta^2)^2\over 1+2S_4/S_2^2} \ge { (1-\theta^2)^2\over2}.$$

At half the Hilbert–Schmidt scale the universal detection probability is at least $9/32$, so 14 independent action probes have miss probability at most $(23/32)^{14}=0.009820<1\%$. For bounded sign vectors the exact variance is $2\{S_4-\sum_i(J^\dagger J)_{ii}^2\}\le S_2^2$, so the same universal bounds apply, although discrete probes can still lie in the kernel with positive probability.

In the four-state certificate the Gaussian relative variance is exactly $1$ and the simulation gives $1.000353$; half-scale detection is $0.779485$. The sign-action relative variance is $0.289570$, reproduced by exhaustive enumeration, with detection probability $0.75$. This is an operator-action result: it does not improve the scalar-query complexity when only bilinear entries are observable, and it does not absorb error from estimating $J$ from a hidden or noisy trajectory.

Fixed-norm action probes remove radial noise

Gaussian directions spend some of their variance on their random length. Let $d\ge2$ be the centered-space dimension, let $u$ be uniform on the sphere of radius $\sqrt d$, and put $Z_u=\lVert Ju\rVert^2$.

Since $\mathbb E[uu^\top]=I$, this remains an unbiased trace estimator: $\mathbb E Z_u=S_2$. The fourth-order spherical moment tensor gives the exact variance

$$\operatorname{Var}(Z_u) ={2\over d+2}\{dS_4-S_2^2\},\qquad {\operatorname{Var}(Z_u)\over S_2^2} \le {d-2\over d+2}.$$

The last inequality again uses $S_4/S_2^2\le1/2$ for a real skew-adjoint matrix. It is sharp for rank two; for $d=2$, $J^\top J$ is a scalar matrix and one fixed-norm action recovers the energy exactly. Paley–Zygmund therefore yields the dimension-dependent guarantee

$$\Pr\!\left(\lVert Ju\rVert\ge\theta\lVert J\rVert_{\rm HS}\right) \ge { (1-\theta^2)^2\over1+\operatorname{Var}(Z_u)/S_2^2} \ge {(1-\theta^2)^2(d+2)\over2d}.$$

At $d=3$ and $\theta=1/2$, the universal detection probability is at least $15/32$. Hence eight independent fixed-norm actions have miss probability at most $(17/32)^8=0.006344<1\%$.

The relative-variance constant is $1/5$, so the same median-of-means theorem below uses blocks of size four and 100 actions in total at $\varepsilon=1/2$, $\delta=0.05$—one quarter of the Gaussian-action certificate. The four-state simulation gives relative variance $0.200079$ against the exact value $0.2$, half-scale detection $0.912885$, zero failures in 4,000 median-of-means trials, and maximum relative error $0.195953$.

This is a rotation-invariant randomized trace estimator for $J^\top J$, in the general framework analyzed by Avron and Toledo (2011). The dimension-dependent improvement comes here from combining the exact spherical moment with the paired spectrum of a skew operator. It retains the same action-oracle and exact-population-operator scope as the Gaussian result; normalizing directions does not address scalar-only access or generator-estimation error.

Orthogonal action blocks sample without replacement

Independent fixed-norm probes can be improved by making the directions orthogonal inside each block. Let $U\in\mathbb R^{d\times k}$ have Haar-uniform orthonormal columns, put $P=UU^\top$, and define the $k$-action estimator

$$Z_k={d\over k}\operatorname{tr}(U^\top J^\top JU) ={d\over k}\sum_{a=1}^k\lVert Ju_a\rVert^2.$$

Haar invariance gives $\mathbb E P=(k/d)I$, so $Z_k$ is unbiased. The second moment tensor of a rank-$k$ Haar projector gives, with $S_2=\operatorname{tr}(J^\top J)$ and $S_4=\operatorname{tr}\{(J^\top J)^2\}$,

$${\operatorname{Var}(Z_k)\over S_2^2} ={2d(d-k)\over k(d-1)(d+2)} \left({S_4\over S_2^2}-{1\over d}\right) \le{(d-k)(d-2)\over k(d-1)(d+2)}.$$

For completeness, orthogonal invariance forces $\mathbb E[P_{ij}P_{\ell m}]=a\delta_{ij}\delta_{\ell m} +b(\delta_{i\ell}\delta_{jm}+\delta_{im}\delta_{j\ell})$. The identities $\operatorname{tr}P=k$ and $\operatorname{tr}(P^2)=k$ determine $b=k(d-k)/\{d(d-1)(d+2)\}$ and $a-k^2/d^2=-2b/d$. Contracting this tensor with $A=J^\top J$ proves $\operatorname{Var}\{\operatorname{tr}(AP)\}=2b[\operatorname{tr}(A^2)-(\operatorname{tr}A)^2/d]$, which is exactly the displayed formula after multiplying by $(d/k)^2$.

The inequality uses the paired nonzero singular values of the real skew operator, and is sharp at rank two. Averaging $k$ independent spherical actions has worst relative variance $(d-2)/\{k(d+2)\}$; orthogonality multiplies it by the finite-population factor $(d-k)/(d-1)$.

Thus a full $d$-frame returns the trace exactly. This exact quadratic specialization complements the broader orthogonal-sampling variance-reduction program of Choromanski, Rowland, Chen and Weller (2019).

The four-state certificate has centered dimension $d=3$. For a Haar two-frame, write its orthogonal complement as the unit vector $v$.

Then $Z_2=(3/2)\{S_2-v^\top J^\top Jv\}$. Every nonzero three-dimensional real skew matrix has $\operatorname{spec}(J^\top J)=(s^2,s^2,0)$, so this is the deterministic interval

$$ {3\over4}S_2\le Z_2\le {3\over2}S_2.$$

Consequently two orthogonal actions alone already give relative error at most $1/2$ for this example; this is stronger than a probability bound. The exact relative variance is $1/20$, and 400,000 simulated Haar frames give $0.050020$, with observed energy ratios in $[0.75,1.5]$. The generic median-of-means construction uses 25 one-frame blocks, hence 50 actions rather than 100 independent spherical actions or 400 Gaussian actions; its 16,000 simulated repetitions have zero failures and maximum relative error $0.259118$.

These gains require a full action oracle and coordinated directions in a common centered-space basis. They do not improve scalar bilinear access, do not make a marginal random-feature guarantee conditional on prescribed features, and do not include error from estimating the population operator $J$.

Exponential confidence from median block means

The polynomial Chebyshev tail can be amplified without assuming any higher moments. Split the squared signals into $k$ independent blocks of size $b$, take the mean inside each block, and let $\widehat S_2$ be the median of those block means. If

$$b\ge {4v\over\varepsilon^2},\qquad k\ge8\log(1/\delta),$$

with $k$ odd, where $v=5$ for Gaussian probes, then

$$\Pr\!\left(\left|{\widehat S_2\over S_2}-1\right| >\varepsilon\right)\le e^{-k/8}\le\delta.$$

Indeed, the variance calculation above and Chebyshev make each block bad with probability at most $1/4$. The median is bad only if at least half the blocks are bad, and Hoeffding's Bernoulli bound gives $e^{-k/8}$. This is the classical repeat-and-take-the-median confidence amplification used by Jerrum, Valiant and Vazirani (1986), applied here to squared skew-resolvent probes.

For Rademacher probes in centered dimension $d$, the sharper variance constant is $v_d=5-8/d$, so the same theorem uses $b\ge4v_d/\varepsilon^2$. At $\varepsilon=1/2$ and $\delta=0.05$, $k=25$ gives the rigorous failure bound $e^{-25/8}=0.0439369$. The Gaussian design uses $b=80$ and 2,000 probes in total.

For the four-state example, $d=3$, so bounded sign probes use $b=38$ and 950 probes. Across 200 Gaussian and 400 sign-probe certificate repetitions, the observed failure frequencies are both zero; maximum relative errors are $0.158228$ and $0.188130$, respectively.

For Gaussian action probes the variance constant is only $v=1$, so the same theorem uses $b=16$ and 400 actions rather than 2,000 bilinear pairs. Across 1,000 certificate repetitions the observed failure frequency is zero and the maximum relative error is $0.208269$.

For the four-state certificate the paired-singular-value bound is attained: $S_4/S_2^2=1/2$ and the exact relative variance is $5$. The exact/simulated fourth moments are $0.013254405/0.013425214$, and the observed half-RMS detection frequency is $0.4933025$ versus the universal $0.09375$ lower bound. These guarantees concern finitely many randomized probes of the exact population operator $J$.

They are not conditional guarantees for prescribed features. The next two results separate the deterministic plug-in map from a finite-path concentration bound for a fully observed chain.

Deterministic stability under generator estimation

For a row-convention irreducible $n$-state generator $Q$, write $\pi Q=0$, $P=\mathbf1\pi$ and

$$G=Q^\#=(Q-P)^{-1}+P.$$

Let $D=\operatorname{diag}\pi$. The canonical Euclidean version of the centered Green–Kubo resolvent and its directional part are

$$\mathcal R=-D^{1/2}GD^{-1/2},\qquad \mathcal J=\operatorname{Anti}\mathcal R.$$

The stationary vector $\sqrt\pi$ is in the kernel, and on its orthogonal complement this is exactly $(-L)^{-1}$ in $L^2(\pi)$. Thus $\|\mathcal J\|_{\rm HS}$ is the basis-independent energy used above.

Let $\widehat Q$ be an estimated irreducible generator and suppose a separate statistical argument supplies the operator-norm guarantee

$$\|\widehat Q-Q\|_2\le\varepsilon.$$

Put $\widehat P=\mathbf1\widehat\pi$, $\widehat G=\widehat Q^\#$, and $\widehat Z=\widehat Q-\widehat P$. The stationary-law identity in the operational direction is

$$\pi-\widehat\pi=\pi(\widehat Q-Q)\widehat G,\qquad \|\pi-\widehat\pi\|_2\le \beta_\pi:=\varepsilon\|\widehat G\|_2.$$

The universal bound above discards the geometry of the observed stationary law. Put $r=\varepsilon\|\widehat G\|_2$. Since $\pi=\widehat\pi+(\pi-\widehat\pi)$, the same identity also gives, when $r<1$,

$$\|\pi-\widehat\pi\|_2 \le {r\|\widehat\pi\|_2\over1-r}.$$

Hence the entirely estimator-centered fixed-point radii

$$\beta_2=\begin{cases} \min\left\{r,{r\|\widehat\pi\|_2\over1-r}\right\},&r<1,\\ r,&r\ge1, \end{cases}$$ $$\beta_i=\min\left\{\beta_2, \varepsilon(\|\widehat\pi\|_2+\beta_2) \|\widehat G e_i\|_2\right\}$$

satisfy $\|\pi-\widehat\pi\|_2\le\beta_2$ and $|\pi_i-\widehat\pi_i|\le\beta_i$. The second formula retains the different sensitivities of the stationary coordinates instead of replacing every column of $\widehat G$ by its operator norm.

Consequently $\|P-\widehat P\|_2\le\sqrt n\,\beta_2$. Define

$$\eta_Z=\varepsilon+\sqrt n\,\beta_2,\qquad \widehat a=\|\widehat Z^{-1}\|_2.$$

If $p_i^-:=\widehat\pi_i-\beta_i>0$ for every $i$ and $\widehat a\eta_Z<1$, Neumann inversion gives the entirely estimator-centered bound

$$\|G-\widehat G\|_2\le B_G:= {\widehat a^2\eta_Z\over1-\widehat a\eta_Z} +\sqrt n\,\beta_2.$$

For the stationary square-root weights, set

$$p_- =\min_i p_i^-,\qquad p_+=\max_i(\widehat\pi_i+\beta_i),$$ $$d_+=\max_i{\beta_i\over\sqrt{p_i^-}+\sqrt{\widehat\pi_i}},\qquad d_-=\max_i{\beta_i\over \sqrt{p_i^-\widehat\pi_i} (\sqrt{p_i^-}+\sqrt{\widehat\pi_i})}.$$

A three-term product expansion, followed by contraction under antisymmetrization, now yields the a posteriori radius

$$\|\widehat{\mathcal J}-\mathcal J\|_{\rm HS}\le\delta_{\rm post},$$ $$\delta_{\rm post}=\sqrt n\left\{ \sqrt{p_+/p_-}\,B_G +{d_+\|\widehat G\|_2\over\sqrt{p_-}} +\sqrt{\widehat\pi_{\max}}\|\widehat G\|_2d_- \right\}.$$

Every quantity on the right is computable from $\widehat Q$ and $\varepsilon$. Therefore, if $S=\|\mathcal J\|_{\rm HS}^2$ and $\widehat S=\|\widehat{\mathcal J}\|_{\rm HS}^2$,

$$\{\max(\sqrt{\widehat S}-\delta_{\rm post},0)\}^2\le S\le (\sqrt{\widehat S}+\delta_{\rm post})^2.$$

This separates deterministic plug-in bias from random-probe error. It is a finite perturbation result, not only a derivative at $\varepsilon=0$. Its group-inverse step is placed against Meyer's finite-chain perturbation theory.

The certificate checks 16 generator-perturbation directions at six scales, for 96 cases total. Actual skew-resolvent errors have fitted orders in $[0.996415,1.003445]$; the a posteriori bounds have orders in $[1.012001,1.016636]$. Relative to the earlier universal choice $\beta_2=r$, the fixed-point stationary radii are $0.504309$–$0.510877$ as large and the final skew radii are $0.602191$–$0.619885$ as large.

The largest observed actual/bound ratios are $0.016276$ for the skew resolvent, $0.858415$ for the stationary law in Euclidean norm, $0.859699$ coordinatewise, and $0.078440$ for the group inverse. The maximum Neumann ratio is $0.044292$, the estimator-centered stationary identity holds within $5.41\times10^{-16}$, and every energy interval contains the population value. The bound is conservative and does not claim sharp constants.

A fully observed path supplies the missing generator radius

Suppose the complete CTMC path is observed, including exact jump times. Write

$$L_i=\int_0^{T_{\rm obs}}\mathbf1_{\{Y_s=i\}}\,ds, \qquad N_{ij}=\#\{i\to j\hbox{ jumps on }[0,T_{\rm obs}]\},$$

and set $\widehat q_{ij}=N_{ij}/L_i$ off the diagonal, with diagonal entries chosen to make row sums zero. Fix before observing the path a finite set $\Theta\subset(0,\infty)$. For each counted transition,

$$\exp\!\left\{\pm\theta(N_{ij}-q_{ij}L_i) -q_{ij}L_i(e^{\pm\theta}-1\mp\theta)\right\}$$

is a mean-one exponential martingale, the multiplicative-intensity construction underlying Aalen's counting-process inference. Markov's inequality applied at both signs and a union bound yield intervals without any known rate cap. Put $x=\log\{2n(n-1)|\Theta|/\delta\}$ and, for $i\ne j$,

$$\ell_{ij}=\max\left\{0,\max_{\theta\in\Theta} {\theta N_{ij}-x\over(e^\theta-1)L_i}\right\},\qquad u_{ij}=\min_{\theta\in\Theta} {x+\theta N_{ij}\over(1-e^{-\theta})L_i}.$$

Then, simultaneously for every directed edge,

$$\Pr\{\ell_{ij}\le q_{ij}\le u_{ij}\text{ for all }i\ne j\}\ge1-\delta.$$

No independence between transition counts is assumed. Define the observable entry radius

$$a_{ij}=\max\{|\widehat q_{ij}-\ell_{ij}|, |u_{ij}-\widehat q_{ij}|\},\qquad s_i=\sum_{j\ne i}a_{ij}.$$

The diagonal error in row $i$ is at most $s_i$. Hence the Frobenius and induced-norm calculations give

$$\|\widehat Q-Q\|_2\le\varepsilon_{\rm path}:=\min\left\{ \sqrt{\sum_{i\ne j}a_{ij}^2+\sum_i s_i^2}, \sqrt{B_1B_\infty}\right\},$$ $$B_\infty=2\max_i s_i,\qquad B_1=\max_j\left\{s_j+\sum_{i\ne j}a_{ij}\right\}.$$

For a small chain the rectangular rate set can instead be propagated exactly. Define the asymmetric off-diagonal error endpoints

$$\underline e_{ij}=\widehat q_{ij}-u_{ij},\qquad \overline e_{ij}=\widehat q_{ij}-\ell_{ij},$$

and, for any choice $e=(e_{ij})_{i\ne j}$, let $E(e)$ have off-diagonal entries $e_{ij}$ and diagonal entries $E_{ii}=-\sum_{j\ne i}e_{ij}$. Then

$$\varepsilon_{\rm box} =\max_{\substack{e_{ij}\in\{\underline e_{ij},\overline e_{ij}\}\\i\ne j}} \|E(e)\|_2 =\max_{\substack{\underline e_{ij}\le e_{ij}\le\overline e_{ij}\\i\ne j}} \|E(e)\|_2.$$

The equality is exact: every point of the box is a convex combination of its vertices, $E(e)$ is affine, and the spectral norm is convex. Thus the same simultaneous confidence event gives $\|\widehat Q-Q\|_2\le\min(\varepsilon_{\rm path},\varepsilon_{\rm box})$.

The enumeration costs $2^{n(n-1)}$ spectral norms, so it is useful only for small state spaces; the displayed Frobenius/induced-norm formula remains the polynomial-time fallback. This step tightens only deterministic propagation of the rate intervals and does not change their statistical coverage.

On the same event, substituting the smaller available radius into the preceding a posteriori theorem gives a $1-\delta$ confidence interval for the population skew energy whenever the estimated generator is irreducible and its displayed positivity and Neumann conditions hold. If some $L_i=0$, this construction supplies no finite radius for that row.

The certificate checks the exponential-martingale identity independently by a tilted Feynman–Kac matrix. On one seeded four-state path with 221,183 jumps and minimum state occupation 20,901.343, the actual generator error is $0.027272$. The 4,096-vertex maximum gives radius $0.113539$, versus $0.147442$ from the analytic rate-cap-free bound and $0.260095$ from the preliminary capped-rate construction.

The actual skew-resolvent error is $0.000755$ against propagated radius $0.677127$; the earlier universal stationary bound gave $1.216804$, and the analytic generator radius now gives $0.946064$ rather than $1.767794$. The fixed-point stationary radius is $0.023883$, versus $0.045214$ from $\|\pi\|_2\le1$. The true energy $0.047000719$ lies in the certified interval $[0,0.799625333]$; the analytic-radius upper endpoint is $1.352930540$.

The example also makes the conservatism transparent: this is a valid finite-path certificate, not a sharp inference method.

The result requires continuous full observation but no rate cap. It is not a theorem for hidden states, discrete sampling, estimated state labels, or misspecified generators, and it does not turn marginal random-probe guarantees into conditional guarantees for prescribed features.

Bounded sign probes: the same finite guarantee, not almost-sure detection

Gaussian coordinates are convenient but not necessary for the finite certificate. Fix any orthonormal basis of the centered space and instead let $X,Y$ have independent Rademacher coordinates, each uniformly distributed on $\{-1,1\}$. For $\Xi_\pm=X^\top JY$, still

$$\mathbb E\Xi_\pm^2=S_2=\|J\|_{\rm HS}^2.$$

The fourth moment has an exact basis-dependent correction. If $q_i=\sum_jJ_{ij}^2$, then direct sign expansion gives

$$\mathbb E\Xi_\pm^4 =3S_2^2+6S_4-12\sum_iq_i^2+4\sum_{i,j}J_{ij}^4 \le 3S_2^2+6S_4\le6S_2^2.$$

The first inequality uses $\sum_{i,j}J_{ij}^4\le\sum_iq_i^2$; the second again uses the paired singular values of a real skew matrix. Therefore the same dimension-free Paley–Zygmund bound, 47-pair one-percent miss certificate, unbiased energy estimator, and relative-variance bound $5/m$ hold for bounded sign probes. This is the bilinear skew-form analogue of Hutchinson's Rademacher trace estimator.

There is a sharper finite-dimensional result. If $d\ge2$ is the dimension of the centered feature space, then $\sum_iq_i=S_2$ and Cauchy–Schwarz gives $\sum_iq_i^2\ge S_2^2/d$. Hence

$$\mathbb E\Xi_\pm^4\le(6-8/d)S_2^2,$$ $$\Pr\!\left(|\Xi_\pm|\ge\theta\|J\|_{\rm HS}\right) \ge\frac{(1-\theta^2)^2}{6-8/d}.$$

For the four-state certificate the centered dimension is $d=3$, so the fourth-moment constant improves from $6$ to $10/3$. At half RMS, 25 independent sign pairs now suffice for a one-percent miss bound, compared with the dimension-free count 47.

The constant $6$ nevertheless cannot be reduced uniformly over dimensions. For $d\ge5$, take normalized sine and cosine vectors $a,b$ on the $d$-cycle and $J_d=(ab^\top-ba^\top)/\sqrt2$. Then $\|J_d\|_{\rm HS}=1$ and direct substitution gives

$$\mathbb E(X^\top J_dY)^4=6-\frac{12}{d}+\frac{6}{d^2}\longrightarrow6.$$

Thus the dimension-free fourth-moment constant is asymptotically sharp, even though every fixed dimension receives a strict improvement. The certificate checks this identity through $d=128$, where the ratio is $5.906616211$.

The almost-sure statement does not survive discretization. For $J=\left(\begin{smallmatrix}0&1\\-1&0\end{smallmatrix}\right)$, $\Xi_\pm=X_1Y_2-X_2Y_1$ vanishes with probability exactly $1/2$ despite $J\ne0$. In the four-state certificate, exhaustive enumeration of all 64 sign pairs gives zero probability $1/4$, half-RMS detection probability $1/2$, and fourth moment $0.005697493374$, exactly matching the formula above; its relative energy variance is $1.579139476$, smaller than the Gaussian value $5$.

Exact enumeration for 100 additional random skew matrices in dimensions two through six verifies the fourth-moment identity with maximum normalized error below $2.2\times10^{-15}$. Thus bounded probes retain the finite detection guarantee but must not be described as almost-sure identifiers. The law depends on the chosen orthonormal basis, while the displayed universal bounds do not.

A sharp quantitative irreversibility bound

The same resolvent gives more than a yes-or-no test. On the centered space write $L=S+A$, where $S^\dagger=S$, $A^\dagger=-A$, and put $B=-S$. Assume $B$ is strictly positive and the energy-normalized skew operator

$$C=B^{-1/2}AB^{-1/2}$$

is bounded. Then $C^\dagger=-C$ and the Green–Kubo resolvent has the exact factorization

$$B^{1/2}(-L)^{-1}B^{1/2}=(I-C)^{-1},$$ $$\operatorname{Sym}(I-C)^{-1}=(I-C^2)^{-1},\qquad \operatorname{Anti}(I-C)^{-1}=C(I-C^2)^{-1}.$$

Let $\eta=\|C\|$ and $u=B^{-1/2}f$. Spectral calculus for the skew-adjoint $C$ gives the sharp comparison

$$\frac1{1+\eta^2}\|u\|^2 \le K^{\rm sym}(f,f)\le\|u\|^2.$$

The upper bound is the reversible Green–Kubo value generated by $S$ alone: a skew perturbation cannot increase the diagonal integrated covariance. For two observables,

$$|K^{\rm anti}(f,g)|\le \rho(\eta) \|B^{-1/2}f\|\,\|B^{-1/2}g\|,\qquad \rho(\eta)=\begin{cases}\eta/(1+\eta^2),&\eta\le1,\\1/2,&\eta\ge1. \end{cases}$$

Thus the energy-normalized directional response is universally at most one half. These are full-space operator bounds, not identification claims for a chosen feature block: its antisymmetric entry may still vanish. The variance-reduction side of this factorization is placed against the nonreversible-sampling analysis of Duncan, Lelièvre and Pavliotis (2016) and the symmetric/antisymmetric Green–Kubo analysis of Pavliotis (2010).

There is a second sharp normalization, and the distinction matters. Let $G=(I-C^2)^{-1}$ be the symmetric normalized Green–Kubo operator. Since $G$ is a function of $C$, the two commute and

$$G^{-1/2}\operatorname{Anti}(I-C)^{-1}G^{-1/2}=C,$$ $$G^{-1/2}(I-C)^{-1}G^{-1/2}=I+C.$$

Consequently the intrinsic sector bounds, normalized by the actual symmetric Green–Kubo form, are

$$|K^{\rm anti}(f,g)|\le\eta \sqrt{K^{\rm sym}(f,f)K^{\rm sym}(g,g)},$$ $$|K(f,g)|\le\sqrt{1+\eta^2} \sqrt{K^{\rm sym}(f,f)K^{\rm sym}(g,g)}.$$

Both constants are sharp in the spectral sense. The universal one-half bound above uses the reversible energy norm $\|B^{-1/2}\cdot\|$; it is not a universal bound after normalization by the reduced, irreversible Green–Kubo variance. This prevents the two natural notions of relative directional response from being conflated.

If features $f_r$ multiply slow operators $A_r$, the first corrector decomposes exactly as

$$\sum_{r,s}K_{rs}A_rA_s =\frac12\sum_{r,s}K^{\rm sym}_{rs}\{A_r,A_s\} +\sum_{r<s}K^{\rm anti}_{rs}[A_r,A_s].$$

Reversal leaves the anticommutator part unchanged and flips the commutator part. A chosen feature family can therefore reveal direction only when its antisymmetric Green–Kubo block is nonzero and the corresponding slow operators do not commute.

The scalar periodic-diffusion specialization

The constant-coefficient example is not special. Let

$$dY_t=b(Y_t)\,dt+\sqrt{2a(Y_t)}\,dW_t\pmod {2\pi},\qquad a>0,$$

have a smooth positive stationary density $p$. Its stationary probability current

$$J=b(y)p(y)-(a(y)p(y))'$$

is constant around the circle, and its generator has the exact decomposition

$$L=p^{-1}\partial_y(ap\,\partial_y)+Jp^{-1}\partial_y=:S+A.$$

In $L^2(p)$, $S$ is self-adjoint and $A$ is skew-adjoint. Consequently the stationary time-reversed diffusion has generator $L^\dagger=S-A$, current $-J$, and drift $b^\dagger=b-2J/p$.

In particular, a regular scalar diffusion on the circle is reversible if and only if $J=0$. This is the periodic specialization of the diffusion time-reversal formula of Haussmann and Pardoux (1986).

The general theorem now gives $K_{L^\dagger}(f,g)=K_L(g,f)$ for every integrable centered pair. If $J\ne0$, the full form cannot be symmetric: under the invertibility assumption above that would make $L$ self-adjoint, contradicting $A\ne0$. Hence some smooth centered pair detects the current, although a particular pair may miss it.

The certificate also uses a genuinely nonconstant example, $p(y)=(1+0.35\cos y)/(2\pi)$, $a(y)=0.8+0.2\sin y$, and $J=0.08$. A flux-form periodic CTMC preserves the self/skew-adjoint decomposition exactly.

For $(\cos y,\sin y+0.25\cos2y)$, its $K_{12}-K_{21}$ converges to a nonzero limit at rate 2.001, while reversing $J$ transposes the full matrix with error below $2.2\times10^{-13}$. The same certificate independently checks the general theorem on a four-state nonreversible chain with a nonuniform invariant law: the reversed Green–Kubo matrix agrees with the transpose to machine precision, the full centered form detects irreversibility, and a scalar feature does not.

On that four-state chain, the normalized resolvent factorization, its symmetric part and its skew part agree with their closed forms within $1.23\times10^{-15}$. Here $\eta=0.535874194$; the symmetric spectrum is $[0.776903497,1]$, attaining the lower bound $1/(1+\eta^2)$, while the skew norm is $0.416322535$, attaining $\eta/(1+\eta^2)$. In the intrinsic symmetric Green–Kubo metric, the measured skew norm is $0.535874194=\eta$ and the full norm is $1.134531248=\sqrt{1+\eta^2}$; both matrix identities hold within $4.5\times10^{-16}$.

The central limit condition sets the order of the correction. A stationary ergodic driver with long memory, whose correlation decays like $t^{-\alpha}$ with $0 < \alpha < 1$, has no finite Green–Kubo integral.

With $\varepsilon$ the fluctuation time scale, the variance of its integrated variance scales like $\varepsilon^\alpha$ in place of $\varepsilon$, and so does the smile correction. The chain below mixes exponentially fast and has the order-$\varepsilon$ correction.

The model

The log price $X_t$ and the variance $v_t$ follow Heston's model with zero correlation. The variance level and the vol-of-vol depend on a hidden regime $y_t$:

$$dX_t = -\tfrac12 v_t\,dt + \sqrt{v_t}\,dW_t, \qquad dv_t = \kappa\,\big(\theta_{y_t} - v_t\big)\,dt + \xi_{y_t}\sqrt{v_t}\,dZ_t .$$

The regime is a three-state chain that cycles mostly one way, $1 \to 2 \to 3 \to 1$. The reverse chain has the same stationary law and the same symmetric Green–Kubo matrix. The variance reverts over $1/\kappa$, half a year, while the regimes last a few weeks.

Parameters: $\kappa = 2$; levels $\theta = (0.02,\ 0.10,\ 0.05)$; vol-of-vol $\xi = (0.2,\ 0.4,\ 1.0)$; $v_0 = 0.04$; $S_0 = 100$; zero interest; one year to expiry. Forward rate 2 and reverse rate 0.1 per year along the cycle, multiplied by a speed factor $m$.

The quantity

European call prices and their Black–Scholes implied volatilities, for a regime drawn from its stationary law. The prices come from the characteristic function of the return,

$$\phi(u) = \mathbb{E}\big[e^{iu(X_T - X_0)} \mid v_0\big],$$

through Lewis's formula:

$$C = S_0 - \frac{\sqrt{S_0K}}{\pi}\int_0^\infty \operatorname{Re}\big[e^{iu\log(S_0/K)}\,\phi(u - \tfrac i2)\big]\,\frac{du}{u^2 + 1/4}.$$

Averaging

When the chain switches infinitely fast, the limit is Heston's model with the stationary averages $\bar\theta$ and $\bar\xi^2$.

The first-order rule

After the Fourier transform in $x$, the switched parameters $\theta$ and $\xi^2$ multiply two operators in the variance:

$$A_\theta = \kappa\,\partial_v, \qquad A_{\xi^2} = \tfrac12\,v\,\partial_{vv} .$$

They do not commute:

$$\big[A_\theta,\,A_{\xi^2}\big] f = \tfrac{\kappa}{2}\big(\partial_v(v f'') - v f'''\big) = \tfrac{\kappa}{2}\,f'' .$$

By the first-order rule, the averaged operator gains

$$\sum_{j,k} K^{\mathrm{sym}}_{jk}\,A_j A_k \;+\; K^{\mathrm{anti}}_{\theta\xi^2}\,\tfrac{\kappa}{2}\,\partial_{vv}, \qquad K^{\mathrm{anti}}_{\theta\xi^2} = \tfrac12\big(K_{\theta\xi^2} - K_{\xi^2\theta}\big).$$

The last term adds a constant diffusion to the variance. Its sign flips when the cycle is reversed, and the symmetric terms do not change.

The intuition is a lag. The variance follows its level with a delay of about $1/\kappa$.

When high vol-of-vol tends to follow a high level, it arrives while the variance is high and shakes it harder. When high vol-of-vol comes first, it acts on a low variance.

Method

For each Fourier frequency the coupled equations for the three regimes are solved on one finite-difference grid in $v$, which gives the numerical solution. The rule is applied on the same grid, as a first-order Duhamel correction to the averaged operator, so both share the discretisation error.

Results

Implied volatilities in percent, with the chain four times faster than the base rates. Each regime then lasts about six weeks.

strikeforward, numericalforward, rulereverse, numericalreverse, ruleaveraged Heston
7023.9823.9623.8823.7824.02
8521.5321.5121.5621.5421.55
10020.6220.6520.7420.8820.62
11521.3221.3021.3621.3721.33
13022.7022.6722.6522.5722.74

The forward cycle has the deeper smile: lower at the money, higher in the wings. The difference is about a tenth of a volatility point. The rule gets its sign and shape, and overstates it at this speed.

As the chain speeds up, the largest implied-volatility error across strikes, in volatility points, is:

SPEED $m$averagedsymmetric part onlyfull ruleforward minus reverse at the money
40.0370.1380.031−0.116
80.0240.0650.009−0.078
160.0140.0310.003−0.047
320.0080.0150.0007−0.026

The full rule's error falls by about four for each doubling, as a second-order remainder should. The effect of the direction falls by about two, as a first-order effect should, and at the fastest speed the rule reproduces it within about 20% at every strike.

Keeping only the symmetric part is worse than no correction. In the forward direction the symmetric and antisymmetric parts of the correction largely cancel, so either one alone misses by more than the average does.

Discrete-cycle certificate: verify_cycle_smile.py.