Robust Inference for Convex Pairwise Difference Estimators

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.

85,832 characters

Robust Inference for Convex Pairwise Difference Estimators



\title{\vspace{-0.0in} Robust Inference for Convex Pairwise Difference Estimators\thanks{This paper was prepared for the Econometric Theory Lecture delivered at the $2025$ International Symposium on Econometric Theory and Applications (SETA), University of Macau (China), June 1--3, 2025. It was also presented at the Econometrics Journal Lecture of the $2024$ $(\text{EC})^2$ Conference (Amsterdam), and the $2025$ Conference in Honor of Bo Honor{\'e}'s $65$th Birthday (Princeton University). We thank the participants at these conferences for their feedback. Cattaneo gratefully acknowledges financial support from the National Science Foundation through grants SES-1947805, DMS-2210561, and SES-2241575. Jansson gratefully acknowledges financial support from the National Science Foundation through grant SES-1947662 and from the Aarhus Center for Econometrics (ACE) funded by the Danish National Research Foundation grant number DNRF186. Nagasawa gratefully acknowledges financial support from the British Academy through grant SRG24$\backslash$241614.}
\bigskip }

\author{Matias D. Cattaneo\thanks{Department of Operations Research and Financial Engineering, Princeton University.} \and
	Michael Jansson\thanks{Department of Economics, UC Berkeley and ACE.} \and
	Kenichi Nagasawa\thanks{Department of Economics, University of Warwick.}}
\maketitle

\begin{abstract}
    This paper develops distribution theory and bootstrap-based inference methods for a broad class of convex pairwise difference estimators. These estimators minimize a kernel-weighted convex-in-parameter function over observation pairs that are similar in terms of certain covariates, where the similarity is governed by a localization (bandwidth) parameter. While classical results establish asymptotic normality under restrictive bandwidth conditions, we show that valid Gaussian and bootstrap-based inference remains possible under substantially weaker assumptions. First, we extend the theory of small bandwidth asymptotics to convex pairwise estimation settings, deriving robust Gaussian approximations even when a smaller than standard bandwidth is used. Second, we employ a debiasing procedure based on generalized jackknifing to enable inference with larger bandwidths, while preserving convexity of the objective function. Third, we construct a novel bootstrap method that adjusts for bandwidth-induced variance distortions, yielding valid inference across a wide range of bandwidth choices. Our proposed inference method enjoys demonstrable more robustness, while retaining the practical appeal of convex pairwise difference estimators.
\end{abstract}

\textit{Keywords:} small bandwidth asymptotics, generalized jackknife, bootstrap, U-process, pairwise comparisons, robust distribution theory.

\thispagestyle{empty}


\thispagestyle{empty}
\clearpage

\setcounter{page}{1}
\pagestyle{plain}

\section{Introduction}

Suppose $\mathbf{z}_1,\dots,\mathbf{z}_n$ is a random sample from the distribution of a random vector $\mathbf{z}$. This paper studies the large-sample properties of the following \textit{convex} pairwise difference estimator:
\begin{equation}\label{eq: Pairwise Difference Estimator}
    \widehat{\boldsymbol{\theta}}_n \in \operatorname*{arg\,min}_{\boldsymbol{\theta}\in\Theta} \binom{n}{2}^{-1} \sum_{i<j} m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) K_{h_n}(\mathbf{w}_i-\mathbf{w}_j), \qquad K_h(\mathbf{u})=\frac{1}{h^d}K\left(\frac{\mathbf{u}}{h}\right),
\end{equation}
where $\Theta \subseteq \mathbb{R}^k$ is a parameter space, $\sum_{i<j}$ denotes $\sum_{j=2}^n \sum_{i=1}^{j-1}$, $(\mathbf{z}_i, \mathbf{z}_j) \mapsto m(\mathbf{z}_i, \mathbf{z}_j; \boldsymbol{\theta})$ is a permutation symmetric function, $K$ is a symmetric, non-negative kernel, $h_n$ is a positive bandwidth (or localization) parameter sequence, $\mathbf{w}$ is a continuously distributed $d$-dimensional subvector of $\mathbf{z}$, and where $\boldsymbol{\theta} \mapsto m(\mathbf{z}_i, \mathbf{z}_j; \boldsymbol{\theta})$ is a \textit{convex} function. Pairwise difference estimation, which relies on local comparisons between observation pairs, has been used to address heterogeneity in nonlinear models. See \citet{Powell_1994_Handbook}, \citet{Honore-Powell_2005_Festschrift}, and \citet*{AradillasLopez-Honore-Powell_2007_IER} for overviews, and Section \ref{Section: Motivating Examples} for three motivating examples.

In contrast to classical extremum estimators, $\widehat{\boldsymbol{\theta}}_n$ is a local $M$-estimator that employs observation pairs $(i, j)$ for which $\mathbf{w}_i$ and $\mathbf{w}_j$ are similar. The bandwidth $h_n$ governs the degree of similarity: When $h_n \to 0$ (as $n \to \infty$), the estimator increasingly focuses on nearly identical-in-$\mathbf{w}$ pairs. In turn, focusing on such pairs is natural in settings where identification can be based on the condition $\mathbf{w}_i \approx \mathbf{w}_j$ (combined with smoothness assumptions). The localization introduces a familiar trade-off for estimation and inference: A smaller $h_n$ reduces bias from dissimilarity between $\mathbf{w}_i$ and $\mathbf{w}_j$, but increases variance due to fewer available usable pairs. As a consequence, the large-sample behavior of $\widehat{\boldsymbol{\theta}}_n$ depends critically on a delicate bias-variance trade-off determined by $h_n$. This paper develops novel inference methods for convex pairwise difference estimators that are demonstrably more robust to bandwidth choice than existing methods.

Under regularity conditions and assuming that
\begin{equation*}
    nh_n^d \to \infty \qquad\text{and}\qquad nh_n^4 \to 0,
\end{equation*}
the pairwise difference estimator is asymptotically linear:
\begin{equation}\label{eq: Asymptotic Linearity}
    \sqrt{n} (\widehat{\boldsymbol{\theta}}_n - \boldsymbol{\theta}_0) = \frac{1}{\sqrt{n}} \sum_{i=1}^n \boldsymbol{\xi}_0(\mathbf{z}_i) + o_\mathbb{P}(1) \rightsquigarrow \mathsf{N}(\mathbf{0}, \mathbb{E}[\boldsymbol{\xi}_0(\mathbf{z})\boldsymbol{\xi}_0(\mathbf{z})']),
\end{equation}
where $\boldsymbol{\theta}_0$ is the estimand and $\boldsymbol{\xi}_0(\cdot)$ is the influence function (whose exact form is given below).
The condition $nh_n^d \to \infty$ lower bounds the level of localization $h_n$ allowed for, while the condition $nh_n^4 \to 0$ upper bounds the level of localization. The purpose of the latter condition is to control a smoothing bias term. The bias condition $nh_n^{4} \to 0$ could be replaced by the weaker condition $nh_n^{2L} \to 0$ if a (higher-order) kernel of order $L>2$ were used, but a higher-order kernel annihilates the convexity of the objective function because higher-order kernels take negative values.

The main results of this paper are obtained by combining three ideas:
\begin{enumerate}
    \item \textit{Small Bandwidth Asymptotics}. Utilizing the framework introduced by \cite*{Cattaneo-Crump-Jansson_2014a_ET}, we establish a more robust Gaussian distributional approximation for the pairwise difference estimator that allows for higher levels of localization by remaining valid even when the condition $nh_n^d \to \infty$ is violated. This generalized distributional approximation shows that, while the localization restriction $nh_n^d \to \infty$ is necessary for establishing asymptotic linearity, a Gaussian approximation can hold under the substantially weaker condition $n^2h_n^d\to\infty$, albeit with a convergence rate and large sample variance that depends explicitly on the level of localization used.

    \item \textit{Debiasing}. Following \cite{Honore-Powell_2005_Festschrift} we debias the pairwise difference estimator using the method of \textit{generalized jackknifing} introduced by \cite{Schucany-Sommers_1977_JASA}. Doing so allows for (larger) bandwidths that violate the bias condition $nh_n^4 \to \infty$. This debiasing approach retains the convexity of the objective function, which is attractive for both theoretical (weaker regularity conditions) and practical (faster computation) reasons. The debiasing procedure combines linearly a collection of convex pairwise difference estimators constructed using different levels of localization. The resulting ensembling-based pairwise difference estimator admits a small bandwidth Gaussian approximation with an associated bias condition of the form $nh_n^{2L} \to 0$, where $L \geq 2$ denotes the order of a certain (equivalent) kernel induced by the debiasing procedure.

    \item \textit{Bootstrapping}. Building on insights in \citet*{Cattaneo-Crump-Jansson_2014b_ET}, we develop a valid bootstrap-based distributional approximation for the debiased pairwise difference estimator rescaling the localization parameter. The nonparametric bootstrap distributional approximation exhibits a mismatch in its asymptotic variance under small bandwidth asymptotics. The mismatch is characterized by a known multiplicative factor involving the localization parameter $h_n$. As a result, bootstrapping the (debiased) pairwise difference estimator with a different localization parameter (namely, $3^{1/d}h_n$ rather than $h_n$) leads to a valid bootstrap-based inference procedure also under small bandwidth asymptotics.
\end{enumerate}

In combination, these three ideas enable us to offer a novel resampling-based inference method for (convex) pairwise difference estimators that are demonstrably more robust to a wider set of choices of the localization parameter $h_n$.

Our theoretical work is carefully developed to retain and leverage convexity of the objective function defining the pairwise difference estimator. This feature not only allows for fast implementation of the estimator and resampling-based methods, but also enables us to proceed under relatively weak conditions when obtaining theoretical results. When developing our theoretical results, we rely heavily on the foundational work of \cite{Hjort-Pollard_1993} and \cite{Pollard_1991_ET}, which we apply to the case of $U$-processes.

This paper is connected to several strands of the literature. Contributions to the pairwise difference estimation literature include \cite{Ahn-Powell_1993_JOE}, \cite*{Ahn-Ichimura-Powell-Ruud_2018_JBES}, \cite{AradillasLopez_2012_JOE}, \cite{Blundell-Powell_2004_REStud}, \cite{Hong-Shum_2010_RESTUD}, \cite{Honore_1992_ECMA}, \cite*{Honore-Kyriazidou-Udry_1997_JOE}, \cite{Honore-Powell_1994_JoE}, \cite{Jochmans_2013_ECTJ}, and \cite{Kyriazidou_1997_ECMA}. The theoretical and practical features of small bandwidth asymptotics, and their connection with resampling methods for inference, are discussed in \citet*{Cattaneo-Crump-Jansson_2010_JASA}, \citet{Cattaneo-Crump-Jansson_2014b_ET}, \cite*{Cattaneo-Jansson-Newey_2018_ET}, \cite{Cattaneo-Jansson_2018_ECMA}, \cite{Matsushita-Otsu_2021_Biometrika}, \cite{Cattaneo-Jansson_2022_ET}, \citet*{Cattaneo-Farrell-Jansson-Masini_2025_JOE}, and references therein. The generalized jackknife has been successfully used for debiasing in density weighted average derivative estimation \citep*{Powell-Stock-Stoker_1989_ECMA}, asymptotically linear pairwise difference estimation \citep{Honore-Powell_2005_Festschrift}, nonlinear semiparametric estimation \citep*{Cattaneo-Crump-Jansson_2013_JASA}, monotone estimation \citep*{Cattaneo-Jansson-Nagasawa_2024_AOS}, and random forest estimation \citep*{Cattaneo-Klusowski-Underwood_2025_wp}, among other settings. \cite{Shao-Tu_2012_Book} give a textbook introduction to jackknifing, bootstrapping, and other resampling methods.

The rest of the paper proceeds as follows. Section \ref{Section: Motivating Examples} introduces the three motivating examples that are used throughout the paper to motivate our work and to illustrate the verification of the high-level assumptions imposed. Section \ref{Section: Distributional Approximation and Bootstrap Inference} present our main theoretical distributional and bootstrap results for robust inference employing convex pairwise difference estimators. The proofs of these results are given in Section \ref{Section: Proofs and Other Technical Results}. Section \ref{Section: Sufficient Conditions for Motivating Examples} showcases how the high-level sufficient conditions imposed in our theoretical developments are verified for the three motivating examples. Section \ref{Section: Conclusion} gives final remarks.

\section{Motivating Examples}\label{Section: Motivating Examples}

We use three examples to motivate and illustrate our work. The first example involves an estimator that can be written in closed form (because it has a quadratic-in-$\boldsymbol{\theta}$ function $m(\mathbf{z}_i, \mathbf{z}_j; \boldsymbol{\theta})$), while the other two examples do not. The second example has a smooth-in-$\boldsymbol{\theta}$ function $m(\mathbf{z}_i, \mathbf{z}_j; \boldsymbol{\theta})$, while the third example does not. All three examples have convex-in-$\boldsymbol{\theta}$ functions $m(\mathbf{z}_i, \mathbf{z}_j; \boldsymbol{\theta})$ and employ the following notation: $\mathbf{z}_i = (y_i, \mathbf{x}_i', \mathbf{w}_i')'$, with $y_i$ a scalar outcome variable, $\mathbf{x}_i$ a $k$-dimensional covariate, and $\mathbf{w}_i$ a $d$-dimensional covariate. For more details on the examples, see \cite{Powell_1994_Handbook}, \cite{Honore-Powell_2005_Festschrift}, and \citet{AradillasLopez-Honore-Powell_2007_IER}.

\subsection{Partially Linear Regression Model}

The partially linear regression model studied here is of the form
\begin{equation*}
    y_i =\mathbf{x}_i'\boldsymbol{\theta}_0 + \gamma_0(\mathbf{w}_i) + \varepsilon_i,
\end{equation*}
where $\boldsymbol{\theta}_0$ is the parameter of interest, $\gamma_0(\cdot)$ is an unknown function, and where $\mathbb{E}[\varepsilon_i|\mathbf{x}_i,\mathbf{w}_i]=0$. Defining $\dot{y}_{i,j}=y_i-y_j$ and $\dot{\mathbf{x}}_{i,j}=\mathbf{x}_i-\mathbf{x}_j$, a pairwise difference estimator of $\boldsymbol{\theta}_0$ can be based on
\begin{equation*}
    m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{1}{2}(\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} )^2.
\end{equation*}
Setting $\Theta=\mathbb{R}^k$, the minimization problem defining the estimator admits a closed form solution (provided that a non-negative kernel function is used), namely
\begin{equation*}
    \widehat{\boldsymbol{\theta}}_n = \left( \sum_{i<j} \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}'K_{h_n}(\mathbf{w}_i-\mathbf{w}_j) \right)^{-1} \sum_{i<j} \dot{\mathbf{x}}_{i,j}\dot{y}_{i,j}K_{h_n}(\mathbf{w}_i-\mathbf{w}_j).
\end{equation*}

\subsection{Partially Linear Logit Model}

The partially linear logit model studied here is of the form
\begin{equation*}
    y_i = \mathbbm{1} \{\mathbf{x}_i'\boldsymbol{\theta}_0 + \gamma_0(\mathbf{w}_i) + \varepsilon_i \geq 0 \},
\end{equation*}
where $\boldsymbol{\theta}_0$ is the parameter of interest, $\gamma_0(\cdot)$ is an unknown function, and where
\begin{equation*}
    \mathbb{P}\big[\varepsilon_i \leq u | \mathbf{x}_i,\mathbf{w}_i \big] = \Lambda(u), \qquad \Lambda(u) = \frac{\exp(u)}{1+\exp(u)}.
\end{equation*}
The parameter $\boldsymbol{\theta}_0$ can be estimated using a pairwise difference estimator with $\Theta=\mathbb{R}^k$ and
\begin{equation*}
	m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\mathbbm{1}\{\dot{y}_{i,j}\neq 0\} \left(y_i\ln\Lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}) + y_j\ln\Lambda(-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta})\right).
