Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
100,976 characters · 17 sections · 113 citation commands
Yurinskii's Coupling for Martingales
\footnotetext[1]{ Department of Operations Research and Financial Engineering, Princeton University } \footnotetext[2]{ Department of Statistics, University of California, Davis } \footnotetext[3]{ Statistical Laboratory, University of Cambridge } \let\thefootnote\relax \footnotetext[1]{ \textsuperscript{*}Corresponding author: [email removed]} } \addtocounter{footnote}{-1}\let\thefootnote\svthefootnote
\setcounter{page}{0}\thispagestyle{empty}
\noindentKeywords: coupling, strong approximation, mixingales, martingales, dependent data, Gaussian mixture approximation, time series, empirical processes, uniform inference, series estimation, local polynomial estimation, central limit theorems.
Yurinskii's coupling yurinskii1978error has proven to be an important theoretical tool for developing non-asymptotic distributional approximations in mathematical statistics and applied probability. For a sum $S$ of $n$ independent zero-mean $d$-dimensional random vectors, this coupling technique constructs (on a suitably enlarged probability space) a zero-mean $d$-dimensional Gaussian vector $T$ which has the same covariance matrix as $S$ and which is close to $S$ in probability, bounding the discrepancy $\|S-T\|$ as a function of $n$, $d$, the choice of norm, and some features of the underlying distribution. See, for example, pollard2002user for a textbook introduction, and csorgo1981strong and Lindvall_1992_Book for background references.
When compared to other coupling approaches, such as the celebrated Hungarian construction komlos1975approximation or Zaitsev's coupling zaitsev1987estimates,zaitsev1987gaussian, Yurinskii's approach stands out for its simplicity, robustness, and wider applicability, while also offering tighter couplings in some applications (see below for more discussion and examples). These features have led many scholars to use Yurinskii's coupling to study the distributional properties of high-dimensional statistical procedures in a variety of settings, often with the end goal of developing uncertainty quantification or hypothesis testing methods. For example, in recent years, Yurinskii's coupling has been used to construct Gaussian approximations for the suprema of empirical processes chernozhukov2014gaussian; to establish distribution theory for non-Donsker stochastic $t$-processes generated in nonparametric series regression belloni2015some; to prove distributional approximations for high-dimensional $\ell_p$-norms biau2015high; to develop distribution theory for vector-valued martingales belloni2018high,li2020uniform; to derive a law of the iterated logarithm for stochastic gradient descent optimization methods anastasiou2019normal; to establish uniform distributional results for nonparametric high-dimensional quantile processes belloni2019conditional; to develop distribution theory for non-Donsker stochastic $t$-processes generated in partitioning-based series regression cattaneo2020large; to deduce Bernstein--von Mises theorems in high-dimensional settings ray2021bernstein; and to develop distribution theory for non-Donsker U-processes based on dyadic network data cattaneo2024uniform. There are also many other early applications of Yurinskii's coupling: dudley1983invariance and dehling1983limit establish invariance principles for Banach space-valued random variables, and lecam1988 and sheehy1992uniform obtain uniform Donsker results for empirical processes, to name just a few.
This paper presents a new Yurinskii coupling which encompasses and improves upon all of the results previously available in the literature, offering four new primary features:
Closest to our work are the recent paper by li2020uniform and the unpublished manuscript by belloni2018high, which both investigated distribution theory for martingale data using Yurinskii's coupling and related methods. Specifically, li2020uniform established a Gaussian $\ell_2$-norm Yurinskii coupling for mixingales and martingales under the assumption that the covariance structure has a minimum eigenvalue bounded away from zero. As formally demonstrated in this paper (see Section (ref)), such eigenvalue assumptions can be prohibitively strong in practically relevant applications. In contrast, our Yurinskii coupling does not impose any restrictions on covariance degeneracy (ref), in addition to offering several other new features not present in li2020uniform, including (ref), (ref), (ref), and applicability to general $\ell_p$-norms. In addition, we correct a slight technical inaccuracy in their proof relating to the derivation of bounds in probability (see Remark (ref)).
belloni2018high did not establish a Yurinskii coupling for martingales, but rather a central limit theorem for smooth functions of high-dimensional martingales using the celebrated second-order Lindeberg method chatterjee2006generalization, explicitly accounting for covariance degeneracy. As a consequence, their result could be leveraged to deduce a Yurinskii coupling for martingales with additional, non-trivial technical work (see \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)} for details). Nevertheless, a Yurinskii coupling derived from belloni2018high would not feature (ref), (ref), (ref), or general $\ell_p$-norms, as our results do. We discuss further the connections between our work and the related literature in the upcoming sections, both when introducing our main theoretical results and when presenting examples and statistical applications.
The most general coupling result of this paper (Theorem (ref)) is presented in Section (ref), where we also specialize it to a slightly weaker yet more user-friendly formulation (Proposition (ref)). Our Yurinskii coupling for approximate martingales is a strict generalization of all previous Yurinskii couplings available in the literature, offering a Gaussian mixture strong approximation for approximate martingale vectors in $\ell_p$-norm, with an improved rate of approximation when the third moments of the data are negligible, making no assumptions on the spectrum of the data covariance matrix. A key technical innovation underlying the proof of Theorem (ref) is that we explicitly account for the possibility that the minimum eigenvalue of the variance may be zero, or that its lower bound may be unknown, with the argument proceeding using a carefully tailored regularization. Establishing a coupling to a Gaussian mixture distribution is achieved by an appropriate conditioning argument, leveraging a conditional version of Strassen's theorem \citetext{chen2020jackknife; monrad1991nearby}, along with some related technical work detailed in \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)}. A third-order coupling is obtained via a modification of a standard smoothing technique for Borel sets from classical versions of Yurinskii's coupling \iftoggle{aos}{ (see Lemma SA.2 in the supplementary material cattaneo2025yurinskiisupplement)}{ (see Lemma (ref) in the appendix)}, enabling improved approximation errors whenever third moments are negligible.
In Proposition (ref), we explicitly tune the parameters of the aforementioned regularization to obtain a simpler, parameter-free version of Yurinskii's coupling for approximate martingales, again offering Gaussian mixture coupling distributions and an improved third-order approximation. This specialization of our main result takes an agnostic approach to potential singularities in the data covariance matrix and, as such, may be improved in specific applications where additional knowledge of the covariance structure is available. Section (ref) also presents some further refinements when additional structure is imposed, deriving Yurinskii couplings for mixingales, martingales, and independent data as Corollaries (ref), (ref), and (ref), respectively. We take the opportunity to discuss and correct in Remark (ref) a technical issue which is often neglected pollard2002user, li2020uniform when using Yurinskii's coupling to derive bounds in probability. Section (ref) presents a stylized example portraying the relevance of our main technical results in the context of canonical factor models, illustrating the importance of each of our new Yurinskii coupling features (ref)-- (ref).
Section (ref) considers a substantive application of our main results: strong approximation of martingale empirical processes. We begin with the motivating example of canonical kernel density estimation, demonstrating how Yurinskii's coupling can be applied, and showing in Lemma (ref) why it is essential that we do not place any conditions on the minimum eigenvalue of the variance matrix (ref). We then present a general-purpose strong approximation for martingale empirical processes in Proposition (ref), combining classical results in the empirical process literature van1996weak with our coupling from Corollary (ref). This statement appears to be the first of its kind for martingale data, and when specialized to independent (and not necessarily identically distributed) data, it is shown to be superior to the best known comparable strong approximation result available in the literature berthet2006revisiting. Our improvement comes from using Yurinskii's coupling for the $\ell_\infty$-norm, where berthet2006revisiting apply Zaitsev's coupling zaitsev1987estimates, zaitsev1987gaussian with the larger $\ell_2$-norm.
Section (ref) further illustrates the applicability of our results through two examples in nonparametric regression estimation. Firstly, we deduce strong approximations for partitioning-based least squares series estimators with time series data, applying Corollary (ref) directly and additionally imposing only a mild mixing condition on the regressors. We show that our Yurinskii coupling for martingale vectors delivers the same distributional approximation rate as the best known result for independent data, and discuss how this can be leveraged to yield a feasible statistical inference procedure. We also show that if the residuals have vanishing conditional third moment, an improved rate of Gaussian approximation can be established. Secondly, we deduce a strong approximation for local polynomial estimators with time series data, using our result on martingale empirical processes (Proposition (ref)) and again imposing a mixing assumption. Appealing to empirical process theory is essential here as, in contrast with series estimators, local polynomials do not possess certain additive separability properties. The bandwidth restrictions we require are relatively mild, and, as far as we know, they have not been improved upon even with independent data.
Section (ref) concludes the paper. Appendix (ref) demonstrates how our coupling results can be used to derive distributional Gaussian approximations (central limit theorems) for possibly high-dimensional martingale vectors (Proposition (ref)). This result complements a recent literature on probability and statistics studying the same problem but with independent data buzun2022strong,lopes2022central, chernozhukov2023nearly,kock2024remark. We also present a version of this result employing a covariance estimator (Proposition (ref)), enabling the construction of valid high-dimensional confidence sets via a Gaussian multiplier bootstrap. \iftoggle{aos}{ Finally we present some further results on applications of our theory to deriving distributional approximations for $\ell_p$-norms of high-dimensional martingale vectors in Appendix (ref). }
All proofs are collected in \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)}, where we also include other technical lemmas of potential independent interest \iftoggle{aos}{ , alongside some further results on distributional approximations for $\ell_p$-norms of high-dimensional martingale vectors }.
We write $\|x\|_p$ for $p\in[1,\infty]$ to denote the $\ell_p$-norm if $x$ is a (possibly random) vector or the induced operator $\ell_p$--$\ell_p$-norm if $x$ is a matrix. For $X$ a real-valued random variable and an Orlicz function $\psi$, we use ${\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_\psi$ to denote the Orlicz $\psi$-norm van1996weak and ${\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_p$ for the $L^p(\ensuremath{\mathbb{P}})$ norm where $p\in [1,\infty]$. For a matrix $M$, we write $\|M\|_{\max}$ for the maximum absolute entry and $\|M\|_\ensuremath{\mathrm{F}}$ for the Frobenius norm. We denote positive semi-definiteness by $M \succeq 0$ and write $I_d$ for the $d \times d$ identity matrix.
For scalar sequences $x_n$ and $y_n$, we write $x_n \lesssim y_n$ if there exists a positive constant $C$ such that $|x_n| \leq C |y_n|$ for sufficiently large $n$. We write $x_n \asymp y_n$ to indicate both $x_n \lesssim y_n$ and $y_n \lesssim x_n$. Similarly, for random variables $X_n$ and $Y_n$, we write $X_n \lesssim_\ensuremath{\mathbb{P}} Y_n$ if for every $\varepsilon > 0$ there exists a positive constant $C$ such that $\ensuremath{\mathbb{P}}(|X_n| \geq C |Y_n|) \leq \varepsilon$, and write $X_n \to_\ensuremath{\mathbb{P}} X$ for limits in probability. For real numbers $a$ and $b$ we use $a \vee b = \max\{a,b\}$. We write $\kappa \in \ensuremath{\mathbb{N}}^d$ for a multi-index, where $d \in \ensuremath{\mathbb{N}} = \{0, 1, 2, \ldots\}$, and define $|\kappa| = \sum_{j=1}^d \kappa_j$, along with $\kappa! = \prod_{j=1}^{d} \kappa_j !$, and $x^\kappa = \prod_{j=1}^d x_j^{\kappa_j}$ for $x \in \ensuremath{\mathbb{R}}^d$.
Since our results concern couplings, some statements must be made on a new or enlarged probability space. We omit the details of this for clarity of notation, but technicalities are handled by the Vorob'ev--Berkes--Philipp Theorem dudley1999uniform.
We begin with our most general result: an $\ell_p$-norm Yurinskii coupling for a sum of vector-valued approximate martingale differences to a Gaussian mixture-distributed random vector. The general result is presented in Theorem (ref), while Proposition (ref) gives a simplified and slightly weaker version which is easier to use in many applications. We then further specialize Proposition (ref) to three scenarios with successively stronger assumptions, namely mixingales, martingales, and independent data, in Corollaries (ref), (ref), and (ref) respectively. In each case we allow for possibly random quadratic variations (cf.\ mixing convergence), thereby establishing Gaussian mixture couplings in the general setting. In Remark (ref) we comment on and correct an often overlooked technicality relating to the derivation of bounds in probability from Yurinskii's coupling. As a first illustration of the power of our generalized $\ell_p$-norm Yurinskii coupling, we present in Section (ref) a simple factor model example relating to all three of the aforementioned scenarios, discussing further how our contributions are related to the existing literature.
This theorem offers four novel contributions to the literature on coupling theory and strong approximation, as discussed in the introduction. Firstly (ref), it allows for approximate vector-valued martingales, with the variables $\tilde X_i$ forming martingale differences with respect to $\ensuremath{\mathcal{H}}_i$ by construction, and $U$ quantifying the associated martingale approximation error. Such martingale approximation techniques for sequences of dependent random vectors are well established and have been used in a range of scenarios: see, for example, wu2004martingale, wu2005nonlinear, dedecker2007weak, zhao2008martingale, peligrad2010conditional, atchade2014martingale, cuny2014martingale, magda2018martingale, and references therein. In Section (ref) we demonstrate how this approximation can be established in practice by restricting our general theorem to the special case of mixingales, while the upcoming example in Section (ref) provides an illustration in the context of auto-regressive factor models.
Secondly (ref), Theorem (ref) allows for the resulting coupling variable $T$ to follow a multivariate Gaussian distribution only conditionally, and thus we offer a useful analog of mixing convergence in the context of strong approximation. To be more precise, the random matrix $\sum_{i=1}^{n} V_i$ is the quadratic variation of the constructed martingale $\sum_{i=1}^n \tilde X_i$, and we approximate it using the $\ensuremath{\mathcal{H}}_0$-measurable random matrix $\Sigma$. This yields the coupling variable $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$, which can alternatively be written as $T=\Sigma^{1/2} Z$ with $Z \sim \ensuremath{\mathcal{N}}(0,I_d)$ independent of $\ensuremath{\mathcal{H}}_0$. The errors in this quadratic variation approximation are accounted for by the terms $\ensuremath{\mathbb{P}}(\Omega \npreceq M)$, $\delta_p(M, \eta)$ and $\varepsilon_p(M, \eta)$, utilizing a regularization argument through the free matrix parameter $M$. If a non-random $\Sigma$ is used, then $T$ is unconditionally Gaussian, and one can take $\ensuremath{\mathcal{H}}_0$ to be the trivial $\sigma$-algebra. As demonstrated in our proof, our approach to establishing a mixing approximation is different from naively taking an unconditional version of Yurinskii's coupling and applying it conditionally on $\ensuremath{\mathcal{H}}_0$, which will not deliver the same coupling as in Theorem (ref) for a few reasons. To begin with, we explicitly indicate in the conditions of Theorem (ref) where conditioning is required. Next, our error of approximation is given unconditionally, involving only marginal expectations and probabilities. Finally, we provide a rigorous account of the construction of the conditionally Gaussian coupling variable $T$ via a conditional version of Strassen's theorem \citetext{chen2020jackknife; monrad1991nearby}. Section (ref) illustrates how a strong approximation akin to mixing convergence can arise when the data forms an exact martingale, and Section (ref) gives a simple example relating to factor modeling in statistics and data science.
As a third contribution to the literature (ref), and of particular importance for applications, Theorem (ref) makes no requirements on the minimum eigenvalue of the quadratic variation of the approximating martingale sequence. Instead, our proof technique employs a careful regularization scheme designed to account for any such exact or approximate rank degeneracy in $\Sigma$. This capability is fundamental in some applications, a fact which we illustrate in Section (ref) by demonstrating the significant improvements in strong approximation errors delivered by Theorem (ref) relative to those obtained using prior results in the literature.
Finally (ref), Theorem (ref) gives a third-order strong approximation alongside the usual second-order version considered in all prior literature. More precisely, we observe that an analog of the term $\beta_{p,2}$ is present in the classical Yurinskii coupling and comes from a Lindeberg telescoping sum argument, replacing random variables by Gaussians with the same mean and variance to match the first and second moments. Whenever the third conditional moments of $\tilde X_i$ are negligible (quantified by $\pi_3$), this moment-matching argument can be extended to third-order terms, giving a new quantity $\beta_{p,3}$. At this level of generality, it is not possible to obtain explicit bounds on $\pi_3$ because we make no assumptions on the relationship between the data $X_i$ and the $\sigma$-algebras $\ensuremath{\mathcal{H}}_i$ (and therefore the variables $\tilde X_i$ resulting from the martingale approximation). However, if $X_1, \ldots, X_n$ form martingale differences with respect to $\ensuremath{\mathcal{H}}_0, \ldots, \ensuremath{\mathcal{H}}_n$, then $\tilde X_i = X_i$ almost surely (see Section (ref)). In this setting, assuming that $\ensuremath{\mathbb{E}} \big[ X_i^\kappa \mid \ensuremath{\mathcal{H}}_{i-1} \big] = 0$ for each multi-index $\kappa$ with $|\kappa| = 3$ (e.g.\ if the data is conditionally symmetrically distributed around zero), then using $\beta_{p,3}$ rather than $\beta_{p,2}$ can give smaller coupling approximation errors in (ref). Such a refinement can be viewed as a strong approximation counterpart to classical Edgeworth expansion methods, and we illustrate this phenomenon in our upcoming applications to nonparametric inference (Section (ref)).
The result in Theorem (ref) is given in a somewhat implicit manner, involving infima over the free parameters $t > 0$ and $M \succeq 0$, and it is not clear how to compute these in general. In the upcoming Proposition (ref), we set $M = \nu^2 I_d$ and approximately optimize over $t > 0$ and $\nu > 0$, resulting in a simplified and slightly weaker version of our main general result. In specific applications, where there is additional knowledge of the quadratic variation structure, other choices of regularization schemes may be more appropriate. Nonetheless, the choice $M = \nu^2 I_d$ leads to arguably the principal result of our work, due to its simplicity and utility in statistical applications. For convenience, define the functions $\phi_p : \{1, 2, \ldots\} \to \ensuremath{\mathbb{R}}$, for $p \in [0, \infty]$, by
With $Z \sim \ensuremath{\mathcal{N}}(0, I_d)$ and $t > 0$, these functions satisfy $\ensuremath{\mathbb{P}}( \|Z\|_p > t ) \leq \ensuremath{\mathbb{E}}[\|Z\|_p] / t \leq \phi_p(d) / t$ (see \iftoggle{aos}{ Lemma SA.4 in the supplementary material cattaneo2025yurinskiisupplement}{ Lemma (ref) in the appendix}).
Proposition (ref) makes clear the potential benefit of a third-order coupling when $\pi_3 = 0$, as in this case the bound features $\beta_{p,3}^{1/4}$ rather than $\beta_{p,2}^{1/3}$. If $\pi_3$ is small but non-zero, an analogous result can easily be derived by adjusting the optimal choices of $t$ and $\nu$, but we omit this for clarity of notation. In applications (see Section (ref)), this reduction of the exponent can provide a significant improvement in terms of the dependence of the bound on the sample size $n$, the dimension $d$, and other problem-specific quantities. When using our results for strong approximation, it is usual to set $p = \infty$ to bound the maximum discrepancy over the entries of a vector (to construct uniform confidence sets, for example). In this setting, we have that $\phi_\infty(d) = \sqrt{2 \log 2d}$ has a sub-Gaussian slow-growing dependence on the dimension. The remaining term depends on $\ensuremath{\mathbb{E}}[\|\Omega\|_2]$ and requires that the matrix $\Sigma$ be a good approximation of $\sum_{i=1}^{n} V_i$, while remaining $\ensuremath{\mathcal{H}}_0$-measurable. In some applications (such as factor modeling; see Section (ref)), it can be shown that the quadratic variation $\sum_{i=1}^n V_i$ remains random and $\ensuremath{\mathcal{H}}_0$-measurable even in large samples, giving a natural choice for $\Sigma$.
In the next few sections, we continue to refine Proposition (ref), presenting a sequence of results with increasingly strict assumptions on the dependence structure of the data $X_i$. These allow us to demonstrate the broad applicability of our main results, providing more explicit bounds in settings which are likely to be of special interest. In particular, we consider mixingales, martingales, and independent data, comparing our derived results with those in the existing literature.
In our first refinement, we provide a natural method for bounding the martingale approximation error term $U$. Suppose that $X_i$ form an $\ell_p$-mixingale in $L^1(\ensuremath{\mathbb{P}})$ in the sense that there exist non-negative $c_1, \ldots, c_n$ and $\zeta_0, \ldots, \zeta_n$ such that for all $1 \leq i \leq n$ and $0 \leq r \leq i$,
and for all $1 \leq i \leq n$ and $0 \leq r \leq n-i$,
These conditions are satisfied, for example, if $X_i$ are integrable strongly $\alpha$-mixing random variables mcleish1975invariance, or if $X_i$ are generated by an auto-regressive or auto-regressive moving average process (see Section (ref)), among many other possibilities bradley2005basic. Then, in the notation of Theorem (ref), we have by Markov's inequality that
with $\zeta = 6 \sum_{i=1}^{n} c_i (\zeta_{i} + \zeta_{n-i+1})$. Combining Proposition (ref) with this martingale error bound yields the following result for mixingales.
The closest antecedent to Corollary (ref) is found in li2020uniform, who also considered Yurinskii's coupling for mixingales. Our result improves on this work in the following manner: it removes any requirements on the minimum eigenvalue of the quadratic variation of the mixingale sequence; it allows for general $\ell_p$-norms with $p\in[1,\infty]$; it establishes a coupling to a multivariate Gaussian mixture distribution in general; and it permits third-order couplings (when $\pi_3=0$). These improvements have important practical implications as demonstrated in Section (ref) and Section (ref), where significantly better coupling approximation errors are demonstrated for a variety of statistical applications. On the technical side, our result is rigorously established using a conditional version of Strassen's theorem, a carefully crafted regularization argument, and a third-order Lindeberg method. Furthermore (Remark (ref)), we clarify a technical issue in li2020uniform surrounding the derivation of valid probability bounds for $\|S-T\|_p$.
Corollary (ref) focused on mixingales for simplicity, but, as previously discussed, any method for constructing a martingale approximation $\tilde X_i$ and bounding the resulting error $U$ could be used instead in Proposition (ref) to derive a similar result.
For our second refinement, suppose that $X_i$ form martingale differences with respect to $\ensuremath{\mathcal{H}}_i$. In this case, $\ensuremath{\mathbb{E}}[X_i \mid \ensuremath{\mathcal{H}}_n] = X_i$ and $\ensuremath{\mathbb{E}}[X_i \mid \ensuremath{\mathcal{H}}_0] = 0$, so $U = 0$, and the martingale approximation error term vanishes. Applying Proposition (ref) in this setting directly yields the following result.
The closest antecedents to Corollary (ref) are belloni2018high and li2020uniform, who also (implicitly or explicitly) considered Yurinskii's coupling for martingales. More specifically, li2020uniform established an explicit $\ell_2$-norm Yurinskii coupling for martingales under a strong assumption on the minimum eigenvalue of the martingale quadratic variation, while belloni2018high established a central limit theorem for vector-valued martingale sequences employing the standard second-order Lindeberg method. As such, their proof could be adapted to deduce a Yurinskii coupling for martingales with the help of a conditional version of Strassen's theorem and some additional nontrivial technical work.
Corollary (ref) improves over this prior work as follows. With respect to li2020uniform, our result establishes an $\ell_p$-norm Gaussian mixture Yurinskii coupling for martingales without any requirements on the minimum eigenvalue of the martingale quadratic variation, and permits a third-order coupling if $\pi_3=0$. The first probability bound (ref) in Corollary (ref) gives the same rate of strong approximation as that in Theorem 1 of li2020uniform when $p=2$, with non-random $\Sigma$, and when the eigenvalues of a normalized version of $\Sigma$ are bounded away from zero. In Section (ref) we demonstrate the crucial importance of removing this eigenvalue lower bound restriction in applications involving nonparametric kernel estimators, while in Section (ref) we demonstrate how the availability of a third-order coupling (ref) can give improved approximation rates in applications involving nonparametric series estimators with conditionally symmetrically distributed residual errors. Finally, our technical work improves on li2020uniform in two respects: (i) we employ a conditional version of Strassen's theorem (see \iftoggle{aos}{ Lemma SA.1 in the supplementary material cattaneo2025yurinskiisupplement}{Lemma (ref) in the appendix}) to appropriately handle the conditioning arguments; and (ii) we deduce valid probability bounds for $\|S-T\|_p$, as the following Remark (ref) makes clear.
Technicalities akin to those outlined in Remark (ref) have been both addressed and neglected alike in the prior literature. pollard2002user apparently misses this subtlety, providing an inaccurate bound in probability based on the Yurinskii coupling. li2020uniform seem to make the same mistake in the proof of their Lemma A2, which invalidates the conclusion of their Theorem 1. In contrast, belloni2015some and belloni2019conditional directly provide bounds in $o_\ensuremath{\mathbb{P}}$ instead of $O_\ensuremath{\mathbb{P}}$, circumventing these issues in a manner similar to our approach involving a diverging sequence $R_n$.
To see how this phenomenon applies to our main results, observe that the second-order martingale coupling given as (ref) in Corollary (ref) implies that for any $R_n \to \infty$,
This bound is comparable to that obtained by li2020uniform with $p=2$, albeit with their formulation missing the $R_n$ correction terms. In Section (ref) we discuss further their (amended) result, in the setting of nonparametric series estimation. Our approach using $p = \infty$ obtains superior distributional approximation rates, alongside exhibiting various other improvements such as the aforementioned third-order coupling.
Turning to the comparison with belloni2018high, our Corollary (ref) again offers the same improvements, with the only exception being that the authors did account for the implications of a possibly vanishing minimum eigenvalue. However, their results exclusively concern high-dimensional central limit theorems for vector-valued martingales, and therefore while their findings could in principle enable the derivation of a result similar to our Corollary (ref), this would require additional technical work on their behalf in multiple ways (see \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)}): (i) a correct application of a conditional version of Strassen's theorem \iftoggle{aos}{ (Lemma SA.1 in the supplementary material cattaneo2025yurinskiisupplement)}{(Lemma (ref) in the appendix)}; (ii) the development of a third-order Borel set smoothing technique and associated $\ell_p$-norm moment control \iftoggle{aos}{ (Lemmas SA.2, SA.3, and SA.4)}{(Lemmas (ref), (ref), and (ref))}; (iii) a careful truncation scheme to account for $\Omega\npreceq0$; and (iv) a valid third-order Lindeberg argument \iftoggle{aos}{ (Lemma SA.8)}{(Lemma (ref))}; among others.
As a final refinement, suppose that $X_i$ are independent and zero-mean conditionally on $\ensuremath{\mathcal{H}}_0$, and take $\ensuremath{\mathcal{H}}_i$ to be the filtration generated by $X_1, \ldots, X_i$ and $\ensuremath{\mathcal{H}}_0$ for $1 \leq i \leq n$. Then, taking $\Sigma = \sum_{i=1}^n V_i$ gives $\Omega = 0$, and hence Corollary (ref) immediately yields the following result.
Taking $\ensuremath{\mathcal{H}}_0$ to be trivial, the first inequality (ref) in Corollary (ref) provides an $\ell_p$-norm approximation analogous to that presented in belloni2019conditional. By further restricting to $p=2$, we recover the original Yurinskii coupling as presented in lecam1988 and pollard2002user. Thus, in the independent data setting, our result improves on prior work as follows: (i) it establishes a coupling to a multivariate Gaussian mixture distribution; and (ii) it permits a third-order coupling if $\pi_3=0$.
In this section, we present a simple statistical example of how our improvements over prior coupling results can have important theoretical and practical implications. Consider the stylized factor model
with random variables $L$ taking values in $\ensuremath{\mathbb{R}}^{d \times m}$, $f_i$ in $\ensuremath{\mathbb{R}}^m$, and $\varepsilon_i$ in $\ensuremath{\mathbb{R}}^d$. We interpret $f_i$ as a latent factor variable and $L$ as a random factor loading, with independent (idiosyncratic) disturbances $(\varepsilon_1, \ldots, \varepsilon_n)$. See fan2020statistical, and references therein, for a textbook review of factor analysis in statistics and econometrics.
We employ the above factor model to give a first illustration of the applicability of our main result Theorem (ref), the user-friendly Proposition (ref), and their specialized Corollaries (ref)--(ref). We consider three different sets of conditions to demonstrate the applicability of each of our corollaries for mixingales, martingales, and independent data, respectively. We assume throughout that each $\varepsilon_i$ is zero-mean and finite variance, and that $(\varepsilon_1, \ldots, \varepsilon_n)$ is independent of $L$ and $(f_1, \ldots, f_n)$. Let $\ensuremath{\mathcal{H}}_i$ be the $\sigma$-algebra generated by $L$, $(f_1, \ldots, f_i)$ and $(\varepsilon_1, \ldots, \varepsilon_i)$, with $\ensuremath{\mathcal{H}}_0$ the $\sigma$-algebra generated by $L$ alone.
This simple application to factor modeling gives a preliminary illustration of the power of our main results, encompassing settings which could not be handled by employing Yurinskii couplings available in the existing literature. Even with independent data, we offer new Yurinskii couplings to Gaussian mixture distributions (due to the presence of the common random factor loading $L$), which could be further improved whenever the factors and residuals possess symmetric (conditional) distributions. Furthermore, our results do not impose any restrictions on the minimum eigenvalue of $\Sigma$, thereby allowing for more general factor structures. These improvements are maintained in the martingale, mixingale, and weakly dependent stationary data settings.
In this section, we demonstrate how our main results can be applied to some more substantive problems in statistics. Having until this point studied only finite-dimensional (albeit potentially high-dimensional) random vectors, we now turn our attention to infinite-dimensional stochastic processes. Specifically, we consider empirical processes of the form \[S(f) = \sum_{i=1}^{n} f(X_i), \qquad f \in \ensuremath{\mathcal{F}},\] with $\ensuremath{\mathcal{F}}$ a problem-specific class of real-valued functions, where for each $f \in \ensuremath{\mathcal{F}}$, the variables $f(X_1), \ldots, f(X_n)$ form martingale differences with respect to an appropriate filtration. We construct (conditionally) Gaussian processes $T(f)$ for which upper bounds on the uniform coupling error $\sup_{f \in \ensuremath{\mathcal{F}}} |S(f) - T(f)|$ are precisely quantified. We control the complexity of $\ensuremath{\mathcal{F}}$ using metric entropy under Orlicz norms.
The novel strong approximation results which we present concern the entire martingale empirical process $(S(f):f \in \ensuremath{\mathcal{F}})$, as opposed to just the scalar supremum of the empirical process, $\sup_{f \in \ensuremath{\mathcal{F}}} |S(f)|$. This distinction has been carefully noted by chernozhukov2014gaussian, who studied Gaussian approximation of empirical process suprema in the independent data setting and wrote (p.\ 1565): “A related but different problem is that of approximating whole empirical processes by a sequence of Gaussian processes in the sup-norm. This problem is more difficult than [approximating the supremum of the empirical process].” Indeed, the results we establish in this section are for strong approximations of entire empirical processes by sequences of Gaussian mixture processes in supremum norm, when the data has a martingale difference structure (cf.\ Corollary (ref)). Our results can be further generalized to approximate martingale empirical processes (including mixingale empirical processes; cf.\ Corollary (ref)), but to reduce notation and the technical burden we do not consider this extension.
We begin with a brief study of a canonical example of an empirical process which is non-Donsker (thus precluding the use of uniform central limit theorems) due to the presence of a function class whose complexity increases with the sample size: the kernel density estimator with i.i.d.\ scalar data. We give an overview of our general strategy for strong approximation of stochastic processes via discretization, and show explicitly in Lemma (ref) how it is crucial that we do not impose lower bounds on the eigenvalues of the discretized covariance matrix. Detailed calculations for this section are relegated to \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)} for conciseness.
Let $X_1, \ldots, X_n$ be i.i.d.\ $\ensuremath{\mathrm{Unif}\hspace*{0.2mm}}[0,1]$, take $K(x) = \frac{1}{\sqrt{2 \pi}} e^{-x^2/2}$ the Gaussian kernel and let $h \in (0,1]$ be a bandwidth. Then, for $a \in (0,1/4]$ and $x \in \ensuremath{\mathcal{X}} = [a, 1-a]$ to avoid boundary issues, the kernel density estimator of the true density function $g(x) = 1$ is
Consider establishing a strong approximation for the process $(\hat g(x)-\ensuremath{\mathbb{E}} [ \hat g(x) ] : x\in\ensuremath{\mathcal{X}})$ which is, upon rescaling, non-Donsker whenever the bandwidth decreases to zero in large samples. To match notation with the upcoming general result for empirical processes, set $f_x(u) = \frac{1}{n} (K_h( u - x) - \ensuremath{\mathbb{E}}[K_h( X_i - x)])$ so $S(x) \vcentcolon= S(f_x) = \hat g(x)-\ensuremath{\mathbb{E}} [ \hat g(x) ]$. The next step is standard: a mesh separates the local oscillations of the processes from the finite-dimensional coupling. For $\delta \in (0,1/2)$, set $N = \left\lfloor 1 + \frac{1 - 2a}{\delta} \right\rfloor$ and $\ensuremath{\mathcal{X}}_\delta = (a + (j-1)\delta : 1 \leq j \leq N)$. Letting $T(x)$ be the approximating stochastic process to be constructed, consider the following decomposition:
Writing $S(\ensuremath{\mathcal{X}}_\delta)$ for $\big(S(x) : x \in \ensuremath{\mathcal{X}}_\delta\big)\in \mathbb{R}^N$, and noting that this is a sum of i.i.d.\ random vectors, we apply Corollary (ref) as $\max_{x \in \ensuremath{\mathcal{X}}_\delta} |S(x) - T(x)| = \| S(\ensuremath{\mathcal{X}}_\delta) - T(\ensuremath{\mathcal{X}}_\delta) \|_\infty$. We thus obtain that, for each $\eta > 0$, there exists a Gaussian vector $T(\ensuremath{\mathcal{X}}_\delta)$ with the same covariance matrix as $S(\ensuremath{\mathcal{X}}_\delta)$ satisfying
assuming that $1/h \geq \log 2 N$. By the Vorob'ev--Berkes--Philipp theorem dudley1999uniform, $T(\ensuremath{\mathcal{X}}_\delta)$ extends to a Gaussian process $T(x)$ defined for all $x \in \ensuremath{\mathcal{X}}$ and with the same covariance structure as $S(x)$.
Next, it is not difficult to show by chaining with the Bernstein--Orlicz and sub-Gaussian norms respectively van1996weak that if $\log(N/h) \lesssim \log n$ and $n h \gtrsim \log n$,
Finally, for any sequence $R_n\to\infty$ (Remark (ref)), the resulting bound on the coupling error is
where the mesh size $\delta$ is then optimized to obtain the tightest possible strong approximation. In particular, since $N \lesssim 1/\delta$, setting $\delta \asymp n^{-1/8} h^{5/8} (\log n)^{-1/8}$ yields
which, after standardization by $\sqrt{n h}$, vanishes whenever $R_n (\log n)^3 / (n h^3) \to 0$. This is a more stringent assumption on the bandwidth $h$ than $(\log n) / (n h) \to 0$ imposed by gine2004kernel and cattaneo2024strong when employing a Hungarian construction komlos1975approximation, or $(\log n)^6 / (n h) \to 0$ imposed by chernozhukov2014gaussian when studying in particular the Kolmogorov--Smirnov distance between the scalar suprema. The difference in side restrictions is a result of the specific assumptions imposed and coupling approaches used; see Section (ref) for related discussion.
The discretization strategy outlined above is at the core of the proof strategy for our upcoming Proposition (ref). Since we will consider martingale empirical processes, our proof will rely on Corollary (ref), which, unlike the martingale Yurinskii coupling established by li2020uniform, does not require a lower bound on the minimum eigenvalue of $\Sigma$. Using the simple kernel density example just discussed, we now demonstrate precisely the crucial importance of removing such eigenvalue conditions. The following Lemma (ref) shows that the discretized covariance matrix $\Sigma = n h\operatorname{Var}[S(\ensuremath{\mathcal{X}}_\delta)]$ has exponentially small eigenvalues, which in turn will negatively affect the strong approximation bound if the li2020uniform coupling were to be used instead of the results in this paper.
Figure (ref) shows how the upper bound in Lemma (ref) captures the behavior of the simulated minimum eigenvalue of $\Sigma$. In particular, the smallest eigenvalue decays exponentially fast in the discretization level $\delta$ and the bandwidth $h$. As seen in the calculations above, the coupling rate depends on $\delta / h$, while the bias will generally depend on $h$, implying that both $\delta$ and $h$ must converge to zero to ensure valid statistical inference. In general, this will lead to $\Sigma$ possessing extremely small eigenvalues, rendering strong approximation approaches such as that of li2020uniform ineffective in such scenarios.
The discussion in this section focuses on the strong approximation of the centered process $\hat g(x)-\ensuremath{\mathbb{E}} [ \hat g(x) ]$. In practice, the goal is often rather to approximate $\hat g(x)- g(x)$. The difference between these is captured by the smoothing bias $\ensuremath{\mathbb{E}} [ \hat g(x) ] - g(x)$, which is straightforward to control with $\sup_{x \in \ensuremath{\mathcal{X}}} \big| \ensuremath{\mathbb{E}} [ \hat g(x) ] - g(x) \big| \lesssim \frac{h}{a} e^{-a^2 / (2 h^2)}$. See Section (ref) for further discussion.
We now give our general result on a strong approximation for martingale empirical processes, obtained by applying the first result (ref) in Corollary (ref) with $p=\infty$ to a discretization of the empirical process, as in Section (ref). We then control the increments in the stochastic processes using chaining with Orlicz norms, but note that other tools are available, including generalized entropy with bracketing geer2000empirical and sequential symmetrization rakhlin2015sequential.
A class of functions is said to be pointwise measurable if it contains a countable subclass which is dense under the pointwise convergence topology. For a finite class $\ensuremath{\mathcal{F}}$, write $\ensuremath{\mathcal{F}}(x) = \big(f(x) : f \in \ensuremath{\mathcal{F}}\big)$. Define the set of Orlicz functions
and, for real-valued $Y$, the Orlicz norm ${\vert\kern-0.25ex\vert\kern-0.25ex\vert} Y {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_\psi = \inf \left\{ C > 0: \ensuremath{\mathbb{E}} \left[ \psi(|Y|/C) \leq 1 \right] \right\}$ as in van1996weak.
Proposition (ref) is given in a rather general form to accommodate a range of different settings and applications. In particular, consider the following well-known Orlicz functions.
For these Orlicz functions and when $\Sigma(f, f') = \ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[S(f), S(f')]$ is non-random, the terms involving $\Sigma$ in (ref) can be controlled by the Orlicz $\psi$-norm term; similarly, $J_2$ is bounded by $J_\psi$. Further, $C_\psi$ can be replaced by a universal constant $C$ which does not depend on the parameter $a$. See Section 2.2 in van1996weak for details. If the conditional third moments of $\ensuremath{\mathcal{F}}_\delta(X_i)$ given $\ensuremath{\mathcal{H}}_{i-1}$ are all zero (if $f$ and $X_i$ are appropriately symmetric, for example), then the second inequality in Corollary (ref) can be applied to obtain a tighter coupling inequality; the details of this are omitted for brevity, and the proof would proceed in exactly the same manner.
In general, however, Proposition (ref) allows for a random covariance function, yielding a coupling to a stochastic process that is Gaussian only conditionally. Such a process can equivalently be formally viewed as a mixture of Gaussian processes, writing $T=\Sigma^{1/2} Z$ with an operator square root and where $Z$ is a Gaussian white noise on $\ensuremath{\mathcal{F}}$ independent of $\ensuremath{\mathcal{H}}_0$. This extension is in contrast with much of the existing strong approximation and empirical process literature, which tends to focus on couplings and weak convergence results with marginally Gaussian processes.
A similar approach was taken by berthet2006revisiting, who used a Gaussian coupling due to zaitsev1987estimates,zaitsev1987gaussian along with a discretization method to obtain strong approximations for empirical processes with independent data. They handled fluctuations in the stochastic processes with uniform $L^2$ covering numbers and bracketing numbers where we opt instead for chaining with Orlicz norms. Our version using the (martingale) Yurinskii coupling can improve upon theirs in approximation rate even for independent data under certain circumstances, as follows. Suppose the setup of Proposition 1 in berthet2006revisiting; that is, $X_1, \ldots, X_n$ are i.i.d.\ and $\sup_{\ensuremath{\mathcal{F}}} \|f\|_\infty \leq M$, with the VC-type assumption $\sup_\ensuremath{\mathbb{Q}} N(\varepsilon, \ensuremath{\mathcal{F}}, d_\ensuremath{\mathbb{Q}}) \leq c_0 \varepsilon^{-\nu_0}$ where $d_\ensuremath{\mathbb{Q}}(f,f')^2 = \ensuremath{\mathbb{E}}_\ensuremath{\mathbb{Q}}\big[(f-f')^2\big]$ for a measure $\ensuremath{\mathbb{Q}}$ on $\ensuremath{\mathcal{X}}$ and $M, c_0, \nu_0$ are constants. Then, using uniform $L^2$ covering numbers rather than Orlicz norm chaining in our Proposition (ref) gives the following. Firstly as $X_i$ are i.i.d.\ we take $\Sigma(f, f') = \ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[S(f), S(f')]$ so $\Omega_\delta = 0$. Let $\ensuremath{\mathcal{F}}_\delta$ be a minimal $\delta$-cover of $(\ensuremath{\mathcal{F}}, d_\ensuremath{\mathbb{P}})$ with cardinality $N_\delta \lesssim \delta^{-\nu_0}$ where $\delta \to 0$. It is not difficult to show that $\beta_\delta \lesssim n \delta^{-\nu_0} \sqrt{\log(1/\delta)}$. Theorem 2.2.8 and Theorem 2.14.1 in van1996weak give
where we used the VC-type property to bound the entropy integral. So by our Proposition (ref), for any sequence $R_n \to \infty$ (see Remark (ref)),
where we minimized over $\delta$ in the last step. berthet2006revisiting achieved
showing that our approach achieves a better approximation rate whenever $\nu_0 > 4/3$. In particular, our method is superior in richer function classes with larger VC-type dimension. For example, if $\ensuremath{\mathcal{F}}$ is smoothly parametrized by $\theta \in \Theta \subseteq \ensuremath{\mathbb{R}}^d$ where $\Theta$ contains an open set, then $\nu_0 > 4/3$ corresponds to $d \geq 2$ and our rate is better as soon as the parameter space is more than one-dimensional. The difference in approximation rate is due to Zaitsev's coupling having better dependence on the sample size but worse dependence on the dimension. In particular, Zaitsev's coupling is stated only in $\ell_2$-norm and hence berthet2006revisiting are compelled to use the inequality $\|\cdot\|_\infty \leq \|\cdot\|_2$ in the coupling step, a bound which is loose when the dimension of the vectors (here on the order of $\delta^{-\nu_0}$) is even moderately large. We use the fact that our version of Yurinskii's coupling applies directly to the supremum norm, giving sharper dependence on the dimension.
In Section (ref) we apply Proposition (ref) to obtain strong approximations for local polynomial estimators in the nonparametric regression setting. In contrast with the series estimators of the upcoming Section (ref), local polynomial estimators are not linearly separable and hence cannot be analyzed directly using the finite-dimensional Corollary (ref).
We illustrate the applicability of our previous strong approximation results with two substantial and classical examples in nonparametric regression estimation. Firstly, we present an analysis of partitioning-based series estimators, in which we can apply the finite-dimensional result of Corollary (ref) directly due to an intrinsic linear separability property. Secondly, we consider local polynomial estimators, this time using the stochastic process formulation in Proposition (ref) due to the presence of a non-linearly separable martingale empirical process.
Partitioning-based least squares methods are essential tools for estimation and inference in nonparametric regression, encompassing splines, piecewise polynomials, compactly supported wavelets and decision trees as special cases. See cattaneo2020large for further details and references throughout this section. We illustrate the usefulness of Corollary (ref) by deriving a Gaussian strong approximation for partitioning series estimators based on multivariate martingale data. Proposition (ref) shows how we achieve the best known rate of strong approximation for independent data by imposing an additional mild $\alpha$-mixing condition to control the time series dependence of the regressors.
Consider the nonparametric regression setup with martingale difference residuals defined by $Y_i = \mu(W_i) + \varepsilon_i$ for $ 1 \leq i \leq n$ where the regressors $W_i$ have compact connected support $\ensuremath{\mathcal{W}} \subseteq \ensuremath{\mathbb{R}}^m$, $\ensuremath{\mathcal{H}}_i$ is the $\sigma$-algebra generated by $(W_1, \ldots, W_{i+1}, \varepsilon_1, \ldots, \varepsilon_i)$, $\ensuremath{\mathbb{E}}[\varepsilon_i \mid \ensuremath{\mathcal{H}}_{i-1}] = 0$ and $\mu: \ensuremath{\mathcal{W}} \to \ensuremath{\mathbb{R}}$ is the estimand. Let $p(w)$ be a $k$-dimensional vector of bounded basis functions on $\ensuremath{\mathcal{W}}$ which are locally supported on a quasi-uniform partition cattaneo2020large. Under minimal regularity conditions, the least-squares partitioning-based series estimator is $\hat\mu(w) = p(w)^{\ensuremath{\mathsf{T}}} \hat H^{-1} \sum_{i=1}^n p(W_i) Y_i$ with $\hat H = \sum_{i=1}^n p(W_i) p(W_i)^\ensuremath{\mathsf{T}}$. The approximation power of the estimator $\hat\mu(w)$ derives from letting $k\to\infty$ as $n\to\infty$. The assumptions made on $p(w)$ are mild enough to accommodate splines, wavelets, piecewise polynomials, and certain types of decision trees. For such a tree, $p(w)$ is comprised of indicator functions over $k$ axis-aligned rectangles forming a partition of $\ensuremath{\mathcal{W}}$ (a Haar basis), provided that the partitions are constructed using independent data (e.g., with sample splitting).
Our goal is to approximate the law of the stochastic process $(\hat\mu(w)-\mu(w):w\in\ensuremath{\mathcal{W}})$, which upon rescaling is typically not asymptotically tight as $k \to \infty$ and thus does not converge weakly. Nevertheless, exploiting the intrinsic linearity of the estimator $\hat\mu(w)$, we can apply Corollary (ref) directly to construct a Gaussian strong approximation. Specifically, we write
where $H= \sum_{i=1}^n \ensuremath{\mathbb{E}}\left[p(W_i) p(W_i)^\ensuremath{\mathsf{T}}\right]$ is the expected outer product matrix, $S = \sum_{i=1}^n p(W_i) \varepsilon_i$ is the score vector, and $\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w) = p(w)^{\ensuremath{\mathsf{T}}} \hat H^{-1}\sum_{i=1}^n p(W_i) \mu(W_i) - \mu(w)$. Imposing some mild time series restrictions and assuming stationarity for simplicity, it is not difficult to show (see \iftoggle{aos}{the supplementary material cattaneo2025yurinskiisupplement}{Appendix (ref)}) that $\|\hat H - H\|_1 \lesssim_\ensuremath{\mathbb{P}} \sqrt{n k}$ and $\sup_{w\in\ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)| \lesssim_\ensuremath{\mathbb{P}} k^{-\gamma}$ for some $\gamma>0$, depending on the specific structure of the basis functions, the dimension $m$ of the regressors, and the smoothness of the regression function $\mu$. Thus, it remains to study the $k$-dimensional zero-mean martingale $S$ by applying Corollary (ref) with $X_i=p(W_i) \varepsilon_i$. Controlling the convergence of the quadratic variation term $\ensuremath{\mathbb{E}}[\|\Omega\|_2]$ also requires some time series dependence assumptions; we impose an $\alpha$-mixing condition on $(W_1, \ldots, W_n)$ for illustration bradley2005basic.
The core of the proof of Proposition (ref) involves applying Corollary (ref) with $S = \sum_{i=1}^n p(W_i) \varepsilon_i$ and $p=\infty$ to construct $T \sim \ensuremath{\mathcal{N}}\big(0, \operatorname{Var}[S]\big)$ such that $\|S - T \|_\infty$ is small, and then setting $G(w) = p(w)^\ensuremath{\mathsf{T}} H^{-1} T$. So long as the bias can be appropriately controlled, this result allows for uniform inference procedures such as uniform confidence bands or shape specification testing. The condition $k^3 / n \to 0$ is the same (up to logs) as that imposed by cattaneo2020large for i.i.d. data, which gives the best known strong approximation rate for this problem. Thus, Proposition (ref) gives the same best approximation rate, without requiring any extra restrictions, for $\alpha$-mixing time series data.
Our results improve substantially on li2020uniform: using the notation of our Corollary (ref), and with any sequence $R_n \to \infty$, a valid (see Remark (ref)) version of their martingale Yurinskii coupling is
where $B_n = \sum_{i=1}^n \ensuremath{\mathbb{E}}[\|X_i\|_2^3]$ and $r_n$ is a term controlling the convergence of the quadratic variation, playing a similar role to our term $\ensuremath{\mathbb{E}}[\|\Omega\|_2]$. Under the assumptions of our Proposition (ref), applying this result with $S = \sum_{i=1}^n p(W_i) \varepsilon_i$ yields a rate no better than $\|S-T\|_2 \lesssim_\ensuremath{\mathbb{P}} (n k)^{1/3} R_n$. As such, they attain a rate of strong approximation no faster than
Hence, for this approach to yield a valid strong approximation, the number of basis functions must satisfy $k^5/n \to 0$, a more restrictive assumption than our $k^3 / n \to 0$ (up to logs). This difference is due to li2020uniform using the $\ell_2$-norm version of Yurinskii's coupling rather than the more recently established $\ell_\infty$-norm version. Further, our approach allows for an improved rate of distributional approximation whenever the residuals have zero conditional third moment.
To illustrate the statistical applicability of Proposition (ref), consider constructing a feasible uniform confidence band for the regression function $\mu$, using standardization and Studentization for statistical power improvements. We assume throughout that the bias is negligible. Proposition (ref) and anti-concentration for Gaussian suprema chernozhukov2014anti yield a distributional approximation for the supremum statistic whenever $k^3(\log n)^6 / n \to 0$, giving
where $\rho(w,w') = \ensuremath{\mathbb{E}}[G(w)G(w')]$. Furthermore, using a Gaussian--Gaussian comparison result chernozhukov2013gaussian and anti-concentration again, it is not difficult to show (see the proof of Proposition (ref)) that with $\ensuremath{\mathbf{W}} = (W_1, \ldots, W_n)$ and $\ensuremath{\mathbf{Y}} = (Y_1, \ldots, Y_n)$,
where $\hat G(w)$ is a zero-mean Gaussian process conditional on $\ensuremath{\mathbf{W}}$ and $\ensuremath{\mathbf{Y}}$ with conditional covariance function $\hat\rho(w,w') =\ensuremath{\mathbb{E}}\big[\hat G(w) \hat G(w') \mid \ensuremath{\mathbf{W}}, \ensuremath{\mathbf{Y}} \big] = p(w)^\ensuremath{\mathsf{T}} \hat H^{-1} \widehat{\operatorname{Var}}[S] \hat H^{-1}p(w')$ for some estimator $\widehat{\operatorname{Var}}[S]$ satisfying $\frac{k (\log n)^2}{n} \big\|\widehat{\operatorname{Var}}[S]-\operatorname{Var}[S]\big\|_2 \to_\ensuremath{\mathbb{P}} 0$. For example, one could use the plug-in estimator $\widehat{\operatorname{Var}}[S]=\sum_{i=1}^n p(W_i) p(W_i)^\ensuremath{\mathsf{T}} \hat{\sigma}^2(W_i)$ where $\hat{\sigma}^2(w)$ satisfies $(\log n)^2 \sup_{w \in \ensuremath{\mathcal{W}}} |\hat{\sigma}^2(w)-\sigma^2(w)| \to_\ensuremath{\mathbb{P}} 0$. This leads to the following feasible and asymptotically valid $100(1-\tau)\%$ uniform confidence band for partitioning-based series estimators based on martingale data.
As a second example application we consider nonparametric regression estimation with martingale data employing local polynomial methods Fan-Gijbels_1996_Book. In contrast with the partitioning-based series methods of Section (ref), local polynomials induce stochastic processes which are not linearly separable, allowing us to showcase the empirical process result given in Proposition (ref).
As before, suppose that $Y_i = \mu(W_i) + \varepsilon_i$ for $ 1 \leq i \leq n$ where $W_i$ has compact connected support $\ensuremath{\mathcal{W}} \subseteq \ensuremath{\mathbb{R}}^m$, $\ensuremath{\mathcal{H}}_i$ is the $\sigma$-algebra generated by $(W_1, \ldots, W_{i+1}, \varepsilon_1, \ldots, \varepsilon_i)$, $\ensuremath{\mathbb{E}}[\varepsilon_i \mid \ensuremath{\mathcal{H}}_{i-1}] = 0$, and $\mu: \ensuremath{\mathcal{W}} \to \ensuremath{\mathbb{R}}$ is the estimand. Let $K$ be a kernel function on $\ensuremath{\mathbb{R}}^m$ and $K_h(w) = h^{-m} K(w/h)$ for some bandwidth $h > 0$. Take $\gamma \geq 0$ and let $k = (m+\gamma)!/(m!\gamma!)$ be the number of monomials up to order $\gamma$. Using multi-index notation, let $p(w)$ be the $k$-dimensional vector collecting the monomials $w^{\kappa}/\kappa!$ for $0 \leq |\kappa| \leq \gamma$, and set $p_h(w) = p(w/h)$. The local polynomial regression estimator of $\mu(w)$ is, with $e_1 = (1, 0, \ldots, 0)^\ensuremath{\mathsf{T}} \in \ensuremath{\mathbb{R}}^k$,
Our goal is again to approximate the distribution of the entire stochastic process, $(\hat{\mu}(w)-\mu(w):w\in\ensuremath{\mathcal{W}})$, which upon rescaling is non-Donsker if $h \to 0$, and decomposes as follows:
where $\hat H(w) = \sum_{i=1}^n K_h(W_i-w) p_h(W_i-w) p_h(W_i-w)^\ensuremath{\mathsf{T}}$, $H(w) = \ensuremath{\mathbb{E}} \big[ \hat H(w) \big]$, $S(w)= \sum_{i=1}^n K_h(W_i-w) p_h(W_i-w) \varepsilon_i$ and $\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w) = e_1^\ensuremath{\mathsf{T}} \hat H(w)^{-1} \sum_{i=1}^n K_h(W_i-w) p_h(W_i-w) \mu(W_i) - \mu(w)$. A key distinctive feature of local polynomial regression is that both $\hat H(w)$ and $S(w)$ are functions of the evaluation point $w\in\ensuremath{\mathcal{W}}$; contrast this with the partitioning-based series estimator discussed in Section (ref), for which neither $\hat H$ nor $S$ depend on $w$. Therefore we use Proposition (ref) to obtain a Gaussian strong approximation for the martingale empirical process directly.
Under some mild regularity conditions, including stationarity for simplicity and an $\alpha$-mixing assumption on the time-dependence of the data, we first show $\sup_{w\in\ensuremath{\mathcal{W}}} \|\hat H(w)-H(w)\|_2 \lesssim_\ensuremath{\mathbb{P}} \sqrt{n h^{-2m}\log n}$. Further, $\sup_{w\in\ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)| \lesssim_\ensuremath{\mathbb{P}} h^\gamma$ provided that the regression function is sufficiently smooth. Thus it remains to analyze the martingale empirical process $\big(e_1^\ensuremath{\mathsf{T}} H(w)^{-1} S(w) : w\in\ensuremath{\mathcal{W}}\big)$ via Proposition (ref) by setting
With this approach, we obtain the following result.
If the residuals further satisfy $\ensuremath{\mathbb{E}} \left[ \varepsilon_i^3 \mid \ensuremath{\mathcal{H}}_{i-1} \right] = 0$, then a third-order Yurinskii coupling delivers an improved rate of strong approximation for Proposition (ref); this is omitted here for brevity. For completeness, the proof of Proposition (ref) verifies that if the regression function $\mu(w)$ is $\gamma$ times continuously differentiable on $\ensuremath{\mathcal{W}}$ then $\sup_w |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)| \lesssim_\ensuremath{\mathbb{P}} h^\gamma$. Further, the assumption that $p(w)$ is a vector of monomials is unnecessary in general; any collection of bounded linearly independent functions which exhibit appropriate approximation power will suffice eggermont2009maximum. As such, we can encompass local splines and wavelets, as well as polynomials, and also choose whether or not to include interactions between the regressor variables. The bandwidth restriction of $n h^{3m} \to \infty$ is analogous to that imposed in Proposition (ref) for partitioning-based series estimators, and as far as we know, has not been improved upon for non-i.i.d.\ data.
Applying an anti-concentration result for Gaussian process suprema, such as Corollary 2.1 in chernozhukov2014anti, allows one to write a Kolmogorov--Smirnov bound comparing the law of $\sup_{w \in \ensuremath{\mathcal{W}}}|\hat\mu(w) - \mu(w)|$ to that of $\sup_{w \in \ensuremath{\mathcal{W}}}|T(w)|$. With an appropriate covariance estimator, we can further replace $T(w)$ by a feasible version $\hat T(w)$ or its Studentized counterpart, enabling procedures for uniform inference analogous to the confidence bands constructed in Section (ref). We omit the details of this to conserve space but note that our assumptions on $W_i$ and $\varepsilon_i$ ensure that Studentization is possible even when the discretized covariance matrix has small eigenvalues (Section (ref)), as we normalize only by the diagonal entries.
In this setting of kernel-based local empirical processes, it is essential that our initial strong approximation result (Corollary (ref)) does not impose a lower bound on the eigenvalues of the variance matrix $\Sigma$. This effect was demonstrated by Lemma (ref) and its surrounding discussion in Section (ref), and as such, the result of li2020uniform is unsuited for this application due to its strong minimum eigenvalue assumption. Finally, for the special case of i.i.d.\ data, chernozhukov2014gaussian achieve better rates for approximating the scalar supremum of the $t$-process in Kolmogorov--Smirnov distance by bypassing the step where we first approximate the entire stochastic process (see Section (ref) for a discussion), while cattaneo2024strong obtain better strong approximations for the entire stochastic process under additional assumptions via a generalization of the celebrated Hungarian construction komlos1975approximation,rio1994local.
In this paper we introduced as our main result a new version of Yurinskii's coupling which strictly generalizes all previously known forms of the result. Our formulation gave a Gaussian mixture coupling for approximate martingale vectors in $\ell_p$-norm where $1 \leq p \leq \infty$, with no restrictions on the minimum eigenvalues of the associated covariance matrices. We further showed how to obtain an improved approximation whenever third moments of the data are negligible. We demonstrated the applicability of this main result by first deriving a user-friendly version, and then specializing it to mixingales, martingales, and independent data, illustrating the benefits with a collection of simple factor models. We then considered the problem of constructing uniform strong approximations for martingale empirical processes, demonstrating how our new Yurinskii coupling can be employed in a stochastic process setting. As substantive illustrative applications of our theory to some well established problems in statistical methodology, we showed how to use our coupling results for both vector-valued and empirical process-valued martingales in developing uniform inference procedures for partitioning-based series estimators and local polynomial models in nonparametric regression. At each stage we addressed issues of feasibility, compared our work with the existing literature, and provided implementable statistical inference procedures.
We thank the Editor, Associate Editor, and several reviewers for their comments, which led to a much improved version of this paper. We also thank Jianqing Fan, Alexander Giessing, Boris Hanin, Michael Jansson, Jason Klusowski, Arun Kumar, Boris Shigida, and Rae Yu for comments.
The authors gratefully acknowledge financial support from the National Science Foundation through grant DMS-2210561, and Cattaneo gratefully acknowledges financial support from the National Science Foundation through grant SES-2241575 and from the National Institute of Health through grant R01 GM072611-16.