The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
100,937 characters
Yurinskii's Coupling for Martingales
\title{Yurinskii's Coupling for Martingales}
\author{
Matias D.\ Cattaneo\textsuperscript{1}
\and Ricardo P.\ Masini\textsuperscript{2}
\and William G.\ Underwood\textsuperscript{3*}
}
\maketitle
\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}
\begin{abstract}
Yurinskii's coupling is a popular theoretical tool for non-asymptotic
distributional analysis in mathematical statistics and applied probability,
offering a Gaussian strong approximation with an explicit error bound under
easily verifiable conditions.
Originally stated in $\ell_2$-norm for sums of independent random vectors, it
has recently been extended both to the $\ell_p$-norm,
for $1 \leq p \leq \infty$, and to vector-valued martingales in
$\ell_2$-norm, under some strong conditions.
We present as our main result a Yurinskii coupling for approximate
martingales in $\ell_p$-norm, under substantially weaker conditions than
those previously imposed. Our formulation further allows for the coupling
variable to follow a more general Gaussian mixture distribution, and we
provide a novel third-order coupling method which gives tighter
approximations in certain settings.
We specialize our main result to mixingales, martingales, and independent
data, and derive uniform Gaussian mixture strong approximations for
martingale empirical processes. Applications to nonparametric
partitioning-based and local polynomial regression procedures are provided,
alongside central limit theorems for high-dimensional martingale vectors.
\end{abstract}
\vspace*{10mm}
\noindent\textbf{Keywords}:
coupling,
strong approximation,
mixingales,
martingales,
dependent data,
Gaussian mixture approximation,
time series,
empirical processes,
uniform inference,
series estimation,
local polynomial estimation,
central limit theorems.
\clearpage
\pagebreak
\tableofcontents
\pagebreak
\section{Introduction}
Yurinskii's coupling \citep{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, \citet[Chapter 10]{pollard2002user} for a textbook
introduction, and \citet{csorgo1981strong} and
\citet{Lindvall_1992_Book} for background references.
When compared to other coupling approaches, such as the celebrated Hungarian
construction \citep{komlos1975approximation} or Zaitsev's coupling
\citep{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
\citep{chernozhukov2014gaussian};
to establish distribution theory for non-Donsker stochastic $t$-processes
generated in nonparametric series regression \citep{belloni2015some};
to prove distributional approximations for
high-dimensional $\ell_p$-norms \citep{biau2015high};
to develop distribution theory for vector-valued martingales
\citep{belloni2018high,li2020uniform};
to derive a law of the iterated logarithm for stochastic gradient descent
optimization methods \citep{anastasiou2019normal};
to establish uniform distributional results for nonparametric
high-dimensional quantile processes \citep{belloni2019conditional};
to develop distribution theory for non-Donsker stochastic
$t$-processes generated in partitioning-based series regression
\citep{cattaneo2020large};
to deduce Bernstein--von Mises theorems in high-dimensional settings
\citep{ray2021bernstein};
and to develop distribution theory for non-Donsker U-processes based on dyadic
network data \citep{cattaneo2024uniform}.
There are also many other early applications of Yurinskii's coupling:
\citet{dudley1983invariance} and \citet{dehling1983limit} establish invariance
principles for Banach space-valued random variables, and \citet{lecam1988} and
\citet{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:
\begin{enumerate}[label=(\roman*),leftmargin=*]
\item
\label{it:contribution_approximate_martingale}
It applies to vector-valued
\textit{approximate martingale} data.
\item
\label{it:contribution_gaussian_mixture}
It allows for a \textit{Gaussian mixture} coupling distribution.
\item
\label{it:contribution_degeneracy}
It imposes \textit{no restrictions on degeneracy} of the
data covariance matrix.
\item
\label{it:contribution_third_order}
It establishes a \textit{third-order} coupling to
improve the approximation in certain situations.
\end{enumerate}
Closest to our work are
the recent paper by \citet{li2020uniform}
and the unpublished manuscript by \citet{belloni2018high},
which both investigated distribution theory for martingale
data using Yurinskii's coupling and related methods. Specifically,
\citet{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{sec:kde}), 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{it:contribution_degeneracy}, in addition to offering several other new
features not present in \citet{li2020uniform}, including
\ref{it:contribution_approximate_martingale},
\ref{it:contribution_gaussian_mixture}, \ref{it:contribution_third_order}, 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{rem:coupling_bounds_probability}).
\citet{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
\citep[see][and references therein]{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
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}}
for details). Nevertheless, a
Yurinskii coupling derived from \citet{belloni2018high} would not feature
\ref{it:contribution_approximate_martingale},
\ref{it:contribution_gaussian_mixture}, \ref{it:contribution_third_order}, 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{thm:sa_dependent})
is presented in Section~\ref{sec:main_results}, where we also specialize it to
a slightly weaker yet more user-friendly formulation
(Proposition~\ref{pro:sa_simplified}). 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{thm:sa_dependent} 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{\citealp[Theorem~B.2]{chen2020jackknife};
\citealp[Theorem~4]{monrad1991nearby}},
along with some related technical work
detailed in
\iftoggle{aos}{the supplementary material
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}}.
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
\citep{cattaneo2025yurinskiisupplement})}{
(see Lemma~\ref{lem:smooth_approximation}
in the appendix)},
enabling improved approximation errors
whenever third moments are negligible.
In Proposition~\ref{pro:sa_simplified}, 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{sec:main_results} also presents some further
refinements when additional structure is imposed, deriving Yurinskii couplings
for mixingales, martingales, and independent data as
Corollaries~\ref{cor:sa_mixingale}, \ref{cor:sa_martingale}, and
\ref{cor:sa_indep}, respectively. We take the opportunity to discuss and correct
in Remark~\ref{rem:coupling_bounds_probability} a technical issue which is
often neglected \citep{pollard2002user, li2020uniform} when using Yurinskii's
coupling to derive bounds in probability. Section~\ref{sec:factor} 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{it:contribution_approximate_martingale}--
\ref{it:contribution_third_order}.
Section~\ref{sec:emp_proc} 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{lem:kde_eigenvalue} why it is essential that we do not place any
conditions on the minimum eigenvalue of the variance matrix
\ref{it:contribution_degeneracy}.
We then present a general-purpose strong
approximation for martingale empirical processes in
Proposition~\ref{pro:emp_proc}, combining classical results in the empirical
process literature \citep{van1996weak} with our coupling from
Corollary~\ref{cor:sa_martingale}. 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 \citep{berthet2006revisiting}. Our improvement
comes from using Yurinskii's coupling for the $\ell_\infty$-norm, where
\citet{berthet2006revisiting} apply Zaitsev's coupling
\citep{zaitsev1987estimates, zaitsev1987gaussian} with the larger
$\ell_2$-norm.
Section~\ref{sec:nonparametric} 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{cor:sa_martingale}
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{pro:emp_proc}) 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{sec:conclusion} concludes the paper.
Appendix~\ref{sec:high_dim_clt} demonstrates
how our coupling results can be used to derive
distributional Gaussian approximations
(central limit theorems)
for possibly high-dimensional martingale vectors
(Proposition~\ref{pro:clt}).
This result complements a recent literature on probability
and statistics studying
the same problem but with independent data
\citep[see][and references
therein]{buzun2022strong,lopes2022central,
chernozhukov2023nearly,kock2024remark}.
We also present a version of this result employing
a covariance estimator (Proposition~\ref{pro:bootstrap}),
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{sec:lp}.
}
All proofs are collected in
\iftoggle{aos}{the supplementary material
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}},
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
}{}.
\subsection{Notation}
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
\citep[Section~2.2]{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 \citep[Theorem~1.1.10]{dudley1999uniform}.
\section{Main results}
\label{sec:main_results}
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{thm:sa_dependent}, while
Proposition~\ref{pro:sa_simplified} gives
a simplified and slightly weaker version which is easier to use in many
applications. We then further specialize
Proposition~\ref{pro:sa_simplified} to
three scenarios with successively stronger assumptions, namely mixingales,
martingales, and independent data, in Corollaries~\ref{cor:sa_mixingale},
\ref{cor:sa_martingale}, and \ref{cor:sa_indep} 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{rem:coupling_bounds_probability} 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{sec:factor} a simple factor model example relating to
all three of
the aforementioned scenarios, discussing further how our contributions are
related to the existing literature.
\begin{theorem}[Strong approximation for vector-valued approximate
martingales]
\label{thm:sa_dependent}
Take a complete probability space with a countably generated
filtration $\ensuremath{\mathcal{H}}_0, \ldots, \ensuremath{\mathcal{H}}_n$
for some $n \geq 1$, supporting the
$\ensuremath{\mathbb{R}}^d$-valued square-integrable random vectors
$X_1, \ldots, X_n$.
Let $S = \sum_{i=1}^n X_i$ and
define
\begin{align*}
\tilde X_i
&=
\sum_{r=1}^n \big(\ensuremath{\mathbb{E}}[X_{r} \mid \ensuremath{\mathcal{H}}_{i}] - \ensuremath{\mathbb{E}}[X_{r} \mid \ensuremath{\mathcal{H}}_{i-1}]\big)
& &\text{and}
&U
&=
\sum_{i=1}^{n} \big(
X_i - \ensuremath{\mathbb{E}}[ X_i \mid \ensuremath{\mathcal{H}}_n]
+ \ensuremath{\mathbb{E}}[ X_i \mid \ensuremath{\mathcal{H}}_0 ]
\big).
\end{align*}
Let $V_i = \operatorname{Var}[\tilde X_i \mid \ensuremath{\mathcal{H}}_{i-1}]$ and
define $\Omega = \sum_{i=1}^n V_i - \Sigma$
where $\Sigma$ is an almost surely positive semi-definite $\ensuremath{\mathcal{H}}_0$-measurable
$d \times d$ random matrix.
Then, for each $\eta > 0$ and $p \in [1,\infty]$,
there exists, on an enlarged probability space,
an $\ensuremath{\mathbb{R}}^d$-valued random vector $T$ with
$T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ such that
\begin{align}
\label{eq:sa_dependent}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > 6\eta\big)
&\leq
\inf_{t>0}
\left\{
2 \ensuremath{\mathbb{P}}\big( \|Z\|_p > t \big)
+ \min\left\{
\frac{\beta_{p,2} t^2}{\eta^3},
\frac{\beta_{p,3} t^3}{\eta^4}
+ \frac{\pi_3 t^3}{\eta^3}
\right\}
\right\} \nonumber \\
&\quad+
\inf_{M \succeq 0}
\Big\{ 2 \ensuremath{\mathbb{P}}\big(\Omega \npreceq M\big) + \delta_p(M,\eta)
+ \varepsilon_p(M, \eta)\Big\}
+\ensuremath{\mathbb{P}}\big(\|U\|_p>\eta\big),
\end{align}
where $Z, Z_1,\dots ,Z_n$ are i.i.d.\ standard Gaussian random variables on
$\ensuremath{\mathbb{R}}^d$ independent of $\ensuremath{\mathcal{H}}_n$,
the second infimum is taken over all positive semi-definite
$d \times d$ non-random matrices $M$,
\begin{align*}
\beta_{p,k}
&=
\sum_{i=1}^n \ensuremath{\mathbb{E}}\left[\| \tilde X_i \|^k_2 \| \tilde X_i \|_p
+ \|V_i^{1/2} Z_i \|^k_2 \|V_i^{1/2} Z_i \|_p \right],
&\pi_3
&=
\sum_{i=1}^{n}
\sum_{|\kappa| = 3}
\ensuremath{\mathbb{E}} \Big[ \big|
\ensuremath{\mathbb{E}} [ \tilde X_i^\kappa \mid \ensuremath{\mathcal{H}}_{i-1} ]
\big| \Big]
\end{align*}
for $k \in \{2, 3\}$, with $\pi_3 = \infty$ if the associated
conditional expectation does not exist, and with
\begin{align*}
\delta_p(M,\eta)
&=
\ensuremath{\mathbb{P}}\left(
\big\|\big((\Sigma +M)^{1/2}- \Sigma^{1/2}\big) Z\big\|_p
\geq \eta
\right), \\
\varepsilon_p(M, \eta)
&=
\ensuremath{\mathbb{P}}\left(\big\| (M - \Omega)^{1/2} Z \big\|_p\geq \eta, \
\Omega \preceq M\right).
\end{align*}
\end{theorem}
This theorem offers four novel contributions to the literature on
coupling theory and strong approximation,
as discussed in the introduction.
Firstly \ref{it:contribution_approximate_martingale}, 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,
\citet{wu2004martingale}, \citet{wu2005nonlinear},
\citet{dedecker2007weak}, \citet{zhao2008martingale},
\citet{peligrad2010conditional}, \citet{atchade2014martingale},
\citet{cuny2014martingale}, \citet{magda2018martingale}, and references
therein. In Section~\ref{sec:mixingales} 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{sec:factor} provides an illustration in the context of
auto-regressive factor models.
Secondly \ref{it:contribution_gaussian_mixture},
Theorem~\ref{thm:sa_dependent} 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{thm:sa_dependent} for a few reasons.
To begin with, we explicitly indicate in the
conditions of Theorem~\ref{thm:sa_dependent} 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{\citealp[Theorem~B.2]{chen2020jackknife};
\citealp[Theorem~4]{monrad1991nearby}}.
Section~\ref{sec:martingales}
illustrates how a strong approximation akin to
mixing convergence can arise when the data
forms an exact martingale, and Section~\ref{sec:factor} gives a simple example
relating to factor modeling in statistics and data science.
As a third contribution to the literature
\ref{it:contribution_degeneracy}, and
of particular importance for applications,
Theorem~\ref{thm:sa_dependent} 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{sec:kde} by
demonstrating the significant improvements in strong approximation errors
delivered by Theorem~\ref{thm:sa_dependent} relative to those obtained using
prior results in the literature.
Finally \ref{it:contribution_third_order},
Theorem \ref{thm:sa_dependent} 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{sec:martingales}).
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 \eqref{eq:sa_dependent}.
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{sec:nonparametric}).
\subsection{User-friendly formulation of the main result}
The result in Theorem~\ref{thm:sa_dependent} 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{pro:sa_simplified}, 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
\begin{align*}
\phi_p(d) =
\begin{cases}
\sqrt{pd^{2/p} } & \text{ if } p \in [1,\infty),\\
\sqrt{2\log 2d} & \text{ if } p =\infty.
\end{cases}
\end{align*}
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
\citep{cattaneo2025yurinskiisupplement}}{
Lemma~\ref{lem:gaussian_pnorm} in the appendix}).
\begin{proposition}[Simplified strong approximation for
vector-valued approximate martingales]
\label{pro:sa_simplified}
Assume the setup and notation of Theorem~\ref{thm:sa_dependent}.
For each $\eta > 0$ and $p \in [1,\infty]$,
there exists a random vector $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ satisfying
\begin{align*}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,2} \phi_p(d)^2}{\eta^3}
\right)^{1/3}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}
+\ensuremath{\mathbb{P}}\left(\|U\|_p>\frac{\eta}{6}\right).
\end{align*}
If further $\pi_3 = 0$, then also
\begin{align*}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,3} \phi_p(d)^3}{\eta^4}
\right)^{1/4}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}
+\ensuremath{\mathbb{P}}\left(\|U\|_p>\frac{\eta}{6}\right).
\end{align*}
\end{proposition}
Proposition~\ref{pro:sa_simplified} 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{sec:series}), 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{sec:factor}), 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{pro:sa_simplified}, 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.
\subsection{Mixingales}\label{sec:mixingales}
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$,
\begin{align}
\label{eq:mixingale_1}
\ensuremath{\mathbb{E}} \left[ \left\|
\ensuremath{\mathbb{E}} \left[ X_i \mid \ensuremath{\mathcal{H}}_{i-r} \right]
\right\|_p \right]
&\leq
c_i \zeta_r,
\end{align}
and for all $1 \leq i \leq n$ and $0 \leq r \leq n-i$,
\begin{align}
\label{eq:mixingale_2}
\ensuremath{\mathbb{E}} \left[ \big\|
X_i - \ensuremath{\mathbb{E}} \big[ X_i \mid \ensuremath{\mathcal{H}}_{i+r} \big]
\big\|_p \right]
&\leq
c_i \zeta_{r+1}.
\end{align}
These conditions are satisfied, for example, if $X_i$ are integrable strongly
$\alpha$-mixing random variables \citep{mcleish1975invariance}, or
if $X_i$ are
generated by an auto-regressive or auto-regressive moving average process (see
Section~\ref{sec:factor}), among many other possibilities
\citep{bradley2005basic}. Then, in the notation of
Theorem~\ref{thm:sa_dependent}, we have by Markov's inequality that
\begin{align*}
\ensuremath{\mathbb{P}} \left( \|U\|_p > \frac{\eta}{6} \right)
&\leq
\frac{6}{\eta}
\sum_{i=1}^{n}
\ensuremath{\mathbb{E}} \left[
\big\|
X_i - \ensuremath{\mathbb{E}} \left[ X_i \mid \ensuremath{\mathcal{H}}_n \right]
\big\|_p
+ \big\|
\ensuremath{\mathbb{E}} \left[ X_i \mid \ensuremath{\mathcal{H}}_0 \right]
\big\|_p
\right]
\leq \frac{\zeta}{\eta},
\end{align*}
with $\zeta = 6 \sum_{i=1}^{n} c_i (\zeta_{i} + \zeta_{n-i+1})$.
Combining Proposition~\ref{pro:sa_simplified} with this
martingale error bound yields the following result for mixingales.
\begin{corollary}[Strong approximation for vector-valued mixingales]
\label{cor:sa_mixingale}
Assume the setup and notation of Theorem~\ref{thm:sa_dependent}, and suppose
that the mixingale conditions \eqref{eq:mixingale_1} and
\eqref{eq:mixingale_2} hold. For each $\eta > 0$ and $p \in
[1,\infty]$ there
exists a random vector $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ satisfying
\begin{align*}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,2} \phi_p(d)^2}{\eta^3}
\right)^{1/3}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}
+ \frac{\zeta}{\eta}.
\end{align*}
If further $\pi_3 = 0$ then
\begin{align*}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,3} \phi_p(d)^3}{\eta^4}
\right)^{1/4}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}
+ \frac{\zeta}{\eta}.
\end{align*}
\end{corollary}
The closest antecedent to Corollary~\ref{cor:sa_mixingale} is found in
\citet[Theorem~4]{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{sec:factor} and Section~\ref{sec:nonparametric},
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{rem:coupling_bounds_probability}), we
clarify a technical issue in \citet{li2020uniform} surrounding the derivation
of valid probability bounds for $\|S-T\|_p$.
Corollary~\ref{cor:sa_mixingale} 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{pro:sa_simplified} to derive a similar result.
\subsection{Martingales}\label{sec:martingales}
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{pro:sa_simplified} in this setting
directly yields the following result.
\begin{corollary}[Strong approximation for vector-valued martingales]
\label{cor:sa_martingale}
With the setup and notation of Theorem~\ref{thm:sa_dependent}, suppose $X_i$
is $\ensuremath{\mathcal{H}}_i$-measurable with $\ensuremath{\mathbb{E}}[X_i \mid \ensuremath{\mathcal{H}}_{i-1}] = 0$ for
$1 \leq i \leq n$. Then, for each $\eta > 0$ and $p \in
[1,\infty]$, there is
a random vector $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ with
\begin{align}
\label{eq:sa_martingale_order_2}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,2} \phi_p(d)^2}{\eta^3}
\right)^{1/3}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}.
\end{align}
If further $\pi_3 = 0$ then
\begin{align}
\label{eq:sa_martingale_order_3}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,3} \phi_p(d)^3}{\eta^4}
\right)^{1/4}
+ 17 \left(
\frac{\ensuremath{\mathbb{E}} \left[ \|\Omega\|_2 \right] \phi_p(d)^2}{\eta^2}
\right)^{1/3}.
\end{align}
\end{corollary}
The closest antecedents to Corollary~\ref{cor:sa_martingale} are
\citet{belloni2018high} and \citet{li2020uniform}, who also (implicitly or
explicitly) considered Yurinskii's coupling for martingales. More
specifically,
\citet[Theorem 1]{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 \citet[Theorem 2.1]{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{cor:sa_martingale} improves over this prior work as follows.
With respect to \citet{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
\eqref{eq:sa_martingale_order_2} in
Corollary~\ref{cor:sa_martingale} gives the
same rate of strong approximation as that in Theorem~1 of
\citet{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{sec:kde} we
demonstrate the crucial importance of removing this eigenvalue lower bound
restriction in applications involving nonparametric kernel
estimators, while in
Section~\ref{sec:series} we demonstrate how the availability of a third-order
coupling \eqref{eq:sa_martingale_order_3} can give improved
approximation rates
in applications involving nonparametric series estimators with conditionally
symmetrically distributed residual errors. Finally, our technical
work improves
on \citet{li2020uniform} in two respects: (i) we employ a conditional version
of Strassen's theorem (see
\iftoggle{aos}{
Lemma~SA.1 in
the supplementary material
\citep{cattaneo2025yurinskiisupplement}}{Lemma~\ref{lem:strassen}
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{rem:coupling_bounds_probability} makes clear.
\begin{remark}[Yurinskii's coupling and bounds in probability]
\label{rem:coupling_bounds_probability}
Given a sequence of random vectors $S_n$, Yurinskii's method provides a
coupling in the following form: for each $n$ and any $\eta > 0$,
there exists
a random vector $T_n$ with $\ensuremath{\mathbb{P}}\big(\|S_n - T_n\| > \eta\big) < r_n(\eta)$,
where $r_n(\eta)$ is the approximation error. Crucially, each coupling
variable $T_n$ is a function of the desired approximation level $\eta$ and,
as such, deducing bounds in probability on $\|S_n - T_n\|$ requires some
extra care. One option is to select a sequence $R_n \to \infty$
and note that
$\ensuremath{\mathbb{P}}\big(\|S_n - T_n\| > r_n^{-1}(1 / R_n)\big) < 1 / R_n \to 0$ and hence
$\|S_n - T_n\| \lesssim_\ensuremath{\mathbb{P}} r_n^{-1}(1 / R_n)$. In this case,
$T_n$ depends on
the choice of $R_n$, which can in turn typically be chosen to diverge slowly
enough to cause no issues in applications.
\end{remark}
Technicalities akin to those outlined in
Remark~\ref{rem:coupling_bounds_probability} have been both addressed and
neglected alike in the prior literature. \citet[Chapter 10.4, Example
16]{pollard2002user} apparently misses this subtlety, providing an
inaccurate bound in probability based on the Yurinskii coupling.
\citet{li2020uniform} seem to make the same mistake in the proof of their
Lemma~A2, which invalidates the conclusion of their Theorem~1. In contrast,
\citet{belloni2015some} and \citet{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 \eqref{eq:sa_martingale_order_2} in
Corollary~\ref{cor:sa_martingale} implies that for any $R_n \to \infty$,
\begin{align*}
\|S - T\|_p
\lesssim_\ensuremath{\mathbb{P}}
\beta_{p,2}^{1/3}
\phi_p(d)^{2/3} R_n
+ \ensuremath{\mathbb{E}}[\|\Omega\|_2]^{1/2}
\phi_p(d) R_n.
\end{align*}
This bound is comparable to that obtained by \citet[Theorem~1]{li2020uniform}
with $p=2$, albeit with their formulation missing the $R_n$ correction terms.
In Section~\ref{sec:series} 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 \citet{belloni2018high}, our
Corollary~\ref{cor:sa_martingale} 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{cor:sa_martingale}, this would require additional
technical work
on their behalf in multiple ways
(see
\iftoggle{aos}{the supplementary material
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}}):
(i) a correct application of a conditional
version of Strassen's theorem
\iftoggle{aos}{
(Lemma~SA.1 in
the supplementary material
\citep{cattaneo2025yurinskiisupplement})}{(Lemma~\ref{lem:strassen}
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{lem:smooth_approximation},
\ref{lem:gaussian_useful}, and \ref{lem:gaussian_pnorm})};
(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{lem:sa_martingale})};
among others.
\subsection{Independence}
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{cor:sa_martingale}
immediately yields the following result.
\begin{corollary}[Strong approximation for sums of independent vectors]
\label{cor:sa_indep}
Assume the setup of Theorem~\ref{thm:sa_dependent},
and suppose $X_i$ are independent given $\ensuremath{\mathcal{H}}_0$,
with $\ensuremath{\mathbb{E}}[X_i \mid \ensuremath{\mathcal{H}}_0] = 0$.
Then, for each $\eta > 0$ and $p \in [1,\infty]$,
with $\Sigma = \sum_{i=1}^n V_i$,
there exists $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ satisfying
\begin{align}
\label{eq:sa_indep_order_2}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,2} \phi_p(d)^2}{\eta^3}
\right)^{1/3}.
\end{align}
If further $\pi_3 = 0$ then
\begin{align*}
\ensuremath{\mathbb{P}}\big(\|S-T\|_p > \eta\big)
&\leq
24 \left(
\frac{\beta_{p,3} \phi_p(d)^3}{\eta^4}
\right)^{1/4}.
\end{align*}
\end{corollary}
Taking $\ensuremath{\mathcal{H}}_0$ to be trivial, the first inequality \eqref{eq:sa_indep_order_2}
in Corollary~\ref{cor:sa_indep} provides an $\ell_p$-norm approximation
analogous to that presented in \cite{belloni2019conditional}. By further
restricting to $p=2$, we recover the original Yurinskii coupling as presented
in \citet[Theorem~1]{lecam1988} and \citet[Theorem~10]{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$.
\subsection{Stylized example: factor modeling}
\label{sec:factor}
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
\begin{align*}
X_i = L f_i + \varepsilon_i, \qquad 1 \leq i \leq n,
\end{align*}
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 \citet{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{thm:sa_dependent}, the
user-friendly Proposition~\ref{pro:sa_simplified}, and their specialized
Corollaries~\ref{cor:sa_mixingale}--\ref{cor:sa_indep}. 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.
\begin{enumerate}[label=(\roman*)]
\item \emph{Independent data}.
Suppose that the factors $(f_1, \ldots,
f_n)$ are independent conditional on $L$ and satisfy
$\ensuremath{\mathbb{E}} [ f_i \mid L ] = 0$.
Then, since $X_i$ are independent conditional on $\ensuremath{\mathcal{H}}_0$ and with
$\ensuremath{\mathbb{E}} [ X_i \mid \ensuremath{\mathcal{H}}_0 ] = \ensuremath{\mathbb{E}} [ L f_i + \varepsilon_i \mid L ] = 0$,
we can apply Corollary~\ref{cor:sa_indep} to $\sum_{i=1}^n X_i$.
In general, we will obtain a coupling variable which has the Gaussian
mixture distribution $T \mid \ensuremath{\mathcal{H}}_0 \sim \ensuremath{\mathcal{N}}(0, \Sigma)$ where
$\Sigma= \sum_{i=1}^n (L\operatorname{Var}[f_i \mid L]L^\ensuremath{\mathsf{T}} +\operatorname{Var}[\varepsilon_i])$.
In the special case where $L$ is non-random
and $\ensuremath{\mathcal{H}}_0$ is trivial, the coupling is Gaussian. Furthermore,
if $f_i\mid L$ and $\varepsilon_i$ are symmetric about zero
and bounded almost surely, then $\pi_3=0$, and the coupling is improved.
\item \emph{Martingales}.
Suppose instead that we assume only a martingale
condition on the latent factor variables so that
$\ensuremath{\mathbb{E}} \left[ f_i \mid L, f_1, \ldots, f_{i-1} \right] = 0$.
Then $\ensuremath{\mathbb{E}} [ X_i \mid \ensuremath{\mathcal{H}}_{i-1} ]
= L\, \ensuremath{\mathbb{E}} \left[ f_i \mid \ensuremath{\mathcal{H}}_{i-1} \right] = 0$
and Corollary~\ref{cor:sa_martingale} is applicable to $\sum_{i=1}^n X_i$.
The preceding comments on Gaussian mixture distributions
and third-order couplings continue to apply.
\item \emph{Mixingales}.
Finally, assume that the factors follow the
auto-regressive model $f_i = A f_{i-1} + u_i$ where
$A \in \ensuremath{\mathbb{R}}^{m \times m}$ is non-random and $(u_1, \ldots, u_n)$ are
zero-mean, independent, and independent of
$(\varepsilon_1, \ldots, \varepsilon_n)$.
Then $\ensuremath{\mathbb{E}} \left[ f_i \mid f_0 \right] = A^i f_0$, so taking
$p \in [1, \infty]$ we see that
$\ensuremath{\mathbb{E}} \big[ \| \ensuremath{\mathbb{E}} [ f_i \mid f_0 ] \|_p \big]
= \ensuremath{\mathbb{E}} \big[ \| A^i f_0 \|_p \big] \leq \|A\|_p^i\,\ensuremath{\mathbb{E}} [ \|f_0\|_p ]$,
and that clearly $f_i - \ensuremath{\mathbb{E}} [ f_i \mid \ensuremath{\mathcal{H}}_n ] = 0$.
Thus, whenever $\|A\|_p < 1$, the geometric sum formula implies that
the mixingale result from Corollary~\ref{cor:sa_mixingale} applies to
$\sum_{i=1}^n X_i$. The conclusions on Gaussian mixture distributions
and third-order couplings parallel the previous cases.
\end{enumerate}
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.
\section{Strong approximation for martingale empirical processes}
\label{sec:emp_proc}
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 \citet{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 \textit{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{cor:sa_martingale}).
Our results can be further generalized to \emph{approximate} martingale
empirical processes (including mixingale empirical processes;
cf.\ Corollary \ref{cor:sa_mixingale}), but
to reduce notation and the technical burden
we do not consider this extension.
\subsection{Motivating example: kernel density estimation}
\label{sec:kde}
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{lem:kde_eigenvalue}
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
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}}
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
\begin{align*}
\hat g(x)
&=
\frac{1}{n}
\sum_{i=1}^{n}
K_h( X_i - x),
\qquad K_h(u) = \frac{1}{h} K\left( \frac{u}{h} \right).
\end{align*}
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:
\begin{align*}
\sup_{x \in \ensuremath{\mathcal{X}}}
\big|S(x) - T(x)\big|
&\leq
\iftoggle{aos}{\!\!}{}
\sup_{|x-x'| \leq \delta}
\iftoggle{aos}{\!\!}{}
\big|S(x) - S(x') \big|
+ \max_{x \in \ensuremath{\mathcal{X}}_\delta}
\big|S(x) - T(x)\big|
+ \iftoggle{aos}{\!\!}{}
\sup_{|x-x'| \leq \delta}
\iftoggle{aos}{\!\!}{}
\big|T(x) - T(x')\big|.
\end{align*}
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{cor:sa_indep} 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
\begin{align*}
\ensuremath{\mathbb{P}}\left(
\|S(\ensuremath{\mathcal{X}}_\delta) - T(\ensuremath{\mathcal{X}}_\delta)\|_\infty > \eta
\right)
&\leq
31 \left(
\frac{N \log 2 N}{\eta^3 n^2 h^2}
\right)^{1/3}
\end{align*}
assuming that $1/h \geq \log 2 N$.
By the Vorob'ev--Berkes--Philipp theorem
\citep[Theorem~1.1.10]{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 \citep[Section~2.2]{van1996weak} that if
$\log(N/h) \lesssim \log n$ and $n h \gtrsim \log n$,
\begin{align*}
\sup_{|x-x'| \leq \delta}
\big\|S(x) - S(x') \big\|_\infty
&\lesssim_\ensuremath{\mathbb{P}}
\delta
\sqrt{\frac{\log n}{n h^3}}, \quad\text{and}\quad
\sup_{|x-x'| \leq \delta}
\big\|T(x) - T(x')\big\|_\infty
\lesssim_\ensuremath{\mathbb{P}}
\delta
\sqrt{\frac{\log n}{n h^3}}.
\end{align*}
Finally, for any sequence $R_n\to\infty$
(Remark~\ref{rem:coupling_bounds_probability}),
the resulting bound on the coupling error is
\begin{align*}
\sup_{x \in \ensuremath{\mathcal{X}}}
\big| S(x) - T(x) \big| \lesssim_\ensuremath{\mathbb{P}}
\left( \frac{N \log 2N}{n^2 h^2} \right)^{1/3} R_n
+ \delta \sqrt{\frac{\log n}{n h^3}},
\end{align*}
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
\begin{align*}
\sup_{x \in \ensuremath{\mathcal{X}}}
\big| S(x) - T(x) \big| \lesssim_\ensuremath{\mathbb{P}}
\left( \frac{(\log n)^3}{n^5 h^7} \right)^{1/8} R_n
\end{align*}
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 \citet{gine2004kernel}
and \citet{cattaneo2024strong} when employing a Hungarian
construction \citep{komlos1975approximation}, or $(\log n)^6 / (n h) \to 0$
imposed by \citet{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{sec:local_poly} for related discussion.
The discretization strategy outlined above is at the core of the
proof strategy
for our upcoming Proposition~\ref{pro:emp_proc}. Since we will consider
martingale empirical processes, our proof will rely on
Corollary~\ref{cor:sa_martingale}, which, unlike the martingale Yurinskii
coupling established by \citet{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{lem:kde_eigenvalue} 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 \citet{li2020uniform} coupling were to be
used instead of the results in this paper.
\begin{lemma}[Minimum eigenvalue of a
kernel density estimator covariance matrix]
\label{lem:kde_eigenvalue}
The minimum eigenvalue of
$\Sigma=n h\operatorname{Var}[S(\ensuremath{\mathcal{X}}_\delta)] \in \ensuremath{\mathbb{R}}^{N \times N}$
satisfies the upper bound
\begin{align*}
\lambda_{\min}(\Sigma)
&\leq
2 e^{-h^2/\delta^2}
+ \frac{h}{\pi a \delta}
e^{-a^2 / h^2}.
\end{align*}
\end{lemma}
Figure~\ref{fig:min_eig} shows how the upper bound in Lemma
\ref{lem:kde_eigenvalue} 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 \citet{li2020uniform}
ineffective in such scenarios.
\begin{figure}[ht]
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{sim_2.pdf}
\caption{Bandwidth $h = 0.03$}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width=\textwidth]{sim_1.pdf}
\caption{Bandwidth $h = 0.01$}
\end{subfigure}
\caption{
Upper bounds on the minimum eigenvalue of the discretized covariance
matrix in kernel density estimation,
with $n=100$ and $a = 0.2$.
Simulated: the kernel density estimator is simulated,
resampling the data $100$ times
to estimate its covariance.
Computing matrix: the minimum eigenvalue of the limiting covariance
matrix $\Sigma$ is computed explicitly.
Upper bound: the bound derived in
Lemma~\ref{lem:kde_eigenvalue}
is shown.
}
\label{fig:min_eig}
\end{figure}
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{sec:nonparametric} for further
discussion.
\subsection{General result for martingale empirical processes}
We now give our general result on a strong approximation for
martingale empirical processes, obtained by applying
the first result \eqref{eq:sa_martingale_order_2} in
Corollary~\ref{cor:sa_martingale} with $p=\infty$
to a discretization of the empirical process,
as in Section~\ref{sec:kde}.
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 \citep{geer2000empirical}
and sequential symmetrization \citep{rakhlin2015sequential}.
A class of functions is said to be \emph{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
\begin{align*}
\Psi
&=
\Big\{
\psi: [0, \infty) \to [0, \infty)
\text{ convex increasing, }
\psi(0) = 0,\
\limsup_{x,y \to \infty} \tfrac{\psi(x) \psi(y)}{\psi(C x y)} < \infty
\text{ for } C > 0
\Big\}
\end{align*}
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 \citet[Section~2.2]{van1996weak}.
\begin{proposition}[Strong approximation for martingale empirical processes]
\label{pro:emp_proc}
Let $X_i$ be random variables for $1 \leq i \leq n$ taking values in a
measurable space $\ensuremath{\mathcal{X}}$, and $\ensuremath{\mathcal{F}}$ be a pointwise measurable class of
functions from $\ensuremath{\mathcal{X}}$ to $\ensuremath{\mathbb{R}}$. Let $\ensuremath{\mathcal{H}}_0, \ldots, \ensuremath{\mathcal{H}}_n$ be a
filtration such
that each $X_i$ is $\ensuremath{\mathcal{H}}_i$-measurable, with $\ensuremath{\mathcal{H}}_0$ the trivial
$\sigma$-algebra, and suppose that $\ensuremath{\mathbb{E}}[f(X_i) \mid \ensuremath{\mathcal{H}}_{i-1}] = 0$ for all
$f \in \ensuremath{\mathcal{F}}$. Define $S(f) = \sum_{i=1}^n f(X_i)$ for $f\in\ensuremath{\mathcal{F}}$ and let
$\Sigma: \ensuremath{\mathcal{F}} \times \ensuremath{\mathcal{F}} \to \ensuremath{\mathbb{R}}$ be an almost surely positive semi-definite
$\ensuremath{\mathcal{H}}_0$-measurable random function. Suppose that for a non-random
metric $d$ on $\ensuremath{\mathcal{F}}$, constant $L$ and $\psi \in \Psi$,
\begin{align}
\label{eq:emp_proc_var}
\Sigma(f,f) - 2\Sigma(f,f') + \Sigma(f',f')
+ {\big\vert\kern-0.35ex\big\vert\kern-0.35ex\big\vert} S(f) - S(f') {\big\vert\kern-0.35ex\big\vert\kern-0.35ex\big\vert}_\psi^2
&\leq L^2 d(f,f')^2 \quad \text{a.s.}
\end{align}
Then for each $\eta > 0$ there is a
process $T(f)$ indexed by $f\in\ensuremath{\mathcal{F}}$
which, conditional on $\ensuremath{\mathcal{H}}_0$, is zero-mean and Gaussian,
satisfying
$\ensuremath{\mathbb{E}}\big[ T(f) T(f') \mid \ensuremath{\mathcal{H}}_0 \big] = \Sigma(f,f')$
for all $f, f' \in \ensuremath{\mathcal{F}}$, and for all $t > 0$ has
\begin{align*}
&\ensuremath{\mathbb{P}}\left(
\sup_{f \in \ensuremath{\mathcal{F}}}
\big| S(f) - T(f) \big|
\geq C_\psi(t + \eta)
\right)
\leq
C_\psi
\inf_{\delta > 0}
\inf_{\ensuremath{\mathcal{F}}_\delta}
\Bigg\{
\frac{\beta_\delta^{1/3} (\log 2 |\ensuremath{\mathcal{F}}_\delta|)^{1/3}}{\eta } \\
&\qquad\quad+
\left(\frac{\sqrt{\log 2 |\ensuremath{\mathcal{F}}_\delta|}
\sqrt{\ensuremath{\mathbb{E}}\left[\|\Omega_\delta\|_2\right]}}{\eta }\right)^{2/3}
+ \psi\left(\frac{t}{L J_\psi(\delta)}\right)^{-1}
+ \exp\left(\frac{-t^2}{L^2 J_2(\delta)^2}\right)
\Bigg\},
\end{align*}
where $\ensuremath{\mathcal{F}}_\delta$ is any finite $\delta$-cover of $(\ensuremath{\mathcal{F}},d)$
and $C_\psi$ is a constant depending only on $\psi$, with
\begin{align*}
\beta_\delta
&= \sum_{i=1}^n
\ensuremath{\mathbb{E}}\left[ \|\ensuremath{\mathcal{F}}_\delta(X_i)\|^2_2\|\ensuremath{\mathcal{F}}_\delta(X_i)\|_\infty
+ \|V_i(\ensuremath{\mathcal{F}}_\delta)^{1/2}Z_i\|^2_2
\|V_i(\ensuremath{\mathcal{F}}_\delta)^{1/2}Z_i\|_\infty \right],
\end{align*}
\vspace*{-6mm}
\begin{align*}
V_i(\ensuremath{\mathcal{F}}_\delta)
&=
\ensuremath{\mathbb{E}}\big[\ensuremath{\mathcal{F}}_\delta(X_i) \ensuremath{\mathcal{F}}_\delta(X_i)^\ensuremath{\mathsf{T}} \mid \ensuremath{\mathcal{H}}_{i-1} \big],
&
\Omega_\delta
&=
\sum_{i=1}^n V_i(\ensuremath{\mathcal{F}}_\delta) - \Sigma(\ensuremath{\mathcal{F}}_\delta), \\
J_\psi(\delta)
&=
\int_0^\delta \psi^{-1}\big( N_\varepsilon \big)
\,\mathrm{d}\varepsilon
+ \delta \psi^{-1} \big( N_\delta^2 \big),
&
J_2(\delta)
&= \int_0^\delta \sqrt{\log N_\varepsilon}
\,\mathrm{d}\varepsilon,
\end{align*}
where $N_\delta = N(\delta, \ensuremath{\mathcal{F}}, d)$
is the $\delta$-covering number of $(\ensuremath{\mathcal{F}}, d)$
and $Z_i$ are i.i.d.\ $\ensuremath{\mathcal{N}}\big(0, I_{|\ensuremath{\mathcal{F}}_\delta|}\big)$
independent of $\ensuremath{\mathcal{H}}_n$.
If $\ensuremath{\mathcal{F}}_\delta$ is a minimal $\delta$-cover
of $(\ensuremath{\mathcal{F}}, d)$, then $|\ensuremath{\mathcal{F}}_\delta| = N_\delta$.
\end{proposition}
Proposition~\ref{pro:emp_proc}
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.
\begin{description}
\item[Polynomial:]
$\psi(x) = x^a$ for $a \geq 2$
has ${\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_2 \leq {\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_\psi$ and
$\sqrt{\log x} \leq \sqrt{a} \psi^{-1}(x)$.
\item[Exponential:]
$\psi(x) = \exp(x^a) - 1$ for $a \in [1,2]$
has ${\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_2 \leq 2{\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_\psi$ and
$\sqrt{\log x} \leq \psi^{-1}(x)$.
\item[Bernstein:]
$\psi(x) = \exp\hspace*{-0.8mm}
\Big(\hspace*{-0.6mm}
\Big(\frac{\sqrt{1+2ax}-1}{a}\Big)^{\hspace*{-0.6mm}2}
\Big)-1$
for $a > 0$ has
${\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_2 \leq (1+a){\vert\kern-0.25ex\vert\kern-0.25ex\vert} X {\vert\kern-0.25ex\vert\kern-0.25ex\vert}_\psi$ and
$\sqrt{\log x}~\leq~\psi^{-1}(x)$.
\end{description}
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 \eqref{eq:emp_proc_var} 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 \citet{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{cor:sa_martingale} 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{pro:emp_proc} 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 \citet{berthet2006revisiting}, who used a
Gaussian coupling due to
\citet{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 \citet{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{pro:emp_proc}
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 \citet{van1996weak} give
\begin{align*}
\ensuremath{\mathbb{E}}\left[
\sup_{d_\ensuremath{\mathbb{P}}(f,f') \leq \delta}
\Big(
|S(f) - S(f')|
+ |T(f) - T(f')|
\Big)
\right]
&\lesssim
\sup_\ensuremath{\mathbb{Q}}
\int_0^\delta
\sqrt{n \log N(\varepsilon, \ensuremath{\mathcal{F}}, d_\ensuremath{\mathbb{Q}})}
\,\mathrm{d}\varepsilon
\iftoggle{aos}{\\ &\lesssim}{\lesssim}
\delta \sqrt{n\log(1/\delta)},
\end{align*}
where we used the VC-type property to bound the entropy integral.
So by our Proposition~\ref{pro:emp_proc},
for any sequence $R_n \to \infty$
(see Remark~\ref{rem:coupling_bounds_probability}),
\begin{align*}
\sup_{f \in \ensuremath{\mathcal{F}}}
\big| S(f) - T(f) \big|
&\lesssim_\ensuremath{\mathbb{P}}
n^{1/3} \delta^{-\nu_0/3}
\sqrt{\log(1/\delta)} R_n
+ \delta \sqrt{n\log(1/\delta)}
\lesssim_\ensuremath{\mathbb{P}}
n^{\frac{2+\nu_0}{6+2\nu_0}}
\sqrt{\log n} R_n,
\end{align*}
where we minimized over $\delta$ in the last step.
\citet[Proposition~1]{berthet2006revisiting} achieved
\begin{align*}
\sup_{f \in \ensuremath{\mathcal{F}}}
\big| S(f) - T(f) \big|
&\lesssim_\ensuremath{\mathbb{P}}
n^{\frac{5\nu_0}{4+10\nu_0}}
(\log n)^{\frac{4+5\nu_0}{4+10\nu_0}},
\end{align*}
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
\citet[Equation~5.3]{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{sec:local_poly} we apply Proposition~\ref{pro:emp_proc} to
obtain strong approximations for local polynomial estimators in the
nonparametric regression setting. In contrast with the series
estimators of the
upcoming Section~\ref{sec:series}, local polynomial estimators are
not linearly
separable and hence cannot be analyzed directly using the finite-dimensional
Corollary~\ref{cor:sa_martingale}.
\section{Applications to nonparametric regression}
\label{sec:nonparametric}
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{cor:sa_martingale} directly
due to an intrinsic linear separability property. Secondly, we consider local
polynomial estimators, this time using the stochastic process
formulation in Proposition~\ref{pro:emp_proc} due to
the presence of a non-linearly separable martingale empirical process.
\subsection{Partitioning-based series estimators}
\label{sec:series}
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 \citet{cattaneo2020large} for further details and references throughout
this section. We illustrate the usefulness of
Corollary~\ref{cor:sa_martingale}
by deriving a Gaussian strong approximation for partitioning series estimators
based on multivariate martingale data. Proposition~\ref{pro:series} 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
\citep[Assumption~2]{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{cor:sa_martingale} directly to
construct a Gaussian
strong approximation. Specifically, we write
\begin{align*}
\hat\mu(w) - \mu(w)
&= p(w)^\ensuremath{\mathsf{T}} H^{-1} S
+ p(w)^\ensuremath{\mathsf{T}} \big(\hat H^{-1} - H^{-1}\big) S
+ \ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w),
\end{align*}
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
\citep{cattaneo2025yurinskiisupplement}}{Appendix~\ref{sec:proofs}})
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{cor:sa_martingale} 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 \citep{bradley2005basic}.
\begin{proposition}[Strong approximation for partitioning series estimators]
\label{pro:series}
Consider the nonparametric regression setup described above
and further assume the following:
\begin{enumerate}[label=(\roman*)]
\item
$(W_i, \varepsilon_i)_{1 \leq i \leq n}$
is strictly stationary.
\item
$W_1, \ldots, W_n$ is $\alpha$-mixing with mixing coefficients
satisfying $\sum_{j=1}^\infty \alpha(j) < \infty$.
\item
$W_i$ has a Lebesgue density on $\ensuremath{\mathcal{W}}$
which is bounded above and away from zero.
\item
$\ensuremath{\mathbb{E}}\big[|\varepsilon_i|^3 \big] < \infty$
and
$\ensuremath{\mathbb{E}}\big[\varepsilon_i^2 \mid \ensuremath{\mathcal{H}}_{i-1}\big]=\sigma^2(W_i)$
is bounded away from zero.
\item
$p(w)$ forms a basis with $k$ features satisfying
Assumptions~2 and~3 in \citet{cattaneo2020large}.
\end{enumerate}
Then, for any sequence $R_n \to \infty$,
there is a zero-mean Gaussian process
$G(w)$ indexed on $\ensuremath{\mathcal{W}}$
with $\operatorname{Var}[G(w)] \asymp\frac{k}{n}$
satisfying
$\ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[G(w), G(w')]
= \ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[p(w)^\ensuremath{\mathsf{T}} H^{-1} S,\, p(w')^\ensuremath{\mathsf{T}} H^{-1} S]$
and
\begin{align*}
\sup_{w \in \ensuremath{\mathcal{W}}}
\left| \hat\mu(w) - \mu(w) - G(w) \right|
&\lesssim_\ensuremath{\mathbb{P}}
\sqrt{\frac{k}{n}}
\left( \frac{k^3 (\log k)^3}{n} \right)^{1/6} R_n
+ \sup_{w \in \ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)|
\end{align*}
assuming the number of basis functions satisfies $k^3 / n \to 0$.
If further $\ensuremath{\mathbb{E}} \left[ \varepsilon_i^3 \mid \ensuremath{\mathcal{H}}_{i-1} \right] = 0$ then
\begin{align*}
\sup_{w \in \ensuremath{\mathcal{W}}}
\left| \hat\mu(w) - \mu(w) - G(w) \right|
&\lesssim_\ensuremath{\mathbb{P}}
\sqrt{\frac{k}{n}}
\left( \frac{k^3 (\log k)^2}{n} \right)^{1/4} R_n
+ \sup_{w \in \ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)|.
\end{align*}
\end{proposition}
The core of the proof of Proposition~\ref{pro:series} involves applying
Corollary~\ref{cor:sa_martingale} 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 \citet{cattaneo2020large} for i.i.d. data, which gives the
best known strong approximation rate for this problem. Thus,
Proposition~\ref{pro:series} gives the same best approximation rate, without
requiring any extra restrictions, for $\alpha$-mixing time series data.
Our results improve substantially on \citet[Theorem~1]{li2020uniform}: using
the notation of our Corollary~\ref{cor:sa_martingale}, and with any sequence
$R_n \to \infty$, a valid (see Remark~\ref{rem:coupling_bounds_probability})
version of their martingale Yurinskii coupling is
\begin{align*}
\|S-T\|_2
\lesssim_\ensuremath{\mathbb{P}}
d^{1/2} r^{1/2}_n
+ (B_n d)^{1/3} R_n,
\end{align*}
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{pro:series}, 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
\begin{align*}
\sup_{w \in \ensuremath{\mathcal{W}}}
\left| \hat\mu(w) - \mu(w) - G(w) \right|
&\lesssim_\ensuremath{\mathbb{P}}
\sqrt{\frac{k}{n}}
\left( \frac{k^5}{n} \right)^{1/6} R_n
+ \sup_{w \in \ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)|.
\end{align*}
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
\citet{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{pro:series}, 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{pro:series} and
anti-concentration for
Gaussian suprema \citep[Corollary~2.1]{chernozhukov2014anti} yield
a distributional approximation for the supremum statistic whenever
$k^3(\log n)^6 / n \to 0$, giving
\begin{align*}
\sup_{t \in \ensuremath{\mathbb{R}}}
\left|
\ensuremath{\mathbb{P}}\left(
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|
\frac{\hat\mu(w)-\mu(w)}{\sqrt{\rho(w,w)}}
\right| \leq t
\right)
-
\ensuremath{\mathbb{P}}\left(
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|
\frac{G(w)}{\sqrt{\rho(w,w)}}
\right| \leq t
\right)
\right|
&\to 0,
\end{align*}
where $\rho(w,w') = \ensuremath{\mathbb{E}}[G(w)G(w')]$. Furthermore, using a Gaussian--Gaussian
comparison result \citep[Lemma~3.1]{chernozhukov2013gaussian} and
anti-concentration again, it is not difficult to show (see the proof of
Proposition~\ref{pro:series}) that with $\ensuremath{\mathbf{W}} = (W_1, \ldots, W_n)$ and
$\ensuremath{\mathbf{Y}} = (Y_1, \ldots, Y_n)$,
\begin{align*}
\sup_{t \in \ensuremath{\mathbb{R}}}
\left|
\ensuremath{\mathbb{P}}\left(
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|
\frac{\hat\mu(w)-\mu(w)}{\sqrt{\hat\rho(w,w)}}
\right| \leq t
\right)
- \ensuremath{\mathbb{P}}\left(
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|
\frac{\hat G(w)}{\sqrt{\hat\rho(w,w)}}
\right| \leq t \biggm| \ensuremath{\mathbf{W}}, \ensuremath{\mathbf{Y}}
\right)
\right|
&\to_\ensuremath{\mathbb{P}} 0,
\end{align*}
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.
\begin{proposition}[Feasible uniform confidence bands for partitioning
series estimators]
\label{pro:series_feasible}
Assume the setup as described above. Then
\begin{align*}
\ensuremath{\mathbb{P}}\Big(
\mu(w) \in
\Big[
\hat\mu(w) \pm \hat q(\tau)
\sqrt{\hat\rho(w,w)}
\Big]
\ \text{for all }
w \in \ensuremath{\mathcal{W}} \Big)
\to 1-\tau,
\end{align*}
where
\vspace*{-2mm}
\begin{align*}
\hat{q}(\tau)
&=
\inf
\left\{
t \in \ensuremath{\mathbb{R}}:
\ensuremath{\mathbb{P}}\left(
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|
\frac{\hat G(w)}{\sqrt{\hat\rho(w,w)}}
\right|
\leq t
\Bigm| \ensuremath{\mathbf{W}}, \ensuremath{\mathbf{Y}}
\right)
\geq \tau
\right\}
\end{align*}
is the conditional quantile of the supremum of the Studentized Gaussian
process. This can be estimated by resampling the conditional law of
$\hat G(w) \mid \ensuremath{\mathbf{W}}, \ensuremath{\mathbf{Y}}$ with a discretization of $w \in \ensuremath{\mathcal{W}}$.
\end{proposition}
\subsection{Local polynomial estimators}
\label{sec:local_poly}
As a second example application we consider nonparametric
regression estimation
with martingale data employing local polynomial methods
\citep{Fan-Gijbels_1996_Book}. In contrast with the partitioning-based series
methods of Section~\ref{sec:series}, local polynomials induce stochastic
processes which are not linearly separable, allowing us to showcase the
empirical process result given in Proposition \ref{pro:emp_proc}.
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$,
\begin{align*}
\hat{\mu}(w)
&=
e_1^\ensuremath{\mathsf{T}}\hat{\beta}(w)
&\text{where} &
&\hat{\beta}(w)
&=
\operatorname*{arg\,min}_{\beta \in \ensuremath{\mathbb{R}}^{k}}
\sum_{i=1}^n
\left(Y_i - p_h(W_i-w)^\ensuremath{\mathsf{T}} \beta \right)^2
K_h(W_i-w).
\end{align*}
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:
\begin{align*}
\hat{\mu}(w)-\mu(w)
&= e_1^\ensuremath{\mathsf{T}} H(w)^{-1} S(w)
+ e_1^\ensuremath{\mathsf{T}} \big(\hat H(w)^{-1} - H(w)^{-1}\big) S(w)
+ \ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)
\end{align*}
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{sec:series}, for which neither $\hat H$ nor $S$ depend on $w$.
Therefore we use Proposition \ref{pro:emp_proc} 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{pro:emp_proc} by setting
\begin{align*}
\ensuremath{\mathcal{F}} = \left\{
(W_i, \varepsilon_i) \mapsto
e_1^\ensuremath{\mathsf{T}} H(w)^{-1}
K_h(W_i-w) p_h(W_i-w) \varepsilon_i
: w \in \ensuremath{\mathcal{W}}
\right\}.
\end{align*}
With this approach, we obtain the following result.
\begin{proposition}[Strong approximation for local polynomial estimators]
\label{pro:local_poly}
Under the nonparametric regression setup described above,
assume further that
\begin{enumerate}[label=(\roman*)]
\item
$(W_i, \varepsilon_i)_{1 \leq i \leq n}$
is strictly stationary.
\item
$(W_i, \varepsilon_i)_{1 \leq i \leq n}$
is $\alpha$-mixing with mixing coefficients
$\alpha(j) \leq e^{-2 j / C_\alpha}$
for some $C_\alpha > 0$.
\item
$W_i$ has a Lebesgue density on $\ensuremath{\mathcal{W}}$
which is bounded above and away from zero.
\item
$\ensuremath{\mathbb{E}}\big[e^{|\varepsilon_i|/C_\varepsilon}\big] < \infty$
for $C_\varepsilon > 0$ and
$\ensuremath{\mathbb{E}}\left[\varepsilon^2_i \mid \ensuremath{\mathcal{H}}_{i-1}\right]=\sigma^2(W_i)$
is bounded away from zero.
\item
$K$ is a non-negative Lipschitz
compactly supported kernel with
$\int K(w) \,\mathrm{d}w < \infty$.
\end{enumerate}
Then for any $R_n \to \infty$,
there is a zero-mean Gaussian process
$T(w)$ on $\ensuremath{\mathcal{W}}$
with $\operatorname{Var}[T(w)] \asymp\frac{1}{n h^m}$
satisfying
$\ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[T(w), T(w')]
= \ensuremath{\mathrm{Cov}\hspace*{0.2mm}}[e_1^\ensuremath{\mathsf{T}} H(w)^{-1} S(w),\, e_1^\ensuremath{\mathsf{T}} H(w')^{-1} S(w')]$
and
\begin{align*}
\sup_{w \in \ensuremath{\mathcal{W}}}
\left|\hat \mu(w) - \mu(w) - T(w) \right|
&\lesssim_\ensuremath{\mathbb{P}}
\frac{R_n}{\sqrt{n h^m}}
\left(
\frac{(\log n)^{m+4}}{n h^{3m}}
\right)^{\frac{1}{2m+6}}
+ \sup_{w \in \ensuremath{\mathcal{W}}} |\ensuremath{\mathrm{Bias}\hspace*{0.2mm}}(w)|,
\end{align*}
provided that the bandwidth sequence satisfies
$n h^{3m} \to \infty$.
\end{proposition}
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{pro:local_poly}; this is omitted here for
brevity. For completeness, the proof of Proposition~\ref{pro:local_poly}
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 \citep{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{pro:series} 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 \citet{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{sec:series}. 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{sec:kde}), 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{cor:sa_martingale}) does not impose a lower bound on the
eigenvalues of the variance matrix $\Sigma$. This effect was demonstrated by
Lemma \ref{lem:kde_eigenvalue} and its surrounding discussion in
Section~\ref{sec:kde}, and as such, the result of \citet{li2020uniform} is
unsuited for this application due to its strong minimum eigenvalue assumption.
Finally, for the special case of i.i.d.\ data,
\citet[Remark~3.1]{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{sec:emp_proc} for a
discussion), while \citet{cattaneo2024strong} obtain better strong
approximations for the entire stochastic process under additional
assumptions via a generalization of the celebrated Hungarian
construction \citep{komlos1975approximation,rio1994local}.
\section{Conclusion}
\label{sec:conclusion}
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.
\section*{Acknowledgments}
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.
\section*{Funding}
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.
\bibliographystyle{hapalike}
\bibliography{refs}
\clearpage