\end{equation*}
The minimization problem defining the estimator does not admit a closed form solution, but (provided that a non-negative kernel function is used) it is convex because $u\mapsto -\ln\Lambda(u)$ is.

\subsection{Partially Linear Tobit Model}

The partially linear censored regression model studied here is of the form
\begin{equation*}
    y_i =\max\{\mathbf{x}_i'\boldsymbol{\theta}_0+\gamma_0(\mathbf{w}_i)+\varepsilon_i,0\},
\end{equation*}
where $\boldsymbol{\theta}_0$ is the parameter of interest, $\gamma_0(\cdot)$ is an unknown function, $\mathbf{x}_i\protect\mathpalette{\protect\independenT}{\perp} \varepsilon_i|\mathbf{w}_i$, and where the conditional distribution of  $\varepsilon_i$ given $\mathbf{w}_i$ admits a Lebesgue density. A pairwise difference estimator of $\boldsymbol{\theta}_0$ can be obtained by setting $\Theta=\mathbb{R}^k$ and employing

\begin{equation*}
    m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = m_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})-\tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol0),
\end{equation*}
where
\begin{equation*}
    \tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) =
    \begin{cases}
        |y_i| - \big( \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} + y_j \big) \operatorname*{sgn}(y_i) & \text{if } \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} \leq -y_j \\
        \big|\dot{y}_{i,j} - \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} \big| & \text{if } -y_j < \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta} < y_i \\
        |y_j| + \big(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}-y_i \big) \operatorname*{sgn}(y_j) & \text{if } y_i \leq \dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}
    \end{cases}.
\end{equation*}
Because $\tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol0)$ does not depend on $\boldsymbol{\theta}$, the presence of $\tilde{m}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol0)$ in $m_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})$ does not affect the minimization problem defining the estimator. Nevertheless, it is theoretically attractive to work with $m_{\mathtt{PLT}}$ rather than $\tilde{m}_{\mathtt{PLT}}$, as doing so allows for weaker regularity conditions for the existence of the expectation of the objective function.

For future reference, we note that $m_{\mathtt{PLT}}$ admits the alternative representation
\begin{equation*}
    m_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) =
    \begin{cases}
        \big|\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}\big|-\big|\dot{y}_{i,j}\big| &\text{ if } y_i>0,y_j>0 \\
        \max\{ y_i-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}, 0\}-y_i &\text{ if } y_i>0,y_j=0 \\
        \max\{ y_j+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}, 0\}-y_j &\text{ if } y_i=0,y_j>0 \\
        0 &\text{ if } y_i=0,y_j=0
    \end{cases}.
\end{equation*}
The function $\boldsymbol{\theta}\mapsto m_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})$ is convex and therefore so is the minimization problem defining the estimator (provided that a non-negative kernel function is used).

\section{Distributional Approximation and Bootstrap Inference}\label{Section: Distributional Approximation and Bootstrap Inference}

As is standard in the literature, we generalize \eqref{eq: Pairwise Difference Estimator} slightly and define our estimator $\widehat{\boldsymbol{\theta}}_n=\widehat{\boldsymbol{\theta}}_n(h_n)$ to be any approximate minimizer of $\widehat{M}_n(\boldsymbol{\theta};h_n)$, where
\begin{equation*}
     \widehat{M}_n(\boldsymbol{\theta};h) = \binom{n}{2}^{-1}\sum_{i<j} m(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) K_h(\mathbf{w}_i-\mathbf{w}_j).
\end{equation*}
To be specific, we require
\begin{equation*}
    \widehat{M}_n(\widehat{\boldsymbol{\theta}}_n(h);h) \leq \inf_{\boldsymbol{\theta}\in\Theta} \widehat{M}_n(\boldsymbol{\theta};h) + o_{\mathbb{P}}( n^{-1} ).
\end{equation*}

The objective function $\widehat{M}_n$ is a sample counterpart of the function $M$ given by
\begin{equation*}
    M(\boldsymbol{\theta};h) = \mathbb{E}\big[\widehat{M}_n(\boldsymbol{\theta};h)\big]
                 = \mathbb{E}\big[ m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})  K_{h}(\mathbf{w}_1-\mathbf{w}_2)\big].
\end{equation*}
Under regularity conditions, this function approximates, as $h\downarrow 0$, a function $M_0$, which (does not depend on $K$ and) admits a unique minimizer, namely the parameter of interest $\boldsymbol{\theta}_0$.

For the purposes of analyzing $\widehat{\boldsymbol{\theta}}_n$ it is convenient to define $\boldsymbol{\theta}_n = \boldsymbol{\theta}(h_n)$, where
\begin{equation*}
    \boldsymbol{\theta}(h) \in \operatorname*{arg\,min}_{\boldsymbol{\theta}\in\Theta} M(\boldsymbol{\theta};h)
\end{equation*}
is interpretable as a (fixed-$h$) ``pseudo'' parameter. With the help of $\boldsymbol{\theta}_n$ we can decompose the estimation error $\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0$ into a (non-stochastic) ``bias'' component $\boldsymbol{\theta}_n-\boldsymbol{\theta}_0$ and a ``noise'' component $\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_n$. Each component can be analyzed separately and in both cases the analysis will leverage convexity.

\subsection{Regularity Conditions}

The following assumption guarantees, among other things, that $\boldsymbol{\theta}_n$ is well defined for large $n$ and that the bias component $\boldsymbol{\theta}_n-\boldsymbol{\theta}_0$ vanishes asymptotically; for details, see Lemma \ref{Lemma: Existence and Convergence of theta(h)} of Section \ref{Section: A Useful Lemma}.

\begin{assumption}\label{Assumption: Convergence of M}
    \begin{enumerate}[(i)]
        \item \label{Assumption: Convergence of M - kernel function} The kernel function $K$ is a symmetric, bounded probability density.

        \item \label{Assumption: Convergence of M - convexity} $\Theta\subseteq\mathbb{R}^k$ is convex, $(\mathbf{z},\bar{\mathbf{z}}) \mapsto m(\mathbf{z},\bar{\mathbf{z}};\boldsymbol{\theta})$ is permutation symmetric, and $\boldsymbol{\theta} \mapsto m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})$ is convex with probability one.

        \item \label{Assumption: Convergence of M - density of w} The distribution of $\mathbf{w}$ admits a Lebesgue density $f_{\mathbf{w}}$, which is bounded and continuous on its support $\mathcal{W}$.

        \item \label{Assumption: Convergence of M - well defined Mn} For each  $\boldsymbol{\theta}\in\Theta$,
        \begin{equation*}
            \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})] + \mathbb{E}\left[\sup_{\mathbf{v}\in\mathcal{W}}\mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{v} ]f_{\mathbf{w}}(\mathbf{v})\right]<\infty
        \end{equation*}
        and (with probability one)
        \begin{equation*}
            \lim_{\mathbf{u}\to \mathbf{0}} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}+\mathbf{u}] = \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}].
        \end{equation*}

        \item \label{Assumption: Convergence of M - well behaved M0} On $\Theta$, the function $M_0$ given by
        \begin{equation*}
            M_0(\boldsymbol{\theta}) = \int_{\mathcal{W}} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})^2d\mathbf{w}
        \end{equation*}
        is uniquely minimized at an interior point $\boldsymbol{\theta}_0$.

    \end{enumerate}
\end{assumption}

The next assumption enables us to analyze the asymptotic properties of the noise component $\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_n$. To accommodate examples (such as the partially linear Tobit model) where $\boldsymbol{\theta} \mapsto m(\mathbf{z}_i,\mathbf{z}_j,\boldsymbol{\theta})$ is not fully differentiable, we assume the existence of derivative-like functions $\mathbf{s}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})\in\mathbb{R}^k$ and $\mathbf{H}(\mathbf{w}_i,\mathbf{w}_j;\boldsymbol{\theta},\mathbf{t})\in\mathbb{R}^{k\times k}$ such that, for any direction $\mathbf{t} \in \mathbb{R}^k$, the (remainder) terms
\begin{equation*}
    r_\mathbf{t}(\boldsymbol{\theta},\tau)=\frac{m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}+\mathbf{t} \tau)-m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})}{\tau} - \mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})' \mathbf{t}
\end{equation*}
and
\begin{equation*}
    R_\mathbf{t}(\boldsymbol{\theta},\tau)=\frac{\mathbb{E}[r_\mathbf{t}(\boldsymbol{\theta},\tau)|\mathbf{w}_1,\mathbf{w}_2]}{\tau} - \frac{1}{2}\mathbf{t}'\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t})\mathbf{t}
