Paper deep dive

Stationary Bias and Extrapolation in Nonlinear Two-Timescale Stochastic Approximation

A. Ch. Madhusudanarao, Rahul Singh

Year: 2026Venue: arXiv preprintArea: cs.LGType: PreprintEmbeddings: 236

Intelligence

Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 91%

Last extracted: 10/11/2026, 4:37:42 AM

Summary

This paper analyzes the stationary bias in nonlinear two-timescale stochastic approximation (SA) driven by an exogenous finite-state Markov chain. The authors derive a first-order bias expansion that remains uniform as the ratio of slow to fast step sizes vanishes. Key findings include the identification of a mixed bias term proportional to ε²/η (where ε is the slow step size and η is the fast step size), which complicates bias reduction via Richardson-Romberg extrapolation. The study verifies these theoretical results using a nonlinear Markov example and applies them to temporal-difference learning with gradient correction (TDC), demonstrating the impact of step-size paths on finite-run extrapolation performance.

Entities (9)

Stochastic Approximation · algorithm · 98%Stationary Bias · concept · 97%Markov Chain · stochastic-process · 96%Two-Timescale Recursions · method · 95%Lyapunov Equation · mathematical-object · 93%Richardson-Romberg Extrapolation · technique · 92%Temporal-Difference Learning · application · 90%Gradient Correction (TDC) · algorithm · 88%Fast-Manifold Coordinates · coordinate-system · 85%

Relation Signals (5)

Stationary Bias → dependson → Step Size Ratio

confidence 94% · The expansion reveals a mixed contribution ε²/η alongside terms linear in each step size.

Temporal-Difference Learning → uses → Two-Timescale Recursions

confidence 92% · We study fixed-policy value estimation using temporal-difference learning with gradient correction (TDC)... We verify localization for temporal-difference learning...

Richardson-Romberg Extrapolation → requires → Path-Matched Weights

confidence 91% · along power-law step-size paths, the bias exponents need not be integers, so Richardson–Romberg extrapolation requires weights matched to the path.

Covariance Equation → regularizedby → Fast-Manifold Coordinates

confidence 89% · Fast-manifold coordinates keep the associated covariance equation regular in this limit.

Stationary Bias → causedby → Nonlinearity

confidence 87% · The connection to covariance comes from nonlinearity: a Taylor expansion of the averaged update around equilibrium contains quadratic terms...

Cypher Suggestions (0)

No Cypher suggestions yet.

Abstract

Abstract:Constant-step stochastic approximation generally has a nonzero stationary mean error that persists under time averaging. This paper studies that error for nonlinear two-timescale recursions driven by an exogenous finite-state Markov chain. Under stated smoothness assumptions and conditions on the stationary distribution, we derive a first-order bias expansion whose error bound remains uniform as the slow step size becomes much smaller than the fast step size. Fast-manifold coordinates keep the associated covariance equation regular in this limit. For fast step $\eta$ and slow step $\varepsilon$, the expansion reveals a mixed contribution $\varepsilon^2/\eta$ alongside terms linear in each step size. This dependence matters for bias reduction: along power-law step-size paths, the bias exponents need not be integers, so Richardson--Romberg extrapolation requires weights matched to the path. An exactly solvable nonlinear Markov example verifies the coefficients. We verify localization for temporal-difference learning and compare finite-run extrapolation at equal update budgets. For finite runs, we bound the initialization error of tail averages on both timescales under an additional coupling assumption. In the special case of additive independent noise, signed third-moment cancellation yields a sharper remainder.

PDF

Open source PDF →Open local PDF →

Trouble viewing inline? Open PDF directly →

Full Text

232,477 characters extracted from source content.

Expand or collapse full text