\end{equation*}
are suitably small for $\boldsymbol{\theta}$ near $\boldsymbol{\theta}_0$, $\tau>0$ near zero, and $\mathbf{w}_1 \approx \mathbf{w}_2$. As further discussed below, functions $\mathbf{s}$ and $\mathbf{H}$ satisfying the following assumption exist (and are relatively easy to find) in each of our motivating examples.

\begin{assumption}\label{Assumption: Asymptotic Distribution}
    \begin{enumerate}[(i)]
        \item \label{Assumption: Asymptotic Distribution - differentiability} For each $\mathbf{t} \in\mathbb{R}^k$, there is some $\delta>0$ such that
        \begin{align*}
            & \mathbb{E}\left[\sup_{ \tau\in (0,\delta),\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta, \mathbf{w}_2 \in\mathcal{W}}
                      \big|\mathbb{E}[ r_\mathbf{t}(\boldsymbol{\theta},\tau) | \mathbf{z}_1,\mathbf{w}_2 ] \big|f_{\mathbf{w}}(\mathbf{w}_2)^2 \right]< \infty,\\
            & \mathbb{E}\left[\sup_{ \tau\in (0,\delta),\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta, \mathbf{w}_2 \in\mathcal{W}}
                      \mathbb{E}[r_\mathbf{t}(\boldsymbol{\theta},\tau)^2 | \mathbf{w}_1,\mathbf{w}_2 ] f_{\mathbf{w}}(\mathbf{w}_2) \right] < \infty,\\
            &\mathbb{E}\left[\sup_{\tau\in (0,\delta),\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta,  \mathbf{w}_2 \in\mathcal{W}}
                     \left|R_\mathbf{t}(\boldsymbol{\theta},\tau)\right| f_{\mathbf{w}}(\mathbf{w}_2) \right]<\infty,
        \end{align*}
        and (with probability one)
        \begin{align*}
            & \lim_{\tau\downarrow0,(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})}
              \mathbb{E}[ r_\mathbf{t}(\boldsymbol{\theta},\tau) | \mathbf{z}_1=\mathbf{z},\mathbf{w}_2 =\mathbf{w}+\mathbf{u} ] = 0, \\
            & \lim_{\tau\downarrow0,(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})}
              \mathbb{E}[ r_\mathbf{t}(\boldsymbol{\theta},\tau)^2 | \mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}+\mathbf{u} ] = 0, \\
            & \lim_{\tau\downarrow0,(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})}
              \mathbb{E}[ R_\mathbf{t}(\boldsymbol{\theta},\tau) | \mathbf{w}_1=\mathbf{w},\mathbf{w}_2 =\mathbf{w}+\mathbf{u} ] = 0.
        \end{align*}

        \item \label{Assumption: Asymptotic Distribution - moment bounds} There is some $\delta>0$ and some function $b$ with
        \begin{equation*}
            \sup_{\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta} \|\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})\|\leq b(\mathbf{z}_1)b(\mathbf{z}_2),
        \end{equation*}
        such that
        \begin{equation*}
            \mathbb{E}[b(\mathbf{z})^4] + \sup_{\mathbf{w}\in\mathcal{W}}\mathbb{E}[b(\mathbf{z})^4|\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})<\infty
        \end{equation*}
        and
        \begin{equation*}
            \mathbb{E}\left[\sup_{\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta, \mathbf{w}_2\in\mathcal{W}}                       \|\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t})\| f_{\mathbf{w}}(\mathbf{w}_2) \right] < \infty \qquad \text{for each } \mathbf{t} \in \mathbb{R}^k.
        \end{equation*}

        \item \label{Assumption: Asymptotic Distribution - variance ingredients} There exist functions $\mathbf{G}_0,\boldsymbol{\xi}_0$, and $\boldsymbol{\Xi}_0$ such that, for each $\mathbf{t}\in\mathbb{R}^k$ (and with probability one),
        \begin{equation*}
            \lim_{(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})} \mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta},\mathbf{t}) f_{\mathbf{w}}(\mathbf{w}) = \mathbf{G}_0(\mathbf{w}),
        \end{equation*}
        \begin{equation*}
            \lim_{(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})} 2\mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w}+\mathbf{u} ] f_{\mathbf{w}}(\mathbf{w}) = \boldsymbol{\xi}_0(\mathbf{z}),
        \end{equation*}
        and
        \begin{equation*}
            \lim_{(\boldsymbol{\theta},\bar{\boldsymbol{\theta}},\mathbf{u})\to (\boldsymbol{\theta}_0,\boldsymbol{\theta}_0,\mathbf{0}) } \mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\bar{\boldsymbol{\theta}})'|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}+\mathbf{u} ] f_{\mathbf{w}}(\mathbf{w})= \boldsymbol{\Xi}_0(\mathbf{w}).
        \end{equation*}

        \item \label{Assumption: Asymptotic Distribution - nonsingularity} $\boldsymbol{\Gamma}_0 = \mathbb{E}[\mathbf{G}_0(\mathbf{w})], \boldsymbol{\Sigma}_0 = \mathbb{E}[\boldsymbol{\xi}_0(\mathbf{z})\boldsymbol{\xi}_0(\mathbf{z})']$, and $\mathbb{E}[\boldsymbol{\Xi}_0(\mathbf{w})]$ are positive definite.
    \end{enumerate}
\end{assumption}

\subsection{Small Bandwidth Asymptotics}

Defining
\begin{equation*}
    \mathbf{V}_n = \mathbf{V}_n(h_n)=\boldsymbol{\Gamma}_0^{-1}\left[ n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(K)\right]\boldsymbol{\Gamma}_0^{-1}, \qquad \boldsymbol{\Delta}_0(K) = \mathbb{E}[ \boldsymbol{\Xi}_0(\mathbf{w})] \int_{\mathbb{R}^d} K^2(\mathbf{u})d\mathbf{u},
\end{equation*}
and letting $\Phi_k$ denote the distribution function of a $k$-dimensional standard Gaussian random vector, we have the following result.

\begin{thm}\label{Theorem: Asymptotic Distribution}
    Suppose Assumptions \ref{Assumption: Convergence of M} and \ref{Assumption: Asymptotic Distribution} hold. If $n^2h_n^d\to\infty$ and if $h_n\to 0$, then
    \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \mathbf{V}_n^{-1/2}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0.
    \end{equation*}
\end{thm}

Under the assumptions of Theorem \ref{Theorem: Asymptotic Distribution}, the convergence rate of $\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_n$ equals
\begin{equation*}
    \rho_n = \sqrt{\min\left(n,\binom{n}{2}h_n^d\right)},
\end{equation*}
the magnitude of $\mathbf{V}_n^{-1/2}$. Provided that the bias is ``small'' in the sense that $\|\boldsymbol{\theta}_n-\boldsymbol{\theta}_0\|=o(\rho_n^{-1})$, Theorem \ref{Theorem: Asymptotic Distribution} therefore encompasses the following three distinct large-sample regimes:
\begin{itemize}
    \item \textit{Asymptotic Linearity}: If $nh_n^d\to\infty$, then $\sqrt{n}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0)$ satisfies \eqref{eq: Asymptotic Linearity}. In particular, $\sqrt{n}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0)$ convergences in law to a mean-zero Gaussian distribution with asymptotic variance
    \begin{equation*}
        \lim_{n\to\infty} n \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1} \boldsymbol{\Sigma}_0 \boldsymbol{\Gamma}_0^{-1}.
    \end{equation*}

    \item \textit{Root-$n$ Consistency without Asymptotic Linearity}: If $nh_n^d\to 2c \in(0,\infty)$, then $\sqrt{n}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0)$ is not asymptotically linear, but convergences in law to a mean-zero Gaussian distribution with asymptotic variance
    \begin{equation*}
        \lim_{n\to\infty} n \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1}\left[ \boldsymbol{\Sigma}_0 + \frac{1}{c} \boldsymbol{\Delta}_0(K)\right]\boldsymbol{\Gamma}_0^{-1}.
    \end{equation*}

    \item \textit{Slower than Root-n Consistency}: If $nh_n^d\to 0$ (but $n^2h_n^d\to\infty$), then $\sqrt{n^2h^d/2}(\widehat{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0)$ converges weekly to a mean-zero Gaussian distribution with asymptotic variance
    \begin{equation*}
        \lim_{n\to\infty} \binom{n}{2}h_n^d \mathbf{V}_n(h_n) = \boldsymbol{\Gamma}_0^{-1}\boldsymbol{\Delta}_0(K)\boldsymbol{\Gamma}_0^{-1}.
    \end{equation*}
\end{itemize}
The small bandwidth component (i.e., the term involving $\boldsymbol{\Delta}_0(K)$) in $\mathbf{V}_n$ captures the additional uncertainty generated from increasing the localization of the observations pairs. Incorporating this component in the approximate variance is key to enabling us to replace the condition $nh_n^d\to\infty$ by the weaker condition $n^2h_n^d\to\infty$ when obtaining a Gaussian approximation. As demonstrated by \cite{Cattaneo-Farrell-Jansson-Masini_2025_JOE}, incorporating the small bandwidth component can furthermore lead to a higher-order corrected distributional approximation even under asymptotic linearity.

\subsection{Debiasing}

In Theorem \ref{Theorem: Asymptotic Distribution}, we centered the estimator $\widehat{\boldsymbol{\theta}}_n=\widehat{\boldsymbol{\theta}}_n(h_n)$ at $\boldsymbol{\theta}_n=\boldsymbol{\theta}(h_n)$ to circumvent bias issues. This section focuses on the bias term $\boldsymbol{\theta}_n-\boldsymbol{\theta}_0$ and introduces an automatic debiasing approach under the assumption that $\boldsymbol{\theta}_n-\boldsymbol{\theta}_0$ can be expanded in even powers of $h_n$. To be specific, we follow \citet[Section 3.3]{Honore-Powell_2005_Festschrift} and dicuss debiasing under the following high-level condition.

\begin{assumption}\label{Assumption: Bias of thetahat}
    For some even $L\geq0$, $\boldsymbol{\theta}(\cdot)$ admits $\mathbf{b}_{2l}\in\mathbb{R}^k$ (for $l=1,\dots, L/2$) such that
    \begin{equation*}
        \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l} + o(h^L)\qquad \text{as } h\downarrow 0.
    \end{equation*}
\end{assumption}

The ease with which Assumption \ref{Assumption: Bias of thetahat} can be verified depends on the magnitude of $L$. For instance, Assumption \ref{Assumption: Convergence of M} implies that Assumption \ref{Assumption: Bias of thetahat} holds with $L=0$. Under additional smoothness conditions and using symmetry of $K$, the following result gives conditions under which Assumption \ref{Assumption: Bias of thetahat} holds with $L=2$. When stating the result, we employ standard multi-index notation: for $\boldsymbol{\alpha} = (\alpha_1,\ldots,\alpha_d)'\in\mathbb{Z}_+^d$, $\mathbf{v} = (v_1,\ldots,v_d)'\in\mathbb{R}^d$, and a sufficiently smooth-in-$\mathbf{v}$ function $f(\mathbf{w},\mathbf{v}),$
\begin{equation*}
    \partial_{\mathbf{v}}^{\boldsymbol{\alpha}} f(\mathbf{w},\mathbf{v}) = \frac{\partial^{|\boldsymbol{\alpha}|}}{\partial v_1^{\alpha_1} \cdots \partial v_d^{\alpha_d}} f(\mathbf{w},\mathbf{v}), \qquad |\boldsymbol{\alpha}| = \sum_{j=1}^d \alpha_j.
\end{equation*}

\begin{prop}\label{Proposition: Bias Expansion L=2}
    Suppose Assumptions \ref{Assumption: Convergence of M}-\ref{Assumption: Asymptotic Distribution} hold and that
    \begin{enumerate}[(i)]
        \item $\int_{\mathbb{R}^d} \|\mathbf{u}\|^2 K(\mathbf{u})d\mathbf{u}<\infty$, and
        \item With probability one, $\mathbf{v} \mapsto \boldsymbol{\psi}(\mathbf{w},\mathbf{v}) = \mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{v}]f_{\mathbf{w}}(\mathbf{v})$ is twice continuously differentiable with $\mathbb{E}[\sup_{\mathbf{v}\in\mathcal{W}}\|\partial_{\mathbf{v}}^{\boldsymbol{\alpha}} \boldsymbol{\psi}(\mathbf{w},\mathbf{v})\|]<\infty$ for all $\boldsymbol{\alpha}\in\mathbb{Z}_+^d$ with $|\boldsymbol{\alpha}|\leq 2$.
    \end{enumerate}
    Then $\boldsymbol{\theta}(\cdot)$ admits a $\mathbf{b}_2\in\mathbb{R}^k$ such that
    \begin{equation*}
        \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 =  \mathbf{b}_2 h^2 + o(h^2) \qquad \text{ as } h\downarrow 0.
    \end{equation*}
\end{prop}

The proof of Proposition \ref{Proposition: Bias Expansion L=2} leverages convexity and may therefore be of independent interest. The convexity argument in question can furthermore be adapted to form the basis of a verification by induction of Assumption \ref{Assumption: Bias of thetahat} with $L>2$. Details are provided in Section \ref{Section: Verifying Assumption 3}, which describes the induction step for general $L$ and states explicit (smoothness) conditions under which Assumption \ref{Assumption: Bias of thetahat} holds with $L=4$.

To describe the debiasing procedure based on generalized jackknifing, we maintain Assumption \ref{Assumption: Bias of thetahat}, define $c_0=1$, and let $\mathbf{c}=(c_0,\dots, c_{L/2})'$ be a vector of (distinct) positive constants such that the vector
\begin{equation*}
    \begin{pmatrix}
        \lambda_0(\mathbf{c}) \\
        \lambda_1(\mathbf{c}) \\
        \vdots \\
        \lambda_{L/2}(\mathbf{c})
    \end{pmatrix}
    =
    \begin{pmatrix}
        1 & 1 &\dots &1 \\
        1 & c_1^2 &\dots &c_{L/2}^2 \\
        \vdots & &\ddots & \\
        1& c_1^L & \dots & c_{L/2}^L
    \end{pmatrix}^{-1}
    \begin{pmatrix}
        1 \\
        0 \\
        \vdots \\
        0
    \end{pmatrix}
\end{equation*}
is well defined. The debiased estimator is
\begin{equation*}
    \widetilde{\boldsymbol{\theta}}_n = \widetilde{\boldsymbol{\theta}}_n(\mathbf{c},h_n) = \sum_{l=0}^{L/2} \lambda_{l}(\mathbf{c}) \widehat{\boldsymbol{\theta}}_n(h_{n,l}), \qquad h_{n,l}=c_lh_n,
\end{equation*}
the construction of which involves solving $L/2+1$ convex optimization problems. As defined, the debiased estimator is a generalization of the original pairwise difference estimator because if $L=0$, then $\mathbf{c}=1=\lambda_0$ and therefore $\widetilde{\boldsymbol{\theta}}_n = \widehat{\boldsymbol{\theta}}_n$.

The next theorem generalizes Theorem \ref{Theorem: Asymptotic Distribution} by establishing the small bandwidth Gaussian approximation for $\widetilde{\boldsymbol{\theta}}_n$. To state the theorem, let
\begin{equation*}
    \bar{\boldsymbol{\theta}}_n = \bar{\boldsymbol{\theta}}_n(\mathbf{c},h_n) = \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \boldsymbol{\theta}(h_{n,l})
\end{equation*}
and
\begin{equation*}
    \bar{\mathbf{V}}_n = \bar{\mathbf{V}}_n(\mathbf{c},h_n) = \boldsymbol{\Gamma}_0^{-1}\left[n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1} h_n^{-d}\boldsymbol{\Delta}_0(\bar{K}) \right]\boldsymbol{\Gamma}_0^{-1}, \qquad \bar{K}(\mathbf{u}) = \bar{K}(\mathbf{u};\mathbf{c}) = \sum_{l=0}^{L/2}\lambda_{l}(\mathbf{c})K_{c_l}(\mathbf{u}).
\end{equation*}
As they should, the expressions have the feature that if $L=0$, then $\bar{\boldsymbol{\theta}}_n = \boldsymbol{\theta}_n$ and $\bar{\mathbf{V}}_n=\mathbf{V}_n$. Another noteworthy feature of the expressions is that debiasing via generalized jackknifing affects the variance $\bar{\mathbf{V}}_n$ only through the kernel shape entering its small bandwidth component.

\begin{thm}\label{Theorem: Generalized Jackknifing}
    Suppose Assumptions \ref{Assumption: Convergence of M} and \ref{Assumption: Asymptotic Distribution} hold. If $n^2h_n^d\to\infty$ and if $h_n\to 0$, then
    \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \bar{\mathbf{V}}_n^{-1/2}(\widetilde{\boldsymbol{\theta}}_n-\bar{\boldsymbol{\theta}}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0.
    \end{equation*}
    As a consequence, if also Assumption \ref{Assumption: Bias of thetahat} holds and if $nh_n^{2L}\to0$, then
    \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}\left[ \bar{\mathbf{V}}_n^{-1/2}(\widetilde{\boldsymbol{\theta}}_n-\boldsymbol{\theta}_0 ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to 0
    \end{equation*}
\end{thm}

The magnitude of $\bar{\mathbf{V}}_n^{-1/2}$ is the same as that of $\mathbf{V}_n^{-1/2}$. As a consequence, with obvious modifications the discussion of $\widehat{\boldsymbol{\theta}}_n$ following Theorem \ref{Theorem: Asymptotic Distribution} applies to $\widetilde{\boldsymbol{\theta}}_n$, the only noteworthy difference being that (by design), the relevant ``small bias'' condition is different (and typically milder) in the case of $\widetilde{\boldsymbol{\theta}}_n$.

It is worth noting that the equivalent kernel $\bar{K}$ is of higher order, even though the debiased estimator $\widetilde{\boldsymbol{\theta}}_n$ only employs estimators constructed using second-order kernels, hereby retaining the desired convexity for implementation. To be specific, if $\int_{\mathbb{R}^d} \|\mathbf{u}\|^{L+2} K(\mathbf{u})d\mathbf{u} < \infty$, then
\begin{equation*}
    \int_{\mathbb{R}^d} \bar{K}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L/2}\lambda_l(\mathbf{c}) \int_{\mathbb{R}^d} K_{c_{l}}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L/2}\lambda_l(\mathbf{c}) = 1
\end{equation*}
and, for $\boldsymbol{\alpha}=(\alpha_1,\dots,\alpha_d)'\in\mathbb{Z}^d_+$ with $0<|\boldsymbol{\alpha}|\leq L+1$,
\begin{equation*}
    \int_{\mathbb{R}^d} \mathbf{u}^{\boldsymbol{\alpha}} \bar{K}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \int_{\mathbb{R}^d} \mathbf{u}^{\boldsymbol{\alpha}} K_{c_l}(\mathbf{u})d\mathbf{u} = \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) c_l^{|\boldsymbol{\alpha}|} \int_{\mathbb{R}^d} \mathbf{v}^{\boldsymbol{\alpha}} K(\mathbf{v})d\mathbf{v} =0
\end{equation*}
where the last equality uses the defining property of $\{\lambda_l(\mathbf{c})\}$ and symmetry of $K$, and where $\mathbf{u}^{\boldsymbol{\alpha}}$ denotes $\prod_{j=1}^d u_j^{\alpha_j}$ for $\mathbf{u}=(u_1,\dots,u_d)'\in\mathbb{R}^d$. In other words, $\bar{K}$ is of order $L+2$.

\subsection{Bootstrapping}\label{Section: Bootstrapping}

To develop feasible inference procedures that do not require (explicit) estimation of $\bar{\mathbf{V}}_n$, we consider nonparametric bootstrap-based approximations to the distribution of $\widetilde{\boldsymbol{\theta}}_n$. Since $\widetilde{\boldsymbol{\theta}}_n=\widehat{\boldsymbol{\theta}}_n$ when $L=0$, results for $\widehat{\boldsymbol{\theta}}_n$ can be extracted by setting $L=0$ in what follows.

Letting $\mathbf{z}^*_{1,n},\dots,\mathbf{z}^*_{n,n}$ denote a random sample from the empirical distribution of $\mathbf{z}_1,\dots,\mathbf{z}_n$, the defining property of $\widehat{\boldsymbol{\theta}}_n^*(h)$, the nonparametric bootstrap analogue of $\widehat{\boldsymbol{\theta}}_n(h)$, is the following:
\begin{equation*}
    \widehat{M}_n^*(\widehat{\boldsymbol{\theta}}_n^*(h);h) \leq \inf_{\boldsymbol{\theta}\in\Theta} \widehat{M}_n^*(\boldsymbol{\theta};h)  + o_{\mathbb{P}}(n^{-1}),
\end{equation*}
where
\begin{equation*}
    \widehat{M}_n^*(\boldsymbol{\theta};h) = \binom{n}{2}^{-1} \sum_{i<j}m(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n};\boldsymbol{\theta}) K_h(\mathbf{w}^*_{i,n}-\mathbf{w}^*_{j,n}).
\end{equation*}
Similarly, the nonparametric bootstrap analogue of $\widetilde{\boldsymbol{\theta}}_n$ is
\begin{equation*}
    \widetilde{\boldsymbol{\theta}}_n^* = \widetilde{\boldsymbol{\theta}}_n^*(\mathbf{c})= \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}^*(h_{n,l}).
\end{equation*}

The following theorem characterizes the large sample properties of $\widetilde{\boldsymbol{\theta}}_n^*-\widehat{\boldsymbol{\theta}}_n^*$, the bootstrap counterpart of $\widehat{\boldsymbol{\theta}}_n^*-\boldsymbol{\theta}_n$. In perfect analogy with the results in \cite{Cattaneo-Crump-Jansson_2014b_ET}, we find that the bootstrap distribution estimator is consistent only when $nh_n^d\to\infty$, but otherwise exhibits a variance inflation making the distributional approximation inconsistent. To state the result, let $\mathbb{P}_n^*[\cdot]$ denote $\mathbb{P}[\cdot | \mathbf{z}_1,\ldots,\mathbf{z}_n]$, let $\to_\mathbb{P}$ denote convergence in probability, and define
\begin{equation*}
    \bar{\mathbf{V}}_n^* = \bar{\mathbf{V}}_n^*(\mathbf{c},h_n) = \boldsymbol{\Gamma}_0^{-1}\left[n^{-1}\boldsymbol{\Sigma}_0 + 3\binom{n}{2}^{-1} h_n^{-d}\boldsymbol{\Delta}_0(\bar{K}) \right]\boldsymbol{\Gamma}_0^{-1}.
\end{equation*}

\begin{thm}\label{Theorem: Bootstrapping}
    Suppose Assumptions \ref{Assumption: Convergence of M}-\ref{Assumption: Asymptotic Distribution} hold and that, for $\boldsymbol{\theta}$ near $\boldsymbol{\theta}_0$, $m(\mathbf{z},\mathbf{z};\boldsymbol{\theta})=0$ and $\mathbf{s}(\mathbf{z},\mathbf{z};\boldsymbol{\theta})=\mathbf{0}$ (with probability one). If $n^2h_n^d\to\infty$ and if $h_n\to0$, then
    \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}^*\left[ \bar{\mathbf{V}}_n^{*-1/2}(\widetilde{\boldsymbol{\theta}}_n^*-\widetilde{\boldsymbol{\theta}}_n ) \leq \mathbf{t} \right] - \Phi_k(\mathbf{t}) \right| \to_\mathbb{P} 0.
    \end{equation*}
\end{thm}