Stationary Bias and Extrapolation in Nonlinear Two-Timescale Stochastic Approximation A.Ch. Madhusudanarao Affiliation: Department of Computer Science and Automation Affiliation: Indian Institute of Science, Bengaluru, India Email: madhusudanar@iisc.ac.in Rahul Singh Affiliation: Laboratoire de Recherche de l’EPITA, Paris, France Email: rahulsingh0188@gmail.com Abstract Constant-step stochastic approximation generally has a nonzero stationary mean error that persists under time averaging. This paper studies that error for nonlinear two-timescale recursions driven by an exogenous finite-state Markov chain. Under stated smoothness assumptions and conditions on the stationary distribution, we derive a first-order bias expansion whose error bound remains uniform as the slow step size becomes much smaller than the fast step size. Fast-manifold coordinates keep the associated covariance equation regular in this limit. For fast step η and slow step ε , the expansion reveals a mixed contribution ε2/η ^2/η alongside terms linear in each step size. This dependence matters for bias reduction: along power-law step-size paths, the bias exponents need not be integers, so Richardson–Romberg extrapolation requires weights matched to the path. An exactly solvable nonlinear Markov example verifies the coefficients. We verify localization for temporal-difference learning and compare finite-run extrapolation at equal update budgets. For finite runs, we bound the initialization error of tail averages on both timescales under an additional coupling assumption. In the special case of additive independent noise, signed third-moment cancellation yields a sharper remainder. 1 Introduction Stochastic approximation (SA) studies iterative algorithms that update their variables using noisy observations. Constant-step SA generally approaches an invariant distribution whose mean differs from the equilibrium of the averaged dynamics. This difference is the stationary bias. Averaging the iterates zkz_k to give z¯N=N−1∑k=1Nzk z_N=N^-1 _k=1^Nz_k reduces fluctuations but retains this stationary bias since under stationarity z¯N=z1E z_N=Ez_1. In two-timescale SA, two coupled sets of variables are updated with different step sizes. Write the fast step as η, the slow step as ε , and their ratio as ρ=ε/ηρ= /η. Our goal is to determine how these two rates jointly affect the stationary bias. The connection to covariance comes from nonlinearity: a Taylor expansion of the averaged update around equilibrium contains quadratic terms whose expectations involve products of coordinate fluctuations. Consequently, the leading curvature contribution to the mean bias depends on the stationary covariance. Correlations between the current iterate and the observation can produce an additional bias contribution, which our analysis also accounts for. To determine this covariance, we expand the stationary balance of second moments near equilibrium. At leading order, contraction under the linearized dynamics balances the fluctuations supplied by the noise. This balance takes the form of a Lyapunov equation, a linear matrix equation for the leading covariance coefficient. For temporally correlated observations, its noise covariance includes correlations across time. Solving and controlling this equation is therefore an intermediate step in computing the nonlinear mean bias. The difficulty is obtaining covariance estimates that remain useful as ρ↓0ρ 0. On the fast clock, where k iterations correspond to time t=kηt=kη, the slow drift is scaled by ρ. A scalar slow mode v˙=−aρv v=-aρ v, a>0a>0, for example, decays as e−aρte^-aρ t. Thus the slow spectral gap, the smallest exponential decay rate among the slow linearized modes, tends to zero with ρ. Bounds for the Lyapunov equation accumulate perturbations over this decay time and can introduce a factor 1/ρ1/ρ, as illustrated by ∫0∞e−2aρtt=(2aρ)−1 _0^∞e^-2aρ t\,dt=(2aρ)^-1. However, the noise and remainder terms in the slow covariance equation also carry small powers of ρ. Retaining these factors prevents an artificial loss of uniformity and allows the covariance estimate to support a uniform bias expansion. We study nonlinear recursions driven by an exogenous finite-state Markov chain: the next input depends on the input history only through its current state, and its transition law is unaffected by the algorithm. We obtain a first-order stationary bias expansion uniform in small positive ρ, meaning that the remainder constant does not diverge as ρ↓0ρ 0. Coordinates relative to the fast equilibrium manifold, the local family of equilibria of the fast averaged dynamics as the slow variable varies, regularize the covariance equation by exposing the scales of its fast and slow blocks. Writing mη,ρm_η,ρ for the stationary mean displacement from equilibrium, we obtain mη,ρ=ηb0+εb1+ε2ηb2+O⁡(η3/2+ε3η2).m_η,ρ=η b_0+ b_1+ ^2ηb_2+O\! (η^3/2+ ^3η^2 ). Here b0,b1,b2b_0,b_1,b_2 are coefficient vectors determined by the update rule and noise law, independent of the step sizes. The mixed term matters for Richardson–Romberg extrapolation, which forms a weighted combination of estimates at different step sizes to cancel leading bias terms. Along the power-law path ε=ηp =η^p, p>1p>1, its exponent is 2p−12p-1, alongside 11 and p. We derive weights adapted to these powers and verify nonzero mixed bias in an exactly solvable nonlinear Markov model. We also bound the initialization error of tail averages, formed after discarding an initial segment of the trajectory, under an additional coupling assumption. This assumption controls the difference between a trajectory from a prescribed initial state and a stationary trajectory constructed on the same probability space. Stationary existence and localization, expressed through moment bounds around the equilibrium, are explicit hypotheses; local stability alone does not establish them. We study fixed-policy value estimation using temporal-difference learning with gradient correction (TDC) [20]. We verify stationary localization through uniform block contraction and compare finite-run bias, variance, and mean-square error for ordinary averaging and extrapolation at equal update budgets. Related work. The ordinary differential equation (ODE) framework for SA is developed in [14, 3]. Decreasing-step two-timescale theory includes [2, 18]; finite-time results treat linear recursions [4], including Markov inputs [13] and nonlinear systems [6], with Markov central limit theory and applications to TDC and related gradient-TD methods in [10]. Invariant-distribution and weak-error analyses, which study errors in expectations of observables, appear for constant-step stochastic gradient methods in [5, 7]. Bias expansions and Richardson–Romberg correction for Markov SA are studied by Huo et al. [11], while Allmeier and Gast [1] compute first-order bias using Lyapunov equations. For nonlinear SA, Huo et al. [12] distinguish bias contributions from memory, nonlinearity, and their interaction. The closest two-timescale work, Kwon et al. [15], studies stationary laws, bias, variance, averaging, and extrapolation for linear recursions. Our work addresses nonlinear curvature and derives a covariance expansion uniform as the step-size ratio vanishes. Extrapolation methods are also developed in [8, 19, 17]. Notation. Norms are Euclidean or induced matrix norms; ⊤ denotes transpose. A matrix is Hurwitz if all its eigenvalues have negative real parts. Constants in uniform bounds are independent of 0<η≤η00<η≤ _0 and 0<ρ≤ρ00<ρ≤ _0, with sufficiently small fixed upper bounds. Unless indicated otherwise, expectations are stationary. 2 Model and Assumptions 2.1 Recursion and Coordinates Consider xk+1 x_k+1 =xk+ηH(Yk,xk,θk), =x_k+η H(Y_k,x_k, _k), (1) θk+1 _k+1 =θk+εF(Yk,xk,θk), = _k+ F(Y_k,x_k, _k), (2) where xk∈ℝdxx_k ^d_x, θk∈ℝdθ _k ^d_θ, η>0η>0, and ε=ηρ =ηρ with 0<ρ<10<ρ<1. Thus xkx_k uses the larger step size and θk _k the smaller one. The uniform results below concern sufficiently small positive η and ρ. The exogenous input (Yk)(Y_k) is a Markov chain on a finite state space Y, with transition matrix P, i.e., ℙ⁡(Yk+1=y′∣ℱk)=P⁡(Yk,y′)P(Y_k+1=y _k)=P(Y_k,y ), where ℱkF_k is the sigma-field generated by the inputs and iterates through time k. Its invariant distribution is μ. For a function a on Y, write Pa(y)=∑y′P⁡(y,y′)a(y′),μa=∑yμ⁡(y)a(y).Pa(y)= _y P(y,y )a(y ), μ a= _yμ(y)a(y). Set z=(x⊤,θ⊤)⊤z=(x ,θ ) and G⁡(y,z)=(H(y,x,θ)⊤,F(y,x,θ)⊤)⊤G(y,z)=(H(y,x,θ) ,F(y,x,θ) ) . The averaged update field is G¯(z)=μG(⋅,z)=(H¯(x,θ)F¯(x,θ)), G(z)=μ G(·,z)= H(x,θ) F(x,θ), where the average is taken over the input while holding z fixed. The joint recursion is therefore Dρ=diag⁡(Idx,ρIdθ),zk+1=zk+ηDρG(Yk,zk). gatheredD_ρ=diag(I_d_x,ρ I_d_θ),\\ z_k+1=z_k+η D_ρG(Y_k,z_k). gathered (3) Fix an equilibrium z⋆=((x⋆)⊤,(θ⋆)⊤)⊤z =((x ) ,(θ ) ) of the averaged field, so G¯(z⋆)=0 G(z )=0. For each θ near θ⋆θ , let λ⁡(θ)λ(θ) denote the corresponding local equilibrium of the fast averaged dynamics: H¯(λ⁡(θ),θ)=0,λ⁡(θ⋆)=x⋆. H(λ(θ),θ)=0, λ(θ )=x . The local uniqueness and C2C^2 regularity of this branch are specified in Assumption A1. The reduced slow drift is g⁡(θ)=F¯(λ⁡(θ),θ),g⁡(θ⋆)=0.g(θ)= F(λ(θ),θ), g(θ )=0. It describes the slow update when the fast variable is at the equilibrium associated with the current slow variable. Define the Jacobian blocks, all evaluated at z⋆z , by A A =∂xH¯(z⋆), = _x H(z ), B B =∂θH¯(z⋆), = _θ H(z ), C C =∂xF¯(z⋆), = _x F(z ), D D =∂θF¯(z⋆). = _θ F(z ). The local stability assumption (Assumption A2) makes A invertible. Differentiating the fast equilibrium identity H¯(λ(θ),θ)=0 H(λ(θ),θ)=0 gives ADλ(θ⋆)+B=0ADλ(θ )+B=0, and the chain rule then gives J=DG¯(z⋆)=(ABCD),Λ=Dλ(θ⋆)=−A−1B,S=Dg(θ⋆)=D+CΛ=D−CA−1B. gatheredJ=D G(z )= pmatrixA&B\\ C&D pmatrix,\\ =Dλ(θ )=-A^-1B,\\ S=Dg(θ )=D+C =D-CA^-1B. gathered (4) Thus Λ describes the first-order movement of the fast equilibrium as θ changes, while S is the Jacobian of the reduced slow drift. To separate fast tracking from slow displacement, define locally u^=x−λ⁡(θ),v=θ−θ⋆. u=x-λ(θ), v=θ-θ . (5) Here u u is the exact tracking error: it measures how far the fast variable is from the equilibrium associated with the current value of θ. Replacing that equilibrium by its tangent approximation x⋆+Λvx + v gives the tangent tracking error u=x−x⋆−Λvu=x-x - v. Equivalently, (uv)=T⁡(z−z⋆),T=(Idx−Λ0Idθ). uv=T(z-z ), T= pmatrixI_d_x&- \\ 0&I_d_θ pmatrix. (6) The two tracking errors satisfy the exact identity u u =u−rλ(v), =u-r_λ(v), rλ(v) r_λ(v) =λ⁡(θ⋆+v)−x⋆−Λv. =λ(θ +v)-x - v. Since λ is C2C^2 near θ⋆θ , its second derivative is bounded on a sufficiently small closed ball around θ⋆θ . Taylor’s theorem therefore gives, for some r>0r>0 and M<∞M<∞, ‖rλ(v)‖≤M2‖v‖2(‖v‖≤r),\|r_λ(v)\|≤ M2\|v\|^2 (\|v\|≤ r), where M bounds ‖D2λ‖\|D^2λ\| on that ball, so that the remainder is locally quadratic. The tangent coordinates also simplify the coupling in the equilibrium linearization. Using AΛ+B=0A +B=0, the linearized averaged recursion becomes uk+1−uk u_k+1-u_k =η⁡[(A−ρΛC)uk−ρΛSvk], =η [(A-ρ C)u_k-ρ Sv_k ], vk+1−vk v_k+1-v_k =ηρ(Cuk+Svk). =ηρ(Cu_k+Sv_k). The identity AΛ+B=0A +B=0 cancels the order-one effect of slow displacement on the fast linearized drift. After factoring out the common step size η, the remaining cross coefficients are −ρΛS-ρ S and ρCρ C, both of order ρ. 2.2 Standing Assumptions Write ek=zk−z⋆e_k=z_k-z for the error process. For a stationary draw (Y,z)∼πη,ρ(Y,z) _η,ρ from a joint invariant law, define e=z−z⋆,mη,ρ=e,Γη,ρ=Cov⁡(e).e=z-z , m_η,ρ=Ee, _η,ρ=Cov(e). (7) Thus e is a stationary random vector, while mη,ρm_η,ρ and Γη,ρ _η,ρ are its deterministic mean and covariance. All moments below are taken under πη,ρ _η,ρ. Assumption A1 (Local smoothness). The averaged field G¯ G is C3C^3 near z⋆z . Near θ⋆θ , the fast equilibrium λ⁡(θ)λ(θ) is locally unique and C2C^2, with λ⁡(θ⋆)=x⋆λ(θ )=x . Assumption A2 (Local stability). The matrices A and S in (4) are Hurwitz. Assumption A3 (Mixing Markov input). For some n0≥1n_0≥ 1, β>0β>0, and probability measure ν on Y, Pn0(y,E)≥βν(E)P^n_0(y,E)≥βν(E) for every y∈y∈ Y and every subset E⊆E Y. Assumption A4 (Stationary fourth-moment localization). There exist η0,ρ0,C>0 _0, _0,C>0 such that, for every 0<η≤η00<η≤ _0 and 0<ρ≤ρ00<ρ≤ _0, the joint process (Yk,zk)(Y_k,z_k) admits an invariant law πη,ρ _η,ρ satisfying πη,ρ‖z−z⋆‖4≤Cη2.E_ _η,ρ\|z-z \|^4≤ Cη^2. (8) The constant C is independent of η and ρ. All stationary results refer to this selected family of invariant laws. Stationary existence and this moment bound do not follow from local Hurwitz stability alone. Proposition 4 gives sufficient conditions, and Appendix H.9 verifies localization for the nonlinear Markov example. Assumption A5 (Update regularity). Each G⁡(y,⋅)G(y,·) is globally C2C^2, with supy,z(‖DzG(y,z)‖+‖Dz2G(y,z)‖)≤K. _y,z (\|D_zG(y,z)\|+\|D_z^2G(y,z)\| )≤ K. (9) These conditions allow updates growing linearly in z. Together with local C3C^3 smoothness and (8), they imply the expected Taylor estimates in Lemma 4. The stationary covariance and mean-bias expansions require no separate slow-coordinate moment or cross-covariance assumption. The stronger block and coordinate conditions in Assumption A7 are used only when invoking the supplementary blockwise fluctuation statements. The stationary iterates need not have bounded support. 3 Geometry, Noise, and Covariance 3.1 The Two Relaxation Clocks On the fast timescale, where one iteration corresponds to time η, the averaged ODE is z˙=DρG¯(z) z=D_ρ G(z). Its linearization at z⋆z , expressed in the tangent coordinates w=T⁡(z−z⋆)w=T(z-z ), is w˙=Lρw w=L_ρw, where Lρ=TDρJT−1=(A−ρΛC−ρΛSρCρS).L_ρ=TD_ρJT^-1= pmatrixA-ρ C&-ρ S\\ ρ C&ρ S pmatrix. (10) Thus etLρe^tL_ρ propagates an initial linearized perturbation over time t; its subscripts below indicate the output and input coordinate groups. The following bounds separate the fast and slow decay rates and control their interaction, providing the linear estimates needed for uniform covariance analysis. Proposition 1 (Two-clock semigroup bounds). Under Assumptions A1–A2, there exist C∗,c,η0,ρ0>0C_*,c, _0, _0>0 such that for 0<ρ≤ρ00<ρ≤ _0 and t≥0t≥ 0, ‖[etLρ]uu‖ \|[e^tL_ρ]_u\| ≤C∗(e−ct+ρe−cρt), ≤ C_*(e^-ct+ρ e^-cρ t), (11) ‖[etLρ]uv‖+‖[etLρ]vu‖ \|[e^tL_ρ]_uv\|+\|[e^tL_ρ]_vu\| ≤C∗ρ(e−ct+e−cρt), ≤ C_*ρ(e^-ct+e^-cρ t), ‖[etLρ]vv‖ \|[e^tL_ρ]_v\| ≤C∗e−cρt. ≤ C_*e^-cρ t. For 0<η≤η00<η≤ _0 and integers k≥0k≥ 0, the same inequalities hold with etLρe^tL_ρ replaced by (I+ηLρ)k(I+η L_ρ)^k, and with e−cte^-ct and e−cρte^-cρ t replaced by e−cηke^-cη k and e−cεke^-c k, respectively, where ε=ηρ =ηρ. Here (I+ηLρ)k(I+η L_ρ)^k propagates the linearized averaged recursion wk+1=(I+ηLρ)wkw_k+1=(I+η L_ρ)w_k through k iterations, corresponding to fast time t=ηkt=η k and slow time ρt=εkρ t= k. The proof of both bounds is given in Appendix B. 3.2 Poisson Representation We identify the effective covariance of the Markov noise, accounting for correlations across time. This supplies the noise term in the Lyapunov equation (14), whose solution determines the leading-order stationary covariance. The latter is used below to quantify how nonlinear drift converts fluctuations into stationary bias. To obtain the noise covariance, we use the Poisson equation introduced next. For fixed z, let G~(y,z)=G⁡(y,z)−G¯(z) G(y,z)=G(y,z)- G(z). We consider the Poisson equation for the Markov transition operator P, with a zero-mean normalization: −P=G~,μ=0.U-PU= G, =0. (12) Under Assumption A3, geometric mixing of the Markov input YkY_k gives the solution ⁡(⋅,z)=∑j≥0PjG~(⋅,z)U(·,z)= _j≥ 0P^j G(·,z) [9, Section 5]. Set ξk=G⁡(Yk,z⋆) _k=G(Y_k,z ) and ζk+1(z)=⁡(Yk+1,z)−(P)(Yk,z) _k+1(z)=U(Y_k+1,z)-(PU)(Y_k,z). Lemma 1 (Long-run covariance). Under Assumption A3, ζk+1(z) _k+1(z) is a martingale difference for fixed z, meaning ⁡[ζk+1(z)∣ℱk]=0E[ _k+1(z) _k]=0. For the stationary input chain Y0∼μY_0 μ, write Kh=⁡[ξ0ξh⊤]K_h=E[ _0 _h ]. Then Q Q =limm→∞m−1[(∑k=0m−1ξk)(∑k=0m−1ξk)⊤] = _m→∞m^-1E [ ( _k=0^m-1 _k ) ( _k=0^m-1 _k ) ] (13) =K0+∑h≥1(Kh+Kh⊤)=⁡[ζ1(z⋆)ζ1(z⋆)⊤]. =K_0+ _h≥ 1(K_h+K_h )=E[ _1(z ) _1(z ) ]. The sum in (13) accumulates the equilibrium noise over m steps. To compute its covariance per step, we use the Poisson equation to write ∑k=0m−1ξk _k=0^m-1 _k =∑k=0m−1ζk+1(z⋆)+⁡(Y0,z⋆)−⁡(Ym,z⋆). = _k=0^m-1 _k+1(z )+U(Y_0,z )-U(Y_m,z ). The first term is a martingale sum. The remaining difference depends only on the initial and final noise states and is bounded uniformly in m, since Y is finite. Its covariance and cross terms vanish after division by m, yielding the last equality in (13); see Appendix C.2. 3.3 Covariance of Stationary Iterates We now approximate Γη,ρ=Cov⁡(z−z⋆) _η,ρ=Cov(z-z ) under the stationary law πη,ρ _η,ρ. This covariance measures fluctuations around the stationary mean and enters the bias expansion in Section 4. Write Jρ=DρJ_ρ=D_ρJ, which is Hurwitz for sufficiently small ρ by Proposition 1. Define Σρ _ρ as the unique solution of the Lyapunov equation JρΣρ+ΣρJρ⊤+DρQDρ=0.J_ρ _ρ+ _ρJ_ρ +D_ρQD_ρ=0. (14) Appendix Lemma 5 derives this equation from the linearized recursion. The next proposition shows that ηΣρη _ρ approximates Γη,ρ _η,ρ, the covariance of the joint iterates zk=(xk⊤,θk⊤)⊤z_k=(x_k , _k ) generated by (3) under the stationary law πη,ρ _η,ρ. Proposition 2 (Uniform covariance approximation). Under Assumptions A1–A5, for all sufficiently small positive η,ρη,ρ, Γη,ρ=ηΣρ+Rη,ρΓ,‖Rη,ρΓ‖≤Cη3/2. _η,ρ=η _ρ+R _η,ρ, \|R _η,ρ\|≤ Cη^3/2. (15) The constant C is independent of η and ρ, so the bound remains valid as ρ↓0ρ 0. Proof idea. The complete proof of Proposition 2 is given in Appendix D.2. Apply the Poisson correction (90) and then the tangent transformation (6); the leading noise then has zero conditional mean. Stationarity then gives a covariance balance, with remainder terms controlled by the moment and regularity assumptions. Crucially, the remainder in the slow covariance equation carries a factor ρ, offsetting its O⁡(1/ρ)O(1/ρ) relaxation time. Solving this balance and removing the correction gives the stated bound; see Appendix D.2. ∎ 4 Stationary Mean Bias We derive the stationary mean bias from the balance G(Yk,zk)=0EG(Y_k,z_k)=0. Theorem 1 combines the covariance approximation of Proposition 2, a second-order Taylor expansion of G¯(z) G(z) about z⋆z , and a Poisson-corrector identity to obtain mη,ρ=ηb(ρ)+O⁡(η3/2)m_η,ρ=η b(ρ)+O(η^3/2). The coefficient separates curvature and direct Markov contributions. Lemma 2 supplies the auxiliary bound ‖mη,ρ‖=O⁡(η)\|m_η,ρ\|=O(η) used in the theorem’s proof. Both bounds are uniform as ρ↓0ρ 0. Define the curvature map and direct Markov response by (M)i=12Tr⁡(D2G¯i(z⋆)M),C(M)_i= 12Tr(D^2 G_i(z )M), rρ=∑y,y′μ⁡(y)P(y,y′)Dz(y′,z⋆)DρG(y,z⋆).r_ρ= _y,y μ(y)P(y,y )D_zU(y ,z )D_ρG(y,z ). (16) Theorem 1 (Fixed-ratio bias). If Assumptions A1–A5 hold, mη,ρ=ηb(ρ)+O⁡(η3/2),b⁡(ρ)=−J−1⁡(Σρ)+rρ, gatheredm_η,ρ=η b(ρ)+O(η^3/2),\\ b(ρ)=-J^-1\C( _ρ)+r_ρ\, gathered (17) with a remainder uniform for 0<ρ≤ρ00<ρ≤ _0. Lemma 2 (A priori mean bound). Under the same assumptions, ‖mη,ρ‖≤Cη\|m_η,ρ\|≤ Cη uniformly. The proof is given in Appendix D.3; Appendix D.4 supplies the detailed identities and remainder bounds. 5 Bias Expansion as the Step-Size Ratio Vanishes. Section 4 leaves the ratio dependence implicit in b⁡(ρ)b(ρ). As ρ↓0ρ 0, the slow relaxation rate vanishes on the fast clock, making the covariance equation singular at ρ=0ρ=0. We rescale the covariance blocks to obtain regular limiting equations and expand b⁡(ρ)b(ρ). Theorem 2 identifies contributions proportional to η, ε , and ε2/η ^2/η, with an O⁡(η3/2+ε3/η2)O(η^3/2+ ^3/η^2) remainder whose implied constant is independent of η and ρ. This expansion also guides bias reduction using runs at different step sizes. For example, at fixed ratio, halving both steps halves the leading bias, so twice the finer-step average minus the coarser-step average cancels that leading stationary bias. Section 6 adapts this cancellation to runs in which the step-size ratio also changes. Theorem 2 (Singular expansion). Under Assumptions A1–A5, uniformly for small positive η,ρη,ρ, mη,ρ m_η,ρ =ηb0+εb1+ε2ηb2+O⁡(η3/2+ε3η2). =η b_0+ b_1+ ^2ηb_2+O\! (η^3/2+ ^3η^2 ). (18) Proof sketch. In tangent coordinates, write TΣρT⊤=(U⁡(ρ)ρR(ρ)ρR(ρ)⊤ρV(ρ)).T _ρT = pmatrixU(ρ)&ρ R(ρ)\\ ρ R(ρ) &ρ V(ρ) pmatrix. Rescaling the cross and slow covariance equations gives a linear system with polynomial coefficients in ρ. At ρ=0ρ=0, it is triangular and invertible by Hurwitz stability of A and S. Its solution is therefore analytic near zero. Expanding Σρ _ρ through order ρ2ρ^2, using the affine dependence of rρr_ρ on ρ, and substituting into Theorem 1 gives the result. The detailed proof and coefficient formulas are in Appendix E.1. ∎ In tangent coordinates, the leading fast covariance, fast–slow cross-covariance, and slow covariance are ηU(ρ)=O⁡(η)η U(ρ)=O(η), ηρR(ρ)=O⁡(ε)ηρ R(ρ)=O( ), and ηρV(ρ)=O⁡(ε)ηρ V(ρ)=O( ), respectively. The absolute O⁡(η3/2)O(η^3/2) covariance remainder can obscure the O⁡(ε)O( ) terms when ρ≲ηρ η. The sharper O⁡(ε)O( ) bounds for the actual slow covariance and fast–slow cross-covariance are additional conditions collected in Assumption A7. They are not needed for Proposition 2 or Theorems 1–2. Proposition 4 gives sufficient conditions under which these stronger bounds hold. 6 Extrapolation Along Step-Size Paths Time averaging reduces random fluctuations but does not remove stationary bias. We reduce this bias by combining averages from several runs at different pairs of fast and slow step sizes. For each run, we choose a fast step η and set the slow step to ε=ηp =η^p, with the same fixed p>1p>1 across runs. Thus, reducing η also determines how ε decreases. This relation between the two step sizes defines the step-size path. The expansion in Section 5 describes how the bias changes along this path, allowing us to choose weights so that the leading stationary bias terms cancel. Proposition 3 gives a three-run combination that cancels the contributions proportional to the fast and slow steps. We also describe how to cancel further terms and bound the remaining bias. Before choosing extrapolation weights, we identify the bias powers along the chosen step-size path and compare them with the remainder. For ε=ηp =η^p, p>1p>1, Theorem 2 gives m⁡(η)=ηb0+ηpb1+η2p−1b2+O⁡(η3/2+η3p−2).m(η)=η b_0+η^pb_1+η^2p-1b_2+O(η^3/2+η^3p-2). (19) A nonzero bias term can be distinguished from the remainder when it decays more slowly as η↓0η 0. Since 3p−2>2p−1>p>13p-2>2p-1>p>1, comparison with the O⁡(η3/2)O(η^3/2) remainder determines which displayed terms satisfy this condition: ηpb1η^pb_1 does so for 1<p<3/21<p<3/2, and η2p−1b2η^2p-1b_2 for 1<p<5/41<p<5/4. Under the additional additive-i.i.d. conditions of Theorem 4, the sharper O⁡(η2)O(η^2) remainder extends the latter range to 1<p<3/21<p<3/2. At fixed ρ, let Mρ(h)M_ρ(h) denote the stationary mean at steps (η,ε)=(h,ρh)(η, )=(h,ρ h). Halving both steps halves the leading bias, so twice the finer-step mean minus the coarser-step mean cancels this contribution: 2Mρ(h/2)−Mρ(h)=z⋆+O⁡(h3/2).2M_ρ(h/2)-M_ρ(h)=z +O(h^3/2). (20) Along ε=ηp =η^p, reducing the fast step must preserve this relation: the corresponding pairs are (ηj,εj)=(2−jh,(2−jh)p)( _j, _j)=(2^-jh,(2^-jh)^p). The bias terms proportional to h and hph^p now decrease by different factors between successive runs. Proposition 3 gives weights for combining the three stationary means at j=0,1,2j=0,1,2 so that both terms cancel, and bounds the remaining bias. Proposition 3 (Three-level path extrapolation). Let M⁡(h)M(h) be the stationary mean at (h,hp)(h,h^p) under Theorem 2. Set q1=2−1q_1=2^-1, qp=2−pq_p=2^-p, and (w0,w1,w2)=(q1qp,−q1−qp,1)(1−q1)(1−qp).(w_0,w_1,w_2)= (q_1q_p,-q_1-q_p,1)(1-q_1)(1-q_p). (21) Then ∑j=02wjM(2−jh)=z⋆+O⁡(hmin⁡2p−1,3/2). _j=0^2w_jM(2^-jh)=z +O\! (h \2p-1,3/2\ ). If the stationary weak remainder is O⁡(h2)O(h^2), replace 3/23/2 by 22. Proof. The weights satisfy ∑jwj=1 _jw_j=1, ∑jwj2−j=0 _jw_j2^-j=0, and ∑jwj2−jp=0 _jw_j2^-jp=0. Apply these identities to (19). ∎ Appendix E gives weights for canceling additional bias terms and explains the limits imposed by the available remainder bounds. Extrapolation can increase sampling variance, so reducing bias does not necessarily improve accuracy at fixed computational cost. An exactly solvable nonlinear Markov example in Appendix H verifies that the mixed bias term can be nonzero and illustrates these extrapolation rules using exact stationary means. 7 Finite-Time Consequences A stationary coupling decomposes the error into initialization, centered fluctuation, and stationary mean. Under Assumptions A1–A5 and the additional block localization in Assumption A7, ‖u‖2=O⁡(η),‖v‖2=O⁡(ε+η2).E\|u\|^2=O(η), \|v\|^2=O( +η^2). (22) The η2η^2 term is squared bias and can dominate the slow variance on thin paths. Assumption A6 (Two-clock coupling). A coupling with a stationary version satisfies ‖zk−zkst‖≤C1+V⁡(z0)(e−cηk+e−cεk),E\|z_k-z_k st\|≤ C\1+V(z_0)\(e^-cη k+e^-c k), (23) where V≥0V≥ 0 has uniformly bounded invariant expectation. Any meeting-time contribution from different input initializations is included in this bound. This additional condition does not follow from local linear stability. Theorem 3 (Tail-average bias and initialization). For z¯K,N=L−1∑k=KN−1zk z_K,N=L^-1 _k=K^N-1z_k, L=N−K>0L=N-K>0, under Assumption A6, ∥K,N∥≤C1+V(z0)[ \| T_K,N\|≤ C\1+V(z_0)\ [ e−cηKmin1,(ηL)−1+e−cεKmin1,(εL)−1], e^-cη K \1,(η L)^-1\+e^-c K \1,( L)^-1\ ], (24) where K,N=z¯K,N−zst T_K,N=E z_K,N-Ez st. If Assumptions A1–A5 also hold, then z¯K,N−z⋆=ηb(ρ)+O⁡(η3/2)+K,N.E z_K,N-z =η b(ρ)+O(η^3/2)+ T_K,N. (25) Proof. Average (23) and use, for h=η,εh=η, , L−1∑k=KN−1e−chk L^-1 _k=K^N-1e^-chk =e−chK1−e−chL⁡(1−e−ch) =e^-chK 1-e^-chLL(1-e^-ch) ≤Ce−chKmin⁡1,(hL)−1. ≤ Ce^-chK \1,(hL)^-1\. Then apply Theorem 1. ∎ Stationary initialization makes K,N=0 T_K,N=0. For finite-run extrapolation, weighted initialization errors must lie below the retained remainder. Appendix I derives conditional blockwise mean-square bounds under additional stationary moment and coordinate L2L^2 coupling assumptions. 8 Application: temporal-difference learning We study fixed-policy value estimation using temporal-difference learning with gradient correction (TDC) [20, 16]. We compare stationary-bias predictions with simulations and evaluate extrapolation at equal update budgets. 8.1 Setup and Recursion Let SkS_k be the state at time k under a fixed policy, Rk+1R_k+1 the reward on the transition to Sk+1S_k+1, and γ∈(0,1)γ∈(0,1) the discount factor. For a fixed scalar feature ϕ⁡(s)φ(s), we approximate the discounted value function by Vθ(s)=ϕ⁡(s)θV_θ(s)=φ(s)θ, with parameter θ∈ℝθ . Write ϕk=ϕ⁡(Sk) _k=φ(S_k) and define δk=Rk+1+γϕk+1θk−ϕkθk. _k=R_k+1+γ _k+1 _k- _k _k. TDC uses a fast auxiliary variable wkw_k and a slow value parameter θk _k [20]: wk+1 w_k+1 =wk+η⁡(δk−ϕkwk)ϕk, =w_k+η( _k- _kw_k) _k, θk+1 _k+1 =θk+ε⁡[δkϕk−γϕk+1ϕkwk]. = _k+ [ _k _k-γ _k+1 _kw_k ]. Both updates use the current parameters and the same observed transition. We use a Markov reward process with state space =0,1,2,3,4 S=\0,1,2,3,4\, representing positions along a five-site chain. From an interior state s∈1,2,3s∈\1,2,3\, the process stays at s with probability 1/21/2 and moves to each neighboring state with probability 1/41/4. At either boundary, an attempted outward move becomes a self-loop: the process stays with probability 3/43/4 and moves to its only neighbor with probability 1/41/4. The reward is Rk+1=Sk+1=4R_k+1=1\S_k+1=4\, so a transition earns reward one exactly when its destination is the rightmost state. The process continues after receiving this reward; state 44 is not terminal. We set ϕ⁡(s)=1+0.05sφ(s)=1+0.05s and γ=0.2γ=0.2. Our estimation target is the TD equilibrium θ⋆≈0.2443164284θ ≈ 0.2443164284, derived in Appendix J.1. The updates are affine, so their curvature contribution vanishes; this example tests the direct Markov contribution to stationary bias. Figure 1: Empirical bias and variability of TDC estimates of θ⋆θ . Top: empirical bias divided by η at fixed ratios ρ=0.1ρ=0.1 and ρ=0.25ρ=0.25, for Markov sampling and independent transition-pair sampling. Dashed lines show the theoretical leading coefficient b0θ+ρb1θb_0^θ+ρ b_1^θ; the dotted line marks the zero stationary bias of the independent-pair control. Bottom: empirical signed bias and standard deviation at equal total update budgets, including burn-in, along ε=ηp =η^p with p=6/5p=6/5 and base step h=0.05h=0.05. All estimates use 256 independent replications. Error bars are pointwise 95% intervals: Student-t intervals for bias and replication-bootstrap intervals for standard deviation. Solid segments connect empirical points. R labels specify the cancelled powers. 8.2 Empirical Bias Estimation and Extrapolation Our experiments have two goals: to compare the bias estimated from simulated trajectories with the theoretical bias expansion, and to assess whether extrapolation improves estimation accuracy at equal total update budgets. Estimating bias from simulated trajectories. Fix the step sizes η and ε=ηρ =ηρ. We generate R=256R=256 independent replications of the TDC recursion. In replication r, we discard the first K updates and average the next L slow iterates: θ^(r)=1L∑ℓ=1LθK+ℓ(r). θ^(r)= 1L _ =1^Lθ^(r)_K+ . The empirical bias is B^η,ρθ=1R∑r=1R(θ^(r)−θ⋆). B^θ_η,ρ= 1R _r=1^R ( θ^(r)-θ ). (26) The target θ⋆θ is the equilibrium of the averaged TD problem, given in (224). Thus, the empirical bias is calculated directly from simulated iterates and the known target. For finite K and L, this estimates the bias of the post-burn average. We interpret it as an estimate of stationary bias when the remaining initialization effect is negligible. Comparison with the theoretical prediction. Theorem 1, specialized to this affine TDC recursion in Appendix J.1, gives, for sufficiently small positive η and ρ, Bη,ρθ:=πη,ρ[θ]−θ⋆=η⁡(b0θ+ρb1θ)+O⁡(η3/2).B^θ_η,ρ:=E_ _η,ρ[θ]-θ =η (b_0^θ+ρ b_1^θ )+O(η^3/2). Consequently, at a fixed ratio ρ, the normalized empirical bias B^η,ρθ/η B^θ_η,ρ/η can be compared with the leading coefficient b0θ+ρb1θb_0^θ+ρ b_1^θ. The top panels of Figure 1 show this comparison at ρ∈0.1,0.25ρ∈\0.1,0.25\. The points and their confidence intervals are obtained from simulations; the dashed horizontal lines mark b0θ+ρb1θb_0^θ+ρ b_1^θ, the predicted limit of the normalized stationary bias Bη,ρθ/ηB^θ_η,ρ/η as η→0η→ 0 at fixed ρ. As a control, we independently sample transition pairs from the same stationary transition distribution. This preserves the averaged TD problem and its target θ⋆θ while removing serial dependence. The stationary bias of this independent-pair recursion is zero. Comparing its empirical bias with that of the Markov recursion therefore examines the effect of temporal dependence. Extrapolation at equal update budgets. We compare ordinary averaging with extrapolated estimates along the step-size path ηj=h2−j,εj=ηjp,h=0.05,p=6/5. _j=h2^-j, _j= _j^p, h=0.05, p=6/5. Within each replication, let θ^j θ_j denote the post-burn average at level j. Ordinary averaging uses level 00. The two-level rule R 1\1\ uses θ^R1=2θ^1−θ^0. θ_R\1\=2 θ_1- θ_0. The three-level rules R 1,2\1,2\ and R 1,p\1,p\ use weights that sum to one and cancel the indicated powers of h. We use total update budgets B∈217,219,221B∈\2^17,2^19,2^21\. A method with m levels allocates N=⌊B/m⌋N= B/m updates to each level, including burn-in. Within a replication, the levels share the same input sequence. Each replication produces one estimate for each method. We compute empirical bias, standard deviation, and RMSE across these independent replication-level estimates. The bottom panels of Figure 1 show empirical bias and standard deviation. RMSE measures their combined contribution to estimation error. Appendix J.1 verifies the assumptions and derives the theoretical bias coefficients. Appendix J.2 specifies the simulation lengths, burn-in, and uncertainty calculations. The nonlinear mixed bias term is illustrated separately in Appendix H. 9 Conclusion For nonlinear two-timescale stochastic approximation, we derive a stationary-bias expansion with terms ηb0+εb1+(ε2/η)b2η b_0+ b_1+( ^2/η)b_2 and uniform remainder control as the step-size ratio vanishes. Along paths ε=ηp =η^p, these contributions have powers 11, p, and 2p−12p-1, motivating Richardson–Romberg extrapolation rules adapted to the chosen path. The TDC experiments illustrate the practical trade-off: cancelling additional bias terms can increase sampling variability, so higher-order cancellation need not improve accuracy at equal update budgets. Sharper Markov remainder bounds and optimal allocation across extrapolation levels remain open questions. AI Use Statement ChatGPT/Codex assisted with paper review, and simulation code. References [1] S. Allmeier and N. Gast (2024) Computing the bias of constant-step stochastic approximation with Markovian noise. Advances in Neural Information Processing Systems 37, p. 137873–137902. Cited by: §1. [2] V. S. Borkar (1997) Stochastic approximation with two time scales. Systems & Control Letters 29 (5), p. 291–294. Cited by: Appendix K, §1. [3] V. S. Borkar (2008) Stochastic approximation: a dynamical systems viewpoint. Texts and Readings in Mathematics, Vol. 48, Hindustan Book Agency. External Links: Document Cited by: §1. [4] G. Dalal, G. Thoppe, B. Szőrényi, and S. Mannor (2018) Finite sample analysis of two-timescale stochastic approximation with applications to reinforcement learning. In Proceedings of the 31st Conference on Learning Theory, Proceedings of Machine Learning Research, Vol. 75, p. 1199–1233. External Links: Link Cited by: Appendix K, §1. [5] A. Dieuleveut, A. Durmus, and F. Bach (2020) Bridging the gap between constant step size SGD and Markov chains. The Annals of Statistics 48 (3), p. 1348–1382. External Links: Document Cited by: §1. [6] T. T. Doan (2023) Nonlinear two-time-scale stochastic approximation: convergence and finite-time performance. IEEE Transactions on Automatic Control 68 (8), p. 4695–4705. External Links: Document Cited by: Appendix K, §1. [7] Y. Feng, T. Gao, L. Li, J. Liu, and Y. Lu (2020) Uniform-in-time weak error analysis for stochastic gradient descent algorithms via diffusion approximation. Communications in Mathematical Sciences 18 (1), p. 163–188. External Links: Document Cited by: §1. [8] N. Frikha and L. Huang (2015) A multi-step richardson–romberg extrapolation method for stochastic approximation. Stochastic Processes and their Applications 125 (11), p. 4066–4101. Cited by: §1. [9] P. W. Glynn and A. Infanger (2024) Solution representations for Poisson’s equation, martingale structure, and the Markov chain central limit theorem. Stochastic Systems 14 (1), p. 47–68. External Links: Document Cited by: §3.2. [10] J. Hu, V. Doshi, and D. Y. Eun (2024) Central limit theorem for two-timescale stochastic approximation with Markovian noise: theory and applications. In Proceedings of the 27th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 238, p. 1477–1485. External Links: Link Cited by: Appendix K, §1. [11] D. Huo, Y. Chen, and Q. Xie (2026) Bias and extrapolation in Markovian linear stochastic approximation with constant step sizes. Mathematics of Operations Research. External Links: Document Cited by: §1. [12] D. L. Huo, Y. Zhang, Y. Chen, and Q. Xie (2024) The collusion of memory and nonlinearity in stochastic approximation with constant stepsize. Advances in Neural Information Processing Systems 37, p. 21699–21762. Cited by: §1. [13] M. Kaledin, E. Moulines, A. Naumov, V. Tadic, and H. Wai (2020) Finite time analysis of linear two-timescale stochastic approximation with Markovian noise. In Conference on Learning Theory, p. 2144–2203. Cited by: Appendix K, §1. [14] H. J. Kushner and G. G. Yin (2003) Stochastic approximation and recursive algorithms and applications. Springer. Cited by: §1. [15] J. Kwon, L. Dotson, Y. Chen, and Q. Xie (2025) Two-timescale linear stochastic approximation: constant stepsizes go a long way. In Proceedings of the 28th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 258, p. 3781–3789. External Links: Link Cited by: Appendix K, §1. [16] V. Mahadevan, C. Chen, S. D. Liu, and S. Zhang (2026) Convergence of two-timescale Markovian stochastic approximations with applications in reinforcement learning. In Proceedings of the 43rd International Conference on Machine Learning, Proceedings of Machine Learning Research, Vol. 306, p. 85598–85668. External Links: Link Cited by: §8. [17] P. Mangold, A. O. Durmus, A. Dieuleveut, S. Samsonov, and E. Moulines (2025) Refined analysis of constant step size federated averaging and federated Richardson–Romberg extrapolation. In Proceedings of the 28th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research, Vol. 258, p. 5023–5031. External Links: Link Cited by: §1. [18] A. Mokkadem and M. Pelletier (2006) Convergence rate and averaging of nonlinear two-time-scale stochastic approximation algorithms. The Annals of Applied Probability 16 (3), p. 1671–1702. External Links: Document Cited by: Appendix K, §1. [19] M. Sheshukova, D. Belomestny, A. Durmus, E. Moulines, A. Naumov, and S. Samsonov (2025) Nonasymptotic analysis of stochastic gradient descent with the Richardson–Romberg extrapolation. In The Thirteenth International Conference on Learning Representations, External Links: Link Cited by: §1. [20] R. S. Sutton, H. R. Maei, D. Precup, S. Bhatnagar, D. Silver, C. Szepesvári, and E. Wiewiora (2009) Fast gradient-descent methods for temporal-difference learning with linear function approximation. In Proceedings of the 26th Annual International Conference on Machine Learning, p. 993–1000. External Links: Document, Link Cited by: §1, §8.1, §8. Supplementary Derivations for Stationary Bias and Extrapolation in Nonlinear Two-Timescale Stochastic Approximation Contents 1 Introduction 2 Model and Assumptions 2.1 Recursion and Coordinates 2.2 Standing Assumptions 3 Geometry, Noise, and Covariance 3.1 The Two Relaxation Clocks 3.2 Poisson Representation 3.3 Covariance of Stationary Iterates 4 Stationary Mean Bias 5 Bias Expansion as the Step-Size Ratio Vanishes. 6 Extrapolation Along Step-Size Paths 7 Finite-Time Consequences 8 Application: temporal-difference learning 8.1 Setup and Recursion 8.2 Empirical Bias Estimation and Extrapolation 9 Conclusion References A Notation, assumptions, and proof roadmap A.1 The Algorithm and Its Two Clocks A.2 Averaged Dynamics and the Fast Equilibrium Manifold A.3 Coordinates Relative to the Fast Equilibrium Manifold A.4 Compact Notation A.5 Standing Assumptions A.6 Quantities of Interest and Step-Size Regimes B ODE geometry of two timescales B.1 Exact Fast-Manifold Dynamics B.2 Linearization at the Equilibrium B.3 Two-Clock Semigroup Bounds C Poisson tools and Taylor estimates C.1 Long-Run Covariance of the Markov Noise C.2 Martingale Decomposition via the Poisson Equation C.3 Taylor Estimates for the Averaged Drift and Poisson Corrector D Uniform stationary covariance and mean bias D.1 Linear Covariance Balance and Its Small-Step Limit D.2 Proof of Proposition : Uniform Covariance Approximation D.3 Proof of Lemma and Theorem D.4 Stationary Mean Bias E The singular ratio limit and path-aware extrapolation E.1 Proof of Theorem E.2 Overview of the Coefficient Expansion E.3 Rescaled Covariance Equations E.4 Derivation of the Covariance Expansion E.5 Bias Coefficients E.6 Visibility along ε=ηp =η^p E.7 Richardson–Romberg Extrapolation along Step-Size Paths F Additive Independent Noise G Proof of localization under uniform block contraction H A nonlinear Markov example with an exact bias formula H.1 Model and Markov Input H.2 Stationary Representations H.3 Exact Second Moments H.4 Stationary Mean H.5 Fixed-Ratio Bias H.6 Singular Expansion and Coefficient Recovery H.7 Extrapolation along Power-law Paths H.8 Agreement with the General Bias Formula H.9 Stationary Localization and Coupling Verification I From the stationary calculation to finite time I.1 Stationary Block Scales I.2 Initialization and Blockwise Mean-Square Bounds J Application details and additional simulations J.1 Application: Linear Temporal-Difference Learning with Correction J.2 Simulation Protocol J.3 Results K Additional related work Appendix A Notation, assumptions, and proof roadmap We first recall the notation and the assumptions from the main paper. The general proofs come next, followed by the additive-noise refinement, a sufficient localization condition, the worked example, finite-time consequences, and the temporal-difference learning application and additional experiments. Location Purpose Appendix B Fast-manifold geometry and the two-clock semigroup bounds. Appendix C Poisson representation, long-run noise covariance, and Taylor estimates. Appendix D Proofs of the uniform covariance approximation and stationary mean bias. Appendix E Singular covariance coefficients and extrapolation along step-size paths. Appendix F The sharper remainder under additive independent noise. Appendix G A sufficient condition for stationary localization. Appendix H Exact formulas and localization and coupling verification for the Markov example. Appendix I Additional block localization and conditional mean-square bounds. Appendix J Temporal-difference learning model, assumption verification, and Monte Carlo protocol for TDC and the exact nonlinear example. Appendix K gives additional related work. Constants in uniform estimates are independent of sufficiently small positive η and ρ, as in the main paper. A.1 The Algorithm and Its Two Clocks The fast variable is xk∈ℝdxx_k ^d_x and the slow variable is θk∈ℝdθ _k ^d_θ. They evolve according to xk+1 x_k+1 =xk+ηH(Yk,xk,θk), =x_k+η H(Y_k,x_k, _k), (27) θk+1 _k+1 =θk+εF(Yk,xk,θk). = _k+ F(Y_k,x_k, _k). (28) Here η>0η>0 is the fast step size, ε>0 >0 is the slow step size, H and F are the two update fields, and YkY_k is the random input. The state space is Euclidean; no projection is applied after either update. A stationary analysis therefore needs a global stability condition (). The random input is an exogenous Markov chain on a finite state space Y with P as its transition matrix. We assume that this chain has a unique stationary distribution μ and mixes geometrically. The ratio ρ=εηρ= η (29) measures how much the slow variable moves during one unit of fast algorithmic time. We work with 0<ρ≤ρ00<ρ≤ _0, where ρ0 _0 is a sufficiently small fixed constant. Equivalently, ε=ηρ =ηρ. The limit η→0η→ 0 controls the size of the stationary fluctuations; the limit ρ→0ρ→ 0 separates the two clocks: fast time tk=kη,slow time sk=kε=ρtk. fast time t_k=kη, time s_k=k =ρ t_k. (30) We distinguish two limiting regimes. In the fixed-ratio regime, η→0η→ 0 while ρ=ε/ηρ= /η remains fixed and positive, so ε=ρη→0 =ρη→ 0. In the singular-ratio regime, both η→0η→ 0 and ρ→0ρ→ 0, so ε=o⁡(η) =o(η) and the two clocks become increasingly separated. A result established separately for each fixed ρ need not hold uniformly in the singular-ratio regime, because its bounds may deteriorate as ρ→0ρ→ 0. The uniform remainder bounds proved below ensure that our stationary expansion also applies in the singular-ratio regime. A.2 Averaged Dynamics and the Fast Equilibrium Manifold This subsection introduces the averaged dynamics and the fast equilibrium manifold, and gives a high-level interpretation of their fast–slow structure. The qualitative discussion of relaxation and tracking motivates the coordinates introduced in the next subsection. The local ODE geometry and quantitative two-clock semigroup bounds are developed in Appendix B. The stationary distribution μ of the exogenous chain defines the averaged update fields H¯(x,θ) H(x,θ) =∑y∈μ⁡(y)H(y,x,θ), = _y∈ Yμ(y)H(y,x,θ), (31) F¯(x,θ) F(x,θ) =∑y∈μ⁡(y)F(y,x,θ). = _y∈ Yμ(y)F(y,x,θ). (32) These fields describe the deterministic motion obtained after averaging the rapid random fluctuations of YkY_k. On the fast time scale t=kηt=kη, the formal averaged ODE is dxdt dxdt =H¯(x,θ), = H(x,θ), (33) dθdt dθdt =ρF¯(x,θ). =ρ\, F(x,θ). (34) When ρ is small, θ barely changes while the fast ODE relaxes. This suggests freezing θ and studying dxdt=H¯(x,θ). dxdt= H(x,θ). (35) Assume that, for every θ in the region of interest, this frozen ODE has a locally unique asymptotically stable equilibrium. Denote it by λ⁡(θ)λ(θ), so H¯(λ(θ),θ)=0. H(λ(θ),θ)=0. (36) The graph x=λ⁡(θ)x=λ(θ) is the fast equilibrium manifold. After the fast transient has died away, the pair (x,θ)(x,θ) should remain close to this graph while θ evolves. Substituting x=λ⁡(θ)x=λ(θ) into the averaged slow update gives the reduced slow field g⁡(θ)=F¯(λ⁡(θ),θ).g(θ)= F(λ(θ),θ). (37) On slow time s=kεs=k , the reduced ODE is dθds=g⁡(θ). dθds=g(θ). (38) Let θ⋆θ be its stable equilibrium, so that g⁡(θ⋆)=0g(θ )=0. The corresponding equilibrium of the full averaged system is x⋆=λ⁡(θ⋆),z⋆=(x⋆θ⋆).x =λ(θ ), z = pmatrixx \\ θ pmatrix. (39) For vector-valued variables, local stability is expressed through the Jacobians A=DxH¯(x⋆,θ⋆),S=Dθg(θ⋆).A=D_x H(x ,θ ), S=D_θg(θ ). (40) We assume that both A and S are Hurwitz: every eigenvalue has strictly negative real part. For initial conditions sufficiently close to the relevant equilibrium, this gives bounds of the form Ce−ctCe^-ct times the initial displacement for the frozen fast ODE, and Ce−csCe^-cs times the initial displacement for the reduced slow ODE, with suitable constants C,c>0C,c>0. Thus both deterministic systems approach their equilibria at an exponential rate, measured on their respective time scales t=kηt=kη and s=kεs=k . We call this exponential approach toward equilibrium exponential relaxation, and use relaxation as shorthand below. Local Hurwitz stability does not ensure an invariant law for the unprojected stochastic recursion; global localization is imposed in . A.3 Coordinates Relative to the Fast Equilibrium Manifold The raw error x−x⋆x-x combines two effects: displacement away from the fast manifold and movement of the manifold caused by θ−θ⋆θ-θ . It is more informative to introduce u^=x−λ⁡(θ),v=θ−θ⋆. u=x-λ(θ), v=θ-θ . (41) These coordinates are initially defined where the local map λ is available. The variable u u measures displacement from the moving fast equilibrium, while v measures slow displacement. Their use throughout the stationary law requires the extension specified in Assumption A7. One update shows why this coordinate is useful. Let x+=x+ηH(Y,x,θ)x^+=x+η H(Y,x,θ) and θ+=θ+εF(Y,x,θ)θ^+=θ+ F(Y,x,θ). When the segment [θ,θ+][θ,θ^+] lies in the domain of the chosen C2C^2 map λ, we have u^+ u^+ =x+−λ⁡(θ+) =x^+-λ(θ^+) =u^+ηH(Y,x,θ)−λ⁡(θ+εF(Y,x,θ))−λ⁡(θ). = u+η H(Y,x,θ)-\λ(θ+ F(Y,x,θ))-λ(θ)\. (42) For λ:ℝdθ→ℝdxλ:R^d_θ ^d_x, the Jacobian with respect to θ has dimensions Dθλ(θ)∈ℝdx×dθ,Dθλ(θ)F(Y,x,θ)∈ℝdx.D_θλ(θ) ^d_x× d_θ, D_θλ(θ)F(Y,x,θ) ^d_x. If λ has a locally Lipschitz derivative (in particular, if λ∈C2λ∈ C^2), Taylor expansion gives λ⁡(θ+εF)−λ⁡(θ)=εDθλ(θ)F+O⁡(ε2‖F‖2).λ(θ+ F)-λ(θ)= D_θλ(θ)F+O( ^2\|F\|^2). Using ε=ηρ =ηρ yields u^+ u^+ =u^+η⁡[H⁡(Y,x,θ)−ρDθλ(θ)F(Y,x,θ)]+O⁡(ε2‖F⁡(Y,x,θ)‖2). = u+η [H(Y,x,θ)-ρ D_θλ(θ)F(Y,x,θ) ]+O( ^2\|F(Y,x,θ)\|^2). (43) The factor ρ=ε/ηρ= /η is now visible in front of the slow forcing of the fast-manifold error. This factor must be retained. A generic bound for the stacked recursion can replace the two relaxation rates by the slow spectral gap, producing a spurious factor 1/ρ1/ρ. In the coordinate (u^,v)( u,v), slow forcing has size ρ and acts over a relaxation time of order 1/ρ1/ρ on the fast clock; the two factors compensate. Near z⋆z , reserve u for the linear approximation of the nonlinear tracking error: u=(x−x⋆)−Dθλ(θ⋆)(θ−θ⋆)and henceu^=u+O⁡(‖θ−θ⋆‖2).u=(x-x )-D_θλ(θ )(θ-θ ) hence u=u+O(\|θ-θ \|^2). (44) To show that the difference between u u and its linear approximation u is O⁡(‖θ−θ⋆‖2)O(\|θ-θ \|^2), let δx=x−x⋆δ x=x-x , v=θ−θ⋆v=θ-θ , and Λ=Dθλ(θ⋆) =D_θλ(θ ). Since x⋆=λ⁡(θ⋆)x =λ(θ ), the C2C^2 Taylor formula is λ⁡(θ⋆+v)=x⋆+Λv+rλ(v),rλ(v)=∫01(1−s)Dθ2λ(θ⋆+sv)[v,v]s.λ(θ +v)=x + v+r_λ(v), r_λ(v)= _0^1(1-s)D_θ^2λ(θ +sv)[v,v]\,ds. A bound ‖Dθ2λ‖≤M\|D_θ^2λ\|≤ M along the segment gives ‖rλ(v)‖≤(M/2)‖v‖2\|r_λ(v)\|≤(M/2)\|v\|^2. Therefore u^=x−λ⁡(θ)=δx−Λv−rλ(v)=u−rλ(v). u=x-λ(θ)=δ x- v-r_λ(v)=u-r_λ(v). This proves equation 44 in vector norm; the minus sign is absorbed by the big-O notation. Differentiability alone would give only u^=u+o⁡(‖v‖) u=u+o(\|v\|). The nonlinear version is useful for localization and transient estimates; the linearized version is convenient when computing covariance and bias coefficients. A.4 Compact Notation Define the stacked state and update z=(xθ),G⁡(y,z)=(H⁡(y,x,θ)F⁡(y,x,θ)),z= pmatrixx\\ θ pmatrix, G(y,z)= pmatrixH(y,x,θ)\\ F(y,x,θ) pmatrix, (45) and define the block scaling matrix Dρ=(Idx00ρIdθ).D_ρ= pmatrixI_d_x&0\\ 0&ρ I_d_θ pmatrix. (46) Here IdxI_d_x and IdθI_d_θ are identity matrices of dimensions dxd_x and dθd_θ. The recursions in equations 27 and 28 are then exactly equivalent to zk+1=zk+ηDρG(Yk,zk).z_k+1=z_k+η D_ρG(Y_k,z_k). (47) The lower block has the smaller step ηρ=εηρ= . A.5 Standing Assumptions The core stationary assumptions are stated in Section 2.2; the additional coupling and block conditions are stated separately. Assumptions A1 and A2. Local smoothness and Hurwitz stability specify the deterministic geometry and the fixed tangent transformation. Assumption A3. The exogenous finite-state chain satisfies the stated Doeblin condition. Assumption A4. A selected family of joint invariant laws satisfies ‖z−z⋆‖4≤Cη2E\|z-z \|^4≤ Cη^2, uniformly in the step sizes. Assumption A5. The update fields have globally bounded first and second derivatives. Lemma 4 derives the expected Taylor and Poisson-corrector bounds from these conditions. Assumption A6. The L1L^1 coupling controls initialization of expected tail averages. The blockwise mean-square statements additionally require the coordinate L2L^2 coupling bounds and the uniform stationary V2V_2 moment displayed in Appendix I. Assumption A7. The stronger raw and centered block moments, tangent cross-covariance bound, and stationary extension and remainder conditions for λ support supplementary blockwise fluctuation statements. They are not assumptions of the core covariance and bias expansions. The local deterministic and Markov assumptions determine the coefficients; fourth-moment localization and update regularity justify the stationary expansions. Additional block conditions and coupling are invoked only for the conclusions that use them. A.6 Quantities of Interest and Step-Size Regimes Suppose that for sufficiently small (η,ρ)(η,ρ) the joint chain (Yk,zk)(Y_k,z_k) admits an invariant law, denoted by πη,ρ _η,ρ. The stationary mean displacement from the averaged equilibrium is mη,ρ=πη,ρ[z]−z⋆.m_η,ρ=E_ _η,ρ[z]-z . (48) The centered stationary covariance is Γη,ρ=πη,ρ[(z−πη,ρz)(z−πη,ρz)⊤]. _η,ρ=E_ _η,ρ [(z-E_ _η,ρz)(z-E_ _η,ρz) ]. (49) The superscript ⊤ denotes matrix transpose. The corresponding mean-square distance decomposes as πη,ρ‖z−z⋆‖2=tr⁡(Γη,ρ)+‖mη,ρ‖2,E_ _η,ρ\|z-z \|^2=tr( _η,ρ)+\|m_η,ρ\|^2, (50) where trtr is the matrix trace and ∥⋅∥\|·\| is the Euclidean norm. Finite-time concentration, stationary covariance, and stationary bias expansions thus address distinct properties of the recursion. For a trajectory started from an arbitrary initial condition, let K≥0K≥ 0 be the burn-in length and let N>KN>K. We average the iterates zK,…,zN−1z_K,…,z_N-1, so the averaging window contains L=N−KL=N-K iterates. The tail average is z¯K,N=1L∑k=KN−1zk. z_K,N= 1L _k=K^N-1z_k. (51) Let ν denote the initial law of the joint chain. Adding and subtracting the invariant mean gives the exact decomposition ν[z¯K,N]−z⋆ _ν[ z_K,N]-z =mη,ρ+K,N(ν), =m_η,ρ+ T_K,N(ν), K,N(ν) T_K,N(ν) :=1L∑k=KN−1(ν[zk]−πη,ρ[z]). := 1L _k=K^N-1 (E_ν[z_k]-E_ _η,ρ[z] ). The first term is the persistent stationary displacement. Under stationary initialization, K,N=0 T_K,N=0 for every K,NK,N; averaging reduces sampling variability but leaves the mean displacement unchanged. The second term is the initialization transient. To bound the initialization transient K,N(ν) T_K,N(ν), we need an estimate of how quickly the expected iterates approach the stationary mean. Theorem 3 bounds this contribution under the two-clock coupling assumption. The ratio ρ determines the form of the persistent displacement. Fixed-ratio analysis holds ρ∈(0,ρ0]ρ∈(0, _0] fixed and lets η↓0η 0 with ε=ρη =ρη, giving a representation mη,ρ=ηb(ρ)+Rη,ρ.m_η,ρ=η b(ρ)+R_η,ρ. Singular-ratio analysis also lets ρ↓0ρ 0. At ρ=0ρ=0 the slow block of DρD_ρ vanishes, the slow equation freezes, and its spectral gap collapses. A remainder bound at each fixed ρ need not be uniform in this limit. Appendix E develops the ratio expansion and explains how its powers determine the extrapolation weights. Appendix B ODE geometry of two timescales This section explains how coordinates relative to the fast equilibrium manifold separate fast tracking from slow evolution. As ρ→0ρ→ 0, slow relaxation takes increasingly long on the fast clock. We therefore keep the two relaxation rates and the factors of ρ in the coupling terms explicit. We first derive the exact coordinate dynamics for both the averaged ODE and the discrete recursion (Subsection B.1). We then linearize at the equilibrium to identify the fast and reduced slow Jacobian blocks (Subsection B.2). This linear structure leads to bounds that retain both relaxation rates, with constants independent of sufficiently small positive ρ (Subsection B.3). These deterministic calculations provide tools for the stationary covariance and bias proofs developed later; stationary existence and localization require additional assumptions. B.1 Exact Fast-Manifold Dynamics We use the clocks (30), the fast equilibrium graph (36), the reduced field (37), and the fast-manifold coordinates (u^,v)( u,v) of (41), and derive the dynamics in these coordinates. To distinguish the exact coordinate u u from its tangent approximation u, use a neutral dummy variable w∈ℝdxw ^d_x and define h⁡(w,v) h(w,v) =H¯(λ⁡(θ⋆+v)+w,θ⋆+v), = H (λ(θ +v)+w,θ +v ), (52) f⁡(w,v) f(w,v) =F¯(λ⁡(θ⋆+v)+w,θ⋆+v), = F (λ(θ +v)+w,θ +v ), (53) ℓ⁡(v) (v) =Dθλ(θ⋆+v)∈ℝdx×dθ. =D_θλ(θ +v) ^d_x× d_θ. (54) Here ℓ⁡(v) (v) is the Jacobian of λ with respect to θ (equivalently, with respect to v); it is not a time derivative. At the actual state w=u^w= u, since x=λ⁡(θ⋆+v)+u^x=λ(θ +v)+ u. The chain rule applied to the fast-clock ODE equations 33 and 34 gives the exact fast-clock equations u^˙=h⁡(u^,v)−ρℓ(v)f(u^,v),v˙=ρf(u^,v). u=h( u,v)-ρ (v)f( u,v), v=ρ f( u,v). (55) The fast-manifold identity (36) becomes h⁡(0,v)=0for every admissible v.h(0,v)=0 every admissible v. (56) Thus a tangential displacement v cannot create an order-one forcing in the normal equation while u^=0 u=0. The only such forcing visible in (55) carries the explicit factor ρ. A one-timescale norm bound for the stacked system obscures these factors. There is an exact discrete counterpart. Abbreviate Hk=H⁡(Yk,xk,θk),Fk=F⁡(Yk,xk,θk).H_k=H(Y_k,x_k, _k), F_k=F(Y_k,x_k, _k). (57) Set Δk=θk+1−θk=ηρFk _k= _k+1- _k=ηρ F_k. Whenever the segment [θk,θk+1][ _k, _k+1] lies in the domain of the chosen C2C^2 extension of λ, applying the fundamental theorem of calculus twice gives λ⁡(θk+Δk)−λ⁡(θk) λ( _k+ _k)-λ( _k) =∫01Dθλ(θk+sΔk)Δks = _0^1D_θλ( _k+s _k) _k\,ds =Dθλ(θk)Δk+∫01∫0sDθ2λ(θk+tΔk)[Δk,Δk]ts =D_θλ( _k) _k+ _0^1 _0^sD_θ^2λ( _k+t _k)[ _k, _k]\,dt\,ds =Dθλ(θk)Δk+∫01(1−s)Dθ2λ(θk+sΔk)[Δk,Δk]ds. =D_θλ( _k) _k+ _0^1(1-s)D_θ^2λ( _k+s _k)[ _k, _k]\,ds. The last equality interchanges the order of integration. Subtracting this identity from xk+1−xk=ηHkx_k+1-x_k=η H_k gives u^k+1−u^k= u_k+1- u_k= ηHk−ηρDθλ(θk)Fk η H_k-ηρ D_θλ( _k)F_k −η2ρ2∫01(1−s)Dθ2λ(θk+sηρFk)[Fk,Fk]ds, -η^2ρ^2 _0^1(1-s)D_θ^2λ( _k+sηρ F_k)[F_k,F_k]\,ds, (58) vk+1−vk= v_k+1-v_k= ηρFk. ηρ F_k. (59) Here Dθ2λ(θ)[a,b]∈ℝdxD_θ^2λ(θ)[a,b] ^d_x is a vector-valued bilinear second derivative. The factor 1/21/2 is already contained in ∫01(1−s)s=1/2 _0^1(1-s)\,ds=1/2. No derivative of FkF_k enters: FkF_k is held fixed as the direction of the one-step increment. This identity involves no asymptotic approximation. Because u^k=xk−λ⁡(θk) u_k=x_k-λ( _k), updating θk _k also changes the reference fast equilibrium λ⁡(θk)λ( _k). The first-order contribution of this change to the tracking error is −ηρDθλ(θk)Fk-ηρ D_θλ( _k)F_k, which has prefactor ηρηρ. The curvature remainder in (58) accounts for the nonlinear part of this change and has prefactor η2ρ2η^2ρ^2. After division by η to express increments per unit of fast time, these prefactors become ρ and ηρ2ηρ^2, respectively. The curvature remainder therefore contains an additional step-size factor as well as two powers of ρ. If the remaining coefficient factors are uniformly controlled, integration against a slow kernel gives the scaling estimates ∫0∞e−cρtρt=1c,∫0∞e−cρtηρ2t=ηρc. _0^∞e^-cρ tρ\,dt= 1c, _0^∞e^-cρ tηρ^2\,dt= ηρc. The preceding integrals illustrate how a factor ρ in a forcing term offsets the factor 1/ρ1/ρ arising from slow relaxation. The stationary covariance proof verifies this cancellation directly in the fixed tangent coordinates (u,v)(u,v). When the covariance equation is partitioned according to these coordinates, its (v,v)(v,v) block corresponds to covariances between components of the slow variable v=θ−θ⋆v=θ-θ . We call this the slow–slow block. The remainder in this block is O⁡(ρη3/2)O(ρη^3/2), as shown in (104). The inverse estimate in Lemma 6 weights this remainder by 1/ρ1/ρ. The factor ρ in the remainder therefore cancels this weight, giving an O⁡(η3/2)O(η^3/2) contribution to the covariance error, uniformly in sufficiently small positive ρ. B.2 Linearization at the Equilibrium Recall the four Jacobian blocks at z⋆z introduced in Subsection 2.1: A=DxH¯(z⋆),B=DθH¯(z⋆),C=DxF¯(z⋆),D=DθF¯(z⋆).A=D_x H(z ), B=D_θ H(z ), C=D_x F(z ), D=D_θ F(z ). (60) Differentiating (36) at θ⋆θ gives AΛ+B=0,Λ=Dλ(θ⋆)=−A−1B.A +B=0, =Dλ(θ )=-A^-1B. (61) The inverse A−1A^-1 is available when the frozen fast equilibrium is locally exponentially stable. The reduced field was defined in equation 37 by g⁡(θ)=F¯(λ⁡(θ),θ)g(θ)= F(λ(θ),θ). Its chain rule is Dθg(θ⋆)=DxF¯(x⋆,θ⋆)Dθλ(θ⋆)+DθF¯(x⋆,θ⋆).D_θg(θ )=D_x F(x ,θ )D_θλ(θ )+D_θ F(x ,θ ). Substituting the Jacobian blocks gives S=Dg(θ⋆)=CΛ+D=D−CA−1B.S=Dg(θ )=C +D=D-CA^-1B. (62) For the local calculation, use the tangent coordinates (u,v)(u,v) of (44). Write δx=x−x⋆δ x=x-x and δθ=θ−θ⋆δθ=θ-θ . Retaining the first-order terms of the fast-clock ODE at z⋆z gives the first-variation system dt(δxδθ)=(ABρCρD)(δxδθ), ddt pmatrixδ x\\ δθ pmatrix= pmatrixA&B\\ ρ C&ρ D pmatrix pmatrixδ x\\ δθ pmatrix, (63) For finite deviations with r=‖(δx,δθ)‖r=\|(δ x,δθ)\|, the omitted remainders are O⁡(r2)O(r^2) in the first row and ρO(r2)ρ O(r^2) in the second. Now substitute δx=u+Λvδ x=u+ v. The lower row gives v˙=ρC(u+Λv)+ρDv=ρCu+ρ⁡(CΛ+D)v=ρCu+ρSv. v=ρ C(u+ v)+ρ Dv=ρ Cu+ρ(C +D)v=ρ Cu+ρ Sv. Because Λ is fixed in this tangent-coordinate calculation, u˙=δx˙−Λv˙ u= δ x- v. Thus u˙ u =A⁡(u+Λv)+Bv−Λ⁡ρC(u+Λv)+ρDv =A(u+ v)+Bv- \ρ C(u+ v)+ρ Dv\ =(A−ρΛC)u+AΛ+B−ρΛ(CΛ+D)v =(A-ρ C)u+\A +B-ρ (C +D)\v =(A−ρΛC)u−ρΛSv. =(A-ρ C)u-ρ Sv. The identities AΛ+B=0A +B=0 and CΛ+D=SC +D=S therefore give dt(uv)=Lρ(uv),Lρ=(A−ρΛC−ρΛSρCρS). ddt pmatrixu\\ v pmatrix=L_ρ pmatrixu\\ v pmatrix, L_ρ= pmatrixA-ρ C&-ρ S\\ ρ C&ρ S pmatrix. (64) Every off-diagonal interaction in these coordinates has a factor ρ. At ρ=0ρ=0, the u-block is the stable matrix A and the v-block is frozen. For small positive ρ, the diagonal slow block ρSρ S is Hurwitz, and the complete slow equation is v˙=ρSv+ρCu v=ρ Sv+ρ Cu. The term ρCuρ Cu describes the coupling from the fast deviation to the slow variable. B.3 Two-Clock Semigroup Bounds We now quantify relaxation in the linearized system: etLρe^tL_ρ maps the initial deviation (u⁡(0),v⁡(0))(u(0),v(0)) to (u⁡(t),v⁡(t))(u(t),v(t)). The proof separates fast and slow invariant subspaces to obtain bounds that retain both relaxation rates and the factors ρ in their interactions, uniformly for small positive ρ, for both the ODE and its Euler discretization. Proposition 1 (Two-clock semigroup bounds; restatement). Under and , there exist constants C∗,c,ρ0>0C_*,c, _0>0, independent of ρ, such that for 0<ρ≤ρ00<ρ≤ _0 and t≥0t≥ 0, ‖[etLρ]uu‖ \|[e^tL_ρ]_u\| ≤C∗(e−ct+ρe−cρt), ≤ C_* (e^-ct+ρ e^-cρ t ), (65) ‖[etLρ]uv‖+‖[etLρ]vu‖ \|[e^tL_ρ]_uv\|+\|[e^tL_ρ]_vu\| ≤C∗ρ(e−ct+e−cρt), ≤ C_*ρ (e^-ct+e^-cρ t ), ‖[etLρ]vv‖ \|[e^tL_ρ]_v\| ≤C∗e−cρt. ≤ C_*e^-cρ t. The same bounds hold for (I+ηLρ)k(I+η L_ρ)^k, with t=ηkt=η k, uniformly for sufficiently small η and integers k≥0k≥ 0. The block subscripts specify the output and input coordinates, in that order. Thus [etLρ]uv[e^tL_ρ]_uv gives the contribution to u⁡(t)u(t) from the initial slow deviation v⁡(0)v(0), while [etLρ]vu[e^tL_ρ]_vu gives the contribution to v⁡(t)v(t) from the initial fast deviation u⁡(0)u(0). The factors e−cte^-ct and e−cρte^-cρ t describe fast and slow relaxation, respectively. Proof. Write the blocks of LρL_ρ as Lρ=(AρℬρCρDρ′),Aρ=A−ρΛC,ℬρ=−ρΛS,Cρ=ρC,Dρ′=ρS.L_ρ= pmatrixA_ρ&B_ρ\\ C_ρ&D_ρ pmatrix, A_ρ=A-ρ C, _ρ=-ρ S, C_ρ=ρ C, D_ρ =ρ S. We seek invariant subspaces of the form u=Xρvu=X_ρv and v=Rρuv=R_ρu: a solution starting in either subspace remains there. These describe the perturbed slow and fast subspaces, respectively. Their invariance equations are AρXρ+ℬρ−Xρ(CρXρ+Dρ′) A_ρX_ρ+B_ρ-X_ρ(C_ρX_ρ+D_ρ ) =0, =0, Cρ+Dρ′Rρ−Rρ(Aρ+ℬρRρ) C_ρ+D_ρ R_ρ-R_ρ(A_ρ+B_ρR_ρ) =0. =0. At ρ=0ρ=0, the derivatives in X and R are respectively X↦AX AX and R↦−RAR -RA, both invertible. The implicit-function theorem therefore gives Xρ,Rρ=O⁡(ρ)X_ρ,R_ρ=O(ρ). The matrix Qρ=(IXρRρI)Q_ρ= pmatrixI&X_ρ\\ R_ρ&I pmatrix is uniformly invertible for small ρ, and Qρ−1LρQρ=diag⁡(Af(ρ),ρSs(ρ)),Af=Aρ+ℬρRρ=A+O⁡(ρ),Ss=S+CXρ=S+O⁡(ρ).Q_ρ^-1L_ρQ_ρ=diag(A_f(ρ),ρ S_s(ρ)), A_f=A_ρ+B_ρR_ρ=A+O(ρ), S_s=S+CX_ρ=S+O(ρ). Fixed quadratic Lyapunov functions for the Hurwitz matrices A and S remain strict Lyapunov functions for these perturbations. Consequently, ‖etAf‖≤Ce−ct,‖eρtSs‖≤Ce−cρt,\|e^tA_f\|≤ Ce^-ct, \|e^ρ tS_s\|≤ Ce^-cρ t, with constants independent of ρ. Finally, writing Qρ=I+EρQ_ρ=I+E_ρ, where EρE_ρ has only off-diagonal blocks of size O⁡(ρ)O(ρ), the Neumann series Qρ−1=I−Eρ+Eρ2−⋯Q_ρ^-1=I-E_ρ+E_ρ^2-·s shows that its diagonal corrections are O⁡(ρ2)O(ρ^2) and its off-diagonal blocks are O⁡(ρ)O(ρ). Multiplication in etLρ=Qρdiag⁡(etAf,eρtSs)Qρ−1e^tL_ρ=Q_ρdiag(e^tA_f,e^ρ tS_s)Q_ρ^-1 then preserves the two clocks and the small cross factors. Multiplying the block matrices proves (65); in fact, the upper-left slow contribution and lower-right fast contribution are both O⁡(ρ2)O(ρ^2). For completeness, the same Lyapunov functions give, uniformly for small η, contraction factors 1−cη1-cη and 1−cηρ1-cηρ for the squared Lyapunov norms of I+ηAfI+η A_f and I+ηρSsI+ηρ S_s, respectively. Indeed, if PfP_f is the fixed fast Lyapunov matrix, then (I+ηAf)⊤Pf(I+ηAf)−Pf=η⁡(Af⊤Pf+PfAf)+η2Af⊤PfAf⪯−cηPf,(I+η A_f) P_f(I+η A_f)-P_f=η(A_f P_f+P_fA_f)+η^2A_f P_fA_f -cη P_f, after decreasing η0 _0; the slow argument is identical with η replaced by ηρηρ. Thus the same block bounds hold for (I+ηLρ)k(I+η L_ρ)^k, with t=ηkt=η k. ∎ The upper equation of equation 64 is u˙=(A−ρΛC)u−ρΛSv. u=(A-ρ C)u-ρ Sv. To isolate the upper-right channel, consider a one-way scalar model with coupling β and stable rates −a,−s-a,-s. Prescribe the uncoupled slow response v⁡(t)=e−ρstv0v(t)=e^-ρ stv_0 and solve u˙=−au+ρβv,a,s>0. u=-au+ρβ v, a,s>0. (66) Then u⁡(t) u(t) =e−atu0+ρβv0∫0te−a⁡(t−r)e−ρsrr =e^-atu_0+ρβ v_0 _0^te^-a(t-r)e^-ρ sr\,dr =e−atu0+ρβv0e−ρst−e−ata−ρs. =e^-atu_0+ρβ v_0 e^-ρ st-e^-ata-ρ s. (67) For small ρ, the denominator is bounded away from zero and the slow tail in u has amplitude O⁡(ρ)O(ρ). The reverse channel is illustrated by a scalar coupling γ in v˙=−ρsv+ργu v=-ρ sv+ργ u. Prescribing the uncoupled fast response u⁡(t)=e−atu0u(t)=e^-atu_0 gives v˙+ρsv=ργu0e−at,dteρstv(t)=ργu0e(ρs−a)t. v+ρ sv=ργ u_0e^-at, ddt\e^ρ stv(t)\=ργ u_0e^(ρ s-a)t. Integrating from zero to t, for a≠ρsa≠ρ s, eρstv(t)−v0=ργu0e(ρs−a)t−1ρs−a.e^ρ stv(t)-v_0=ργ u_0 e^(ρ s-a)t-1ρ s-a. Multiplication by e−ρste^-ρ st yields v⁡(t)=e−ρstv0+ργu0e−ρst−e−ata−ρs,v(t)=e^-ρ stv_0+ργ u_0 e^-ρ st-e^-ata-ρ s, (68) which again has cross amplitude O⁡(ρ)O(ρ) for sufficiently small ρ. At the exceptional resonant value a=ρsa=ρ s, the continuous limiting formula is v⁡(t)=e−atv0+ργu0t.v(t)=e^-at\v_0+ργ u_0t\. The same denominator caveat applies to equation 67; at resonance its convolution is te−atte^-at. One may choose ρ0<a/s _0<a/s to avoid resonance throughout the small-ratio range. These one-way formulas illustrate the two cross bounds in equation 65. The slow exponential integrates to O⁡(1/ρ)O(1/ρ) on the fast clock: ∫0∞e−cρtt=1cρ. _0^∞e^-cρ t\,dt= 1cρ. (69) If a forcing term has size O⁡(ρ)O(ρ) with a constant independent of ρ, integration against this slow exponential gives an O⁡(1)O(1) bound. Discarding the factor ρ before integration can therefore introduce an unnecessary factor 1/ρ1/ρ in the estimate. Returning to iteration number k, t=kηt=kη changes the two exponentials into e−cηkande−cρηk=e−cεk.e^-cη k e^-cρη k=e^-c k. (70) These linear bounds identify the fast and slow relaxation rates relevant to initialization effects. Quantitative initialization bounds for the nonlinear stochastic recursion are developed in Appendix I under the two-clock coupling conditions. Stationary bias calculations start from an invariant law and do not require an initialization estimate. Appendix C Poisson tools and Taylor estimates We first compute the long-run covariance Q of the noise ξk=G⁡(Yk,z⋆) _k=G(Y_k,z ), obtained by holding the algorithmic state at the equilibrium z⋆z while the Markov chain continues to evolve (Subsection C.1). Because successive Markov inputs are correlated, the accumulated noise fluctuations depend on correlations across time, which Q records. This matrix supplies the noise term in the Lyapunov equation for the leading stationary covariance of the iterates. That covariance enters the quadratic Taylor terms of the averaged drift and therefore determines the curvature contribution to the stationary mean bias. Next, for a fixed algorithmic state z, we use the Poisson equation to decompose the centered Markov update into a martingale difference and a telescoping correction, whose sum reduces to initial and final terms (Subsection C.2). We then derive Taylor estimates for the nonlinear remainder terms and for changes in the Poisson solution as the algorithmic state evolves (Subsection C.3). Appendix D uses these tools to establish the stationary covariance and mean-bias expansions. C.1 Long-Run Covariance of the Markov Noise Recall the compact recursion zk+1=zk+ηDρG(Yk,zk),z_k+1=z_k+η D_ρG(Y_k,z_k), where z=(x,θ)z=(x,θ), Dρ=diag⁡(I,ρI)D_ρ=diag(I,ρ I), and ρ=ε/ηρ= /η. Average the update over the stationary law μ of the exogenous Markov chain: G¯(z)=∑yμ⁡(y)G(y,z), define G~(y,z)=G⁡(y,z)−G¯(z). G(z)= _yμ(y)G(y,z), define G(y,z)=G(y,z)- G(z). (71) At the averaged equilibrium z⋆z , we have G¯(z⋆)=0 G(z )=0. The update noise frozen at that point is therefore ξk=G⁡(Yk,z⋆)=G~(Yk,z⋆). _k=G(Y_k,z )= G(Y_k,z ). (72) For each fixed z, G~(⋅,z) G(·,z) is centered under μ. In Sections C.1 and C.2, expectations refer to the stationary Markov path: Y0∼μY_0 μ, followed by transitions with kernel P. In particular, ξk=0E _k=0. This convention specifies the noise covariance without requiring stationary initialization of the algorithm. With Markov data the variance of a long sum contains all cross-time covariances, so the relevant matrix is Q=limm→∞1m[(∑k=0m−1ξk)(∑k=0m−1ξk)⊤].Q= _m→∞ 1mE [ ( _k=0^m-1 _k ) ( _k=0^m-1 _k )^\! ]. (73) Thus Q is the covariance accumulated per step. Write Sm=∑k=0m−1ξkS_m= _k=0^m-1 _k and Kh=⁡[ξ0ξh⊤]K_h=E[ _0 _h ]. Since the sum is centered, Cov⁡(Sm)=⁡[SmSm⊤]=∑i=0m−1∑j=0m−1⁡[ξiξj⊤].Cov(S_m)=E[S_mS_m ]= _i=0^m-1 _j=0^m-1E[ _i _j ]. Stationarity makes each diagonal term equal to K0K_0. For a positive lag h, the m−hm-h pairs with j=i+hj=i+h each contribute KhK_h, whereas the m−hm-h reversed pairs each contribute Kh⊤K_h . Thus the exact finite-sum identity is 1m[SmSm⊤]=K0+∑h=1m−1(1−hm)(Kh+Kh⊤). 1mE[S_mS_m ]=K_0+ _h=1^m-1 (1- hm )(K_h+K_h ). Suppose ∑h≥1‖Kh‖<∞ _h≥ 1\|K_h\|<∞. For every fixed h, the factor h<m(1−h/m)1_\h<m\(1-h/m) tends to one and lies in [0,1][0,1]. The norm of the corresponding summand is bounded by 2‖Kh‖2\|K_h\|, a summable sequence. The dominated convergence theorem for series, applied entrywise with summable dominating sequence 2‖Kh‖2\|K_h\|, therefore gives Q=K0+∑h=1∞(Kh+Kh⊤).Q=K_0+ _h=1^∞ (K_h+K_h ). (74) The summability condition follows from Assumption . Indeed, on the finite state space the centered function ξ⁡(y)=G⁡(y,z⋆)ξ(y)=G(y,z ) is bounded. Geometric mixing gives constants C<∞C<∞ and r∈(0,1)r∈(0,1) such that supy‖(Phξ)(y)‖≤Crhmaxy‖ξ⁡(y)‖. _y\|(P^hξ)(y)\|≤ Cr^h _y\|ξ(y)\|. Conditioning on Y0Y_0 then gives Kh=∑yμ⁡(y)ξ(y)(Phξ)(y)⊤,‖Kh‖≤Crh(maxy⁡‖ξ⁡(y)‖)2.K_h= _yμ(y)ξ(y)(P^hξ)(y) , \|K_h\|≤ Cr^h ( _y\|ξ(y)\| )^2. Consequently, (74) holds here. Stationarity alone would not suffice: for the stationary alternating chain Yk+1=−YkY_k+1=-Y_k, with equally likely initial signs and ξk=Yk _k=Y_k, the partial sums are bounded and Q=0Q=0, but Kh=(−1)hK_h=(-1)^h makes the displayed infinite series divergent. That periodic chain is excluded by geometric mixing. By contrast, ⁡[m−1∑k<mξkξk⊤]=K0E[m^-1 _k<m _k _k ]=K_0 omits every off-diagonal term, so Q generally differs from the one-step covariance. C.2 Martingale Decomposition via the Poisson Equation For a function a of the Markov state, define (Pa)(y)=∑y′P⁡(y,y′)a(y′).(Pa)(y)= _y P(y,y )a(y ). This is the conditional mean of a⁡(Yk+1)a(Y_k+1) given Yk=yY_k=y. For each fixed value of z consider the following equation ⁡(y,z)−(P)(y,z)=G~(y,z),∑yμ⁡(y)(y,z)=0.U(y,z)-(PU)(y,z)= G(y,z), _yμ(y)U(y,z)=0. (75) This is the Poisson equation. The second condition chooses one solution by forcing its stationary mean to be zero. For a finite-state geometrically ergodic chain, the solution is ⁡(⋅,z)=∑j=0∞PjG~(⋅,z).U(\,·\,,z)= _j=0^∞P^j G(\,·\,,z). (76) Define the one-step innovation ζk+1(z)=⁡(Yk+1,z)−(P)(Yk,z). _k+1(z)=U(Y_k+1,z)-(PU)(Y_k,z). (77) Let ℱkF_k contain the iterate and the complete Markov history up to time k. Exogeneity means that, given YkY_k, the next Markov state does not depend on the rest of ℱkF_k. Therefore ⁡[ζk+1(z)∣ℱk]=0.E[ _k+1(z) _k]=0. (78) Indeed, conditional on YkY_k, the expected value of ⁡(Yk+1,z)U(Y_k+1,z) is exactly (P)(Yk,z)(PU)(Y_k,z). The Poisson equation can now be rearranged as G~(Yk,z)=⁡(Yk,z)−⁡(Yk+1,z)+ζk+1(z). G(Y_k,z)=U(Y_k,z)-U(Y_k+1,z)+ _k+1(z). (79) Summing this identity makes the first two terms cancel across time: ∑k=0m−1G~(Yk,z)=⁡(Y0,z)−⁡(Ym,z)+∑k=0m−1ζk+1(z). _k=0^m-1 G(Y_k,z)=U(Y_0,z)-U(Y_m,z)+ _k=0^m-1 _k+1(z). (80) The endpoint difference remains bounded in mean square. The martingale sum has root-mean-square size O⁡(m)O( m). Fix z=z⋆z=z and write Mm=∑k=0m−1ζk+1(z⋆),Bm=(Y0,z⋆)−(Ym,z⋆),Qζ=⁡[ζ1(z⋆)ζ1(z⋆)⊤]. gatheredM_m= _k=0^m-1 _k+1(z ), B_m=U(Y_0,z )-U(Y_m,z ),\\ Q_ζ=E[ _1(z ) _1(z ) ]. gathered Thus Sm=Mm+BmS_m=M_m+B_m. For i<ji<j, the vector ζi(z⋆) _i(z ) is measurable with respect to ℱj−1F_j-1, so ⁡[ζi(z⋆)ζj(z⋆)⊤]=⁡[ζi(z⋆)[ζj(z⋆)⊤∣ℱj−1]]=0.E[ _i(z ) _j(z ) ]=E\! [ _i(z )E[ _j(z ) _j-1] ]=0. The expectation with the indices reversed is also zero: ⁡[ζj(z⋆)ζi(z⋆)⊤]=(⁡[ζi(z⋆)ζj(z⋆)⊤])⊤=0.E[ _j(z ) _i(z ) ]= (E[ _i(z ) _j(z ) ] ) =0. Since the Markov input chain is stationary (Y0∼μY_0 μ), the pair (Yk,Yk+1)(Y_k,Y_k+1) has the same distribution for every k. At the fixed state z⋆z , each innovation is the same function of this pair. Consequently, ⁡[ζk+1(z⋆)ζk+1(z⋆)⊤]=Qζfor every k≥0.E[ _k+1(z ) _k+1(z ) ]=Q_ζ every k≥ 0. Expanding the square therefore yields [MmMm⊤]=mQζ,∥Mm∥2=mtrQζ,‖Mm‖L2=mtr⁡Qζ, gatheredE[M_mM_m ]=mQ_ζ, \|M_m\|^2=mtrQ_ζ,\\ \|M_m\|_L^2= mtrQ_ζ, gathered where ‖X‖L2=(‖X‖2)1/2\|X\|_L^2=(E\|X\|^2)^1/2. This proves the O⁡(m)O( m) assertion in L2L^2. Since (a⊤Mm)2=ma⊤QζaE(a M_m)^2=m\,a Q_ζa, the scale is nondegenerate in a direction a if and only if a⊤Qζa>0a Q_ζa>0. The correction BmB_m consists of the values of the Poisson solution at the initial and final Markov states, Y0Y_0 and YmY_m. We now show that its contribution to the normalized covariance m−1[SmSm⊤]m^-1E[S_mS_m ], including its cross terms with MmM_m, tends to zero as m→∞m→∞. Put CU=2‖⁡(Y0,z⋆)‖L2<∞C_U=2\|U(Y_0,z )\|_L^2<∞. Stationarity and the triangle inequality give ‖Bm‖L2≤CU\|B_m\|_L^2≤ C_U, uniformly in m. Keeping every cross term gives 1m[SmSm⊤]−Qζ=1m[MmBm⊤+BmMm⊤+BmBm⊤]. 1mE[S_mS_m ]-Q_ζ= 1mE[M_mB_m +B_mM_m +B_mB_m ]. The inequality ‖xy⊤‖≤‖x‖‖y‖\|xy \|≤\|x\|\|y\| and Cauchy–Schwarz give ‖1m[SmSm⊤]−Qζ‖≤2CUtr⁡Qζm+CU2m⟶ 0. \| 1mE[S_mS_m ]-Q_ζ \|≤ 2C_U trQ_ζ m+ C_U^2m\ \ 0. We have proved the following. Lemma 3 (Poisson representation). Under , the innovations (77) are martingale differences, and the long-run covariance (73) satisfies Q=⁡[ζ1(z⋆)ζ1(z⋆)⊤].Q=E [ _1(z ) _1(z ) ]. (81) Example: A Symmetric Two-State Chain. Let Yk∈−1,1Y_k∈\-1,1\ be a stationary symmetric two-state Markov chain, with ⁡[Yk+1∣Yk]=qYk,−1<q<1.E[Y_k+1 Y_k]=qY_k, -1<q<1. (82) For an integer h≥0h≥ 0, the Markov property and the tower property of conditional expectation give ⁡[Yk+h+1∣Yk]=⁡[⁡[Yk+h+1∣ℱk+h]∣Yk]=q[Yk+h∣Yk].E[Y_k+h+1 Y_k]=E\! [E[Y_k+h+1 _k+h] Y_k ]=qE[Y_k+h Y_k]. Starting at h=0h=0, induction yields ⁡[Yk+h∣Yk]=qhYkE[Y_k+h Y_k]=q^hY_k, hence ⁡[YkYk+h]=⁡Yk[Yk+h∣Yk]=qhYk2=qh.E[Y_kY_k+h]=E\Y_kE[Y_k+h Y_k]\=q^hEY_k^2=q^h. Since Yk=0EY_k=0, this is also the autocovariance. For negative integer lags, stationarity gives q|h|q^|h|. Fix ω∈ℝdω ^d, d=dx+dθd=d_x+d_θ, and take G~(y,z⋆)=ωy G(y,z )=ω y. The frozen noise is ξk=ωYk _k=ω Y_k, and the Poisson solution is ⁡(y,z⋆)=ωy1−q,U(y,z )= ω y1-q, (83) because (P)(y,z⋆)=qωy/(1−q)(PU)(y,z )=qω y/(1-q). Hence ζk+1(z⋆)=ω⁡(Yk+1−qYk)1−q. _k+1(z )= ω(Y_k+1-qY_k)1-q. (84) Since Yk2=1Y_k^2=1 and ⁡[YkYk+1]=qE[Y_kY_k+1]=q, its covariance is ⁡[ζk+1(z⋆)ζk+1(z⋆)⊤] [ _k+1(z ) _k+1(z ) ] =1−2q2+q2(1−q)2ωω⊤=1+q1−qωω⊤. = 1-2q^2+q^2(1-q)^2ω = 1+q1-qω . (85) The autocovariance sum gives the same answer: Q=(1+2∑h=1∞qh)ωω⊤=1+q1−qωω⊤.Q= (1+2 _h=1^∞q^h )ω = 1+q1-qω . (86) Here K0=ωω⊤K_0=ω for every q, whereas Q depends strongly on q. C.3 Taylor Estimates for the Averaged Drift and Poisson Corrector The stationary covariance and mean calculations require control of two expansions: the averaged drift around the equilibrium, and the Poisson corrector across one algorithmic update. We also need to control the error from evaluating the corrector’s response coefficient at the equilibrium. Write ek=zk−z⋆,Δzk=zk+1−zk=ηDρG(Yk,zk),e_k=z_k-z , z_k=z_k+1-z_k=η D_ρG(Y_k,z_k), and let U be the normalized solution of (12). For the averaged drift, RG¯,k(2)R^(2)_ G,k is the error left after retaining the linear Taylor term, while RG¯,k(3)R^(3)_ G,k is the error left after also retaining the quadratic term. Since G¯(z⋆)=0 G(z )=0 and J=DG¯(z⋆)J=D G(z ), define RG¯,k(2) R^(2)_ G,k =G¯(zk)−Jek, = G(z_k)-Je_k, RG¯,k(3) R^(3)_ G,k =G¯(zk)−Jek−12D2G¯(z⋆)[ek,ek]. = G(z_k)-Je_k- 12D^2 G(z )[e_k,e_k]. Here D2G¯(z⋆)[ek,ek]D^2 G(z )[e_k,e_k] denotes the vector obtained by applying the Hessian of each component of G¯ G to (ek,ek)(e_k,e_k). Near the equilibrium, these remainders satisfy ‖RG¯,k(2)‖=O⁡(‖ek‖2),‖RG¯,k(3)‖=O⁡(‖ek‖3).\|R^(2)_ G,k\|=O(\|e_k\|^2), \|R^(3)_ G,k\|=O(\|e_k\|^3). Thus the superscripts indicate the local remainder orders. For the Poisson corrector, we expand ⁡(Yk+1,⋅)U(Y_k+1,·) from zkz_k to zk+1z_k+1, holding the Markov argument Yk+1Y_k+1 fixed. The error in its linear Taylor approximation is R,k R_U,k =⁡(Yk+1,zk+1)−⁡(Yk+1,zk) =U(Y_k+1,z_k+1)-U(Y_k+1,z_k) −Dz(Yk+1,zk)Δzk. -D_zU(Y_k+1,z_k) z_k. This remainder is quadratic in the one-step displacement Δzk z_k under the bounded second-derivative condition. The linear corrector increment contains the response coefficient Dz(Yk+1,zk)DρG(Yk,zk).D_zU(Y_k+1,z_k)D_ρG(Y_k,z_k). To compute its leading contribution, we evaluate both factors at z⋆z . Define the resulting coefficient error by Δ,k _U,k =Dz(Yk+1,zk)DρG(Yk,zk) =D_zU(Y_k+1,z_k)D_ρG(Y_k,z_k) −Dz(Yk+1,z⋆)DρG(Yk,z⋆). -D_zU(Y_k+1,z )D_ρG(Y_k,z ). The Markov arguments retain their original values; only the iterate argument is replaced by the equilibrium. Consequently, the corrector increment has the exact decomposition ⁡(Yk+1,zk+1)−⁡(Yk+1,zk) (Y_k+1,z_k+1)-U(Y_k+1,z_k) =ηDz(Yk+1,z⋆)DρG(Yk,z⋆)+ηΔ,k+R,k. =η D_zU(Y_k+1,z )D_ρG(Y_k,z )+η _U,k+R_U,k. This identity separates the equilibrium response, the error from freezing its coefficient, and the Taylor remainder over one step. The following lemma derives the expectation bounds below from the regularity and stationary localization assumptions. Expectations are taken along a stationary trajectory of the joint process (Yk,zk)(Y_k,z_k), and the constant C is independent of sufficiently small positive η and ρ: ‖RG¯,k(2)‖ \|R^(2)_ G,k\| ≤Cη, ≤ Cη, ‖RG¯,k(3)‖ \|R^(3)_ G,k\| ≤Cη3/2, ≤ Cη^3/2, ‖R,k‖ \|R_U,k\| ≤Cη2, ≤ Cη^2, ‖Δ,k‖ \| _U,k\| ≤Cη, ≤ C η, ∫01 _0^1 ‖Dz(Yk+1,zk+sΔzk)‖‖DρG(Yk,zk)‖ds≤C. \|D_zU(Y_k+1,z_k+s z_k)\|\,\|D_ρG(Y_k,z_k)\|\,ds≤ C. (87) ∥e∥j≤Cηj/2,j=2,3.E\|e\|^j≤ Cη^j/2, j=2,3. (88) Lemma 4 (Taylor estimates from fourth-moment localization). Suppose that the state space of the exogenous chain is finite and its transition kernel satisfies Assumption A3. Suppose that G⁡(y,⋅)∈C2(ℝdx+dθ)G(y,·)∈ C^2(R^d_x+d_θ) has uniformly bounded first and second derivatives, that G¯ G is C3C^3 on a neighborhood of z⋆z , and that a stationary law satisfies ‖ek‖4≤Cη2E\|e_k\|^4≤ Cη^2 uniformly for 0<ρ≤ρ0≤10<ρ≤ _0≤ 1. Then every estimate in (87) and (88) follows with a constant independent of η and ρ. The coordinate-change estimates (215) also follow if the chosen extension of λ has bounded second derivative on the segments joining θ⋆θ to the stationary values of θ. Proof. All norms below are Euclidean norms or the induced multilinear norms. For 1≤p≤41≤ p≤ 4, Hölder’s inequality gives ‖ek‖p≤(‖ek‖4)p/4≤Cηp/2.E\|e_k\|^p≤(E\|e_k\|^4)^p/4≤ Cη^p/2. The operator I−PI-P is not invertible on the full space of functions. It is invertible on the subspace of μ-centered functions under Assumption A3. The normalization μ=0 =0 in (12) selects precisely this inverse. Here are the details. Centering and the mixing bound. For a function a on Y, define (Πa)(y):=μa,‖a‖∞:=maxy∈⁡‖a⁡(y)‖.( a)(y):=μ a, \|a\|_∞:= _y∈ Y\|a(y)\|. Thus Π projects onto the constant functions, and ker⁡Π is the centered subspace. Invariance of μ gives Π2=Π,PΠ=ΠP=Π. ^2= , P = P= . We may decrease the minorization constant, if necessary, so that 0<β<10<β<1. The Doeblin condition then gives the decomposition Pn0(y,⋅)=βν(⋅)+(1−β)R(y,⋅)P^n_0(y,·)=βν(·)+(1-β)R(y,·) for a Markov kernel R. For probability measures α,α′α,α on Y, we use the total variation distance ‖α−α′‖TV:=supE⊆|α⁡(E)−α′(E)|=12∑y∈|α⁡(y)−α′(y)|.\|α-α \|_ TV:= _E Y|α(E)-α (E)|= 12 _y∈ Y|α(y)-α (y)|. The common βνβν term cancels when comparing two probability distributions, and a Markov kernel does not increase total variation distance. Therefore ‖αPn0−α′Pn0‖TV≤(1−β)‖α−α′‖TV.\|α P^n_0-α P^n_0\|_ TV≤(1-β)\|α-α \|_ TV. Iterating this inequality, using μP=μ P=μ, and treating the remaining fewer than n0n_0 steps by nonexpansiveness yields ‖Pn(y,⋅)−μ‖TV≤(1−β)⌊n/n0⌋.\|P^n(y,·)-μ\|_ TV≤(1-β) n/n_0 . Consequently, for scalar, vector, or derivative-tensor-valued a, ‖(Pn−Π)a‖∞≤2(1−β)⌊n/n0⌋‖a‖∞.\|(P^n- )a\|_∞≤ 2(1-β) n/n_0 \|a\|_∞. For vector or tensor values this follows by the triangle inequality applied to the finite signed sum over Y. The bounded Poisson solution operator. Define the fixed linear operator ℛP:=∑n=0∞(Pn−Π). R_P:= _n=0^∞(P^n- ). The preceding bound gives convergence in operator norm and the explicit estimate ‖ℛP‖∞→∞≤2∑n=0∞(1−β)⌊n/n0⌋=2n0β=:CP.\| R_P\|_∞→∞≤ 2 _n=0^∞(1-β) n/n_0 = 2n_0β=:C_P. For the partial sum through N, telescoping gives (I−P)∑n=0N(Pn−Π)=I−PN+1.(I-P) _n=0^N(P^n- )=I-P^N+1. Taking the limit, and using the centering identities, gives (I−P)ℛP=ℛP(I−P)=I−Π,ΠℛP=ℛPΠ=0.(I-P) R_P= R_P(I-P)=I- , R_P= R_P =0. Hence, for every centered h, u=ℛPhu= R_Ph satisfies (I−P)u=h(I-P)u=h and μu=0μ u=0. It is unique: if (I−P)v=0(I-P)v=0 and μv=0μ v=0, then v=Pnv=P^nv for every n, while the mixing bound forces ‖Pnv‖∞→0\|P^nv\|_∞→ 0. Equivalently, ℛP=(I−P+Π)−1−Π. R_P=(I-P+ )^-1- . Indeed, the preceding identities show that Π+ℛP + R_P is the two-sided inverse of I−P+ΠI-P+ . In particular, on a centered right-hand side, ℛPh=(I−P+Π)−1h R_Ph=(I-P+ )^-1h. Derivatives of the Poisson corrector. We now differentiate the Poisson corrector ⁡(y,z)U(y,z) with respect to z, keeping y fixed. Our goal is to derive uniform bounds on its first and second derivatives from the corresponding derivative bounds on the update field G. The equation in (12) is (I−P)(⋅,z)=G~(⋅,z),G~(⋅,z)=(I−Π)G(⋅,z),μ(⋅,z)=0.(I-P)U(·,z)= G(·,z), G(·,z)=(I- )G(·,z), (·,z)=0. Thus ⁡(⋅,z)=ℛPG~(⋅,z).U(·,z)= R_P G(·,z). Both P and μ are independent of z, so ℛP R_P is a fixed finite-dimensional matrix. Differentiating its finite matrix-vector product requires no interchange of an infinite series with a derivative. For j=1,2j=1,2, Dzj(⋅,z)=ℛP(DzjG(⋅,z)−μDzjG(⋅,z)).D_z^jU(·,z)= R_P (D_z^jG(·,z)-μ D_z^jG(·,z) ). The operator acts on the Markov-state index of these derivative tensors. Assumption A5 gives supy,z‖DzjG(y,z)‖≤K _y,z\|D_z^jG(y,z)\|≤ K, and hence supz‖DzjG(⋅,z)−μDzjG(⋅,z)‖∞≤2K,supy,z‖Dzj(y,z)‖≤2CPK. _z\|D_z^jG(·,z)-μ D_z^jG(·,z)\|_∞≤ 2K, _y,z\|D_z^jU(y,z)\|≤ 2C_PK. This is the precise meaning of commuting with differentiation. Growth bounds for the update field and Poisson corrector. We show that the update field G⁡(y,z)G(y,z) and the Poisson corrector ⁡(y,z)U(y,z) grow at most linearly with the distance ‖z−z⋆‖\|z-z \|, uniformly in the Markov state y. These growth bounds also control the one-step increment Δzk=ηDρG(Yk,zk) z_k=η D_ρG(Y_k,z_k). Set MG:=maxy⁡‖G⁡(y,z⋆)‖<∞M_G:= _y\|G(y,z )\|<∞. Since G¯(z⋆)=0 G(z )=0, we have G~(⋅,z⋆)=G⁡(⋅,z⋆) G(·,z )=G(·,z ) and therefore maxy⁡‖⁡(y,z⋆)‖≤CPMG _y\|U(y,z )\|≤ C_PM_G. Writing e=z−z⋆e=z-z , the fundamental theorem of calculus along the segment z⋆+sez +se, 0≤s≤10≤ s≤ 1, gives ‖G⁡(y,z)‖ \|G(y,z)\| ≤MG+K‖e‖, ≤ M_G+K\|e\|, ‖⁡(y,z)‖ \|U(y,z)\| ≤CPMG+2CPK‖e‖. ≤ C_PM_G+2C_PK\|e\|. Adding these inequalities proves ‖G⁡(y,z)‖+‖⁡(y,z)‖≤C⁡(1+‖z−z⋆‖).\|G(y,z)\|+\|U(y,z)\|≤ C(1+\|z-z \|). Finally, ‖Dρ‖≤1\|D_ρ\|≤ 1 for 0<ρ≤10<ρ≤ 1, so the exact update satisfies ‖Δzk‖=η‖DρG(Yk,zk)‖≤η⁡(MG+K‖ek‖)≤Cη(1+‖ek‖).\| z_k\|=η\|D_ρG(Y_k,z_k)\|≤η (M_G+K\|e_k\| )≤ Cη(1+\|e_k\|). These pointwise bounds hold with constants independent of η and ρ. We now combine them with the stationary fourth-moment bound to derive the expectation estimates stated above. Assumption A5 bounds the second derivatives of G⁡(y,z)G(y,z) uniformly over all y and z. The averaged field G¯ G inherits this bound, so Taylor’s theorem gives ‖RG¯,k(2)‖≤C‖ek‖2.\|R^(2)_ G,k\|≤ C\|e_k\|^2. Choose δ>0δ>0 such that the closed ball of radius δ about z⋆z lies inside the C3C^3 neighborhood. On this ball, ‖RG¯,k(3)‖≤C‖ek‖3\|R^(3)_ G,k\|≤ C\|e_k\|^3; outside it, the global Hessian bound still gives ‖RG¯,k(3)‖≤C‖ek‖2\|R^(3)_ G,k\|≤ C\|e_k\|^2. Consequently, ‖RG¯,k(3)‖≤C‖ek‖3+Cδ−2‖ek‖4≤Cη3/2.E\|R^(3)_ G,k\|≤ CE\|e_k\|^3+Cδ^-2E\|e_k\|^4≤ Cη^3/2. For each realized value of Yk+1Y_k+1, Taylor’s formula along the segment from zkz_k to zk+Δzkz_k+ z_k gives the exact remainder R,k=∫01(1−s)Dz2(Yk+1,zk+sΔzk)[Δzk,Δzk]s.R_U,k= _0^1(1-s)D_z^2U(Y_k+1,z_k+s z_k)[ z_k, z_k]\,ds. Since the second derivative of U is uniformly bounded, ‖R,k‖≤C‖Δzk‖2\|R_U,k\|≤ C\| z_k\|^2 for every realization. Consequently, the preceding linear-growth estimate and the fourth-moment assumption imply ‖Δzk‖2 \| z_k\|^2 ≤Cη2(1+‖ek‖)2 ≤ Cη^2E(1+\|e_k\|)^2 ≤Cη2(1+‖ek‖2)≤Cη2(1+η)≤Cη2. ≤ Cη^2(1+E\|e_k\|^2)≤ Cη^2(1+η)≤ Cη^2. The last inequality uses sufficiently small η, and all constants are uniform in ρ. Because the Taylor bound holds for every realization, dependence between Yk+1Y_k+1 and the iterates does not require an independence argument. For the frozen-coefficient error, add and subtract Dz(Yk+1,zk)DρG(Yk,z⋆)D_zU(Y_k+1,z_k)D_ρG(Y_k,z ). The bounded first derivative of G, the bounded first and second derivatives of U, and the bounded values G⁡(y,z⋆)G(y,z ) imply ‖Δ,k‖≤C‖ek‖,‖Δ,k‖≤Cη.\| _U,k\|≤ C\|e_k\|, \| _U,k\|≤ C η. For every s∈[0,1]s∈[0,1], the bounded first derivative of U, the linear growth of G, and ‖Dρ‖≤1\|D_ρ\|≤ 1 give ‖Dz(Yk+1,zk+sΔzk)‖‖DρG(Yk,zk)‖ \|D_zU(Y_k+1,z_k+s z_k)\|\,\|D_ρG(Y_k,z_k)\| ≤C⁡(1+‖ek‖). ≤ C(1+\|e_k\|). The derivative bound applies throughout the update segment, so no bounded-support assumption on that segment is needed. Integrating over an interval of length one and taking expectations yields ∫01‖Dz(Yk+1,zk+sΔzk)‖‖DρG(Yk,zk)‖s _0^1\|D_zU(Y_k+1,z_k+s z_k)\|\|D_ρG(Y_k,z_k)\|\,ds ≤C⁡(1+‖ek‖)≤C. ≤ C(1+E\|e_k\|)≤ C. By Hölder’s inequality and the stationary fourth-moment bound (8), ‖ek‖≤(‖ek‖4)1/4≤Cη.E\|e_k\|≤ (E\|e_k\|^4 )^1/4≤ C η. Finally, suppose that the second derivative of the chosen extension of λ is uniformly bounded along the segments joining θ⋆θ to stationary values of θ, as specified in the last part of the lemma. Taylor’s theorem gives ‖rλ(vk)‖≤C‖vk‖2.\|r_λ(v_k)\|≤ C\|v_k\|^2. Since vkv_k is a block of eke_k, Hölder’s inequality and the fourth-moment bound give ‖rλ(vk)‖ \|r_λ(v_k)\| ≤C‖vk‖2≤C(‖ek‖4)1/2≤Cη, ≤ CE\|v_k\|^2≤ C (E\|e_k\|^4 )^1/2≤ Cη, ‖rλ(vk)‖2 \|r_λ(v_k)\|^2 ≤C‖vk‖4≤C‖ek‖4≤Cη2. ≤ CE\|v_k\|^4≤ CE\|e_k\|^4≤ Cη^2. ∎ Appendix D Uniform stationary covariance and mean bias This section derives the stationary covariance and mean-bias expansions for the nonlinear recursion. Subsection D.1 studies the linearized, Poisson-corrected recursion and identifies the leading covariance coefficient through a Lyapunov equation. Subsection D.2 then proves that this coefficient approximates the stationary covariance of the original nonlinear recursion, with an error bound that remains uniform as ρ↓0ρ 0, establishing Proposition 2. Finally, Subsection D.4 first establishes the a priori O⁡(η)O(η) mean bound, then combines the covariance approximation with the stationary drift balance and Poisson-corrector identity to obtain the mean-bias expansion in Theorem 1. D.1 Linear Covariance Balance and Its Small-Step Limit Our goal is to identify the leading stationary covariance coefficient from the linear approximation near equilibrium. We first apply a Poisson correction to express the Markov noise as a martingale innovation plus a remainder. Linearizing the drift and evaluating the innovation at equilibrium then gives a recursion whose stationary covariance balances linear contraction against noise input. The innovation has zero conditional mean, so its cross terms with the current state vanish in this balance. Rescaling the covariance by η and taking η↓0η 0 at fixed ρ>0ρ>0 yields a Lyapunov equation for the leading coefficient. Subsection D.2 proves that the stationary covariance of the nonlinear recursion is ηΣρ+O⁡(η3/2)η _ρ+O(η^3/2), with an error constant independent of ρ for sufficiently small positive η and ρ. Let J=DG¯(z⋆)=(ABCD),Jρ=DρJ=(ABρCρD).J=D G(z )= pmatrixA&B\\ C&D pmatrix, J_ρ=D_ρJ= pmatrixA&B\\ ρ C&ρ D pmatrix. (89) Recall that ek=zk−z⋆e_k=z_k-z . Using the Poisson decomposition, we define the corrected error e^k=ek+ηDρ(Yk,zk). e_k=e_k+η D_ρU(Y_k,z_k). (90) Substitute the decomposition (79) into the recursion and use the definition of e^k+1 e_k+1. Before making any approximation, this gives e^k+1= e_k+1= e^k+ηDρG¯(zk)+ηDρζk+1(zk) e_k+η D_ρ G(z_k)+η D_ρ _k+1(z_k) +ηDρ⁡(Yk+1,zk+1)−⁡(Yk+1,zk). +η D_ρ \U(Y_k+1,z_k+1)-U(Y_k+1,z_k) \. (91) By the recursion (3) and the growth bound in (96), ‖zk+1−zk‖≤Cη(1+‖ek‖).\|z_k+1-z_k\|≤ Cη(1+\|e_k\|). The fourth-moment localization (8) therefore gives ‖zk+1−zk‖L2≤Cη\|z_k+1-z_k\|_L^2≤ Cη. Using the Lipschitz bound for U in (96), we obtain ‖ηDρ⁡(Yk+1,zk+1)−⁡(Yk+1,zk)‖L2≤Cη‖zk+1−zk‖L2≤Cη2. \|η D_ρ \U(Y_k+1,z_k+1)-U(Y_k+1,z_k) \ \|_L^2≤ Cη\|z_k+1-z_k\|_L^2≤ Cη^2. Thus the last line of (91) is OL2(η2)O_L^2(η^2), uniformly in sufficiently small positive η and ρ. Linearize G¯ G, replace eke_k by e^k e_k at leading order, and freeze the innovation at z⋆z . The leading recursion is then e^k+1=(I+ηJρ)e^k+ηDρζk+1(z⋆). e_k+1=(I+η J_ρ) e_k+η D_ρ _k+1(z ). (92) This equation defines the auxiliary linear recursion studied in Lemma 5. For fixed sufficiently small ρ>0ρ>0, the lemma shows that its stationary covariance, divided by η, converges as η↓0η 0 to Σρ _ρ, the unique solution of the associated Lyapunov equation. This identifies the leading covariance coefficient for the linear model. In forming this model, we have omitted the nonlinear drift remainder, the drift correction arising from replacing eke_k by e^k e_k, the change of the Poisson corrector over one update, and the variation of the innovation away from equilibrium. Subsection D.2 controls these contributions to establish the covariance approximation for the original nonlinear recursion uniformly as ρ↓0ρ 0. Lemma 5 (Small-step covariance of the linear recursion). Fix a sufficiently small ρ>0ρ>0 such that JρJ_ρ is Hurwitz. For sufficiently small η>0η>0, let Γ^η,ρlin ^lin_η,ρ be the stationary covariance of the linear recursion (92). Then limη↓0Γ^η,ρlinη=Σρ, _η 0 ^lin_η,ρη= _ρ, where Σρ _ρ is the unique solution of JρΣρ+ΣρJρ⊤+DρQDρ=0.J_ρ _ρ+ _ρJ_ρ +D_ρQD_ρ=0. (93) Proof. Write Aη,ρ=I+ηJρA_η,ρ=I+η J_ρ and ρ=DρQDρQ_ρ=D_ρQD_ρ. Since JρJ_ρ is Hurwitz, Aη,ρA_η,ρ is stable for all sufficiently small η. The frozen innovations are bounded because the Markov state space is finite. With the Markov chain in stationarity, stability therefore gives a stationary solution determined by past innovations and having finite covariance. By Lemma 1, the covariance of ζk+1(z⋆) _k+1(z ) is Q. The current state of the linear recursion is measurable with respect to the past, while ζk+1(z⋆) _k+1(z ) has conditional mean zero. Thus their cross-covariance is zero, and stationarity gives the exact balance Γ^η,ρlin=Aη,ρΓ^η,ρlinAη,ρ⊤+η2ρ. ^lin_η,ρ=A_η,ρ ^lin_η,ρA_η,ρ +η^2Q_ρ. Set Ση,ρ=Γ^η,ρlin/η _η,ρ= ^lin_η,ρ/η. Expanding the balance and dividing by η2η^2 yields JρΣη,ρ+Ση,ρJρ⊤+ηJρΣη,ρJρ⊤+ρ=0.J_ρ _η,ρ+ _η,ρJ_ρ +η J_ρ _η,ρJ_ρ +Q_ρ=0. Define the linear map ℒρ(X)=JρX+XJρ⊤L_ρ(X)=J_ρX+XJ_ρ . It is invertible: if ℒρ(X)=0L_ρ(X)=0, then etJρXetJρ⊤e^tJ_ρXe^tJ_ρ is constant in t and tends to zero as t→∞t→∞, so X=0X=0. Since this is a linear map on a finite-dimensional space, injectivity implies invertibility. Consequently, for sufficiently small η, the map X⟼ℒρ(X)+ηJρXJρ⊤X _ρ(X)+η J_ρXJ_ρ is also invertible, and its inverse converges to ℒρ−1L_ρ^-1. Hence Ση,ρ _η,ρ converges to the unique solution of (93). ∎ D.2 Proof of Proposition 2: Uniform Covariance Approximation Under the assumptions of Proposition 2, we prove the following estimate for the stationary covariance Γη,ρ _η,ρ of the process generated by the nonlinear recursion (3): Γη,ρ=ηΣρ+ℛη,ρΓ,‖ℛη,ρΓ‖≤Cη3/2. _η,ρ=η _ρ+R _η,ρ, \|R _η,ρ\|≤ Cη^3/2. (94) All constants below are independent of 0<η≤η00<η≤ _0 and 0<ρ≤ρ0≤10<ρ≤ _0≤ 1. The sharper slow and cross-moment bounds are not needed for this absolute covariance approximation. The covariance error will satisfy a Lyapunov equation whose right-hand side is the remainder in the stationary covariance balance. The next lemma bounds the solution in terms of the blocks of this remainder. The slow block enters the bound with a factor 1/ρ1/ρ. In the proof of Proposition 2, this factor is compensated by the O⁡(ρη3/2)O(ρη^3/2) bound on the slow remainder, yielding a covariance error uniform as ρ↓0ρ 0. Lemma 6 (Lyapunov inverse with scaled slow forcing). Let Lρ=(A−ρΛC−ρΛSρCρS),L_ρ= pmatrixA-ρ C&-ρ S\\ ρ C&ρ S pmatrix, where A and S are Hurwitz. For any symmetric matrix FforceF^force, write Fforce=(FuuforceFuvforce(Fuvforce)⊤Fvvforce),F^force= pmatrixF^force_u&F^force_uv\\ (F^force_uv) &F^force_v pmatrix, using the same (u,v)(u,v) block partition as LρL_ρ. After reducing ρ0>0 _0>0 if necessary, the symmetric solution X of LρX+XLρ⊤=−FforceL_ρX+XL_ρ =-F^force satisfies ‖X‖≤C∗(‖Fuuforce‖+‖Fuvforce‖+ρ−1‖Fvvforce‖),0<ρ≤ρ0.\|X\|≤ C_* (\|F^force_u\|+\|F^force_uv\|+ρ^-1\|F^force_v\| ), 0<ρ≤ _0. (95) Proof. Write X=(URR⊤V)X= pmatrixU&R\\ R &V pmatrix and divide the vv block equation by ρ. The resulting linear map acting on (U,R,V)(U,R,V) has coefficients continuous at ρ=0ρ=0. At that point it is (U,R,V)⟼(AU+UA⊤,AR,CR+R⊤C⊤+SV+VS⊤).(U,R,V) (AU+UA ,\;AR,\;CR+R C +SV+VS ). This map is block triangular. Its diagonal operators are the Lyapunov map associated with A, left multiplication by A, and the Lyapunov map associated with S. Each is invertible because A,SA,S are Hurwitz. Continuity of inversion in finite dimensions therefore gives a uniformly bounded inverse for sufficiently small ρ. Its right-hand side is −(Fuuforce,Fuvforce,Fvvforce/ρ)-(F^force_u,F^force_uv,F^force_v/ρ), proving (95). Stability of LρL_ρ also gives uniqueness for positive small ρ. ∎ Proof of Proposition 2. All expectations are taken under the joint invariant law. Throughout the proof, ‖X‖2=(‖X‖2)1/2\|X\|_2=(E\|X\|^2)^1/2 for a random vector X, and C denotes a constant independent of sufficiently small η and ρ. We first obtain an exact recursion whose leading noise is a martingale innovation. We then derive its stationary covariance balance and estimate the remainder in that balance. The decisive estimate is that the slow–slow remainder carries a factor ρ: this compensates for the factor 1/ρ1/ρ in Lemma 6. Finally, we transfer the result back to the original coordinates. Step 1: Express the recursion in corrected tangent coordinates. The tangent transformation exposes the factors of ρ in the slow dynamics, while the Poisson correction isolates noise with zero conditional mean. Define T=(I−Λ0I),Bρ=TDρ=(I−ρΛ0ρI),Lρ=TDρJT−1.T= pmatrixI&- \\ 0&I pmatrix, B_ρ=TD_ρ= pmatrixI&-ρ \\ 0&ρ I pmatrix, L_ρ=TD_ρJT^-1. For the fixed finite-state Doeblin kernel, the Poisson inverse is bounded on the space of μ-centered functions. Applying this inverse to the first and second derivatives of G~=G−G¯ G=G- G shows that the Poisson corrector U has bounded first and second derivatives, uniformly in the Markov state and z. Its values at z⋆z are bounded because the state space is finite. Together with the derivative bounds on G, this gives ‖G⁡(y,z)‖+‖⁡(y,z)‖ \|G(y,z)\|+\|U(y,z)\| ≤C⁡(1+‖z−z⋆‖), ≤ C(1+\|z-z \|), ‖⁡(y,z′)−⁡(y,z)‖ \|U(y,z )-U(y,z)\| ≤C‖z′−z‖, ≤ C\|z -z\|, ‖ζk+1(z)−ζk+1(z⋆)‖ \| _k+1(z)- _k+1(z )\| ≤C‖z−z⋆‖. ≤ C\|z-z \|. (96) For the adapted iterate, the innovation satisfies ⁡[ζk+1(zk)∣ℱk]=0E[ _k+1(z_k) _k]=0; this uses no independence between YkY_k and zkz_k. Set wk=Tek,w^k=wk+ηBρ(Yk,zk),RG¯,k(2)=G¯(zk)−Jek.w_k=Te_k, w_k=w_k+η B_ρU(Y_k,z_k), R^(2)_ G,k= G(z_k)-Je_k. Applying T to the exact Poisson-corrected recursion (91) gives w^k+1=w^k+η⁡ak+nk+1+k+1,ak=Lρw^k,nk+1=Bρζk+1(zk), w_k+1= w_k+η\a_k+n_k+1+ r_k+1\, a_k=L_ρ w_k, n_k+1=B_ρ _k+1(z_k), (97) where k+1= r_k+1= BρRG¯,k(2)−ηLρBρ(Yk,zk) B_ρR^(2)_ G,k-η L_ρB_ρU(Y_k,z_k) +Bρ⁡(Yk+1,zk+1)−⁡(Yk+1,zk). +B_ρ\U(Y_k+1,z_k+1)-U(Y_k+1,z_k)\. (98) The three terms account for the nonlinear drift, the correction to the linear drift, and the change of the Poisson corrector over one update. In particular, k+1 r_k+1 can depend on Yk+1Y_k+1 and need not be ℱkF_k-measurable. Step 2: Derive the stationary covariance balance. The corrected process is stationary because w^k w_k is a fixed function of the stationary pair (Yk,zk)(Y_k,z_k). Let qk=w^k−w^k,P^=⁡[qkqk⊤].q_k= w_k-E w_k, P=E[q_kq_k ]. Stationarity makes the mean increment in (97) zero. Moreover, qkq_k and aka_k are ℱkF_k-measurable, whereas nk+1n_k+1 has conditional mean zero. Hence ⁡[qknk+1⊤]=0,⁡[aknk+1⊤]=0.E[q_kn_k+1 ]=0, [a_kn_k+1 ]=0. To identify the leading noise covariance, recall that the Y-marginal of the joint invariant law is μ. Thus (Yk,Yk+1)(Y_k,Y_k+1) has law μ⁡(dy)P(y,dy′)μ(dy)P(y,dy ), and ⁡[ζk+1(z⋆)ζk+1(z⋆)⊤]=Q.E[ _k+1(z ) _k+1(z ) ]=Q. Define the error from evaluating the innovation at zkz_k instead of z⋆z by ΔQη,ρ=⁡[ζk+1(zk)ζk+1(zk)⊤]−Q. Q_η,ρ=E[ _k+1(z_k) _k+1(z_k) ]-Q. Expanding ⁡[qk+1qk+1⊤−qkqk⊤]=0E[q_k+1q_k+1 -q_kq_k ]=0 and dividing by η now yields the exact equation LρP^+P^Lρ⊤+ηBρQBρ⊤+ℛη,ρ=0,L_ρ P+ PL_ρ +η B_ρQB_ρ +R_η,ρ=0, (99) where, suppressing time indices on a,n,a,n, r, ℛη,ρ= _η,ρ= ⁡[qk⊤+qk⊤]+ηBρΔQη,ρBρ⊤ [q_k r + rq_k ]+η B_ρ Q_η,ρB_ρ +η[(a+)(a+)⊤+n⊤+n⊤]. + [(a+ r)(a+ r) +n r + rn ]. (100) The three groups are the state–remainder cross terms, the error in the innovation covariance, and the remaining products of increments. The terms n⊤n r and n⊤ rn must be retained because the remainder need not be predictable. It remains to bound these three groups with their slow-coordinate factors intact. Step 3: Bound the remainder and retain its slow factor. The raw fourth-moment assumption gives ‖ek‖2≤Cη,‖ek‖2‖2≤Cη.\|e_k\|_2≤ C η, \|\|e_k\|^2 \|_2≤ Cη. (101) Bounded second derivatives of G¯ G and G¯(z⋆)=0 G(z )=0 imply ‖RG¯,k(2)‖≤C‖ek‖2\|R^(2)_ G,k\|≤ C\|e_k\|^2. Also, the update and (96) give ‖⁡(Yk+1,zk+1)−⁡(Yk+1,zk)‖≤Cη(1+‖ek‖).\|U(Y_k+1,z_k+1)-U(Y_k+1,z_k)\|≤ Cη(1+\|e_k\|). The upper block rows of BρB_ρ, LρL_ρ, and LρBρL_ρB_ρ are uniformly bounded, and their lower block rows are O⁡(ρ)O(ρ). Applying these facts to (97) and (98) yields ‖w^k‖2 \| w_k\|_2 ≤Cη, ≤ C η, ‖k+1u‖2 \| r_k+1^u\|_2 ≤Cη, ≤ Cη, ‖k+1v‖2 \| r_k+1^v\|_2 ≤Cρη, ≤ Cρη, ‖aku‖2 \|a_k^u\|_2 ≤Cη, ≤ C η, ‖akv‖2 \|a_k^v\|_2 ≤Cρη, ≤ Cρ η, ‖nk+1u‖2 \|n_k+1^u\|_2 ≤C, ≤ C, ‖nk+1v‖2 \|n_k+1^v\|_2 ≤Cρ. ≤ Cρ. (102) In particular, ‖qku‖2,‖qkv‖2≤Cη\|q_k^u\|_2,\|q_k^v\|_2≤ C η. For the innovation covariance, the Lipschitz bound in (96) and boundedness of the frozen innovation give ‖ΔQη,ρ‖≤C(‖ek‖+‖ek‖2)≤Cη.\| Q_η,ρ\|≤ CE(\|e_k\|+\|e_k\|^2)≤ C η. (103) We can now estimate each group in (100) using ‖⁡[XY⊤]‖≤‖X‖2‖Y‖2\|E[XY ]\|≤\|X\|_2\|Y\|_2. The first group has uu and uvuv blocks of order η3/2η^3/2. Its vv block is of order ρη3/2ρη^3/2 because each product contains v r^v. The second group has blocks of orders η3/2η^3/2, ρη3/2ρη^3/2, and ρ2η3/2ρ^2η^3/2, respectively. The last group has blocks of orders η2η^2, ρη2ρη^2, and ρ2η2ρ^2η^2. Consequently, for η,ρ≤1η,ρ≤ 1, ‖(ℛη,ρ)uu‖+‖(ℛη,ρ)uv‖ \|(R_η,ρ)_u\|+\|(R_η,ρ)_uv\| ≤Cη3/2, ≤ Cη^3/2, (104) ‖(ℛη,ρ)vv‖ \|(R_η,ρ)_v\| ≤Cρη3/2. ≤ Cρη^3/2. The slow factor comes from the update coefficients. Only the total raw-error moments in (101) were used; no sharper slow-coordinate or cross-moment bound is needed. Step 4: Solve the covariance equation and return to the original variables. The leading covariance in tangent coordinates, Σ~ρ=TΣρT⊤ _ρ=T _ρT , satisfies LρΣ~ρ+Σ~ρLρ⊤+BρQBρ⊤=0.L_ρ _ρ+ _ρL_ρ +B_ρQB_ρ =0. Subtracting η times this equation from (99) shows that P^−ηΣ~ρ P-η _ρ solves a Lyapunov equation with forcing ℛη,ρR_η,ρ. Lemma 6 and (104) therefore give ‖P^−ηΣ~ρ‖ \| P-η _ρ\| ≤C∗(‖(ℛη,ρ)uu‖+‖(ℛη,ρ)uv‖+ρ−1‖(ℛη,ρ)vv‖) ≤ C_* (\|(R_η,ρ)_u\|+\|(R_η,ρ)_uv\|+ρ^-1\|(R_η,ρ)_v\| ) (105) ≤Cη3/2. ≤ Cη^3/2. This is where the retained factor ρ gives a bound uniform as ρ↓0ρ 0. Finally, the Poisson correction changes the covariance by at most the same order. Indeed, write hk=Bρ(Yk,zk),w^k=wk+ηhk.h_k=B_ρU(Y_k,z_k), w_k=w_k+η h_k. The preceding bounds give ‖wk‖2≤Cη\|w_k\|_2≤ C η and ‖hk‖2≤C\|h_k\|_2≤ C. Expanding the centered covariances yields ‖Cov⁡(w^k)−Cov⁡(wk)‖ \|Cov( w_k)-Cov(w_k)\| ≤2η‖wk‖2‖hk‖2+η2‖hk‖22≤Cη3/2. ≤ 2η\|w_k\|_2\|h_k\|_2+η^2\|h_k\|_2^2≤ Cη^3/2. Since Cov⁡(wk)=TΓη,ρT⊤Cov(w_k)=T _η,ρT and T,T−1T,T^-1 are fixed, combining this estimate with (105) gives ‖Γη,ρ−ηΣρ‖≤Cη3/2.\| _η,ρ-η _ρ\|≤ Cη^3/2. This proves (94) and Proposition 2. The argument uses neither a stationary mean-bias expansion nor an assumed covariance approximation. ∎ Theorem 4 sharpens the remainder to O⁡(η2)O(η^2) under its additional additive-i.i.d. hypotheses. An absolute O⁡(η3/2)O(η^3/2) covariance remainder does not establish the sharper O⁡(ε)O( ) slow and cross-covariance bounds on all step-size paths. Those bounds are separate conditions in Assumption A7 for the supplementary fluctuation statements. D.3 Proof of Lemma 2 and Theorem 1 Proof of the lemma and theorem. Stationarity and invertibility of DρD_ρ imply G(Yk,zk)=0EG(Y_k,z_k)=0. Exogeneity and the Poisson equation give the exact moving-corrector identity G~(Yk,zk)=⁡[⁡(Yk+1,zk+1)−⁡(Yk+1,zk)].E G(Y_k,z_k)=E[U(Y_k+1,z_k+1)-U(Y_k+1,z_k)]. The derivative and moment bounds first bound this expression by CηCη. The first-order drift expansion then gives Jmη,ρ=O⁡(η)Jm_η,ρ=O(η), proving the lemma because detJ=detAdetS≠0 J= A S≠ 0. This inverts the fixed matrix J, avoiding the 1/ρ1/ρ loss from Jρ−1J_ρ^-1. Taylor expansion of the corrector, followed by freezing its coefficient at z⋆z , gives G~(Yk,zk)=ηrρ+O⁡(η3/2).E G(Y_k,z_k)=η r_ρ+O(η^3/2). By Proposition 2 and the mean bound, G¯(zk)=Jmη,ρ+η(Σρ)+O⁡(η3/2).E G(z_k)=Jm_η,ρ+ ( _ρ)+O(η^3/2). (106) The mean-square contribution is O⁡(η2)O(η^2). Combining the two expressions with the zero expected drift proves the theorem. Appendix D.4 supplies the full identities and remainder bounds. ∎ D.4 Stationary Mean Bias We now use the covariance approximation from Proposition 2 to determine the stationary mean displacement mη,ρm_η,ρ, proving Theorem 1. The proof starts from the fact that the expected update vanishes under the invariant law. The Poisson equation relates this stationary balance to the dependence between the iterate and the Markov input. A first-order Taylor expansion then gives the preliminary bound mη,ρ=O⁡(η)m_η,ρ=O(η). With this bound in hand, we expand the averaged drift to second order, so that its expectation is expressed through the mean displacement and the stationary covariance. Substituting the covariance approximation and combining the result with the Poisson identity yields the bias expansion, including both the curvature contribution and the contribution from Markov dependence. Recall from Appendix D that ek=zk−z⋆e_k=z_k-z . Expectations without a subscript are under the joint invariant law πη,ρ _η,ρ of (Yk,zk)(Y_k,z_k). Set mη,ρ=ek,Γη,ρ=Cov⁡(ek).m_η,ρ=Ee_k, _η,ρ=Cov(e_k). The Y-marginal of πη,ρ _η,ρ is μ; denote its z-marginal, obtained by integrating out YkY_k, by πη,ρzπ^z_η,ρ. Recall that G¯ G and G~ G are defined in (71). Although G~(⋅,z) G(·,z) has zero μ-mean for each fixed z, G~(Yk,zk)E G(Y_k,z_k) need not vanish under the joint invariant law because YkY_k and zkz_k may be dependent. Equivalently, G¯(zk) G(z_k) is generally not ⁡[G⁡(Yk,zk)∣zk]E[G(Y_k,z_k) z_k]. The Jacobian J=DG¯(z⋆)J=D G(z ), defined in (89), is nonsingular because detJ=detAdetS≠0 J= A\, S≠ 0. Theorem 1 (Fixed-ratio bias; restatement). Under , , , and , for 0<ρ≤ρ00<ρ≤ _0, mη,ρ=ηb(ρ)+O(η3/2),b(ρ)=−J−1(Σρ)+rρ. m_η,ρ=η b(ρ)+O(η^3/2), b(ρ)=-J^-1\C( _ρ)+r_ρ\. (107) with a remainder uniform in ρ. Proof of Theorem 1. Proof. For M∈ℝ(dx+dθ)×(dx+dθ)M ^(d_x+d_θ)×(d_x+d_θ), define (M)i=12Tr(D2G¯i(z⋆)M),i=1,…,dx+dθ.C(M)_i= 12Tr\! (D^2 G_i(z )M ), i=1,…,d_x+d_θ. Define the direct Markov response by rρ=∑y,y′∈μ⁡(y)P(y,y′)Dz(y′,z⋆)DρG(y,z⋆).r_ρ= _y,y ∈ Yμ(y)P(y,y )D_zU(y ,z )D_ρG(y,z ). (108) The weights μ⁡(y)P(y,y′)μ(y)P(y,y ) give the stationary one-step transition law, and U is the Poisson corrector from (75). Stationarity of zk+1=zk+ηDρG(Yk,zk)z_k+1=z_k+η D_ρG(Y_k,z_k), and invertibility of DρD_ρ for ρ>0ρ>0, give G(Yk,zk)=0,G¯(zk)=−G~(Yk,zk).EG(Y_k,z_k)=0, G(z_k)=-E G(Y_k,z_k). (109) The Poisson equation, exogeneity, and stationarity yield G~(Yk,zk) G(Y_k,z_k) =⁡(Yk,zk)−⁡(Yk+1,zk) =E\U(Y_k,z_k)-U(Y_k+1,z_k)\ =⁡(Yk+1,zk+1)−⁡(Yk+1,zk). =E\U(Y_k+1,z_k+1)-U(Y_k+1,z_k)\. (110) Taylor expansion across the update zk+1−zkz_k+1-z_k, followed by the remainder bounds in (87), gives G~(Yk,zk)=ηrρ+O⁡(η3/2).E G(Y_k,z_k)=η r_ρ+O(η^3/2). (111) Consequently, (109) implies G¯(zk)=−ηrρ+O⁡(η3/2).E G(z_k)=-η r_ρ+O(η^3/2). (112) The localization moment bounds give ‖ek‖2=O⁡(η)E\|e_k\|^2=O(η). The moving-corrector identity and the mean-value bound in (87) give G~(Yk,zk)=O⁡(η)E G(Y_k,z_k)=O(η). Also, ‖RG¯,k(2)‖=O⁡(η)E\|R^(2)_ G,k\|=O(η). Taking expectations in G¯(zk)=Jek+RG¯,k(2) G(z_k)=Je_k+R^(2)_ G,k and using (109) therefore yields Jmη,ρ=−G~(Yk,zk)−RG¯,k(2)=O⁡(η).Jm_η,ρ=-E G(Y_k,z_k)-ER^(2)_ G,k=O(η). Since J is nonsingular, mη,ρ=O⁡(η)m_η,ρ=O(η). By Proposition 2, Γη,ρ=ηΣρ+O⁡(η3/2), _η,ρ=η _ρ+O(η^3/2), (113) where Σρ _ρ solves (93). Taylor expansion of the averaged drift gives G¯(zk)=Jek+12D2G¯(z⋆)[ek,ek]+RG¯,k(3). G(z_k)=Je_k+ 12D^2 G(z )[e_k,e_k]+R^(3)_ G,k. Using ⁡[ekek⊤]=Γη,ρ+mη,ρmη,ρ⊤E[e_ke_k ]= _η,ρ+m_η,ρm_η,ρ , the covariance expansion, mη,ρ=O⁡(η)m_η,ρ=O(η), and the third-order remainder bound yields G¯(zk)=Jmη,ρ+η(Σρ)+O⁡(η3/2).E G(z_k)=Jm_η,ρ+ ( _ρ)+O(η^3/2). (114) Equating (112) and (114), then multiplying by J−1J^-1, proves the stated expansion. All remainder bounds above are uniform for 0<ρ≤ρ0≤10<ρ≤ _0≤ 1: the Taylor estimates and Proposition 2 are uniform, ‖Dρ‖≤1\|D_ρ\|≤ 1, and both the Poisson solution operator and J−1J^-1 are independent of ρ. The inverse of DρD_ρ is used only to deduce the exact identity G(Yk,zk)=0EG(Y_k,z_k)=0, so it introduces no factor 1/ρ1/ρ into the error estimates. ∎ Appendix E The singular ratio limit and path-aware extrapolation Appendix D.4 established the stationary mean expansion mη,ρ=ηb(ρ)+O⁡(η3/2),m_η,ρ=η b(ρ)+O(η^3/2), with a remainder uniform for 0<ρ≤ρ00<ρ≤ _0. We now expand b⁡(ρ)b(ρ) as ρ↓0ρ 0 to obtain the expansion stated in Theorem 2. The coefficient vectors b0,b1,b2b_0,b_1,b_2, defined in (125)–(127), multiply η, ε , and ε2/η ^2/η, respectively, and are independent of the step sizes. To derive these coefficients, we expand the covariance coefficient Σρ _ρ and the Markov response rρr_ρ. For the covariance calculation, we use the tangent coordinates u=x−x⋆−Λ⁡(θ−θ⋆),v=θ−θ⋆,u=x-x - (θ-θ ), v=θ-θ , where u is the fast tracking error and v is the slow displacement. In these coordinates, the leading covariance coefficient is Σ~ρ=TΣρT⊤ _ρ=T _ρT . We partition it as Σ~ρ=(Σ~ρuuΣ~ρuv(Σ~ρuv)⊤Σ~ρvv). _ρ= pmatrix _ρ^u& _ρ^uv\\ ( _ρ^uv) & _ρ^v pmatrix. The uu block describes fluctuations within the fast coordinates, the vv block describes fluctuations within the slow coordinates, and the uvuv block describes their cross-covariance. For ρ>0ρ>0, define the rescaled matrices U⁡(ρ)=Σ~ρuu,R⁡(ρ)=ρ−1Σ~ρuv,V⁡(ρ)=ρ−1Σ~ρvv.U(ρ)= _ρ^u, R(ρ)=ρ^-1 _ρ^uv, V(ρ)=ρ^-1 _ρ^v. Thus we leave the fast block unchanged and divide the cross and slow blocks by ρ, writing Σ~ρ=(U⁡(ρ)ρR(ρ)ρR(ρ)⊤ρV(ρ)). _ρ= pmatrixU(ρ)&ρ R(ρ)\\ ρ R(ρ) &ρ V(ρ) pmatrix. Substituting this representation into the transformed Lyapunov equation, we divide its uvuv equation by ρ and its vv equation by ρ2ρ^2. The resulting equations remain informative at ρ=0ρ=0 and determine finite limits for U⁡(ρ)U(ρ), R⁡(ρ)R(ρ), and V⁡(ρ)V(ρ). This allows us to compute the covariance expansion needed for b0,b1,b2b_0,b_1,b_2. E.1 Proof of Theorem 2 Proof. Let Bρ=TDρB_ρ=TD_ρ and Σ~ρ=TΣρT⊤ _ρ=T _ρT . The transformed Lyapunov equation is LρΣ~ρ+Σ~ρLρ⊤+BρQBρ⊤=0.L_ρ _ρ+ _ρL_ρ +B_ρQB_ρ =0. (115) Partition Q according to (x,θ)(x,θ) and seek Σ~ρ=(U⁡(ρ)ρR(ρ)ρR(ρ)⊤ρV(ρ)). _ρ= pmatrixU(ρ)&ρ R(ρ)\\ ρ R(ρ) &ρ V(ρ) pmatrix. (116) After dividing the cross and slow equations by ρ and ρ2ρ^2, respectively, one obtains the regular system 0= 0= (A−ρΛC)U+U(A−ρΛC)⊤ (A-ρ C)U+U(A-ρ C) −ρ2ΛSR⊤−ρ2RS⊤Λ⊤ -ρ^2 SR -ρ^2RS +Qxx−ρ⁡(ΛQθx+QxθΛ⊤)+ρ2ΛQθθΛ⊤, +Q_x-ρ( Q_θ x+Q_xθ )+ρ^2 Q_θ , (117) 0= 0= (A−ρΛC)R−ρΛSV+UC⊤ (A-ρ C)R-ρ SV+UC +ρRS⊤+Qxθ−ρΛQθθ, +ρ RS +Q_xθ-ρ Q_θ, (118) 0= 0= CR+SV+R⊤C⊤+VS⊤+Qθθ. CR+SV+R C +VS +Q_θ. (119) At zero the system is triangular: AU0+U0A⊤+Qxx AU_0+U_0A +Q_x =0, =0, (120) AR0+U0C⊤+Qxθ AR_0+U_0C +Q_xθ =0, =0, (121) SV0+V0S⊤+CR0+R0⊤C⊤+Qθθ SV_0+V_0S +CR_0+R_0 C +Q_θ =0. =0. (122) Hurwitz stability makes this finite-dimensional linear system invertible. Its coefficients are polynomial in ρ, so its inverse and solution are analytic near zero. Successive coefficient matching therefore gives Σρ=Σ0+ρΣ1+ρ2Σ2+O⁡(ρ3),Σ0=diag⁡(U0,0). _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3), _0=diag(U_0,0). (123) The corrector response is affine: rρ=r0+ρr1,r0=[Dx(Y1,z⋆)ξ0x],r1=⁡[Dθ(Y1,z⋆)ξ0θ]. gatheredr_ρ=r_0+ρ r_1, r_0=E[D_xU(Y_1,z ) _0^x],\\ r_1=E[D_θU(Y_1,z ) _0^θ]. gathered (124) Consequently, b0 b_0 =−J−1⁡(Σ0)+r0, =-J^-1\C( _0)+r_0\, (125) b1 b_1 =−J−1⁡(Σ1)+r1, =-J^-1\C( _1)+r_1\, (126) b2 b_2 =−J−1(Σ2). =-J^-1C( _2). (127) Insert these expressions into Theorem 1 and use ηρ=εηρ= . Appendix E details the coefficient equations. ∎ E.2 Overview of the Coefficient Expansion The starting point is the bias coefficient obtained in Appendix D.4: b⁡(ρ)=−J−1⁡(Σρ)+rρ.b(ρ)=-J^-1\C( _ρ)+r_ρ\. Here J=DG¯(z⋆)J=D G(z ) and the curvature map C are independent of ρ. Thus, to expand b⁡(ρ)b(ρ), we need to determine how the leading covariance coefficient Σρ _ρ and the Markov response rρr_ρ depend on the step-size ratio. The dependence of rρr_ρ is explicit. Its defining expression (108) contains Dρ=diag⁡(I,ρI)D_ρ=diag(I,ρ I) linearly, so rρ=r0+ρr1,r_ρ=r_0+ρ r_1, where r0r_0 and r1r_1 are defined in (124). This identity is exact; there are no higher-order powers of ρ in the direct Markov response. The covariance coefficient requires a separate calculation. As shown below, the rescaled covariance equations (134)–(136) form a finite-dimensional linear system that is invertible at ρ=0ρ=0. Its coefficients depend polynomially on ρ, so its solution admits a Taylor expansion near zero. Returning from the rescaled tangent coordinates to the original coordinates gives Σρ=Σ0+ρΣ1+ρ2Σ2+O⁡(ρ3). _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3). The matrices Σ0,Σ1,Σ2 _0, _1, _2 are independent of the step sizes. The following subsections derive the equations that determine them. Substituting these expansions into the formula for b⁡(ρ)b(ρ) and using the linearity of C gives b⁡(ρ)=b0+ρb1+ρ2b2+O⁡(ρ3),b(ρ)=b_0+ρ b_1+ρ^2b_2+O(ρ^3), with b0 b_0 =−J−1⁡(Σ0)+r0, =-J^-1\C( _0)+r_0\, b1 b_1 =−J−1⁡(Σ1)+r1, =-J^-1\C( _1)+r_1\, b2 b_2 =−J−1(Σ2). =-J^-1C( _2). These are the coefficient vectors defined in (125)–(127). The first two coefficients contain both covariance and direct Markov-response terms. The quadratic coefficient b2b_2 contains only the covariance contribution because rρr_ρ has no quadratic term. To express the result in the original step sizes, multiply by η and use ρ=ε/ηρ= /η. The identities ηρ=εηρ= and ηρ2=ε2/ηρ^2= ^2/η yield ηb(ρ)=ηb0+εb1+ε2ηb2+O⁡(ε3η2).η b(ρ)=η b_0+ b_1+ ^2ηb_2+O\! ( ^3η^2 ). (128) Theorem 1 has already established mη,ρ=ηb(ρ)+O⁡(η3/2)m_η,ρ=η b(ρ)+O(η^3/2) with a remainder uniform in ρ. Combining the two expansions therefore gives mη,ρ m_η,ρ =ηb0+εb1+ε2ηb2 =η b_0+ b_1+ ^2ηb_2 +O⁡(η3/2+ε3η2), +O\! (η^3/2+ ^3η^2 ), which is the expansion in Theorem 2. The two remainder terms have different origins: O⁡(η3/2)O(η^3/2) comes from the stationary mean approximation in Theorem 1, whereas O⁡(ε3/η2)O( ^3/η^2) comes from truncating the expansion of b⁡(ρ)b(ρ) after its quadratic term. E.3 Rescaled Covariance Equations At ρ=0ρ=0 both slow drift and slow noise vanish. We therefore pass to fast-manifold coordinates and divide out the known powers of ρ before taking the limit. With J, Λ , and S as in (89), (61), and (62), define T=(I−Λ0I),(uv)=T(x−x⋆θ−θ⋆).T= pmatrixI&- \\ 0&I pmatrix, pmatrixu\\ v pmatrix=T pmatrixx-x \\ θ-θ pmatrix. (129) so that (u,v)(u,v) are the tangent coordinates (44). Let Dρ=diag⁡(I,ρI)D_ρ=diag(I,ρ I). In the new coordinates, the linear drift and noise scaling are Lρ=TDρJT−1 L_ρ=TD_ρJT^-1 =(A−ρΛC−ρΛSρCρS), = pmatrixA-ρ C&-ρ S\\ ρ C&ρ S pmatrix, (130) Bρ=TDρ B_ρ=TD_ρ =(I−ρΛ0ρI). = pmatrixI&-ρ \\ 0&ρ I pmatrix. (131) The leading covariance coefficient from Appendix D becomes Σ~ρ=TΣρT⊤ _ρ=T _ρT and solves LρΣ~ρ+Σ~ρLρ⊤+BρQBρ⊤=0.L_ρ _ρ+ _ρL_ρ +B_ρQB_ρ =0. (132) Partition the long-run noise covariance as Q=(QxxQxθQθxQθθ).Q= pmatrixQ_x&Q_xθ\\ Q_θ x&Q_θ pmatrix. The uu, uvuv, and vv blocks of BρQBρ⊤B_ρQB_ρ are respectively O⁡(1)O(1), O⁡(ρ)O(ρ), and O⁡(ρ2)O(ρ^2). Together with the different fast and slow relaxation rates, these noise scales motivate the following rescaling of the covariance coefficient. For each ρ>0ρ>0, write Σ~ρ=(U⁡(ρ)ρR(ρ)ρR(ρ)⊤ρV(ρ)). _ρ= pmatrixU(ρ)&ρ R(ρ)\\ ρ R(ρ) &ρ V(ρ) pmatrix. (133) For positive ρ, this is an exact change of variables. The assertion that U⁡(ρ)U(ρ), R⁡(ρ)R(ρ), and V⁡(ρ)V(ρ) remain bounded as ρ↓0ρ 0 is justified by the invertibility argument in Appendix E.4. Substitute this representation into (132). Leave the uu equation unchanged, divide the uvuv equation by ρ, and divide the vv equation by ρ2ρ^2. For ρ>0ρ>0, these operations give the equivalent system: 0= 0= (A−ρΛC)U+U(A−ρΛC)⊤−ρ2ΛSR⊤−ρ2RS⊤Λ⊤ (A-ρ C)U+U(A-ρ C) -ρ^2 SR -ρ^2RS +Qxx−ρ⁡(ΛQθx+QxθΛ⊤)+ρ2ΛQθθΛ⊤, +Q_x-ρ( Q_θ x+Q_xθ )+ρ^2 Q_θ , (134) 0= 0= (A−ρΛC)R−ρΛSV+UC⊤+ρRS⊤+Qxθ−ρΛQθθ, (A-ρ C)R-ρ SV+UC +ρ RS +Q_xθ-ρ Q_θ, (135) 0= 0= CR+SV+R⊤C⊤+VS⊤+Qθθ. CR+SV+R C +VS +Q_θ. (136) We can now set ρ=0ρ=0. The resulting equations form a triangular chain: solve first for the fast covariance U0U_0, then the cross-covariance R0R_0, and finally the slow covariance V0V_0: AU0+U0A⊤+Qxx AU_0+U_0A +Q_x =0, =0, (137) AR0+U0C⊤+Qxθ AR_0+U_0C +Q_xθ =0, =0, (138) SV0+V0S⊤+CR0+R0⊤C⊤+Qθθ SV_0+V_0S +CR_0+R_0 C +Q_θ =0. =0. (139) Stability of A and S gives uniqueness at every step. Eliminating R0R_0 makes the slow-noise interpretation explicit. The middle equation gives R0=−A−1(U0C⊤+Qxθ)R_0=-A^-1(U_0C +Q_xθ). Substituting this into the last equation gives SV0+V0S⊤+Qred=0,SV_0+V_0S +Q_red=0, with Qred=(−CA−1I)Q(−CA−1I)⊤.Q_red= pmatrix-CA^-1&I pmatrixQ pmatrix-CA^-1&I pmatrix^\! . (140) The slow system receives both direct noise and fast noise transmitted by the equilibrium response −CA−1-CA^-1. Multiplication by η gives the leading covariance contributions: [ηΣ~ρ]uu=O⁡(η),[ηΣ~ρ]uv=O⁡(ε),[ηΣ~ρ]vv=O⁡(ε).[η _ρ]_u=O(η), [η _ρ]_uv=O( ), [η _ρ]_v=O( ). (141) These are orders for the leading covariance approximation. Proposition 2 bounds the difference between the actual covariance and this approximation by O⁡(η3/2)O(η^3/2). Consequently, under Assumptions A1–A5, it gives ‖Cov⁡(v)‖+|Cov⁡(u,v)|=O⁡(ε+η3/2).\|Cov(v)\|+\|Cov(u,v)\|=O( +η^3/2). When ρ≲ηρ η, the approximation error can be comparable to or larger than the leading O⁡(ε)O( ) terms. The stronger bounds ‖Cov⁡(v)‖=O⁡(ε),‖Cov⁡(u,v)‖=O⁡(ε),\|Cov(v)\|=O( ), \|Cov(u,v)\|=O( ), uniformly over all sufficiently small positive η and ρ, are additional conditions collected in Assumption A7. Extrapolation needs terms through ρ2ρ^2. U=U0+ρU1+ρ2U2+O⁡(ρ3),R=R0+ρR1+O⁡(ρ2),V=V0+ρV1+O⁡(ρ2).U=U_0+ρ U_1+ρ^2U_2+O(ρ^3), R=R_0+ρ R_1+O(ρ^2), V=V_0+ρ V_1+O(ρ^2). Differentiate (134)–(136) to obtain the coefficients by successive Lyapunov or Sylvester solves. Transforming back gives Σρ=Σ0+ρΣ1+ρ2Σ2+O⁡(ρ3). _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3). (142) The coefficients follow from Σρ=T−1Σ~ρT−⁣⊤ _ρ=T^-1 _ρT^- ; in particular, Σ0=diag⁡(U0,0) _0=diag(U_0,0). E.4 Derivation of the Covariance Expansion We now justify the expansion Σρ=Σ0+ρΣ1+ρ2Σ2+O⁡(ρ3) _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3) and derive equations for its coefficients. The argument first establishes that the rescaled covariance equations have a solution analytic in ρ near zero. We then compute its coefficients and transform back to the original coordinates. Existence of the Expansion. Let mS^m denote the space of real symmetric m×m× m matrices, and set =dx×ℝdx×dθ×dθ.X=S^d_x×R^d_x× d_θ×S^d_θ. The rescaled equations (134)–(136) can be written as ⁡(ρ)(U,R,V)+ℱ⁡(ρ)=0. A(ρ)(U,R,V)+ F(ρ)=0. Here ⁡(ρ):→ A(ρ):X contains all terms linear in U,R,VU,R,V, while ℱ⁡(ρ) F(ρ) consists of the terms involving the blocks of Q. In fixed bases, both have polynomial entries in ρ. To prove invertibility at zero, consider the homogeneous system (0)(U,R,V)=0 A(0)(U,R,V)=0: AU+UA⊤ AU+UA =0, =0, AR+UC⊤ AR+UC =0, =0, SV+VS⊤+CR+R⊤C⊤ SV+VS +CR+R C =0. =0. Since A is Hurwitz, its Lyapunov equation in the first line has only the solution U=0U=0. The second line then gives R=0R=0 because A is invertible. Finally, Hurwitz stability of S gives V=0V=0 from the third line. Thus ⁡(0) A(0) has trivial kernel. It maps the finite-dimensional space X into itself and is therefore invertible. Its determinant remains nonzero on a neighborhood of zero. Because its entries are polynomial in ρ, the inverse (ρ)−1 A(ρ)^-1 is analytic on that neighborhood. Consequently, the solution is analytic and is given by (U⁡(ρ),R⁡(ρ),V⁡(ρ))=−(ρ)−1ℱ(ρ).(U(ρ),R(ρ),V(ρ))=- A(ρ)^-1 F(ρ). For sufficiently small positive ρ, assemble these matrices as Σ~ρ=(U⁡(ρ)ρR(ρ)ρR(ρ)⊤ρV(ρ)). _ρ= pmatrixU(ρ)&ρ R(ρ)\\ ρ R(ρ) &ρ V(ρ) pmatrix. The equivalence of the rescaled equations and (132) shows that this matrix solves the transformed Lyapunov equation. Proposition 1 ensures that LρL_ρ is Hurwitz for sufficiently small positive ρ, so that equation has a unique solution. Hence the assembled matrix is precisely TΣρT⊤T _ρT . This establishes the expansion for the leading covariance coefficient. Calculation of the Coefficients To obtain the assembled covariance through order ρ2ρ^2, we need U through order ρ2ρ^2 and R,VR,V through order ρ, since the latter two blocks already carry a factor ρ. Write U⁡(ρ) U(ρ) =U0+ρU1+ρ2U2+O⁡(ρ3), =U_0+ρ U_1+ρ^2U_2+O(ρ^3), R⁡(ρ) R(ρ) =R0+ρR1+O⁡(ρ2), =R_0+ρ R_1+O(ρ^2), V⁡(ρ) V(ρ) =V0+ρV1+O⁡(ρ2). =V_0+ρ V_1+O(ρ^2). These are power-series coefficients; in particular, U2=U′(0)/2U_2=U (0)/2. First compute U0,R0,V0U_0,R_0,V_0 from the successive equations (137)–(139). Next, match the coefficients of ρ in the three rescaled equations. This gives, in the order U1,R1,V1U_1,R_1,V_1, AU1+U1A⊤ AU_1+U_1A =ΛCU0+U0C⊤Λ⊤+ΛQθx+QxθΛ⊤, = CU_0+U_0C + Q_θ x+Q_xθ , (143) AR1 AR_1 =ΛCR0+ΛSV0−U1C⊤−R0S⊤+ΛQθθ, = CR_0+ SV_0-U_1C -R_0S + Q_θ, (144) SV1+V1S⊤ SV_1+V_1S =−CR1−R1⊤C⊤. =-CR_1-R_1 C . (145) Matching the coefficient of ρ2ρ^2 in (134) then gives AU2+U2A⊤= AU_2+U_2A = ΛCU1+U1C⊤Λ⊤ CU_1+U_1C (146) +ΛSR0⊤+R0S⊤Λ⊤−ΛQθθΛ⊤. + SR_0 +R_0S - Q_θ . Each right-hand side contains only known matrices at the time it is used. The equations for U1,V1,U2U_1,V_1,U_2 have unique solutions because the corresponding Lyapunov operators are invertible; the equation for R1R_1 is uniquely solved by multiplication by A−1A^-1. Return to the Original Coordinates. Define Σ~0=(U0000),Σ~1=(U1R0R0⊤V0),Σ~2=(U2R1R1⊤V1). _0= pmatrixU_0&0\\ 0&0 pmatrix, _1= pmatrixU_1&R_0\\ R_0 &V_0 pmatrix, _2= pmatrixU_2&R_1\\ R_1 &V_1 pmatrix. The assembled expansion is Σ~ρ=Σ~0+ρΣ~1+ρ2Σ~2+O⁡(ρ3). _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3). Using T from (129), set Σj=T−1Σ~j(T−1)⊤,j=0,1,2. _j=T^-1 _j(T^-1) , j=0,1,2. Since T is independent of ρ, transforming back yields Σρ=Σ0+ρΣ1+ρ2Σ2+O⁡(ρ3),Σ0=diag⁡(U0,0). _ρ= _0+ρ _1+ρ^2 _2+O(ρ^3), _0=diag(U_0,0). After choosing ρ0 _0 so that [0,ρ0][0, _0] lies inside the analyticity neighborhood, the derivatives needed for these Taylor remainders are uniformly bounded. Thus the final remainder is bounded by Cρ3Cρ^3 for 0≤ρ≤ρ00≤ρ≤ _0, where the value at zero denotes the analytic extension. Substitution into (125)–(127) then gives the three bias coefficients. E.5 Bias Coefficients The preceding subsection determines the covariance coefficients Σ0,Σ1,Σ2 _0, _1, _2. We now use them to compute the bias coefficients b0,b1,b2b_0,b_1,b_2 and establish the stationary mean expansion in the two step sizes. The starting point is the formula from Theorem 1: b⁡(ρ)=−J−1⁡(Σρ)+rρ.b(ρ)=-J^-1\C( _ρ)+r_ρ\. First consider the direct Markov response rρr_ρ. Let (Y0,Y1)(Y_0,Y_1) have the stationary one-step distribution ℙ⁡(Y0=y,Y1=y′)=μ⁡(y)P(y,y′)P(Y_0=y,Y_1=y )=μ(y)P(y,y ), and write ξ0=G⁡(Y0,z⋆)=(ξ0xξ0θ). _0=G(Y_0,z )= pmatrix _0^x\\ _0^θ pmatrix. By (108), rρ=⁡[Dz(Y1,z⋆)Dρξ0].r_ρ=E\! [D_zU(Y_1,z )D_ρ _0 ]. The matrix DρD_ρ leaves the fast component unchanged and multiplies the slow component by ρ. Since the other quantities in this expectation do not depend on ρ, we obtain the exact decomposition rρ r_ρ =r0+ρr1, =r_0+ρ r_1, (147) r0 r_0 =⁡[Dx(Y1,z⋆)ξ0x], =E[D_xU(Y_1,z ) _0^x], r1 r_1 =⁡[Dθ(Y1,z⋆)ξ0θ]. =E[D_θU(Y_1,z ) _0^θ]. Here DxD_xU and DθD_θU are derivatives of the full vector-valued corrector with respect to the fast and slow variables, respectively. Next, insert this decomposition and the covariance expansion (142) into the formula for b⁡(ρ)b(ρ). Because C is linear and J−1J^-1 is fixed, b⁡(ρ)=b0+ρb1+ρ2b2+O⁡(ρ3),b(ρ)=b_0+ρ b_1+ρ^2b_2+O(ρ^3), (148) where b0 b_0 =−J−1⁡(Σ0)+r0, =-J^-1\C( _0)+r_0\, (149) b1 b_1 =−J−1⁡(Σ1)+r1, =-J^-1\C( _1)+r_1\, (150) b2 b_2 =−J−1(Σ2). =-J^-1C( _2). (151) These vectors are independent of both step sizes. There is no quadratic term in rρr_ρ, so b2b_2 has no direct Markov-response term. It can nevertheless depend on temporal correlations through the covariance coefficient Σ2 _2. To pass from the coefficient expansion to the stationary mean, we combine two uniform bounds. Theorem 1 gives ‖mη,ρ−ηb(ρ)‖≤C1η3/2,\|m_η,ρ-η b(ρ)\|≤ C_1η^3/2, while (148) gives ‖b⁡(ρ)−b0−ρb1−ρ2b2‖≤C2ρ3.\|b(ρ)-b_0-ρ b_1-ρ^2b_2\|≤ C_2ρ^3. The constants are independent of sufficiently small positive η,ρη,ρ. Multiplying the second bound by η and adding the first yields ‖mη,ρ−ηb0−ηρb1−ηρ2b2‖≤C1η3/2+C2ηρ3. \|m_η,ρ-η b_0-ηρ b_1-ηρ^2b_2 \|≤ C_1η^3/2+C_2ηρ^3. Using ηρ=εηρ= , ηρ2=ε2/ηρ^2= ^2/η, and ηρ3=ε3/η2ηρ^3= ^3/η^2 proves Theorem 2, restated below. Theorem 2 (Singular expansion; restatement). Under , , , and , there exist η0,ρ0>0 _0, _0>0 such that, for 0<η≤η00<η≤ _0 and 0<ρ≤ρ00<ρ≤ _0, with ε=ηρ =ηρ, mη,ρ m_η,ρ =ηb0+εb1+ε2ηb2 =η b_0+ b_1+ ^2ηb_2 (152) +O⁡(η3/2+ε3η2), +O\! (η^3/2+ ^3η^2 ), where the implied constant is independent of η and ρ. The factor ε2/η ^2/η comes from the quadratic term ηρ2ηρ^2 in the ratio expansion. Although it contains a negative power of η when written in the two step sizes, it satisfies 0<ε2η=ηρ2≤ηρ02.0< ^2η=ηρ^2≤η _0^2. Thus the mixed contribution tends to zero uniformly as η↓0η 0 in the stated parameter range. Along paths with ρ↓0ρ 0, ε2/ηε=ρ⟶0, ^2/η =ρ 0, so its step-size factor is o⁡(ε)o( ). Whether this contribution can be distinguished from the remainder is considered in the next subsection. E.6 Visibility along ε=ηp =η^p The expansion in the preceding subsection displays three bias terms, but some may be smaller than the approximation error. We now compare their sizes along a prescribed relation between the two step sizes. This identifies which terms can be distinguished from the proved remainder and determines the powers relevant to extrapolation in the next subsection. Fix p>1p>1 and choose the step-size path ε=ηp. =η^p. (153) Then ρ=ε/η=ηp−1↓0ρ= /η=η^p-1 0 as η↓0η 0. Define m⁡(η):=mη,ηp−1m(η):=m_η,η^p-1. Substitution into (152) gives m⁡(η) m(η) =ηb0+ηpb1+η2p−1b2 =η b_0+η^pb_1+η^2p-1b_2 (154) +O⁡(η3/2+η3p−2). +O(η^3/2+η^3p-2). Both contributions to the remainder must be retained: O⁡(η3/2)O(η^3/2) comes from the stationary mean approximation, and O⁡(η3p−2)O(η^3p-2) comes from truncating the ratio expansion. A nonzero bias term is resolved by this expansion when, after subtracting the preceding terms, the remaining error is of smaller order than that term. For positive powers of η, a smaller exponent means slower decay as η↓0η 0. Since 1<p<2p−1<3p−2,1<p<2p-1<3p-2, the error O⁡(η3p−2)O(η^3p-2) from omitting higher-order terms in the expansion of b⁡(ρ)b(ρ) is of smaller order than every displayed bias term. The comparison with η3/2η^3/2 therefore determines which terms are resolved. For nonzero corresponding coefficients, ηb0η b_0 is resolved for every p>1p>1, while ηpb1 η^pb_1 :1<p<32, : 1<p< 32, η2p−1b2 η^2p-1b_2 :1<p<54. : 1<p< 54. For example, subtracting the first two contributions and dividing by the mixed-term scale gives m⁡(η)−ηb0−ηpb1η2p−1=b2+O⁡(η5/2−2p+ηp−1). m(η)-η b_0-η^pb_1η^2p-1=b_2+O\! (η^5/2-2p+η^p-1 ). Both error terms tend to zero when 1<p<5/41<p<5/4. At p=5/4p=5/4, the first error bound is only O⁡(1)O(1), so the stated expansion does not guarantee convergence of this normalized expression to b2b_2. Under the additional assumptions of Theorem 4, the stationary mean remainder improves to O⁡(η2)O(η^2). The path expansion then has remainder O⁡(η2+η3p−2)O(η^2+η^3p-2). The same comparison extends the range for a nonzero b1b_1 term to 1<p<21<p<2 and for a nonzero b2b_2 term to 1<p<3/21<p<3/2. These endpoints are strict: when a term and the allowed remainder have the same order, the bound alone does not separate them. More generally, a term ρjρ^j in the analytic expansion of b⁡(ρ)b(ρ) contributes ηρj=η1+j⁡(p−1)ηρ^j=η^1+j(p-1) to ηb(ρ)η b(ρ). Thus the ratio expansion generates the exponents γj=1+j(p−1),j=0,1,2,…:1,p, 2p−1, 3p−2,…. _j=1+j(p-1), j=0,1,2,…: 1,\ p,\ 2p-1,\ 3p-2,…. (155) This sequence describes the powers in ηb(ρ)η b(ρ); expanding it further does not reduce the stationary approximation error. If a higher-order stationary expansion establishes additional terms, such as a term proportional to η2η^2, their exponents must also be included when choosing extrapolation weights. Contributions with the same exponent are combined into a single coefficient along the path. E.7 Richardson–Romberg Extrapolation along Step-Size Paths The preceding subsection identifies the bias powers along a chosen step-size path. We now combine stationary means at several step sizes so that selected bias terms cancel while the equilibrium z⋆z is preserved. We begin with a fixed step-size ratio, then consider the path ε=hp =h^p, and finally describe the weights for cancelling an arbitrary finite set of bias powers. First fix ρ and let M⁡(h)M(h) denote the stationary mean at fast step h and slow step ρhρ h. Theorem 1 gives M⁡(h)=z⋆+hb(ρ)+O⁡(h3/2)M(h)=z +hb(ρ)+O(h^3/2). Halving both step sizes preserves ρ, so the same coefficient b⁡(ρ)b(ρ) appears at the finer level. The weights −1-1 and 22 preserve the equilibrium and cancel the term linear in h, giving 2M(h/2)−M⁡(h)=z⋆+O⁡(h3/2).2M(h/2)-M(h)=z +O(h^3/2). (156) Now fix p>1p>1 and use M⁡(h)M(h) for the stationary mean at fast step h and slow step hph^p. Refinement must preserve this relation, so the step-size pairs are (hj,εj)=(2−jh,hjp)(h_j, _j)=(2^-jh,h_j^p). The preceding two-level rule cancels the term proportional to h but generally leaves the term proportional to hph^p. To cancel both powers while preserving z⋆z , combine three levels with weights satisfying ∑j=02wj=1,∑j=02wj2−j=0,∑j=02wj2−jp=0. _j=0^2w_j=1, _j=0^2w_j2^-j=0, _j=0^2w_j2^-jp=0. (157) The first constraint preserves the common equilibrium. The second and third cancel the h and hph^p terms, whose coefficients at level j are multiplied by 2−j2^-j and 2−jp2^-jp, respectively. Put q1=2−1q_1=2^-1 and qp=2−pq_p=2^-p. Solving these three linear equations gives w0=q1qp(1−q1)(1−qp),w1=−q1+qp(1−q1)(1−qp),w2=1(1−q1)(1−qp). w_0= q_1q_p(1-q_1)(1-q_p), w_1=- q_1+q_p(1-q_1)(1-q_p), w_2= 1(1-q_1)(1-q_p). (158) Applying these weights to (154) removes its first two bias terms. The remaining stationary bias is O⁡(hmin⁡2p−1, 3/2)O(h \2p-1,\,3/2\). The exponent reflects the surviving mixed term and the stationary approximation error. The additional error of order h3p−2h^3p-2 is absorbed into this bound because 3p−2>2p−13p-2>2p-1 for p>1p>1. When the sharper uniform weak remainder is available, as under the additional assumptions of Theorem 4, replace 3/23/2 by 22. More generally, suppose we wish to cancel the distinct positive bias powers γ1,…,γM _1,…, _M. Using M+1M+1 levels of the same geometric refinement, the weights must satisfy ∑j=0Mwj=1,∑j=0Mwj2−jγℓ=0,ℓ=1,…,M. _j=0^Mw_j=1, _j=0^Mw_j2^-j _ =0, =1,…,M. (159) As before, the first equation preserves the equilibrium, and each remaining equation cancels one selected power. These conditions can equivalently be expressed through the polynomial W⁡(t)=∑j=0Mwjtj=∏ℓ=1Mt−2−γℓ1−2−γℓ.W(t)= _j=0^Mw_jt^j= _ =1^M t-2^- _ 1-2^- _ . Each numerator factor gives W⁡(2−γℓ)=0W(2^- _ )=0, while the denominator normalizes the product so that W⁡(1)=1W(1)=1. Expanding this polynomial in powers of t gives the weights w0,…,wMw_0,…,w_M. For example, adding γ3=2p−1 _3=2p-1 to the powers 11 and p gives a four-level rule that also cancels the term containing b2b_2. The resulting bias order remains limited by the remainder established for the stationary expansion; the weight equations alone do not justify a sharper remainder. For finite runs, initialization error must also be controlled at the retained remainder scale to preserve the stated bias order. Appendix F Additive Independent Noise For additive i.i.d. noise, we can sharpen the covariance and mean remainders to O⁡(η2)O(η^2). The key step is an O⁡(η2)O(η^2) bound on centered third moments, obtained from cancellations in their stationary balance equation. We use this bound first to refine the covariance and then to obtain the mean expansion. The theorem below states the required smoothness and moment assumptions explicitly, including localization around the selected equilibrium z⋆z . Theorem 4 (O⁡(η2)O(η^2) remainder under additive i.i.d. noise). Consider zk+1=zk+ηDρG¯(zk)+ζk+1,z_k+1=z_k+η D_ρ\ G(z_k)+ _k+1\, (160) where the innovations are identically distributed, independent of one another and of the past, centered, and have a finite sixth moment. Assume that G¯(z⋆)=0 G(z )=0, that the matrices A and S in (4) are Hurwitz, and that G¯∈C4(ℝdx+dθ) G∈ C^4(R^d_x+d_θ) has bounded derivatives of orders one through four. For every sufficiently small positive η,ρη,ρ, suppose that a stationary law exists and satisfies ∥z−z⋆∥2j≤Cjηj,j=1,2,3,E\|z-z \|^2j≤ C_jη^j, j=1,2,3, (161) with constants independent of the two step sizes. Let Qζ=⁡[ζ1ζ1⊤]Q_ζ=E[ _1 _1 ], and let Σρ _ρ solve (14) with Q=QζQ=Q_ζ. Then, uniformly for 0<ρ≤ρ00<ρ≤ _0, Cov⁡(z−z⋆) (z-z ) =ηΣρ+O⁡(η2), =η _ρ+O(η^2), (162) m:=⁡(z−z⋆) m:=E(z-z ) =−ηJ−1(Σρ)+O⁡(η2). =-η J^-1C( _ρ)+O(η^2). (163) Proof. We first bound the mean and center the recursion. We then establish an inverse estimate that is uniform as ρ↓0ρ 0 and use it to control the third moments and the covariance, in that order. The mean expansion follows from a final Taylor expansion of the stationary drift equation. All constants below are independent of sufficiently small positive η and ρ. Step 1: Mean Control and Centering. Stationary increment balance, together with the centered innovations and the invertibility of DρD_ρ for ρ>0ρ>0, gives G¯(z)=0E G(z)=0. Taylor expansion at z⋆z and (161) give 0=Jm+O⁡(‖z−z⋆‖2)=Jm+O⁡(η).0=Jm+O(E\|z-z \|^2)=Jm+O(η). Since J is fixed and invertible, m=O⁡(η)m=O(η) uniformly in ρ. This preliminary bound does not use the bias expansion to be proved. Put q=T⁡(z−z⋆−m)q=T(z-z -m), where T is defined in (6), and set Bρ=TDρB_ρ=TD_ρ. Define the nonlinear remainder by rm(q)=BρG¯(z⋆+m+T−1q)−Lρq.r_m(q)=B_ρ G(z +m+T^-1q)-L_ρq. The centered recursion is then exactly qk+1=qk+η⁡Lρqk+rm(qk)+Bρζk+1,q_k+1=q_k+η\L_ρq_k+r_m(q_k)+B_ρ _k+1\, (164) with q=0Eq=0 and rm(q)=0Er_m(q)=0. The localization bounds and m=O⁡(η)m=O(η) imply the centered moment bounds below; Taylor expansion at z⋆+mz +m gives the remainder estimates: ‖q‖2j \|q\|^2j ≤Cjηj,j=1,2,3, ≤ C_jη^j, j=1,2,3, (165) ‖rmu(q)‖ \|r_m^u(q)\| ≤C⁡(‖q‖2+η‖q‖+η), ≤ C(\|q\|^2+η\|q\|+η), ‖rmv(q)‖ \|r_m^v(q)\| ≤Cρ(‖q‖2+η‖q‖+η). ≤ Cρ(\|q\|^2+η\|q\|+η). (166) Indeed, the linearization at z⋆+mz +m differs from J by O⁡(η)O(η). Subtracting the expectation of the quadratic remainder introduces a constant of order η. The factor ρ in the second estimate comes from the lower block row of BρB_ρ. We will also use a linear-growth bound for the full centered drift d⁡(q)=Lρq+rm(q)d(q)=L_ρq+r_m(q). Global Lipschitz continuity gives ‖du(q)‖≤C⁡(‖q‖+η),‖dv(q)‖≤Cρ(‖q‖+η),d(q)=0.\|d^u(q)\|≤ C(\|q\|+ η), \|d^v(q)\|≤ Cρ(\|q\|+ η), (q)=0. Step 2: A Uniform Inverse for the Moment Equations. Both the covariance and third-moment equations involve the following operator, with s=2s=2 and s=3s=3, respectively: s,ρ=∑a=1sI⊗(a−1)⊗Lρ⊗I⊗(s−a). K_s,ρ= _a=1^sI (a-1) L_ρ I (s-a). Its inverse need not be uniformly bounded on arbitrary tensors. The estimate we need retains the factor ρ in the block whose indices are all slow coordinates. Specifically, if a tensor R satisfies ∥R∥≤aη,∥Rv⋯v∥≤ρaη,\|R\|≤ a_η, \|R_v·s v\|≤ρ a_η, then ‖s,ρ−1R‖≤Csaη.\| K_s,ρ^-1R\|≤ C_sa_η. (167) Here Rv⋯vR_v·s v denotes the block in which every tensor index is a slow coordinate. To prove the estimate, use the change of basis from the proof of Proposition : Qρ−1LρQρ=diag⁡(Af(ρ),ρSs(ρ)),Qρ=I+O⁡(ρ),Q_ρ^-1L_ρQ_ρ=diag(A_f(ρ),ρ S_s(ρ)), Q_ρ=I+O(ρ), where Af(ρ)A_f(ρ) and Ss(ρ)S_s(ρ) are uniformly Hurwitz. The all-slow block of the transformed forcing (Qρ−1)⊗sR(Q_ρ^-1) sR is still O⁡(ρaη)O(ρ a_η). Indeed, any contribution from a block containing a fast index includes at least one fast-to-slow entry of Qρ−1Q_ρ^-1, which is O⁡(ρ)O(ρ). In the transformed operator, a tensor block with at least one fast index has a uniformly bounded inverse: the integral of its product semigroup contains a uniformly decaying fast factor. The all-slow inverse is O⁡(1/ρ)O(1/ρ), but its forcing is O⁡(ρaη)O(ρ a_η). Thus every solution block is O⁡(aη)O(a_η). Transforming back proves (167). Step 3: The Centered Third-Moment Bound. We now derive the bound on signed third moments that sharpens the covariance remainder. Expanding the stationary third-moment identity using (164) and dividing by η gives an equation of the form 0=3,ρ(q⊗3)+R3.0= K_3,ρE(q 3)+R_3. We identify R3R_3 explicitly to show where independence and centering are used. Write q=qkq=q_k, q′=qk+1q =q_k+1, d=d⁡(q)d=d(q), n=Bρζk+1n=B_ρ _k+1, and h=d+nh=d+n, so that q′=q+ηhq =q+η h. For fixed indices i,j,ℓi,j, , let ∑cyc _cyc denote the sum over (i,j,ℓ)(i,j, ), (j,ℓ,i)(j, ,i), and (ℓ,i,j)( ,i,j), retaining all three placements even when indices coincide. Expanding the product and using stationarity gives the exact identity 0= 0= ∑cyc⁡[hiqjqℓ]+η∑cyc⁡[hihjqℓ]+η2[hihjhℓ]. _cycE[h_iq_jq_ ]+η _cycE[h_ih_jq_ ]+η^2E[h_ih_jh_ ]. The innovation n is centered and independent of q, and hence of d⁡(q)d(q), so every term containing exactly one innovation vanishes. Writing Vij=⁡[ninj]V_ij=E[n_in_j], the terms containing exactly two innovations also vanish: ⁡[ninjqℓ]=Vijqℓ=0,⁡[ninjdℓ]=Vijdℓ=0.E[n_in_jq_ ]=V_ijEq_ =0, [n_in_jd_ ]=V_ijEd_ =0. These cancellations use both q=0Eq=0 and d=0Ed=0. The term containing three innovations remains; no symmetry of the innovation distribution is required. Set Mijℓ=⁡[qiqjqℓ]M_ij =E[q_iq_jq_ ] and r=rm(q)r=r_m(q). Substituting d=Lρq+rd=L_ρq+r in the first-order term identifies 3,ρM K_3,ρM, leaving (R3)ijℓ= (R_3)_ij = ∑cyc⁡[riqjqℓ]+η∑cyc⁡[didjqℓ] _cycE[r_iq_jq_ ]+η _cycE[d_id_jq_ ] +η2[didjdℓ]+η2[ninjnℓ]. +η^2E[d_id_jd_ ]+η^2E[n_in_jn_ ]. The remainder estimate in (166) bounds the first group by C⁡‖q‖4+η‖q‖3+η‖q‖2=O⁡(η2).C\E\|q\|^4+ \|q\|^3+ \|q\|^2\=O(η^2). The linear-growth bound on d and the centered moment bounds give orders η5/2η^5/2 and η7/2η^7/2 for the second and third groups, respectively. The finite third moment of the innovations gives O⁡(η2)O(η^2) for the fourth group. To apply Step 2, we also need the all-slow block. Each term in the first group then contains a factor ρ from rvr^v. The remaining three groups contain at least ρ2ρ^2, ρ3ρ^3, and ρ3ρ^3, respectively, from their slow drift or innovation coordinates. Thus ‖R3‖≤Cη2\|R_3\|≤ Cη^2 and ‖(R3)vvv‖≤Cρη2\|(R_3)_v\|≤ Cρη^2. By (167), maxi,j,ℓ⁡|⁡(qiqjqℓ)|≤Cη2. _i,j, |E(q_iq_jq_ )|≤ Cη^2. (168) Step 4: Covariance Refinement. We next use the signed third-moment bound in the stationary covariance equation. Let Pc=⁡(qq⊤)P_c=E(q ). Independence and centering of the innovations give the exact balance 0= 0= LρPc+PcLρ⊤+ηBρQζBρ⊤ L_ρP_c+P_cL_ρ +η B_ρQ_ζB_ρ +⁡qrm(q)⊤+rm(q)q⊤+ηd⁡(q)d(q)⊤. +E\qr_m(q) +r_m(q)q \+ \d(q)d(q) \. (169) Expand rmr_m through its quadratic term at the stationary mean z⋆+mz +m. In the cross terms with q, the change in the linear coefficient contributes O⁡(η)‖q‖2=O⁡(η2)O(η)E\|q\|^2=O(η^2). The quadratic term contracts with the signed third moments in (168) and is therefore O⁡(η2)O(η^2). The cubic remainder contributes O⁡(‖q‖4)=O⁡(η2)O(E\|q\|^4)=O(η^2). All constants subtracted when centering these terms disappear because q=0Eq=0. Finally, the discretization term satisfies η‖d⁡(q)‖2=O⁡(η2) \|d(q)\|^2=O(η^2). Let R2R_2 denote the sum of the nonlinear cross terms and the discretization term in this balance. Its slow–slow block retains a factor ρ in each nonlinear term and ρ2ρ^2 in the discretization term. Consequently, ‖R2‖≤Cη2,‖(R2)vv‖≤Cρη2.\|R_2\|≤ Cη^2, \|(R_2)_v\|≤ Cρη^2. Subtract the Lyapunov equation for ηΣ~ρ=ηTΣρT⊤η _ρ=η T _ρT and apply (167) with s=2s=2. This yields Pc=ηΣ~ρ+O⁡(η2).P_c=η _ρ+O(η^2). (170) Transforming back with the fixed matrix T proves (162). Since m=O⁡(η)m=O(η), adding mm⊤m also gives the raw second moment needed for the mean expansion: ⁡(z−z⋆)(z−z⋆)⊤=T−1PcT−⁣⊤+mm⊤=ηΣρ+O⁡(η2).E\(z-z )(z-z ) \=T^-1P_cT^- +m =η _ρ+O(η^2). (171) Step 5: Mean Expansion. It remains to substitute the moment estimates into G¯(z)=0E G(z)=0. Expand G¯ G through degree three at z⋆z . Every raw signed third moment is O⁡(η2)O(η^2): using z−z⋆=T−1q+mz-z =T^-1q+m, it is a linear combination of the centered third moments in (168), products of m with entries of PcP_c, and products of three entries of m. The fourth-order remainder is bounded by C‖z−z⋆‖4=O⁡(η2)CE\|z-z \|^4=O(η^2), since the fourth derivatives are bounded. Using (171) for the quadratic term therefore gives 0=Jm+η(Σρ)+O⁡(η2).0=Jm+ ( _ρ)+O(η^2). Multiplication by the fixed matrix J−1J^-1 proves (163), with the required uniform remainder. ∎ Appendix G Proof of localization under uniform block contraction Proposition 4 (A sufficient localization condition). Suppose Assumption A3 holds, each G⁡(y,⋅)G(y,·) is globally C2C^2 with (9), and there exist a,d>0a,d>0, b,c≥0b,c≥ 0, bc<adbc<ad, such that, uniformly in the noise state y and all arguments, ⟨x−x′,H⁡(y,x,θ)−H⁡(y,x′,θ)⟩ x-x ,H(y,x,θ)-H(y,x ,θ) ≤−a‖x−x′‖2, ≤-a\|x-x \|^2, ⟨θ−θ′,F⁡(y,x,θ)−F⁡(y,x,θ′)⟩ θ-θ ,F(y,x,θ)-F(y,x,θ ) ≤−d‖θ−θ′‖2, ≤-d\|θ-θ \|^2, ‖H⁡(y,x,θ)−H⁡(y,x,θ′)‖ \|H(y,x,θ)-H(y,x,θ )\| ≤b‖θ−θ′‖, ≤ b\|θ-θ \|, ‖F⁡(y,x,θ)−F⁡(y,x′,θ)‖ \|F(y,x,θ)-F(y,x ,θ)\| ≤c‖x−x′‖. ≤ c\|x-x \|. (172) Then G¯ G has a unique equilibrium and a globally defined C2C^2 fast equilibrium map. For sufficiently small positive η,ρη,ρ, the joint recursion has a unique invariant law satisfying Assumptions A4 and A7, with ‖z−z⋆‖≤Cη\|Ez-z \|≤ Cη uniformly. These conditions verify localization; the remaining assumptions of the bias theorems remain in force. Proof of Proposition 4. All coordinate norms are Euclidean, and matrix norms are the corresponding operator norms. Constants may depend on the updates, dimensions, and fixed Markov kernel, but not on sufficiently small positive η,ρη,ρ. We write ε=ρη =ρη and take η0,ρ0≤1 _0, _0≤ 1. 1. Averaged Equilibrium and Deterministic Comparison. Averaging preserves (172). For each fixed θ, global Lipschitz continuity and strong dissipativity make x↦x+hH¯(x,θ)x x+h H(x,θ) a contraction for sufficiently small h>0h>0. Its unique fixed point defines λ⁡(θ)λ(θ). For x=λ⁡(θ)x=λ(θ), x′=λ⁡(θ′)x =λ(θ ), δx=x−x′δ x=x-x , and δθ=θ−θ′δθ=θ-θ , the equilibrium identities and the averaged block bounds give 0 0 =⟨δx,H¯(x,θ)−H¯(x′,θ′)⟩ = δ x, H(x,θ)- H(x ,θ ) ≤−a‖δx‖2+b‖δx‖‖δθ‖. ≤-a\|δ x\|^2+b\|δ x\|\,\|δθ\|. Consequently, ‖λ⁡(θ)−λ⁡(θ′)‖≤ba‖θ−θ′‖.\|λ(θ)-λ(θ )\|≤ ba\|θ-θ \|. Strong dissipativity also makes DxH¯D_x H everywhere nonsingular. The implicit function theorem and uniqueness of the fast equilibrium therefore give a globally defined C2C^2 map λ. For the reduced drift g⁡(θ)=F¯(λ⁡(θ),θ)g(θ)= F(λ(θ),θ), the same block bounds yield ⟨δθ,g⁡(θ)−g⁡(θ′)⟩ δθ,g(θ)-g(θ ) ≤−d‖δθ‖2+c‖δθ‖‖λ⁡(θ)−λ⁡(θ′)‖ ≤-d\|δθ\|^2+c\|δθ\|\,\|λ(θ)-λ(θ )\| ≤−(d−cba)‖δθ‖2. ≤- (d- cba )\|δθ\|^2. Thus g is strongly dissipative with constant κg=d−cb/a>0 _g=d-cb/a>0. If LθL_θ bounds ‖DθF¯‖\|D_θ F\|, then ‖g⁡(θ)−g⁡(θ′)‖≤(Lθ+cba)‖θ−θ′‖.\|g(θ)-g(θ )\|≤ (L_θ+ cba )\|θ-θ \|. Writing Lg=Lθ+cb/aL_g=L_θ+cb/a, we obtain ‖θ+hg(θ)−θ′−hg(θ′)‖2≤(1−2κgh+Lg2h2)‖δθ‖2.\|θ+hg(θ)-θ -hg(θ )\|^2≤(1-2 _gh+L_g^2h^2)\|δθ\|^2. For sufficiently small h>0h>0, this is a contraction. Banach’s fixed-point theorem gives a unique zero θ⋆θ of g, and z⋆=(λ⁡(θ⋆),θ⋆)z =(λ(θ ),θ ) is the unique equilibrium of G¯ G. We next compare two trajectories driven by the same input. Choose a0∈(0,a)a_0∈(0,a) and d0∈(0,d)d_0∈(0,d) with bc<a0d0bc<a_0d_0. If necessary, enlarge b,cb,c slightly so that they are positive and still satisfy this inequality. Let L be a uniform bound for the derivatives of both update fields, and write Δx=x−x′ x=x-x and Δθ=θ−θ′ θ=θ-θ . For a common noise state y, set ΔHx=H⁡(y,x,θ)−H⁡(y,x′,θ) H_x=H(y,x,θ)-H(y,x ,θ). Then ‖Δx+ηΔHx‖2≤(1−2aη+L2η2)‖Δx‖2≤(1−a0η)2‖Δx‖2\| x+η H_x\|^2≤(1-2aη+L^2η^2)\| x\|^2≤(1-a_0η)^2\| x\|^2 for sufficiently small η. Taking square roots and applying the cross-Lipschitz bound gives ‖Δxk+1‖≤(1−a0η)‖Δxk‖+bη‖Δθk‖.\| x_k+1\|≤(1-a_0η)\| x_k\|+bη\| _k\|. The analogous calculation for the slow update gives ‖Δθk+1‖≤cε‖Δxk‖+(1−d0ε)‖Δθk‖.\| _k+1\|≤ c \| x_k\|+(1-d_0 )\| _k\|. Since ε≤η ≤η, a single sufficiently small η0 _0, independent of ρ, ensures both bounds and nonnegative diagonal coefficients. Hence (‖Δxk+1‖Δθk+1‖)≤M(‖Δxk‖Δθk‖),M=(1−a0ηbηcε1−d0ε). pmatrix\| x_k+1\|\\ \| _k+1\| pmatrix≤ M pmatrix\| x_k\|\\ \| _k\| pmatrix, M= pmatrix1-a_0η&bη\\ c &1-d_0 pmatrix. (173) The comparison is componentwise, and M is nonnegative. Choose t with c/d0<t<a0/bc/d_0<t<a_0/b and define ‖(x,θ)‖t=max⁡‖x‖,‖θ‖/t.\|(x,θ)\|_t= \\|x\|,\|θ\|/t\. The corresponding weighted row sums of M are 1−(a0−bt)η1-(a_0-bt)η and 1−(d0−c/t)ε1-(d_0-c/t) . Thus every update map has Lipschitz constant at most 1−γε1-γ , where γ=min⁡a0−bt,d0−c/t>0γ= \a_0-bt,d_0-c/t\>0 is independent of the steps. 2. A Common Invariant Set and the Stationary Law. Let BH=maxξ⁡‖H⁡(ξ,0,0)‖,BF=maxξ⁡‖F⁡(ξ,0,0)‖,B_H= _ξ\|H(ξ,0,0)\|, B_F= _ξ\|F(ξ,0,0)\|, and choose R large enough that (a0−bt)R≥BH,(d0t−c)R≥BF.(a_0-bt)R≥ B_H, (d_0t-c)R≥ B_F. The comparison above shows that the convex compact set K=(x,θ):∥x∥≤R,∥θ∥≤tRK=\(x,θ):\|x\|≤ R,\ \|θ\|≤ tR\ is invariant under every noise-state update, uniformly in both step sizes. Starting at zero at time −N-N along a two-sided stationary noise sequence and letting N→∞N→∞ gives a limit Zk=(Xk,Θk)∈KZ_k=(X_k, _k)∈ K. Contraction ensures convergence; the limit is measurable with respect to σ(Yi:i≤k−1)σ(Y_i:i≤ k-1) and is stationary jointly with YkY_k. Any invariant law must be supported on K. Comparison with a trajectory started in K gives distt⁡(zk,K)≤(1−γε)kdistt⁡(z0,K).dist_t(z_k,K)≤(1-γ )^kdist_t(z_0,K). Under an invariant law, the left-hand side has the same distribution for every k, while the right-hand side tends to zero almost surely. The Doeblin condition allows two copies of the noise to couple in finite time almost surely. Their iterates contract thereafter, which proves uniqueness of the joint invariant law. 3. Influence of a Past Input. To obtain separate fast and slow moment bounds, we estimate the effect of changing one past input. Write M=I−ηKρ,Kρ=(a0−b−ρcρd0).M=I-η K_ρ, K_ρ= pmatrixa_0&-b\\ -ρ c&ρ d_0 pmatrix. Its eigenvalues are ℓ±=a0+ρd0±(a0−ρd0)2+4ρbc2. _±= a_0+ρ d_0± (a_0-ρ d_0)^2+4ρ bc2. For sufficiently small ρ0 _0, uniformly over 0<ρ≤ρ00<ρ≤ _0, ℓ+≍1,ℓ−≍ρ,ℓ+−ℓ−≥a0/2,ℓ+−a0=O⁡(ρ). _+ 1, _- ρ, _+- _-≥ a_0/2, _+-a_0=O(ρ). Decrease η0 _0 so that 1−ηℓ±∈[0,1]1-η _±∈[0,1]. The spectral formula Mk=(1−ηℓ+)kKρ−ℓ−Iℓ+−ℓ−+(1−ηℓ−)kℓ+I−Kρℓ+−ℓ−M^k=(1-η _+)^k K_ρ- _-I _+- _-+(1-η _-)^k _+I-K_ρ _+- _- then gives the entrywise bound 0≤Mk≤C(e−κηk+ρe−κεke−κεkρe−κεke−κεk)0≤ M^k≤ C pmatrixe^-κη k+ρ e^-κ k&e^-κ k\\ ρ e^-κ k&e^-κ k pmatrix (174) for some fixed κ>0κ>0. Since the update fields are bounded on K, changing one noise input changes the next iterate by at most C⁡(η,ε)C(η, ) componentwise. Applying (174) for the subsequent k updates bounds its effects on X0X_0 and Θ0 _0 by, respectively, Ak=C⁡(ηe−κηk+εe−κεk),Bk=Cεe−κεk.A_k=C (η e^-κη k+ e^-κ k ), B_k=C e^-κ k. (175) These bounds hold for arbitrary input sequences because every update preserves K. Telescoping therefore bounds the effect of several input changes by the sum of their individual bounds. 4. Centered Moments under Markov Dependence. We combine the influence bounds with a coupling of the Markov inputs. The Doeblin condition allows chains started in any two states to be coupled so that ℙ⁡(Yℓ≠Yℓ′)≤CPqℓ,ℓ≥0,P(Y_ ≠ Y _ )≤ C_Pq , ≥ 0, (176) for fixed q∈(0,1)q∈(0,1) and CP<∞C_P<∞. Couple the endpoints of successive blocks of length n0n_0 using the common minorizing measure, and use identical transitions after coupling. Sampling intermediate paths conditionally on their endpoints gives (176) at all integer times. Set ℱj=σ(Yi:i≤j)F_j=σ(Y_i:i≤ j), j≤−1j≤-1, and define DjX=⁡[X0∣ℱj]−⁡[X0∣ℱj−1],DjΘ=⁡[Θ0∣ℱj]−⁡[Θ0∣ℱj−1].D_j^X=E[X_0 _j]-E[X_0 _j-1], D_j =E[ _0 _j]-E[ _0 _j-1]. Conditionally on ℱj−1F_j-1, ZjZ_j is fixed. Comparing two possible values of YjY_j by coupling their future chains and telescoping the resulting input changes gives, for j=−k−1j=-k-1, ‖D−k−1X‖≤CP∑ℓ=0kqℓAk−ℓ,‖D−k−1Θ‖≤CP∑ℓ=0kqℓBk−ℓ.\|D_-k-1^X\|≤ C_P _ =0^kq A_k- , \|D_-k-1 \|≤ C_P _ =0^kq B_k- . Decrease η0 _0 so that qeκη0<1qe^κ _0<1. The convolutions are then bounded by constant multiples of AkA_k and BkB_k, uniformly in both steps. Enlarging the constant in (175), we obtain ‖D−k−1X‖≤Ak,‖D−k−1Θ‖≤Bk.\|D_-k-1^X\|≤ A_k, \|D_-k-1 \|≤ B_k. Geometric summation gives ∑k≥0Ak2≤Cη,∑k≥0Bk2≤Cε,∑k≥0AkBk≤Cε. _k≥ 0A_k^2≤ Cη, _k≥ 0B_k^2≤ C , _k≥ 0A_kB_k≤ C . (177) The remote past tail of the stationary Doeblin chain is trivial. Reverse martingale convergence and the measurability of Z0Z_0 with respect to ℱ−1F_-1 therefore give the L2L^2 expansions X0−X0=∑j=−∞−1DjX,Θ0−Θ0=∑j=−∞−1DjΘ.X_0-EX_0= _j=-∞^-1D_j^X, _0-E _0= _j=-∞^-1D_j . For each scalar coordinate, conditional Hoeffding’s inequality applied to finite partial sums gives a sub-Gaussian moment generating function with variance proxy bounded by ∑Ak2Σ A_k^2 or ∑Bk2Σ B_k^2. Passing to the limit and integrating the tail bound yields, for each fixed integer j≥1j≥ 1, ‖X0−X0‖2j≤Cjηj,‖Θ0−Θ0‖2j≤Cjεj.E\|X_0-EX_0\|^2j≤ C_jη^j, \| _0-E _0\|^2j≤ C_j ^j. (178) The passage to Euclidean norms changes only dimension-dependent constants. Orthogonality of martingale increments at distinct times also gives ‖Cov⁡(X0,Θ0)‖≤∑k≥0⁡[‖D−k−1X‖‖D−k−1Θ‖]≤Cε.\|Cov(X_0, _0)\|≤ _k≥ 0E [\|D_-k-1^X\|\|D_-k-1 \| ]≤ C . (179) 5. Location of the Stationary Mean. We now locate the stationary mean relative to z⋆z , using stationarity and the Poisson equation. Let z¯=Z0 z=EZ_0 and G~(ξ,z)=G⁡(ξ,z)−G¯(z) G(ξ,z)=G(ξ,z)- G(z). The equation (I−P)(⋅,z)=G~(⋅,z),∑ξμ⁡(ξ)(ξ,z)=0(I-P)V(·,z)= G(·,z), _ξμ(ξ)V(ξ,z)=0 has a solution with bounded first z-derivative on K, since the inverse of I−PI-P on the μ-centered subspace is a fixed finite-dimensional operator. Exogeneity and joint stationarity give G~(Y0,Z0) G(Y_0,Z_0) =⁡(Y0,Z0)−⁡(Y1,Z0) =E\V(Y_0,Z_0)-V(Y_1,Z_0)\ =⁡(Y1,Z1)−⁡(Y1,Z0)=O⁡(η). =E\V(Y_1,Z_1)-V(Y_1,Z_0)\=O(η). The last bound uses convexity of K and ‖Z1−Z0‖≤Cη\|Z_1-Z_0\|≤ Cη on K. Stationary increment balance gives G(Y0,Z0)=0EG(Y_0,Z_0)=0. Taylor expansion at z¯ z, bounded second derivatives, and (178) yield ‖G¯(Z0)−G¯(z¯)‖≤C‖Z0−z¯‖2≤Cη.\|E G(Z_0)- G( z)\|≤ CE\|Z_0- z\|^2≤ Cη. Consequently, ‖G¯(z¯)‖≤Cη\| G( z)\|≤ Cη. For any w=(x,θ)w=(x,θ), set r=‖x−x⋆‖r=\|x-x \| and s=‖θ−θ⋆‖s=\|θ-θ \|. The averaged block bounds give ar≤bs+‖H¯(w)‖,ds≤cr+‖F¯(w)‖.ar≤ bs+\| H(w)\|, ds≤ cr+\| F(w)\|. Since ad−bc>0ad-bc>0, these inequalities imply ‖w−z⋆‖≤C‖G¯(w)‖.\|w-z \|≤ C\| G(w)\|. Applying this with w=z¯w= z proves ‖Z0−z⋆‖≤Cη.\|EZ_0-z \|≤ Cη. (180) 6. Localization in the Exact and Tangent Coordinates. Combining (178) and (180), for j=1,2,3j=1,2,3, ‖X0−x⋆‖2j \|X_0-x \|^2j ≤Cjηj, ≤ C_jη^j, ‖Θ0−θ⋆‖2j \| _0-θ \|^2j ≤Cj(εj+η2j). ≤ C_j( ^j+η^2j). For u^=X0−λ⁡(Θ0) u=X_0-λ( _0) and v=Θ0−θ⋆v= _0-θ , the Lipschitz bound on λ gives ‖u^‖≤‖X0−x⋆‖+(b/a)‖Θ0−θ⋆‖.\| u\|≤\|X_0-x \|+(b/a)\| _0-θ \|. Thus ‖u^‖2j≤CjηjE\| u\|^2j≤ C_jη^j. The inequality ‖W−W‖p≤2p‖W‖p,p≥1,E\|W-EW\|^p≤ 2^pE\|W\|^p, p≥ 1, gives the centered fast tracking bounds. The centered slow bounds are already contained in (178). In tangent coordinates, u=X0−x⋆−Λvu=X_0-x - v, with Λ=Dλ(θ⋆) =Dλ(θ ). Hence Cov⁡(u,v) (u,v) =Cov⁡(X0,Θ0)−ΛCov⁡(Θ0,Θ0), =Cov(X_0, _0)- ( _0, _0), ‖Cov⁡(u,v)‖ \|Cov(u,v)\| ≤Cε, ≤ C , by (179) and the centered slow second-moment bound. The raw fourth moments and ε≤η≤1 ≤η≤ 1 also give ‖Z0−z⋆‖4≤Cη2,E\|Z_0-z \|^4≤ Cη^2, which verifies Assumption A4. To control the coordinate remainder, differentiate H¯(λ(θ),θ)=0 H(λ(θ),θ)=0 twice. For h1,h2∈ℝdθh_1,h_2 ^d_θ, D2λ(θ)[h1,h2]= D^2λ(θ)[h_1,h_2]= −[DxH¯(λ⁡(θ),θ)]−1 -[D_x H(λ(θ),θ)]^-1 ×D2H¯(λ⁡(θ),θ)[(Dλ(θ)h1,h1),(Dλ(θ)h2,h2)]. × D^2 H(λ(θ),θ)[(Dλ(θ)h_1,h_1),(Dλ(θ)h_2,h_2)]. Strong dissipativity gives ‖[DxH¯(λ⁡(θ),θ)]−1‖≤a−1\|[D_x H(λ(θ),θ)]^-1\|≤ a^-1, and ‖Dλ(θ)‖≤b/a\|Dλ(θ)\|≤ b/a. Together with (9), these bounds make D2λD^2λ globally bounded. Therefore ‖rλ(v)‖≤C‖v‖2\|r_λ(v)\|≤ C\|v\|^2, and the raw fourth-moment bound yields ‖rλ(v)‖≤Cη,‖rλ(v)‖2≤Cη2.E\|r_λ(v)\|≤ Cη, \|r_λ(v)\|^2≤ Cη^2. The global C2C^2 map supplies the required extension, completing the verification of Assumption A7. ∎ Appendix H A nonlinear Markov example with an exact bias formula A shared two-state Markov input drives fast and slow linear filters, and a nonlinear readout converts their variance and cross-covariance into stationary bias. The exact mean contains a nonzero ε2/η ^2/η term and allows a direct comparison of the extrapolation rules in Section 6. H.1 Model and Markov Input Let ξk∈−1,+1 _k∈\-1,+1\ be a stationary symmetric Markov chain with transition matrix, in the state order (−1,+1)(-1,+1), Pq=((1+q)/2(1−q)/2(1−q)/2(1+q)/2),−1<q<1.P_q= pmatrix(1+q)/2&(1-q)/2\\ (1-q)/2&(1+q)/2 pmatrix, -1<q<1. (181) The invariant distribution is uniform, and ⁡(ξk+1∣ξk)=qξkE( _k+1 _k)=q _k. Consequently, ξk=0,⁡(ξkξk+h)=q|h|.E _k=0, ( _k _k+h)=q^|h|. (182) The case q=0q=0 gives independent symmetric signs. The scalar long-run covariance is cq=∑h∈ℤ⁡(ξ0ξh)=1+2∑h=1∞qh=1+q1−q.c_q= _h E( _0 _h)=1+2 _h=1^∞q^h= 1+q1-q. (183) Consider Xk+1 X_k+1 =(1−η)Xk+ηξk+1, =(1-η)X_k+η _k+1, (184) Yk+1 Y_k+1 =(1−ε)Yk+εξk+1, =(1- )Y_k+ _k+1, (185) Zk+1 Z_k+1 =Zk+ε⁡−Zk+ϕ⁡(Xk,Yk). =Z_k+ \-Z_k+φ(X_k,Y_k)\. (186) The fast coordinate is X, and the slow coordinates are (Y,Z)(Y,Z). The filters X and Y share the same input. We take 0<ε≤η≤10< ≤η≤ 1 and choose a smooth cutoff χ:ℝ→[0,1]χ:R→[0,1] equal to one on [−1,1][-1,1] and zero outside (−2,2)(-2,2). Define ϕ⁡(x,y)=χ⁡(x)χ(y)(x2+xy).φ(x,y)=χ(x)χ(y)(x^2+xy). (187) Thus ϕφ and all its derivatives are globally bounded. The stationary filters lie in [−1,1]2[-1,1]^2, as shown below, so the cutoff leaves their readout equal to X2+XYX^2+XY. The averaged drift is G¯(x,y,z)=(−x,−y,−z+ϕ⁡(x,y))⊤. G(x,y,z)=(-x,-y,-z+φ(x,y)) . Its unique equilibrium is the origin. The fast equilibrium map is λ≡0λ≡ 0, and the reduced slow drift is (−y,−z)⊤(-y,-z) . H.2 Stationary Representations Set a=1−η,d=1−ε.a=1-η, d=1- . (188) Backward iteration gives the unique stationary causal filters Xk=η∑i=0∞aiξk−i,Yk=ε∑j=0∞djξk−j.X_k=η _i=0^∞a^i _k-i, Y_k= _j=0^∞d^j _k-j. (189) Their weights are nonnegative and sum to one, so |Xk|,|Yk|≤1|X_k|,|Y_k|≤ 1 almost surely. Their means are zero by (182). The stationary slow readout is Zk=ε∑j≥0djϕ(Xk−1−j,Yk−1−j).Z_k= _j≥ 0d^jφ(X_k-1-j,Y_k-1-j). Since |ϕ⁡(x,y)|≤2|φ(x,y)|≤ 2 on [−1,1]2[-1,1]^2, this also gives |Zk|≤2|Z_k|≤ 2. H.3 Exact Second Moments Define the correlations with the current driver by rX=⁡(Xkξk),rY=⁡(Ykξk).r_X=E(X_k _k), r_Y=E(Y_k _k). (190) The Markov property gives ⁡(Xk−1ξk)=qrXE(X_k-1 _k)=qr_X and ⁡(Yk−1ξk)=qrYE(Y_k-1 _k)=qr_Y. Hence rX=aqrX+ηr_X=aqr_X+η and rY=dqrY+εr_Y=dqr_Y+ , so rX=η1−aq,rY=ε1−dq.r_X= η1-aq, r_Y= 1-dq. (191) Squaring the fast recursion and taking stationary expectations yields (1−a2)X2=2aηqrX+η2.(1-a^2)EX^2=2aη qr_X+η^2. Substituting (191) gives X2=η21−a21+aq1−aq. EX^2= η^21-a^2 1+aq1-aq. (192) The same calculation for the slow filter gives Y2=ε21−d21+dq1−dq. EY^2= ^21-d^2 1+dq1-dq. (193) Multiplying the two filter recursions and taking expectations gives (1−ad)(XY) (1-ad)E(XY) =aεqrX+dηqrY+ηε =a qr_X+dη qr_Y+η =ηε1−adq2(1−aq)(1−dq). =η 1-adq^2(1-aq)(1-dq). Therefore ⁡(XY)=ηε(1−adq2)(1−ad)(1−aq)(1−dq). E(XY)= η (1-adq^2)(1-ad)(1-aq)(1-dq). (194) For independent observations, q=0q=0, this reduces to ⁡(XY)=ηεη+ε−ηε.E(XY)= η η+ -η . H.4 Stationary Mean Stationarity of (186) gives Z=ϕ(X,Y)EZ=Eφ(X,Y). Since the cutoff is inactive on the stationary support, Z=X2+⁡(XY). EZ=EX^2+E(XY). (195) The full mean displacement is (0,0,Z)(0,0,EZ). The two quadratic terms contribute the fast variance and the fast–slow cross-covariance, respectively. H.5 Fixed-Ratio Bias Set ε=ρη =ρη and first let η→0η→ 0 at fixed ρ. The exact fast variance is X2=η2−η1+q−qη1−q+qη.EX^2= η2-η 1+q-qη1-q+qη. (196) Taylor expansion at η=0η=0 gives X2=cq2η+O⁡(η2).EX^2= c_q2η+O(η^2). (197) Similarly, substituting a=1−ηa=1-η and d=1−ρηd=1-ρη into (194) gives ⁡(XY)=ρη1−q2+q2(1+ρ)η−q2ρη2(1+ρ−ρη)(1−q+qη)(1−q+qρη).E(XY)=ρη\, 1-q^2+q^2(1+ρ)η-q^2ρη^2(1+ρ-ρη)(1-q+qη)(1-q+qρη). (198) For fixed q∈(−1,1)q∈(-1,1), the rational factor is smooth near (η,ρ)=(0,0)(η,ρ)=(0,0), so ⁡(XY)=cqηρ1+ρ+O⁡(η2ρ).E(XY)=c_qη ρ1+ρ+O(η^2ρ). Combining the two moments yields Z=cqη12+ρ1+ρ+O⁡(η2),EZ=c_qη \ 12+ ρ1+ρ \+O(η^2), (199) uniformly for sufficiently small positive ρ. Thus the leading fixed-ratio coefficient is bZ(ρ)=cq12+ρ1+ρ.b_Z(ρ)=c_q \ 12+ ρ1+ρ \. (200) H.6 Singular Expansion and Coefficient Recovery Using ρ1+ρ=ρ−ρ2+ρ31+ρ ρ1+ρ=ρ-ρ^2+ ρ^31+ρ (201) in (199) gives Z=cq2η+cqηρ−cqηρ2+O⁡(ηρ3+η2).EZ= c_q2η+c_qηρ-c_qηρ^2+O(ηρ^3+η^2). (202) Returning to the two step sizes, Z=cq2η+cqε−cqε2η+O⁡(η2+ε3η2). EZ= c_q2η+c_q -c_q ^2η+O\! (η^2+ ^3η^2 ). (203) The vector coefficients are therefore b0=(00cq/2),b1=(00cq),b2=(00−cq).b_0= pmatrix0\\ 0\\ c_q/2 pmatrix, b_1= pmatrix0\\ 0\\ c_q pmatrix, b_2= pmatrix0\\ 0\\ -c_q pmatrix. (204) In particular, the mixed coefficient is nonzero for every fixed q∈(−1,1)q∈(-1,1). Write biZb_i^Z for the third component of bib_i. The coefficients can also be recovered algebraically through b^0 b_0 =Zη=b0Z+O⁡(ρ+η), = EZη=b_0^Z+O(ρ+η), b^1 b_1 =Z−ηb0Zε=b1Z+O⁡(ρ+ηρ), = EZ-η b_0^Z =b_1^Z+O\! (ρ+ ηρ ), b^2 b_2 =Z−ηb0Z−εb1Zε2/η=b2Z+O⁡(ρ+ηρ2). = EZ-η b_0^Z- b_1^Z ^2/η=b_2^Z+O\! (ρ+ ηρ^2 ). (205) Along ρ=ηp−1ρ=η^p-1, these bounds guarantee convergence of all three quantities for 1<p<3/21<p<3/2. The upper threshold is sharp when the regular η2η^2 coefficient is nonzero. H.7 Extrapolation along Power-law Paths For fixed p>1p>1, setting ε=ηp =η^p in (203) gives Z=cq2η+cqηp−cqη2p−1+O⁡(η3p−2+η2).EZ= c_q2η+c_qη^p-c_qη^2p-1+O(η^3p-2+η^2). (206) For 1<p<3/21<p<3/2, the mixed term precedes the regular quadratic remainder. More generally, analyticity of the exact formula in (η,ρ)(η,ρ) gives path powers k+j⁡(p−1)k+j(p-1), with k≥1k≥ 1 and j≥0j≥ 0. Contributions with the same exponent must be combined. Applying the extrapolation weights from Appendix E to the expansion above gives the following cancellation properties. For 1<p<3/21<p<3/2, cancelling the integer powers 1,2\1,2\ leaves the ηpη^p term, whereas cancelling 1,p\1,p\ leaves bias of order η2p−1η^2p-1. Cancelling 1,p,2p−1\1,p,2p-1\ gives the bound O⁡(ηmin⁡3p−2,2)O(η \3p-2,2\). This last exponent is attained when the corresponding combined coefficient is nonzero; additional cancellation can give a faster rate. H.8 Agreement with the General Bias Formula To apply (107), reindex the driver as sk=ξk+1s_k= _k+1. The update then has the form G⁡(sk,zk)G(s_k,z_k), where, for z=(x,y,z3)z=(x,y,z_3), G⁡(s,z)=(−x+s−y+s−z3+ϕ⁡(x,y)),G¯(z)=(−x−y−z3+ϕ⁡(x,y)).G(s,z)= pmatrix-x+s\\ -y+s\\ -z_3+φ(x,y) pmatrix, G(z)= pmatrix-x\\ -y\\ -z_3+φ(x,y) pmatrix. At the origin, J=−I3,Dρ=diag⁡(1,ρ,ρ).J=-I_3, D_ρ=diag(1,ρ,ρ). (207) Let ℓ=(1,1,0)⊤ =(1,1,0) . The centered update is ℓs s, and Pqs=qsP_qs=qs, so the centered Poisson solution is ⁡(s,z)=ℓs1−q.U(s,z)= s1-q. (208) It is independent of z, hence rρ=0.r_ρ=0. (209) The temporal dependence enters through the long-run covariance Q=cqℓℓ⊤=cq(110110000).Q=c_q =c_q pmatrix1&1&0\\ 1&1&0\\ 0&0&0 pmatrix. (210) The nonzero entries of the Lyapunov solution satisfy −2(Σρ)11+cq -2( _ρ)_11+c_q =0, =0, −(1+ρ)(Σρ)12+ρcq -(1+ρ)( _ρ)_12+ρ c_q =0, =0, −2ρ(Σρ)22+ρ2cq -2ρ( _ρ)_22+ρ^2c_q =0. =0. Thus Σρ=cq(1/2ρ/(1+ρ)0ρ/(1+ρ)ρ/20000). _ρ=c_q pmatrix1/2&ρ/(1+ρ)&0\\ ρ/(1+ρ)&ρ/2&0\\ 0&0&0 pmatrix. (211) The only nonzero Hessian of the averaged drift at the origin is D2G¯3(0)=(210100000).D^2 G_3(0)= pmatrix2&1&0\\ 1&0&0\\ 0&0&0 pmatrix. Its contraction with the covariance gives ⁡(Σρ)=cq(0012+ρ1+ρ).C( _ρ)=c_q pmatrix0\\ 0\\ 12+ ρ1+ρ pmatrix. (212) Since −J−1=I3-J^-1=I_3 and rρ=0r_ρ=0, the general formula reproduces (200). If the filters instead use independent centered drivers, their cross-covariance vanishes and the leading bias comes only from the fast variance. The mixing and coupling estimates are uniform when q ranges over a compact subset of (−1,1)(-1,1). Mixing slows as q approaches either endpoint: q↑1q 1 also makes cqc_q diverge, while q↓−1q -1 approaches periodic alternation. The exact rational mean formulas remain regular at q=−1q=-1, although the uniform mixing conditions used in the general theory no longer hold there. H.9 Stationary Localization and Coupling Verification For this model, u^=X u=X, v=(Y,Z)v=(Y,Z), and λ≡0λ≡ 0. We verify the fourth-moment condition in Assumption A4 and the additional block conditions in Assumption A7: πη,ρ‖u^‖2j _ _η,ρ u ^2j ≤Cηj, ≤ Cη^j, πη,ρ‖v‖2j _ _η,ρ v ^2j ≤Cεj+Cη2j, ≤ C ^j+Cη^2j, πη,ρ‖u^−u^‖2j _ _η,ρ u-E u ^2j ≤Cηj, ≤ Cη^j, πη,ρ‖v−v‖2j _ _η,ρ v-Ev ^2j ≤Cεj,j=1,2,3, ≤ C ^j, j=1,2,3, (213) Fix −1<q<1-1<q<1 and 0<ε≤η≤η0≤10< ≤η≤ _0≤ 1. Constants may depend on q, but not on the step sizes. Independent-Input Representation and Sensitivity. Put r=|q|r=|q| and choose σ∈−1,1σ∈\-1,1\ with q=rσq=rσ, taking σ=1σ=1 when q=0q=0. Let Wk=(Bk,ωk)W_k=(B_k, _k), k∈ℤk , be independent pairs, with Bk∼Bernoulli⁡(r)B_k (r) and ωk _k an independent symmetric sign. The recursion ξk=Bkσξk−1+(1−Bk)ωk _k=B_kσ _k-1+(1-B_k) _k defines a stationary causal chain: almost surely, a reset Bk=0B_k=0 occurs at a finite time in the past. Its invariant distribution is uniform and its conditional mean is qξk−1q _k-1, so its transition matrix is PqP_q. The stationary representations above give |Xk|,|Yk|≤1|X_k|,|Y_k|≤ 1 and |Zk|≤2|Z_k|≤ 2. Replace W0W_0 by an independent copy, keep all other inputs fixed, and denote the resulting processes by primes. The drivers agree after their first common reset at a positive time. For p∈2,4,6p∈\2,4,6\, set g0=1g_0=1 and gj=rj/pg_j=r^j/p for j≥1j≥ 1. Then ‖ξj−ξj′‖Lp≤2gj\| _j-ξ _j\|_L^p≤ 2g_j. Write δj,pA=‖Aj−Aj′‖Lpδ^A_j,p=\|A_j-A _j\|_L^p and let ∗* denote convolution on the nonnegative integers. With aj=η(1−η)j,dj=ε(1−ε)j,a_j=η(1-η)^j, d_j= (1- )^j, the filter representations and Minkowski’s inequality give δ⋅,pX≤2a∗g,δ⋅,pY≤2d∗g.δ^X_·,p≤ 2a*g, δ^Y_·,p≤ 2d*g. On the invariant square, ϕ⁡(x,y)=x2+xyφ(x,y)=x^2+xy and |ϕ⁡(x,y)−ϕ⁡(x′,y′)|≤3|x−x′|+|y−y′|.|φ(x,y)-φ(x ,y )|≤ 3|x-x |+|y-y |. Thus δϕ≤3δX+δYδ^φ≤ 3δ^X+δ^Y and δZ≤⁡(d∗δϕ)δ^Z (d*δ^φ), where (h)0=0(Sh)_0=0 and (h)j=hj−1(Sh)_j=h_j-1 for j≥1j≥ 1. This shift accounts for the use of (Xk,Yk)(X_k,Y_k) in the update of Zk+1Z_k+1 and preserves every sequence norm used below. Since ‖a‖ℓ1=‖d‖ℓ1=1,‖a‖ℓ22=η2−η,‖d‖ℓ22=ε2−ε,‖g‖ℓ1=11−r1/p,\|a\|_ ^1=\|d\|_ ^1=1, \|a\|_ ^2^2= η2-η, \|d\|_ ^2^2= 2- , \|g\|_ ^1= 11-r^1/p, Young’s convolution inequality yields ‖δ⋅,pX‖ℓ2≤Cpη,‖δ⋅,pY‖ℓ2≤Cpε,‖δ⋅,pϕ‖ℓ1≤Cp,‖δ⋅,pZ‖ℓ2≤Cpε.\|δ^X_·,p\|_ ^2≤ C_p η, \|δ^Y_·,p\|_ ^2≤ C_p , \|δ^φ_·,p\|_ ^1≤ C_p, \|δ^Z_·,p\|_ ^2≤ C_p . Stationary Moments and Cross-Covariance. Let ℋi=σ(Wj:j≤i)H_i=σ(W_j:j≤ i) and iA0=⁡(A0∣ℋi)−⁡(A0∣ℋi−1)P_iA_0=E(A_0 _i)-E(A_0 _i-1). Each stationary coordinate is bounded and causal, and the independent-input tail sigma-field is trivial. Consequently, A0−A0=∑i≤0iA0in Lp.A_0-EA_0= _i≤ 0P_iA_0 L^p. Replacing WiW_i by an independent copy and using conditional Jensen’s inequality gives ‖−jA0‖Lp≤δj,pA\|P_-jA_0\|_L^p≤δ^A_j,p. Orthogonality for p=2p=2, and the martingale square-function inequality followed by Minkowski’s inequality for p=4,6p=4,6, imply ‖A0−A0‖Lp≤Cp(∑j≥0‖−jA0‖Lp2)1/2≤Cp‖δ⋅,pA‖ℓ2.\|A_0-EA_0\|_L^p≤ C_p ( _j≥ 0\|P_-jA_0\|_L^p^2 )^1/2≤ C_p\|δ^A_·,p\|_ ^2. The centered ppth moments are therefore O⁡(ηp/2)O(η^p/2) for X and O⁡(εp/2)O( ^p/2) for Y,ZY,Z. Moreover, X=Y=0EX=EY=0 and Z=X2+⁡(XY)=O⁡(η),EZ=EX^2+E(XY)=O(η), by stationarity, Cauchy–Schwarz, and ε≤η ≤η. Combining these estimates with |A|p≤2p−1|A−A|p+|A|pE|A|^p≤ 2^p-1\E|A-EA|^p+|EA|^p\ proves (213). The sensitivity bounds also give ‖δ⋅,2X‖ℓ1≤C,‖δ⋅,2Y‖ℓ∞≤2‖d‖ℓ∞‖g‖ℓ1≤Cε,‖δ⋅,2Z‖ℓ∞≤‖d‖ℓ∞‖δ⋅,2ϕ‖ℓ1≤Cε.\|δ^X_·,2\|_ ^1≤ C, \|δ^Y_·,2\|_ ^∞≤ 2\|d\|_ ^∞\|g\|_ ^1≤ C , \|δ^Z_·,2\|_ ^∞≤\|d\|_ ^∞\|δ^φ_·,2\|_ ^1≤ C . Orthogonality of projections at distinct times therefore yields, for A=Y,ZA=Y,Z, |Cov⁡(X,A)|≤∑j≥0δj,2Xδj,2A≤‖δ⋅,2X‖ℓ1‖δ⋅,2A‖ℓ∞≤Cε.|Cov(X,A)|≤ _j≥ 0δ^X_j,2δ^A_j,2≤\|δ^X_·,2\|_ ^1\|δ^A_·,2\|_ ^∞≤ C . Since λ=0λ=0 globally, u=u^=Xu= u=X and rλ=0r_λ=0. This verifies the remaining covariance, extension, and coordinate-remainder conditions in Assumption A7. The raw fourth moments give ‖(X,Y,Z)‖4=O⁡(η2),E\|(X,Y,Z)\|^4=O(η^2), verifying Assumption A4. The cutoff makes all derivatives of ϕφ bounded, and the other drift terms are linear. Hence the updates satisfy the global derivative bounds, and Proposition 2 applies to this example. Coupling from a Prescribed Initial State. The cutoff makes ϕφ globally Lipschitz; let LϕL_φ be a Lipschitz constant. Couple a trajectory from any deterministic (X0,Y0,Z0)(X_0,Y_0,Z_0) and initial observed sign s0=ξ1s_0= _1 to the stationary trajectory. Use common reset variables (Bi,ωi)(B_i, _i) from time i=2i=2 onward. The initial observed signs may differ, but the drivers agree after their first common reset. Thus ‖ξi−ξist‖L2≤2r(i−1)/2,i≥1,\| _i- _i st\|_L^2≤ 2r^(i-1)/2, i≥ 1, with the bound at i=1i=1 equal to 22, including when r=0r=0. Put DX,k=‖Xk−Xkst‖L2D_X,k=\|X_k-X_k st\|_L^2 and define DY,k,DZ,kD_Y,k,D_Z,k similarly. The filter equations give DX,k D_X,k ≤(1−η)kDX,0+2η∑i=1k(1−η)k−ir(i−1)/2, ≤(1-η)^kD_X,0+2η _i=1^k(1-η)^k-ir^(i-1)/2, DY,k D_Y,k ≤(1−ε)kDY,0+2ε∑i=1k(1−ε)k−ir(i−1)/2. ≤(1- )^kD_Y,0+2 _i=1^k(1- )^k-ir^(i-1)/2. Choose η0 _0 so that 1−η0>r1- _0> r. Geometric summation gives DX,k≤Ce−ηk(DX,0+η),DY,k≤Ce−εk(DY,0+ε).D_X,k≤ Ce^-η k(D_X,0+η), D_Y,k≤ Ce^- k(D_Y,0+ ). The Lipschitz bound also gives DZ,k≤(1−ε)kDZ,0+εLϕ∑j=0k−1(1−ε)k−1−j(DX,j+DY,j).D_Z,k≤(1- )^kD_Z,0+ L_φ _j=0^k-1(1- )^k-1-j(D_X,j+D_Y,j). For 0<ρ≤ρ0<10<ρ≤ _0<1, the fast-to-slow convolution is bounded by Cρ(e−cηk+e−cεk)Cρ(e^-cη k+e^-c k). The slow-to-slow convolution is bounded by Ce−cεkCe^-c k after decreasing c>0c>0 to absorb the factor kεk . Writing Fk=e−cηkF_k=e^-cη k and Sk=e−cεkS_k=e^-c k, we obtain DX,k D_X,k ≤CFkAu, ≤ CF_kA_u, DY,k+DZ,k D_Y,k+D_Z,k ≤C⁡ρ⁡(Fk+Sk)Au+SkAv, ≤ C\ρ(F_k+S_k)A_u+S_kA_v\, Au A_u =DX,0+η,Av=DY,0+DZ,0+ε. =D_X,0+η, A_v=D_Y,0+D_Z,0+ . The bounded stationary coordinates ensure Au2+Av2≤C⁡(1+‖(X0,Y0,Z0)‖2)A_u^2+A_v^2≤ C(1+\|(X_0,Y_0,Z_0)\|^2). Because λ=0λ=0 globally, these imply the exact-coordinate block L2L^2 coupling bounds with V2(z)=‖z‖2V_2(z)=\|z\|^2. Cauchy–Schwarz gives Assumption A6 with V⁡(z)=‖z‖V(z)=\|z\|; both size functions have uniformly bounded stationary expectations. Appendix I From the stationary calculation to finite time The stationary block scales and mean-square bounds below use Assumptions A1–A5 together with the additional localization condition stated next. Constants are independent of 0<η≤η00<η≤ _0 and 0<ρ≤ρ0≤10<ρ≤ _0≤ 1, with ε=ηρ =ηρ. Stationary expectations are taken under the selected joint invariant law πη,ρ _η,ρ. Assumption A7 (Additional block localization). For the selected invariant laws in Assumption A4, the local map λ has a C2C^2 extension to a neighborhood containing the stationary slow states. In the resulting exact coordinates u^=x−λ⁡(θ) u=x-λ(θ) and v=θ−θ⋆v=θ-θ , suppose that, for j=1,2,3j=1,2,3, ‖u^‖2j \| u\|^2j ≤Cηj, ≤ Cη^j, ‖v‖2j \|v\|^2j ≤C⁡(εj+η2j), ≤ C( ^j+η^2j), ‖u^−u^‖2j \| u-E u\|^2j ≤Cηj, ≤ Cη^j, ‖v−v‖2j \|v-Ev\|^2j ≤Cεj. ≤ C ^j. (214) In tangent coordinates, ‖Cov⁡(u,v)‖≤Cε\|Cov(u,v)\|≤ C . The coordinate remainder rλ(v)=λ⁡(θ⋆+v)−x⋆−Λvr_λ(v)=λ(θ +v)-x - v satisfies ‖rλ(v)‖≤Cη,‖rλ(v)‖2≤Cη2.E\|r_λ(v)\|≤ Cη, \|r_λ(v)\|^2≤ Cη^2. (215) All constants are independent of sufficiently small positive η,ρη,ρ. Appendices G and H.9 verify this assumption under their respective conditions. The main covariance and bias expansions do not require it. I.1 Stationary Block Scales The exact and tangent tracking errors satisfy u=u^+rλ(v)u= u+r_λ(v). By (214) and (215), ‖u−u‖2≤2‖u^−u^‖2+2‖rλ−rλ‖2=O⁡(η).E\|u-Eu\|^2≤ 2E\| u-E u\|^2+2E\|r_λ-Er_λ\|^2=O(η). The reverse comparison follows from u^=u−rλ(v) u=u-r_λ(v), and the slow coordinate is unchanged. Thus ‖u−u‖2 u-Eu ^2 =O⁡(η), =O(η), ‖v−v‖2 v-Ev ^2 =O⁡(ε), =O( ), (216) ‖⁡(u,v)‖ (u,v) =O⁡(η). =O(η). (217) For the mean estimate, Lemma 2 gives ⁡(z−z⋆)=O⁡(η)E(z-z )=O(η), and (u,v)=T⁡(z−z⋆)(u,v)=T(z-z ) with fixed T. The expected-remainder bound also gives u^=O⁡(η)E u=O(η). If D2λD^2λ is bounded on the segments from θ⋆θ to stationary slow states, Taylor’s formula gives ‖rλ(v)‖≤C‖v‖2\|r_λ(v)\|≤ C\|v\|^2; the second and fourth moments in (214) then imply (215). Using ‖W‖2=‖W−W‖2+‖W‖2E\|W\|^2=E\|W-EW\|^2+\|EW\|^2, we obtain ‖u‖2 u ^2 =O⁡(η), =O(η), ‖v‖2 v ^2 =O⁡(ε+η2). =O( +η^2). (218) The same fast bound holds for u u. The slow η2η^2 term comes from squared bias and may dominate the variance when ε≪η2 η^2. I.2 Initialization and Blockwise Mean-Square Bounds Let zkstz_k^st be a stationary version coupled to a trajectory from a prescribed state z0z_0. Assumption A6 requires ⁡‖zk−zkst‖≤C1+V⁡(z0)(e−cηk+e−cεk),E z_k-z_k^st ≤ C\1+V(z_0)\ (e^-cη k+e^-c k ), (219) where V≥0V≥ 0 and supη,ρπη,ρV(z)<∞ _η,ρE_ _η,ρV(z)<∞. Expectations include the stationary starting state and any coupling of differently initialized Markov inputs. Theorem 3 gives the corresponding tail-average initialization bound. For blockwise mean-square control, assume that the chosen extension of λ is defined at the slow states of both coupled trajectories. Require, in addition, the following exact-coordinate L2L^2 bounds, with Fk=e−cηkF_k=e^-cη k and Sk=e−cεkS_k=e^-c k: ‖u^k−u^kst‖L2 \| u_k- u_k^st\|_L^2 ≤C⁡(Fk+ρSk)Au+ρ⁡(Fk+Sk)Av, ≤ C\(F_k+ρ S_k)A_u+ρ(F_k+S_k)A_v\, ‖vk−vkst‖L2 \|v_k-v_k^st\|_L^2 ≤C⁡ρ⁡(Fk+Sk)Au+(Sk+ρFk)Av, ≤ C\ρ(F_k+S_k)A_u+(S_k+ρ F_k)A_v\, Au2+Av2 A_u^2+A_v^2 ≤C⁡1+V2(z0). ≤ C\1+V_2(z_0)\. Here Au,Av≥0A_u,A_v≥ 0 bound the corresponding L2L^2 initial discrepancies, V2≥0V_2≥ 0, and supη,ρπη,ρV2(z)<∞ _η,ρE_ _η,ρV_2(z)<∞. These bounds must include any contribution from different initial input states; they do not follow from the L1L^1 condition (219). Appendix H.9 verifies them for the nonlinear Markov example. Squaring the block bounds, using (a+b)2≤2a2+2b2(a+b)^2≤ 2a^2+2b^2, and adjusting the positive decay constant give ‖u^k−u^kst‖2 \| u_k- u_k^st\|^2 ≤C⁡1+V2(z0)(e−cηk+ρ2e−cεk), ≤ C\1+V_2(z_0)\(e^-cη k+ρ^2e^-c k), ‖vk−vkst‖2 \|v_k-v_k^st\|^2 ≤C⁡1+V2(z0)(ρ2e−cηk+e−cεk). ≤ C\1+V_2(z_0)\(ρ^2e^-cη k+e^-c k). Combining these estimates with (218) and ‖a+b‖2≤2‖a‖2+2‖b‖2\|a+b\|^2≤ 2\|a\|^2+2\|b\|^2 yields ‖u^k‖2 u_k ^2 ≤Cη+C⁡1+V2(z0)(e−cηk+ρ2e−cεk), ≤ Cη+C\1+V_2(z_0)\ (e^-cη k+ρ^2e^-c k ), (220) ‖vk‖2 v_k ^2 ≤C⁡(ε+η2)+C⁡1+V2(z0)(ρ2e−cηk+e−cεk). ≤ C( +η^2)+C\1+V_2(z_0)\ (ρ^2e^-cη k+e^-c k ). (221) Appendix J Application details and additional simulations J.1 Application: Linear Temporal-Difference Learning with Correction We use temporal-difference learning with correction (TDC) to approximate the discounted value function of a fixed Markov reward process. The state space is 0,1,2,3,4\0,1,2,3,4\. From an interior state, the process stays put with probability 1/21/2 and moves to either neighbor with probability 1/41/4; at a boundary, the missing outward move is replaced by a self-loop. The invariant distribution is uniform. Set γ=0.2,ϕ(s)=1+0.05s,Rk+1=Sk+1=4,Vθ(s)=ϕ(s)θ.γ=0.2, φ(s)=1+0.05s, R_k+1=1\S_k+1=4\, V_θ(s)=φ(s)θ. For the temporal-difference error δk=Rk+1+γϕ(Sk+1)θk−ϕ⁡(Sk)θk, _k=R_k+1+γφ(S_k+1) _k-φ(S_k) _k, the fast auxiliary variable w and the slow value parameter θ follow wk+1 w_k+1 =wk+η⁡[δk−ϕ⁡(Sk)wk]ϕ(Sk), =w_k+η\,[ _k-φ(S_k)w_k]φ(S_k), (222) θk+1 _k+1 =θk+ε[δkϕ(Sk)−γϕ(Sk+1)ϕ(Sk)wk],ε=ηρ. = _k+ \,[ _kφ(S_k)-γφ(S_k+1)φ(S_k)w_k], =ηρ. Both updates use the old parameters (wk,θk)(w_k, _k) and the same observed transition Yk=(Sk,Sk+1)Y_k=(S_k,S_k+1). Averaged Objective and Equilibrium. For a transition y=(s,s′)y=(s,s ), define Cy=ϕ(s)2,Ay=ϕ(s)ϕ(s)−γϕ(s′),Dy=γϕ(s)ϕ(s′),by=ϕ(s)s′=4.C_y=φ(s)^2, A_y=φ(s)\φ(s)-γφ(s )\, D_y=γφ(s)φ(s ), b_y=φ(s)1\s =4\. Thus Ay=Cy−DyA_y=C_y-D_y and, in the coordinate order z=(w,θ)z=(w,θ), G⁡(y,z)=Jyz+by(11),Jy=(−Cy−Ay−Dy−Ay).G(y,z)=J_yz+b_y 11, J_y= pmatrix-C_y&-A_y\\ -D_y&-A_y pmatrix. (223) Let ν⁡(s,s′)=P⁡(s,s′)/5ν(s,s )=P(s,s )/5 on the 13 admissible transitions, and use bars for ν-averages. Direct summation gives C¯=1.215,A¯=0.9721,b¯=0.2375. C=1.215, A=0.9721, b=0.2375. The fast equilibrium is λ⁡(θ)=(b¯−A¯θ)/C¯λ(θ)=( b- Aθ)/ C. Substituting it into the averaged slow update gives g⁡(θ)=A¯C¯(b¯−A¯θ)=−ℒ′(θ),ℒ⁡(θ)=(b¯−A¯θ)22C¯.g(θ)= A C( b- Aθ)=-L (θ), (θ)= ( b- Aθ)^22 C. Here ℒL is one half of the mean-squared projected Bellman error for this scalar feature. Its unique minimizer and the joint equilibrium are θ⋆=b¯A¯≈0.2443164284,z⋆=(0,θ⋆).θ = b A≈ 0.2443164284, z =(0,θ ). (224) The fast and reduced slow Jacobians are −C¯<0- C<0 and S=−A¯2/C¯<0S=- A^2/ C<0. Verification of the Assumptions. The transition-pair driver is exogenous, with kernel K((i,j),(u,v))=u=jP(j,v).K((i,j),(u,v))=1\u=j\P(j,v). Although this kernel has zeros, it satisfies a finite-step Doeblin condition. Indeed, P4(j,u)≥1/256P^4(j,u)≥ 1/256 for every pair of states, so K5((i,j),(u,v))=P4(j,u)P(u,v)≥5256ν(u,v).K^5((i,j),(u,v))=P^4(j,u)P(u,v)≥ 5256ν(u,v). The updates are affine with globally bounded derivatives. Since 1≤ϕ⁡(s)≤1.21≤φ(s)≤ 1.2, the global block-contraction constants in Proposition 4 can be chosen as a=1,b=1.2,c=0.288,d=0.76,bc=0.3456<ad=0.76.a=1, b=1.2, c=0.288, d=0.76, bc=0.3456<ad=0.76. These follow from Cy≥1C_y≥ 1, 0<Ay≤1.20<A_y≤ 1.2, Dy≤0.288D_y≤ 0.288, and Ay≥0.76A_y≥ 0.76. Together with the global affine fast equilibrium and the two negative averaged Jacobians, Proposition 4 verifies Assumptions A1–A5, including existence, uniqueness, and uniform fourth-moment localization of the stationary law for sufficiently small steps. For the IID control, transitions are drawn independently from the same distribution ν. This preserves the averaged objective and equilibrium, while removing serial dependence; these IID pairs are not required to form a contiguous state trajectory. Theoretical Bias Coefficients. The Poisson equation for the finite-state kernel K determines the centered corrector ⁡(y,z)U(y,z). Because G is affine, D2G¯=0D^2 G=0, and the leading bias contains only the direct Markov response: rρ=∑y,y′ν⁡(y)K(y,y′)Dz(y′,z⋆)DρG(y,z⋆),b⁡(ρ)=−J−1rρ=b0+ρb1.r_ρ= _y,y ν(y)K(y,y )D_zU(y ,z )D_ρG(y,z ), b(ρ)=-J^-1r_ρ=b_0+ρ b_1. For this instance, b0θ≈0.00624921b_0^θ≈ 0.00624921 and b1θ≈−0.40166550b_1^θ≈-0.40166550. Consequently b2=0b_2=0 for this application. This example tests the direct contribution of Markov dependence to stationary bias. The nonlinear curvature contribution is illustrated by the exact model in Appendix H. In the IID control, stationarity and independence give G¯(z)=0 G(Ez)=0; therefore the stationary mean equals z⋆z exactly throughout the verified small-step regime. Optional Analytic Characterization of Stationary Bias. The affine structure also permits a separate analytic characterization of the stationary mean. This characterization is not needed to compute the empirical bias in (26). Write zk=(wk,θk)⊤z_k=(w_k, _k) and define ξy=G⁡(y,z⋆),My=I2+ηDρJy,Dρ=diag⁡(1,ρ). _y=G(y,z ), M_y=I_2+η D_ρJ_y, D_ρ=diag(1,ρ). The error recursion is zk+1−z⋆=MYk(zk−z⋆)+ηDρξYk.z_k+1-z =M_Y_k(z_k-z )+η D_ρ _Y_k. For each of the 13 admissible transitions y, let hy=πη,ρ[(zk−z⋆)Yk=y].h_y=E_ _η,ρ [(z_k-z )1\Y_k=y\ ]. Stationarity and the transition-pair kernel K give hy′=∑yK⁡(y,y′)Myhy+ην(y)Dρξy.h_y = _yK(y,y ) \M_yh_y+ην(y)D_ρ _y \. (225) These are 26 scalar linear equations. Their solution determines the stationary bias through Bη,ρθ=(01)∑yhy.B^θ_η,ρ= pmatrix0&1 pmatrix _yh_y. This identity retains the full dependence on the step sizes. It provides an additional theoretical characterization for the affine example; the empirical estimates are obtained from simulated trajectories. J.2 Simulation Protocol We simulate the TDC model in (222) and the nonlinear model in Appendix H. All data are synthetic. The TDC target is θ⋆θ in (224); the nonlinear model has target Z⋆=0Z =0. These targets are specified by the averaged problems. Each setting uses R=256R=256 independent replications. Within a replication, extrapolation levels share the same input sequence. For a given method, let ϑ^(r) ^(r) denote its estimate in replication r: either a post-burn time average or a weighted combination of such averages. With target ϑ⋆ , define ϑ^¯=1R∑r=1Rϑ^(r). = 1R _r=1^R ^(r). We compute the empirical bias, sample standard deviation, and RMSE as b b =ϑ^¯−ϑ⋆, = - , s^2 s^2 =1R−1∑r=1R(ϑ^(r)−ϑ^¯)2, = 1R-1 _r=1^R ( ^(r)- )^2, RMSE RMSE =1R∑r=1R(ϑ^(r)−ϑ⋆)21/2. = \ 1R _r=1^R ( ^(r)- )^2 \^1/2. All three quantities are computed from the simulated replication-level estimates. For TDC, an approximate pointwise 95% confidence interval for the bias of the finite-run estimator is b^±tR−1,0.975s^R, b\ ±\ t_R-1,0.975 s R, where tR−1,0.975t_R-1,0.975 is the corresponding Student-t quantile. For normalized bias, both the estimate and the interval endpoints are divided by η. Standard-deviation and RMSE intervals, and the nonlinear model’s bias intervals, use the percentile bootstrap with 2,000 resamples of entire replications. RMSE ratios use paired resampling of methods. All intervals are pointwise 95% intervals. Dependent iterations within a trajectory are not treated as independent observations. For a finite burn-in, empirical bias includes any remaining initialization effect. The confidence intervals quantify Monte Carlo uncertainty; they do not account for this systematic transient error. Stationary-bias interpretation therefore also requires adequate burn-in. Extrapolation Rules and Update Budgets. The base fast step is h, and level j uses ηj=h2−j _j=h2^-j and εj=ηjp _j= _j^p, with p=6/5p=6/5. The methods are ordinary averaging and extrapolation targeting =1G=\1\, 1,2\1,2\, or 1,p\1,p\. For m=1+||m=1+|G| levels, the weights solve ∑j=0m−1wj=1,∑j=0m−1wj2−jγ=0(γ∈). _j=0^m-1w_j=1, _j=0^m-1w_j2^-jγ=0 (γ ). We apply the weights to post-burn tail averages, rather than to individual terminal iterates. A recursion update advances one complete state at one level. With total budget B, each of m levels receives N=⌊B/m⌋N= B/m updates, including its burn-in; at most m−1m-1 updates are unused. Shared random inputs do not reduce this accounting. Equal update counts are not a claim of equal wall-clock cost. Exact-model Runs. We use q∈0,0.5q∈\0,0.5\ and h∈0.1,0.05,0.025,0.0125h∈\0.1,0.05,0.025,0.0125\, starting the iterates at zero and the sign chain in its uniform invariant distribution. Level j discards Kj=⌈20/εj⌉K_j= 20/ _j updates. The fixed-ratio study sets ρ=0.25ρ=0.25 and uses 240,000 updates per replication. The path-bias study uses 2,400,000 updates per level, whereas the separate equal-budget comparison fixes B=240,000B=240,000 total updates per method. The longer path-bias runs were chosen after a pilot showed that the smallest corrected bias was close to Monte Carlo uncertainty; the step sizes, coefficients, and seeds were unchanged. These long runs are not used as the equal-budget comparison. TDC Runs. We initialize w0=θ0=0w_0= _0=0 and draw S0S_0 uniformly. The fixed-ratio study uses ρ∈0.1,0.25ρ∈\0.1,0.25\ and η∈0.1,0.05,0.025,0.0125η∈\0.1,0.05,0.025,0.0125\, with 256 independent replications for each Markov and IID transition-pair setting. Each replication retains 2202^20 updates after discarding K=⌈80/ε⌉K= 80/ updates. The IID control draws independent pairs from ν; it does not replace the Markov trajectory by overlapping pairs of IID states. The TDC budget study uses ηj=h2−j _j=h2^-j and εj=ηjp _j= _j^p, with p=6/5p=6/5, h=0.05h=0.05, and total update budgets B∈217,219,221B∈\2^17,2^19,2^21\. We compare ordinary averaging with the extrapolation rules targeting =1G=\1\, 1,2\1,2\, and 1,p\1,p\ defined above. For m levels, every level runs for N=⌊B/m⌋N= B/m updates and discards the common burn-in K=⌈80/minj⁡εj⌉K= 80/ _j _j . All levels share a Markov trajectory and average the same post-burn iteration indices. The total cost is mNmN, including all mKmK discarded updates. The error target is θ⋆θ in (224). Bias intervals use the replication-level Student-t construction; RMSE intervals and ratios use 2,000 bootstrap resamples, paired across methods for ratios. As a simulation-based burn-in diagnostic, we compare K=64,000K=64,000 and K=128,000K=128,000 at ρ=0.1ρ=0.1 and η=0.0125η=0.0125, keeping the retained length and replication seeds fixed. We compute the change in the estimated mean and a paired confidence interval across replications. Changing the burn-in also changes the retained stochastic window, so this comparison is a sensitivity diagnostic rather than an exact measurement of initialization error. J.3 Results All reported bias estimates are computed from simulated post-burn averages across 256 independent replications, subtracting the corresponding equilibrium target. The simulations do not use a stationary linear-system solution or an exact stationary-bias formula. The protocol and uncertainty calculations are described in Appendix J.2. Temporal-Difference Learning. Figure 1 reports the empirical normalized bias under Markov and independent-pair sampling, together with the bias and standard deviation of the equal-budget estimators. Figure 2 reports the corresponding RMSE relative to θ⋆θ . This comparison accounts jointly for residual bias and sampling variability when the total number of updates, including burn-in across all extrapolation levels, is fixed. Figure 2: TDC estimation error at equal total update budgets. Points show empirical RMSE relative to θ⋆θ for ordinary averaging and three extrapolation rules. The step sizes are ηj=h2−j _j=h2^-j and εj=ηjp _j= _j^p, with h=0.05h=0.05 and p=6/5p=6/5. Each budget includes burn-in and all updates across levels. Error bars are pointwise 95% bootstrap intervals based on 256 independent replications and 2,000 bootstrap resamples. Solid segments connect empirical points. R labels specify the cancelled powers. Nonlinear Toy Example. We simulate the model of Appendix H, whose equilibrium target is Z⋆=0Z =0. Figure 3 compares empirical normalized bias with the theoretical leading coefficient at a fixed step-size ratio. Figure 4 separates the long-run bias study from the comparison of RMSE at equal total update budgets. The latter comparison includes the sampling variability introduced by combining estimates from several levels. Figure 3: Empirical bias of the nonlinear toy model at ρ=ε/η=0.25ρ= /η=0.25, for q=0q=0 and q=0.5q=0.5. Points show estimated bias relative to Z⋆=0Z =0, divided by η. Dashed lines show the theoretical leading coefficient cq(1/2+ρ/(1+ρ))c_q (1/2+ρ/(1+ρ) ), where cq=(1+q)/(1−q)c_q=(1+q)/(1-q). Error bars are pointwise 95% bootstrap intervals based on 256 independent replications and 2,000 bootstrap resamples. Each replication uses 240,000 updates, including burn-in. Figure 4: Empirical extrapolation results for the nonlinear toy model along ε=η6/5 =η^6/5. Columns correspond to q=0q=0 and q=0.5q=0.5. Top: estimated signed bias from runs using 2,400,000 updates per level, including burn-in. The vertical scale is logarithmic away from zero and linear for |bias|≤10−5|bias|≤ 10^-5. Bottom: empirical RMSE at a total budget of 240,000 updates per method, including burn-in across all levels. The top and bottom rows therefore represent distinct experiments with different averaging lengths. All error bars are pointwise 95% bootstrap intervals based on 256 independent replications and 2,000 bootstrap resamples. Solid segments connect empirical points. R labels specify the cancelled powers, with p=6/5p=6/5. Appendix K Additional related work Classical two-timescale stochastic approximation with decreasing steps uses singular ODEs, averaging, and weak convergence [2, 18]. Finite-sample and limit results cover linear Markov recursions, nonlinear convergence rates, lock-in, and controlled-Markov central limit theory [4, 13, 6, 10]. Their primary target is convergence as the steps vanish over time, rather than the stationary mean displacement produced by two fixed steps. Kaledin et al. [13] give finite-time mean-square bounds for linear Markov recursions with constant, piecewise-constant, and diminishing step sizes. Our distinction is the nonlinear stationary mean-bias coefficient expansion, rather than constant-step analysis alone. The closest constant-step comparison is the linear TTSA analysis of Kwon et al. [15], which treats stationarity, finite-time behavior, mean and covariance, and joint extrapolation. Linear tracking error can already contain a scale such as ε2/η ^2/η. We address the nonlinear question of how stationary covariance moves the invariant mean. Expanding the covariance response in ρ produces the ε2/η ^2/η correction, and curvature converts it into mean bias. Fast–slow cross-covariance is one source, but not the only possible source. Linear recursions with centered additive noise often have invariant mean exactly at the ODE equilibrium. A quadratic readout provides a simple mechanism for converting singular covariance into visible mean bias.

Instant research discovery