Because $\bar{\mathbf{V}}_n^{-1}\bar{\mathbf{V}}_n^*\to \mathbf{I}_k$ if and only if $nh_n^d\to\infty$ (where $\mathbf{I}_k$ denotes the $k$-dimensional identity matrix), under the assumptions of Theorem \ref{Theorem: Bootstrapping}
\begin{equation*}
    \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^*\left[ \widetilde{\boldsymbol{\theta}}_n^*-\widetilde{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \bar{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] \right| \to_\mathbb{P} 0
\end{equation*}
if and only if $nh_n^d\to\infty$. In particular, if $\liminf_{n\to\infty} nh_n^d<\infty$, then the nonparametric bootstrap is inconsistent, albeit conservative in the sense that the (approximate) variance under the bootstrap distribution is larger than the (approximate) variance of the asymptotic distribution: $\bar{\mathbf{V}}_n^* > \bar{\mathbf{V}}_n$ in a positive definite sense.

The variance inflation problem associated with the nonparametric bootstrap under the small bandwidth regime can be easily fixed by appropriately rescaling the bandwidth used for the bootstrap implementation of the pairwise estimator: employing
\begin{equation*}
    \breve{\boldsymbol{\theta}}_n^* = \breve{\boldsymbol{\theta}}_n^*(\mathbf{c},h_n)= \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}^*(3^{1/d}h_{n,l}).
\end{equation*}
and centering its distribution at
\begin{equation*}
    \breve{\boldsymbol{\theta}}_n = \breve{\boldsymbol{\theta}}_n(\mathbf{c},h_n)= \sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \widehat{\boldsymbol{\theta}}_{n}(3^{1/d}h_{n,l}).
\end{equation*}
automatically adjusts the bootstrap variance, leading to a consistent distributional approximation. Indeed, the following result is an immediate consequence of Theorems \ref{Theorem: Generalized Jackknifing} and \ref{Theorem: Bootstrapping}.

\begin{coro}\label{Corollary: Bootstrapping}
   If the assumptions of Theorem \ref{Theorem: Bootstrapping} hold, then
   \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^*\left[ \breve{\boldsymbol{\theta}}_n^* - \breve{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \bar{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] \right| \to_\mathbb{P} 0.
    \end{equation*}
    As a consequence, if also Assumption \ref{Assumption: Bias of thetahat} holds and if $nh_n^{2L}\to0$, then
    \begin{equation*}
        \sup_{\mathbf{t}\in\mathbb{R}^k} \left| \mathbb{P}_n^*\left[ \breve{\boldsymbol{\theta}}_n^* - \breve{\boldsymbol{\theta}}_n \leq \mathbf{t} \right] - \mathbb{P}\left[ \widetilde{\boldsymbol{\theta}}_n - \boldsymbol{\theta}_{0} \leq \mathbf{t} \right] \right| \to_\mathbb{P} 0.
    \end{equation*}
\end{coro}

The statement of Corollary \ref{Corollary: Bootstrapping} emphasizes the rate-adaptive nature of the consistency property enjoyed by the bootstrap distributional approximation. The result has immediate implications for robust inference. For example, letting $\alpha \in (0,1)$, $\mathbf{a}\in\mathbb{R}^k$ be a fixed vector, and using the ``percentile method'' \citep[in the terminology of][]{vanderVaart_1998_Book}, the (nominal) level $1 - \alpha$ bootstrap confidence interval for $\mathbf{a}'\boldsymbol{\theta}_0$ is
\begin{equation*}
    \breve{\mathsf{CI}}^*_n(1-\alpha) = \left[\mathbf{a}'\widetilde{\boldsymbol{\theta}}_n - \breve{q}^*_{1-\alpha/2,n} ~ , ~ \mathbf{a}'\widetilde{\boldsymbol{\theta}}_n - \breve{q}^*_{\alpha/2,n} \right], \qquad \breve{q}^*_{t,n} = \inf\left\{ q \in \mathbb{R} : \mathbb{P}_n^* [\mathbf{a}'\breve{\boldsymbol{\theta}}_n^* - \mathbf{a}'\breve{\boldsymbol{\theta}}_n \leq q] \geq t \right\}.
\end{equation*}
If Assumptions \ref{Assumption: Convergence of M}-\ref{Assumption: Bias of thetahat} hold and if $n^2h_n^d\to\infty$ and $nh_n^{2L}\to0$, then
\begin{equation*}
    \lim_{n \to \infty} \mathbb{P}\left[\mathbf{a}'\boldsymbol{\theta}_0 \in \breve{\mathsf{CI}}^*_n(1 - \alpha)\right] = 1 - \alpha.
\end{equation*}

\section{Proofs and Other Technical Results}\label{Section: Proofs and Other Technical Results}

\subsection{A Useful Lemma}\label{Section: A Useful Lemma}

The following lemma is used in the proofs of Theorems \ref{Theorem: Asymptotic Distribution}-\ref{Theorem: Bootstrapping}.

\begin{lem}\label{Lemma: Existence and Convergence of theta(h)}
    Suppose that Assumption \ref{Assumption: Convergence of M} holds. Then $\operatorname*{arg\,min}_{\boldsymbol{\theta}\in\Theta}M(\boldsymbol{\theta};h)$ is non-empty for $h>0$ near zero and
    \begin{equation*}
        \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 = o(1) \qquad \text{as } h \downarrow 0.
    \end{equation*}
    If also Assumption \ref{Assumption: Asymptotic Distribution} holds, then $\mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}(h) ) K_h(\mathbf{w}_1-\mathbf{w}_2)]=\mathbf{0}$ for $h>0$ near zero.
\end{lem}

\begin{proof}
    For every $\boldsymbol{\theta} \in \Theta$,
    \begin{align*}
        M(\boldsymbol{\theta};h) &= \int_\mathcal{W} \int_{\mathbb{R}^d} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{w}-\mathbf{u} h] f_\mathbf{w}(\mathbf{w}) f_\mathbf{w}(\mathbf{w} - \mathbf{u} h) K(\mathbf{u}) d\mathbf{u} d\mathbf{w} \\
        &\to \int_\mathcal{W} \mathbb{E}[m(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta})|\mathbf{w}_1=\mathbf{w}_2=\mathbf{w}] f_\mathbf{w}(\mathbf{w})^2 d\mathbf{w} = M_0(\boldsymbol{\theta}) \qquad \text{as } h \downarrow 0,
    \end{align*}
    the convergence being uniform on compact subsets of $\Theta$ because $\boldsymbol{\theta} \mapsto M(\boldsymbol{\theta};h)$ is convex \citep[e.g.,][Lemma 1]{Hjort-Pollard_1993}.

    Take any $\epsilon > 0$ with $\Theta^\epsilon_0 = \{\boldsymbol{\theta} \in \mathbb{R}^k: \|\boldsymbol{\theta}-\boldsymbol{\theta}_0\| \leq \epsilon\} \subseteq \Theta$. By the preceding paragraph,
    \begin{equation*}
        \sup_{\boldsymbol{\theta} \in \Theta^\epsilon_0} |M(\boldsymbol{\theta};h) - M_0(\boldsymbol{\theta})| \leq \frac{1}{2} \left(\inf_{\boldsymbol{\theta}:\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|=\epsilon} M_0(\boldsymbol{\theta}) - M_0(\boldsymbol{\theta}_0)\right)
    \end{equation*}
    for $h>0$ near zero. For any such $h$ and any $\boldsymbol{\theta} \in \Theta \setminus \Theta^\epsilon_0$, we have
    \begin{equation*}
        \eta = \frac{\epsilon}{\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|} \in (0,1),
    \end{equation*}
    and therefore, by convexity of $\boldsymbol{\theta} \mapsto M(\boldsymbol{\theta};h)$,
    \begin{equation*}
        M(\eta \boldsymbol{\theta} + (1-\eta) \boldsymbol{\theta}_0;h) \leq \eta M(\boldsymbol{\theta};h) + (1-\eta) M(\boldsymbol{\theta}_0;h),
    \end{equation*}
    which rearranges as
    \begin{equation*}
        M(\boldsymbol{\theta};h) - M(\boldsymbol{\theta}_0;h)
        \geq \frac{1}{\eta} [M(\eta \boldsymbol{\theta} + (1-\eta) \boldsymbol{\theta}_0;h) - M(\boldsymbol{\theta}_0;h)] \geq 0.
    \end{equation*}
    As a consequence,
    \begin{equation*}
        \inf_{\boldsymbol{\theta} \in \Theta}M(\boldsymbol{\theta};h) = \inf_{\boldsymbol{\theta} \in \Theta^\epsilon_0}M(\boldsymbol{\theta};h) = \min_{\boldsymbol{\theta} \in \Theta^\epsilon_0}M(\boldsymbol{\theta};h),
    \end{equation*}
    where the last equality uses continuity of $\boldsymbol{\theta} \mapsto M(\boldsymbol{\theta};h)$ and compactness of $\Theta^\epsilon_0$.

    The above argument shows in particular that $\boldsymbol{\theta}(h) \in \Theta^\epsilon_0$ for $h>0$ near zero.

    If also Assumption \ref{Assumption: Asymptotic Distribution} holds, then, for $\boldsymbol{\theta}$ near $\boldsymbol{\theta}_0$, $h>0$ near zero, and for any $\mathbf{t} \in \mathbb{R}^k$,
    \begin{align*}
        \left| \frac{M(\boldsymbol{\theta}+\mathbf{t}\tau;h)-M(\boldsymbol{\theta};h)}{\tau} - \mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}) K_h(\mathbf{w}_1-\mathbf{w}_2)]' \mathbf{t} \right| &\leq \left| \mathbb{E}[r_\mathbf{t}(\boldsymbol{\theta},\tau) K_h(\mathbf{w}_1-\mathbf{w}_2)] \right| \\
        &\to 0 \qquad \text{as } \tau \downarrow 0,
    \end{align*}
    implying that for $h>0$ near zero, $\boldsymbol{\theta} \mapsto M(\boldsymbol{\theta};h)$ is (directionally) differentiable near $\boldsymbol{\theta}_0$, the directional derivative $\mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}) K_h(\mathbf{w}_1-\mathbf{w}_2)]' \mathbf{t}$ being zero when $\boldsymbol{\theta} = \boldsymbol{\theta}(h)$ because $\boldsymbol{\theta}(h)$ minimizes $M(\boldsymbol{\theta};h).$
\end{proof}

\subsection{Proof of Theorems \ref{Theorem: Asymptotic Distribution} and \ref{Theorem: Generalized Jackknifing}}

Theorem \ref{Theorem: Asymptotic Distribution} can be obtained from Theorem \ref{Theorem: Generalized Jackknifing} by setting $L=0$ in the latter, so it suffices to prove Theorem \ref{Theorem: Generalized Jackknifing}. To do so, for $l \in\{0,\dots,L/2\}$, let $\widehat{\boldsymbol{\theta}}_{n,l}=\widehat{\boldsymbol{\theta}}_n(h_{n,l}),\boldsymbol{\theta}_{n,l}=\boldsymbol{\theta}(h_{n,l}),$ and
\begin{equation*}
    \widehat{\mathbf{U}}_{n,l} = \binom{n}{2}^{-1} \sum_{i<j}  \mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j), \qquad \mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)= \mathbf{s}_{n,l}(\mathbf{z}_i,\mathbf{z}_j) - \mathbb{E}[\mathbf{s}_{n,l}(\mathbf{z}_1,\mathbf{z}_2)],
\end{equation*}
where
\begin{equation*}
    \mathbf{s}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)=\mathbf{s}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}_{n,l})K_{h_{n,l}}(\mathbf{w}_i-\mathbf{w}_j).
\end{equation*}
By Lemma \ref{Lemma: Existence and Convergence of theta(h)}, $\lim_{n \to \infty}\boldsymbol{\theta}_{n,l}=\boldsymbol{\theta}_0$ and $\mathbb{E}[\mathbf{s}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)]=0$ for large $n$.

Suppose that
\begin{equation}\label{eq: Hjort-Pollard stochastic expansion}
    \widehat{\boldsymbol{\theta}}_{n,l}-\boldsymbol{\theta}_{n,l} = -\boldsymbol{\Gamma}_0^{-1}\widehat{\mathbf{U}}_{n,l} + o_{\mathbb{P}}(\rho_n^{-1}) \quad \text{for} \quad l \in\{0,\dots,L/2\}.
\end{equation}
Then
\begin{equation*}
    \widetilde{\boldsymbol{\theta}}_n-\bar{\boldsymbol{\theta}}_n = -\boldsymbol{\Gamma}_0^{-1}\widetilde{\mathbf{U}}_n + o_{\mathbb{P}}(\rho_n^{-1}),
\end{equation*}
where
\begin{equation*}
    \widetilde{\mathbf{U}}_n = \sum_{l=0}^{L/2} \lambda_l(\mathbf{c})\widehat{\mathbf{U}}_{n,l}
\end{equation*}
satisfies
\begin{equation}\label{eq: U-statistic asymptotic normality}
    \left[ n^{-1}\boldsymbol{\Sigma}_0 + \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(\bar{K})\right]^{-1/2} \widetilde{\mathbf{U}}_n \leadsto \mathcal{N} (\mathbf{0}_{k \times 1},\mathbf{I}_k)
\end{equation}
because, letting $\sum_i$ denote $\sum_{i=1}^n$,
\begin{equation*}
    \widetilde{\mathbf{U}}_n = \widetilde{\mathbf{L}}_n + \widetilde{\mathbf{W}}_n,
\end{equation*}
where
\begin{equation*}
    \widetilde{\mathbf{L}}_n = n^{-1} \sum_i \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i), \qquad \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i)=\sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \mathbb{E}[\mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j)|\mathbf{z}_i] \qquad (j \neq i),
\end{equation*}
and
\begin{equation*}
    \widetilde{\mathbf{W}}_n = \binom{n}{2}^{-1} \sum_{i<j} \widetilde{\boldsymbol{\omega}}_n(\mathbf{z}_i,\mathbf{z}_j), \qquad \widetilde{\boldsymbol{\omega}}_n(\mathbf{z}_i,\mathbf{z}_j)=\sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \big[\mathbf{s}^{\mu}_{n,l}(\mathbf{z}_i,\mathbf{z}_j) - \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i) - \widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_j)\big],
\end{equation*}
satisfy
\begin{equation*}
    \left(\begin{array}{c} \sqrt{n}\widetilde{\mathbf{L}}_n \\  \sqrt{\binom{n}{2}h_n^d} \widetilde{\mathbf{W}}_n \end{array}\right) \leadsto \mathcal{N}\left( \left[\begin{array}{c} \mathbf{0}_{k \times 1} \\ \mathbf{0}_{k \times 1} \end{array}\right], \left[\begin{array}{cc} \boldsymbol{\Sigma}_0 &\mathbf{0}_{k \times k} \\ \mathbf{0}_{k \times k} &\boldsymbol{\Delta}_0(\bar{K})  \end{array}\right] \right),
\end{equation*}
as can be shown by means of the Cram{\'e}r-Wold device and the central limit theorem of \cite{Heyde-Brown_1970}, the latter being applicable because it follows from routine calculations that for every $\boldsymbol{\mu}_1,\boldsymbol{\mu}_2 \in \mathbb{R}^k$, we have
\begin{equation*}
    \varsigma_n^2 = \sum_i\mathbb{V}[g_{i,n}] = \boldsymbol{\mu}_1'\boldsymbol{\Sigma}_0\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2'\boldsymbol{\Delta}_0(\bar{K})\boldsymbol{\mu}_2 + o(1),
\end{equation*}
\begin{equation*}
    \sum_i \mathbb{E}[g_{i,n}^4] = o(1),
\end{equation*}
and
\begin{equation*}
    \mathbb{V}\left[\sum_i\varsigma_{i,n}^2-\varsigma_n^2\right] = o(1),\qquad \varsigma_{i,n}^2=\mathbb{V}[g_{i,n}|\mathbf{z}_1,\dots,\mathbf{z}_{i-1}],
\end{equation*}
where
\begin{equation*}
    g_{i,n} = g_{i,n}(\boldsymbol{\mu}) = \frac{2}{\sqrt{n}} \boldsymbol{\mu}_1'\widetilde{\boldsymbol{\ell}}_n(\mathbf{z}_i) + \sqrt{\binom{n}{2}^{-1}h_n^d}\sum_{j=1}^{i-1}\boldsymbol{\mu}_2'\widetilde{\boldsymbol{\omega}}_{n}(\mathbf{z}_i,\mathbf{z}_j).
\end{equation*}
The proof of Theorem \ref{Theorem: Generalized Jackknifing} can therefore be completed by verifying \eqref{eq: Hjort-Pollard stochastic expansion}.

To do so, we leverage convexity. For any $l \in \{0,\dots,L/2\}$ and any $\mathbf{t} \in \mathbb{R}^k$, it can be shown that
\begin{equation*}
    \lim_{\tau \downarrow 0,h \downarrow 0,\boldsymbol{\theta} \to \boldsymbol{\theta}_0} \mathbb{E} \big[ \mathbb{E}[ r_{\mathbf{t}}(\boldsymbol{\theta},\tau) K_h(\mathbf{w}_1-\mathbf{w}_2) | \mathbf{z}_1 ] ^2 \big] = 0
\end{equation*}
and
\begin{equation*}
    \lim_{\tau \downarrow 0,h \downarrow 0,\boldsymbol{\theta} \to \boldsymbol{\theta}_0} h^d \mathbb{E}[ r_{\mathbf{t}}(\boldsymbol{\theta},\tau)^2 K_h(\mathbf{w}_1-\mathbf{w}_2)^2 ] = 0,
\end{equation*}
and it therefore follows from a Hoeffding decomposition that
\begin{align*}
	&\rho_n^2\Big[\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \Big] \\
    &=  \rho_n^2 [M(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - M(\boldsymbol{\theta}_{n,l};h_{n,l})] + \mathbf{t}'\rho_n\widehat{\mathbf{U}}_{n,l} + o_\mathbb{P}(1).
\end{align*}
Moreover,
\begin{equation*}
    \rho_n^2 [M(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - M(\boldsymbol{\theta}_{n,l};h_{n,l}) ] \to \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t},
\end{equation*}
and proceeding as in the proof of \eqref{eq: U-statistic asymptotic normality} it can be shown that $\rho_n\widehat{\mathbf{U}}_{n,l} = O_\mathbb{P}(1)$. Because $\boldsymbol{\Gamma}_0$ is positive definite and because $\mathbf{t} \mapsto \widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l})$ is convex (almost surely), the corollary following \citet[Lemma 2]{Hjort-Pollard_1993} implies that \eqref{eq: Hjort-Pollard stochastic expansion} holds.

\subsection{Proof of Theorem \ref{Theorem: Bootstrapping}}

The proof of Theorem \ref{Theorem: Bootstrapping} is a natural bootstrap analog of the proof of Theorem \ref{Theorem: Generalized Jackknifing}.

For $l \in\{0,\dots,L/2\}$, let $\widehat{\boldsymbol{\theta}}^*_{n,l}=\widehat{\boldsymbol{\theta}}^*_n(h_{n,l})$ and
\begin{equation*}
    \widehat{\mathbf{U}}^*_{n,l} = \binom{n}{2}^{-1} \sum_{i<j}  \mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n}), \qquad \mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n})= \mathbf{s}_{n,l}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n}) - \mathbb{E}^*_n[\mathbf{s}_{n,l}(\mathbf{z}^*_{1,n},\mathbf{z}^*_{2,n})],
\end{equation*}
where $\mathbb{E}^*_n[\cdot]$ denotes $\mathbb{E}[\cdot|\mathbf{z}_1,\dots,\mathbf{z}_n]$.

It suffices to show that
\begin{equation}\label{eq: Hjort-Pollard stochastic expansion - bootstrap}
    \widehat{\boldsymbol{\theta}}^*_{n,l}-\boldsymbol{\theta}_{n,l} = -\boldsymbol{\Gamma}_0^{-1}\big(\widehat{\mathbf{U}}^*_{n,l} + \widehat{\mathbf{U}}_{n,l} \big) + o_{\mathbb{P}}(\rho_n^{-1}) \quad \text{for} \quad l \in\{0,\dots,L/2\}
\end{equation}
and that
\begin{equation}\label{eq: U-statistic asymptotic normality - bootstrap}
    \left[ n^{-1}\boldsymbol{\Sigma}_0 + 3 \binom{n}{2}^{-1}h_n^{-d}\boldsymbol{\Delta}_0(K)\right]^{-1/2} \widetilde{\mathbf{U}}^*_n \leadsto_\mathbb{P} \mathcal{N} (\mathbf{0}_{k \times 1},\mathbf{I}_k),
\end{equation}
where
\begin{equation*}
    \widetilde{\mathbf{U}}^*_n = \sum_{l=0}^{L/2} \lambda_l(\mathbf{c})\widehat{\mathbf{U}}^*_{n,l},
\end{equation*}
and where $\leadsto_\mathbb{P}$ denotes weak convergence in probability.

For every $\mathbf{t} \in \mathbb{R}^k$, using a Hoeffding decomposition and the fact that $m(\mathbf{z},\mathbf{z};\boldsymbol{\theta}_{n,l})=0$ and $\mathbf{s}(\mathbf{z},\mathbf{z};\boldsymbol{\theta}_{n,l})=\boldsymbol0$ for large $n$, we have
\begin{align*}
	&\rho_n^2\Big[\widehat{M}^*_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}^*_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \Big] \\
    &= \rho_n^2 (1+o(1)) [\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) ] + \mathbf{t}'\rho_n\widehat{\mathbf{U}}^*_{n,l} + o_\mathbb{P}(1),
\end{align*}
where it can be shown that $\rho_n\widehat{\mathbf{U}}^*_{n,l}=O_\mathbb{P}(1)$ and where it follows from the proof of \eqref{eq: Hjort-Pollard stochastic expansion} that
\begin{equation*}
    \rho_n^2 \Big[\widehat{M}_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \Big] = \frac{1}{2}\mathbf{t}'\boldsymbol{\Gamma}_0\mathbf{t} + \mathbf{t}'\rho_n\widehat{\mathbf{U}}_{n,l} + o_\mathbb{P}(1),
\end{equation*}
where $\rho_n\widehat{\mathbf{U}}_{n,l}=O_\mathbb{P}(1)$. In other words,
\begin{align*}
    &\rho_n^2\Big[\widehat{M}^*_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l}) - \widehat{M}^*_n(\boldsymbol{\theta}_{n,l};h_{n,l}) \Big] \\
    &= \frac{1}{2}\mathbf{t}'\boldsymbol{\Gamma}_0\mathbf{t} + \mathbf{t}'\rho_n\big(\widehat{\mathbf{U}}^*_{n,l} + \widehat{\mathbf{U}}_{n,l}\big) + o_\mathbb{P}(1) \qquad \text{for every } \mathbf{t} \in \mathbb{R}^k.
\end{align*}
Because $\boldsymbol{\Gamma}_0$ is positive definite and because $\mathbf{t} \mapsto \widehat{M}^*_n(\boldsymbol{\theta}_{n,l}+\mathbf{t} \rho_n^{-1};h_{n,l})$ is convex (almost surely),  the corollary following \citet[Lemma 2]{Hjort-Pollard_1993} implies that \eqref{eq: Hjort-Pollard stochastic expansion - bootstrap} holds.

To prove \eqref{eq: U-statistic asymptotic normality - bootstrap}, we begin by decomposing $\widetilde{\mathbf{U}}^*_n$ as
\begin{equation*}
    \widetilde{\mathbf{U}}^*_n = \widetilde{\mathbf{L}}^*_n + \widetilde{\mathbf{W}}^*_n,
\end{equation*}
where
\begin{equation*}
    \widetilde{\mathbf{L}}^*_n = n^{-1} \sum_i \widetilde{\boldsymbol{\ell}}^*_n(\mathbf{z}^*_{i,n}), \qquad \widetilde{\boldsymbol{\ell}}^*_n(\mathbf{z}^*_{i,n})=\sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \mathbb{E}^*_n[\mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n})|\mathbf{z}^*_{i,n}] \qquad (j \neq i),
\end{equation*}
and
\begin{equation*}
    \widetilde{\mathbf{W}}^*_n = \binom{n}{2}^{-1} \sum_{i<j} \widetilde{\boldsymbol{\omega}}^*_n(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n}), \qquad \widetilde{\boldsymbol{\omega}}^*_n(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n})=\sum_{l=0}^{L/2} \lambda_l(\mathbf{c}) \big[\mathbf{s}^{\mu,*}_{n,l}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n}) - \widetilde{\boldsymbol{\ell}}^*_n(\mathbf{z}^*_{i,n}) - \widetilde{\boldsymbol{\ell}}^*_n(\mathbf{z}^*_{j,n})\big].
\end{equation*}
Defining
\begin{equation*}
    \pi_n = \frac{\sqrt{nh_n^d}}{1+\sqrt{nh_n^d}},
\end{equation*}
routine calculations can be used to show that for every $\boldsymbol{\mu}_1,\boldsymbol{\mu}_2 \in \mathbb{R}^k$, we have
\begin{equation*}
    \widehat{\varsigma}_n^2 = \sum_i\mathbb{V}^*_n[g^*_{i,n}] = \boldsymbol{\mu}_1'\big[\pi_n^2\boldsymbol{\Sigma}_0 + 4(1-\pi_n)^2\boldsymbol{\Delta}_0(\bar{K})\big]\boldsymbol{\mu}_1 + \boldsymbol{\mu}_2'\boldsymbol{\Delta}_0(\bar{K})\boldsymbol{\mu}_2 + o_\mathbb{P}(1),
\end{equation*}
\begin{equation*}
    \sum_i \mathbb{E}^*_n[g^{*4}_{i,n}] = o_\mathbb{P}(1),
\end{equation*}
and
\begin{equation*}
    \mathbb{V}\left[\sum_i\widehat{\varsigma}_{i,n}^2-\widehat{\varsigma}_n^2\right] = o_\mathbb{P}(1), \qquad \widehat{\varsigma}_{i,n}^2=\mathbb{V}^*_n[g^*_{i,n}|\mathbf{z}^*_{1,n},\dots,\mathbf{z}^*_{i-1,n}],
\end{equation*}
where $\mathbb{V}^*_n[\cdot]$ denotes $\mathbb{V}[\cdot|\mathbf{z}_1,\dots,\mathbf{z}_n]$ and where
\begin{equation*}
    g^*_{i,n} = g^*_{i,n}(\boldsymbol{\mu}) = \frac{\pi_n}{\sqrt{n}} 2\boldsymbol{\mu}_1'\widetilde{\boldsymbol{\ell}}^*_n(\mathbf{z}^*_{i,n}) + \sqrt{\binom{n}{2}^{-1}h_n^d}\sum_{j=1}^{i-1}\boldsymbol{\mu}_2'\widetilde{\boldsymbol{\omega}}^*_{n}(\mathbf{z}^*_{i,n},\mathbf{z}^*_{j,n}).
\end{equation*}
The Cram{\'e}r-Wold device and the central limit theorem of \cite{Heyde-Brown_1970} therefore imply that if $\pi_n \to \pi_0 \in [0,1]$, then
\begin{equation*}
    \left(\begin{array}{c} \sqrt{n}\pi_n\widetilde{\mathbf{L}}^*_n \\  \sqrt{\binom{n}{2}h_n^d} \widetilde{\mathbf{W}}^*_n \end{array}\right) \leadsto_\mathbb{P} \mathcal{N}\left( \left[\begin{array}{c} \mathbf{0}_{k \times 1} \\ \mathbf{0}_{k \times 1} \end{array}\right], \left[\begin{array}{cc} \pi_0^2\boldsymbol{\Sigma}_0 + 4(1-\pi_0)^2\boldsymbol{\Delta}_0(\bar{K}) &\mathbf{0}_{k \times k} \\ \mathbf{0}_{k \times k} &\boldsymbol{\Delta}_0(\bar{K})  \end{array}\right] \right).
\end{equation*}
Whether or not $\pi_n$ is convergent, the result \eqref{eq: U-statistic asymptotic normality - bootstrap} can be obtained from the preceding display by arguing along subsequences (if necessary).

\subsection{Verifying Assumption \ref{Assumption: Bias of thetahat}}\label{Section: Verifying Assumption 3}

It follows from Lemma \ref{Lemma: Existence and Convergence of theta(h)} that if Assumption \ref{Assumption: Convergence of M} holds, then so does Assumption \ref{Assumption: Bias of thetahat} with $L=0$. This observation provides the base case for an induction argument. To describe the induction step, suppose that for some even $L \geq 0$, we have
\begin{equation*}
    \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l} + o(h^L) \qquad \text{ as } h\downarrow 0.
\end{equation*}
Suppose also that, for every $\mathbf{t} \in \mathbb{R}^k$ and some $\boldsymbol{\beta}_{L+2} \in \mathbb{R}^k$, we have
\begin{align*}
    &h^{-2(L+2)} \left[ M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l} + \mathbf{t} h^{L+2};h \right) - M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l}; h  \right)\right] \\
    &= \mathbf{t}' \boldsymbol{\beta}_{L+2} + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \qquad \text{as } h \downarrow 0.
\end{align*}
Then, the corollary following \citet[Lemma 2]{Hjort-Pollard_1993} implies that
\begin{align*}
    h^{-(L+2)} \left( \boldsymbol{\theta}(h) - \boldsymbol{\theta}_0 - \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l} \right) &= \operatorname*{arg\,min}_{\mathbf{t} \in \mathbb{R}^k} M \left(\boldsymbol{\theta}_0 + \sum_{l=1}^{L/2} \mathbf{b}_{2l} h^{2l} + \mathbf{t} h^{L+2};h \right) \\
    &= - \boldsymbol{\Gamma}_0^{-1} \boldsymbol{\beta}_{L+2} + o(1) \qquad \text{ as } h\downarrow 0;
\end{align*}
that is, defining $\mathbf{b}_{L+2} = - \boldsymbol{\Gamma}_0^{-1} \boldsymbol{\beta}_{L+2}$, we have
\begin{equation*}
    \boldsymbol{\theta}(h) -\boldsymbol{\theta}_0 = \sum_{l=1}^{(L+2)/2} \mathbf{b}_{2l} h^{2l} + o(h^{L+2}) \qquad \text{ as } h\downarrow 0.
\end{equation*}

To complete the proof of Proposition \ref{Proposition: Bias Expansion L=2}, it therefore suffices to note that (for every $\mathbf{t} \in \mathbb{R}^k$ and) under the assumptions of the proposition, we have
\begin{equation*}
    h^{-4} \left[ M (\boldsymbol{\theta}_0 + \mathbf{t} h^2;h ) - M (\boldsymbol{\theta}_0 ; h )\right]
    = \mathbf{t}' \boldsymbol{\beta}_2 + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \qquad \text{as } h \downarrow 0,
\end{equation*}
where
\begin{equation*}
    \boldsymbol{\beta}_2 = \frac{1}{2}\sum_{i=1}^d \int_\mathcal{W} \frac{\partial^2 \boldsymbol{\psi}(\mathbf{w},\mathbf{v})}{\partial v_i^2}\Big\vert_{\mathbf{v}=\mathbf{w}} f_{\mathbf{w}}(\mathbf{w})d\mathbf{w} \int_{\mathbb{R}^k} u_i^2 K(\mathbf{u})d\mathbf{u}.
\end{equation*}

Similarly, if in addition to the assumptions of Proposition \ref{Proposition: Bias Expansion L=2} it is assumed that for every $\mathbf{t} \in \mathbb{R}^k$ and for some $\boldsymbol{\beta}_4 \in \mathbb{R}^k$, we have
\begin{equation*}
    h^{-8} \left[ M (\boldsymbol{\theta}_0 + \mathbf{b}_2 h^2 + \mathbf{t} h^4;h ) - M (\boldsymbol{\theta}_0 + \mathbf{b}_2 h^2; h )\right] = \mathbf{t}' \boldsymbol{\beta}_4 + \frac{1}{2} \mathbf{t}' \boldsymbol{\Gamma}_0 \mathbf{t} + o(1) \qquad \text{as } h \downarrow 0,
\end{equation*}
then Assumption \ref{Assumption: Bias of thetahat} holds with $L=4$. One set of sufficient conditions for this to occur is that Assumptions \ref{Assumption: Convergence of M}-\ref{Assumption: Asymptotic Distribution} hold and that, for every $\mathbf{t} \in \mathbb{R}^k$, the following are satisfied (with probability one):

\begin{enumerate}[(i)]
    \item $\int_{\mathbb{R}^d} \|\mathbf{u}\|^4 K(\mathbf{u})d\mathbf{u}<\infty$.

    \item $\mathbf{v} \mapsto \boldsymbol{\psi}(\mathbf{w},\mathbf{v}) = \mathbb{E}[\mathbf{s}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{w}_1=\mathbf{w},\mathbf{w}_2=\mathbf{v}]f_{\mathbf{w}}(\mathbf{v})$ is four times continuously differentiable with $\mathbb{E}[\sup_{\mathbf{v}\in\mathcal{W}}\|\partial_{\mathbf{v}}^{\boldsymbol{\alpha}} \boldsymbol{\psi}(\mathbf{w},\mathbf{v})\|]<\infty$ for all $\boldsymbol{\alpha}\in\mathbb{Z}_+^d$ with $|\boldsymbol{\alpha}|\leq 4$.

    \item $f_\mathbf{w}$ is twice continuously differentiable and $\mathbf{v} \mapsto \mathbf{H}(\mathbf{w},\mathbf{v};\boldsymbol{\theta}_0,\mathbf{t})$ is twice continuously differentiable with $\mathbb{E}[\sup_{\mathbf{v}\in\mathcal{W}}\|\partial_{\mathbf{v}}^{\boldsymbol{\alpha}} \mathbf{H}(\mathbf{w},\mathbf{v};\boldsymbol{\theta}_0,\mathbf{t})f_\mathbf{w}(\mathbf{v})\|]<\infty$ for all $\boldsymbol{\alpha}\in\mathbb{Z}_+^d$ with $|\boldsymbol{\alpha}|\leq 2$.

    \item For some function $\dot\mathbf{H}(\mathbf{w},\mathbf{v};\boldsymbol{\theta},\mathbf{t})\in\mathbb{R}^{k\times k}$, $\mathbf{v} \mapsto \dot\mathbf{H}(\mathbf{w},\mathbf{v};\boldsymbol{\theta}_0,\mathbf{t})$ is continuous,
    \begin{equation*}
        \mathbb{E}\left[\sup_{\mathbf{v}\in\mathcal{W}}\|\dot \mathbf{H}(\mathbf{w},\mathbf{v};\boldsymbol{\theta}_0,\mathbf{t})f_\mathbf{w}(\mathbf{v})\|\right]<\infty,
    \end{equation*}
    \begin{equation*}
        \lim_{\tau\downarrow0,(\boldsymbol{\theta},\mathbf{u})\to(\boldsymbol{\theta}_0,\mathbf{0})} \left\| \frac{\mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta}+\tau\mathbf{t},\mathbf{t})-\mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta},\mathbf{t})}{\tau}-\dot\mathbf{H}(\mathbf{w},\mathbf{w}+\mathbf{u};\boldsymbol{\theta},\mathbf{t}) \right\|= 0,
    \end{equation*}
    and, for some $\delta>0$,
    \begin{equation*}
        \mathbb{E}\left[\sup_{ \tau\in (0,\delta),\|\boldsymbol{\theta}-\boldsymbol{\theta}_0\|<\delta, \mathbf{w}_2 \in\mathcal{W}} \left\| \frac{\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta}+\tau\mathbf{t},\mathbf{t})-\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t})}{\tau}-\dot\mathbf{H}(\mathbf{w}_1,\mathbf{w}_2;\boldsymbol{\theta},\mathbf{t}) \right\|\right]< \infty.
    \end{equation*}
\end{enumerate}

\section{Sufficient Conditions for Motivating Examples}\label{Section: Sufficient Conditions for Motivating Examples}

To demonstrate the plausibility of Assumptions \ref{Assumption: Convergence of M} and \ref{Assumption: Asymptotic Distribution}, we revisit the examples of Section \ref{Section: Motivating Examples}. In each example, Assumptions \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - convexity} holds and Assumption \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - density of w} is fairly primitive, so we focus on giving primitive sufficient conditions for Assumptions \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - well defined Mn}-\ref{Assumption: Convergence of M - well behaved M0} and \ref{Assumption: Asymptotic Distribution}.

\subsection{Partially Linear Regression Model}
We take $\mathbf{s}=\mathbf{s}_{\mathtt{PLR}}$ and $\mathbf{H}=\mathbf{H}_{\mathtt{PLR}}$, where
\begin{equation*}
    \mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{\partial}{\partial \boldsymbol{\theta}} m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\dot{\mathbf{x}}_{i,j} (\dot{y}_{i,j}-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta})
\end{equation*}
and
\begin{align*}
    \mathbf{H}_{\mathtt{PLR}}(\mathbf{w}_i,\mathbf{w}_j)
    &= \frac{\partial^2}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}'} \mathbb{E}[ m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j]
     = \frac{\partial}{\partial \boldsymbol{\theta}'} \mathbb{E}[ s_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] \\
    &= \mathbb{E}[\dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' | \mathbf{w}_i,\mathbf{w}_j],
\end{align*}
the latter depending on neither $\boldsymbol{\theta}$ nor $\mathbf{t}$ (because $m_{\mathtt{PLR}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})$ is quadratic in $\boldsymbol{\theta}$).

Under mild conditions, Assumptions \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - well defined Mn}-\ref{Assumption: Convergence of M - well behaved M0} and \ref{Assumption: Asymptotic Distribution} hold with
\begin{equation*}
    \boldsymbol{\xi}_0(\mathbf{z}) = 2\mathbb{E}[\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
\begin{equation*}
    \boldsymbol{\Xi}_0(\mathbf{w}) = \mathbb{E}[\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLR}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
and
\begin{equation*}
    \mathbf{G}_0(\mathbf{w}) = \mathbf{H}_{\mathtt{PLR}}(\mathbf{w},\mathbf{w}) f_{\mathbf{w}}(\mathbf{w}).
\end{equation*}
For instance, it suffices to set $b(\mathbf{z}) = (1+\|\mathbf{x}\|) (1+ |\varepsilon| + |\gamma_0(\mathbf{w})| +\|\mathbf{x}\|)$ and to assume that
\begin{enumerate}[(i)]
    \item The functions $\mathbf{w}\mapsto \gamma_0(\mathbf{w})$, $\mathbf{w}\mapsto \mathbb{E}[\mathbf{x}|\mathbf{w}]$, $\mathbf{w}\mapsto \mathbb{E}[\mathbf{x}\mathbf{x}'|\mathbf{w}]$, $\mathbf{w}\mapsto \mathbb{E}[\varepsilon^2|\mathbf{w}]$, $\mathbf{w}\mapsto \mathbb{E}[\mathbf{x}\varepsilon^2|\mathbf{w}]$, and $\mathbf{w}\mapsto \mathbb{E}[\mathbf{x}\mathbf{x}'\varepsilon^2|\mathbf{w}]$ are continuous on $\mathcal{W}$.

    \item $\mathbb{E}[(1+\|\mathbf{x}\|^4)\varepsilon^4] + \sup_{\mathbf{w}\in\mathcal{W}}\mathbb{E}[(1+\|\mathbf{x}\|^4)\varepsilon^4 |\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w})<\infty$ and
    \begin{equation*}
        \mathbb{E}[(1+\|\mathbf{x}\|^4)\gamma_0(\mathbf{w})^4 + \|\mathbf{x}\|^8] + \sup_{\mathbf{w}\in\mathcal{W}}\mathbb{E}[(1+\|\mathbf{x}\|^4)\gamma_0(\mathbf{w})^4 + \|\mathbf{x}\|^8|\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})<\infty.
    \end{equation*}

    \item With probability one, $\mathbb{V}[\mathbf{x}|\mathbf{w}]$ is positive definite and $\mathbb{V}[\varepsilon|\mathbf{x},\mathbf{w}]>0$.
\end{enumerate}

\subsection{Partially Linear Logit Model}
We take $\mathbf{s}=\mathbf{s}_{\mathtt{PLL}}$ and $\mathbf{H}=\mathbf{H}_{\mathtt{PLL}}$, where
\begin{equation*}
    \mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \frac{\partial}{\partial \boldsymbol{\theta}} m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = -\dot{\mathbf{x}}_{i,j} (y_i - \Lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}))\mathbbm{1}\{\dot{y}_{i,j}\neq0\}
\end{equation*}
and
\begin{align*}
    \mathbf{H}_{\mathtt{PLL}}(\mathbf{w}_i,\mathbf{w}_j;\boldsymbol{\theta}) &= \frac{\partial^2}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}'} \mathbb{E}[ m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] = \frac{\partial}{\partial \boldsymbol{\theta}'} \mathbb{E}[ s_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) | \mathbf{w}_i,\mathbf{w}_j] \\
    &= \mathbb{E}[ \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \lambda(\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta})\mathbbm{1}\{\dot{y}_{i,j}\neq0\}| \mathbf{w}_i,\mathbf{w}_j],
\end{align*}
where the latter does not depend on $\mathbf{t}$ (because $m_{\mathtt{PLL}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})$ is twice differentiable in $\boldsymbol{\theta}$).

Under mild conditions, Assumptions \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - well defined Mn}-\ref{Assumption: Convergence of M - well behaved M0} and \ref{Assumption: Asymptotic Distribution} hold with
\begin{equation*}
    \boldsymbol{\xi}_0(\mathbf{z}) = 2\mathbb{E}[\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
\begin{equation*}
    \boldsymbol{\Xi}_0(\mathbf{w}) = \mathbb{E}[\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
and
\begin{equation*}
    \mathbf{G}_0(\mathbf{w}) = \mathbf{H}_{\mathtt{PLL}}(\mathbf{w},\mathbf{w};\boldsymbol{\theta}_0)f_{\mathbf{w}}(\mathbf{w}).
\end{equation*}
For instance, it suffices to set $b(\mathbf{z}) = 1+\|\mathbf{x}\|$ and to assume that, for some $\delta>0$,
\begin{enumerate}[(i)]
    \item The function $\mathbf{w}\mapsto \gamma_0(\mathbf{w})$ is continuous on $\mathcal{W}$. Also, the conditional distribution of $\mathbf{x}$ given $\mathbf{w}$ admits a density $f_{\mathbf{x}|\mathbf{w}}$ with respect to some measure $\rho$ such that $\mathbf{w}\mapsto f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}|\mathbf{w})$ is continuous on $\mathcal{W}$ (with probability one) and
    \begin{equation*}
        \int_{\mathbb{R}^k} (1+ \|\mathbf{x}\|^2) \sup_{\|\mathbf{u}\|\leq\delta} f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}| \mathbf{w} + \mathbf{u})d\rho(\mathbf{x})<\infty \qquad \text{for every } \mathbf{w}\in\mathcal{W}.
    \end{equation*}

    \item $\mathbb{E}[\|\mathbf{x}\|^4] + \sup_{\mathbf{w}\in\mathcal{W}}\mathbb{E}[\|\mathbf{x}\|^4 |\mathbf{w}]f_{\mathbf{w}}(\mathbf{w})<\infty.$

    \item With probability one, $\mathbb{V}[\mathbf{x}|\mathbf{w}]$ is positive definite.
\end{enumerate}

\subsection{Partially Linear Tobit Model}
We take $\mathbf{s}=\mathbf{s}_{\mathtt{PLT}}$ and $\mathbf{H}=\mathbf{H}_{\mathtt{PLT}}$, where
\begin{equation*}
    \mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) = \dot{\mathbf{x}}_{i,j} \left(\mathbbm{1}\{y_j>\max(y_i-\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta},0)\} - \mathbbm{1}\{y_i>\max(y_j+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta},0)\}\right)
\end{equation*}
and
\begin{align*}
    \mathbf{H}_{\mathtt{PLT}}(\mathbf{w}_i,\mathbf{w}_j;\boldsymbol{\theta},\mathbf{t})
    &= \mathbb{E}[ \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \left(\mathbbm{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}>0\}+\mathbbm{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}=0,\dot{\mathbf{x}}_{i,j}'\mathbf{t}\geq0\}\right)\eta_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta})| \mathbf{w}_i,\mathbf{w}_j]  \\
    &\quad + \mathbb{E}[ \dot{\mathbf{x}}_{i,j}\dot{\mathbf{x}}_{i,j}' \left(\mathbbm{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}<0\}+\mathbbm{1}\{\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}=0,\dot{\mathbf{x}}_{i,j}'\mathbf{t}<0\}\right)\eta_{\mathtt{PLT}}(\mathbf{z}_j,\mathbf{z}_i;\boldsymbol{\theta})| \mathbf{w}_i,\mathbf{w}_j],
\end{align*}
with
\begin{align*}
    \eta_{\mathtt{PLT}}(\mathbf{z}_i,\mathbf{z}_j;\boldsymbol{\theta}) &= 2\int_0^{\infty} f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_i'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_i)+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}|\mathbf{w}_i) f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_j'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_j)|\mathbf{w}_j)d\varepsilon  \\
    &\quad + f_{\varepsilon|\mathbf{w}}(-\mathbf{x}_i'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_i)+\dot{\mathbf{x}}_{i,j}'\boldsymbol{\theta}|\mathbf{w}_i) \int_{-\infty}^0  f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}_j'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w}_j)|\mathbf{w}_j)d\varepsilon.
\end{align*}

Under mild conditions, Assumptions \ref{Assumption: Convergence of M}\ref{Assumption: Convergence of M - well defined Mn}-\ref{Assumption: Convergence of M - well behaved M0} and \ref{Assumption: Asymptotic Distribution} hold with
\begin{equation*}
    \boldsymbol{\xi}_0(\mathbf{z}) = 2\mathbb{E}[\mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)|\mathbf{z}_1=\mathbf{z},\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
\begin{equation*}
    \boldsymbol{\Xi}_0(\mathbf{w}) = \mathbb{E}[\mathbf{s}_{\mathtt{PLT}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)\mathbf{s}_{\mathtt{PLL}}(\mathbf{z}_1,\mathbf{z}_2;\boldsymbol{\theta}_0)'|\mathbf{w}_1=\mathbf{w}_2=\mathbf{w} ] f_{\mathbf{w}}(\mathbf{w}),
\end{equation*}
and
\begin{equation*}
    \mathbf{G}_0(\mathbf{w}) = \mathbf{H}_{\mathtt{PLT}}(\mathbf{w},\mathbf{w};\boldsymbol{\theta}_0,\boldsymbol{\theta}_0)f_{\mathbf{w}}(\mathbf{w}).
\end{equation*}

For instance, it suffices to set $b(\mathbf{z}) = 1+\|\mathbf{x}\|$ and to assume that, for some $\delta>0$,
\begin{enumerate}[(i)]
    \item The function $\mathbf{w}\mapsto \gamma_0(\mathbf{w})$ is continuous on $\mathcal{W}$. Also, the conditional distribution of $\mathbf{x}$ given $\mathbf{w}$ admits a a density $f_{\mathbf{x}|\mathbf{w}}$ with respect to some measure $\rho$ such that $\mathbf{w}\mapsto f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}|\mathbf{w})$ is continuous on $\mathcal{W}$ (with probability one) and
    \begin{equation*}
        \int_{\mathbb{R}^k} (1+ \|\mathbf{x}\|^2) \sup_{\|\mathbf{u}\|\leq\delta} f_{\mathbf{x}|\mathbf{w}}(\mathbf{x}| \mathbf{w} + \mathbf{u})d\rho(\mathbf{x})<\infty \qquad \text{for every } \mathbf{w}\in\mathcal{W}.
    \end{equation*}
    In addition, the function $(\varepsilon,\mathbf{w})\mapsto f_{\varepsilon|\mathbf{w}}(\varepsilon|\mathbf{w})$ is continuous and bounded and the function
    \begin{equation*}
        (\mathbf{x},\mathbf{w})\mapsto \int_\mathbb{R} \sup_{|u|+\|\mathbf{u}\|\leq\delta}f_{\varepsilon|\mathbf{w}}(\varepsilon-\mathbf{x}'\boldsymbol{\theta}_0-\gamma_0(\mathbf{w})+u|\mathbf{w}+\mathbf{u})d\varepsilon
    \end{equation*}
    is bounded.

    \item $\mathbb{E}[\|\mathbf{x}\|^4] + \sup_{\mathbf{w}\in\mathcal{W}}(1+\mathbb{E}[\|\mathbf{x}\|^4 |\mathbf{w}])f_{\mathbf{w}}(\mathbf{w})<\infty.$

    \item With probability one, $\mathbb{V}[\mathbf{x}|\mathbf{w}]$ is positive definite.
\end{enumerate}

\section{Conclusion}\label{Section: Conclusion}

This paper has developed bandwidth robust distribution theory and bootstrap-based inference procedures for a broad class of convex pairwise difference estimators. Our theoretical work is based on small bandwidth asymptotics and carefully leverages convexity. The theory is illustrated by means of three prominent examples. In addition to expanding the scope of small bandwidth asymptotics, our results lay the groundwork for several promising avenues of future research. First, our methods could be generalized to develop bandwidth selection based on higher-order stochastic expansions. Second, they could be expanded to allow for pairwise difference estimators based on generated regressors, a class of estimators that sometimes arises in the context of control function and related econometric methods. Third, when the objective function is smooth, plug-in variance estimation could be developed as an alternative to bootstrap inference. Finally, our current results do not cover settings where the objective function is sufficiently non-smooth to result in non-Gaussian distributional approximations. We plan to investigate these research directions in upcoming work.

\bibliography{CJN_2026_ET--bib}
\bibliographystyle{ecta}