Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.
Debiased Machine Learning U-Statistics
bibunit[ectabib]
\begin{abstract}
We propose a method to debias estimators based on U-statistics with Machine Learning (ML) first-steps. Standard plug-in estimators often suffer from regularization and model-selection biases, producing invalid inferences. We show that Debiased Machine Learning (DML) estimators can be constructed within a U-statistics framework to correct these biases while preserving desirable statistical properties. The approach delivers simple, robust estimators with provable asymptotic normality and good finite-sample performance. We apply our method to three problems: inference on Inequality of Opportunity (IOp) using the Gini coefficient of ML-predicted incomes given circumstances, inference on predictive accuracy via the Area Under the Curve (AUC), and inference on linear models with ML-based sample-selection corrections. Using European survey data, we present the first debiased estimates of income IOp. In our empirical application, commonly employed ML-based plug-in estimators systematically underestimate IOp, while our debiased estimators are robust across ML methods.\newline\newlineJEL\ Classification: C13; C14;
C21; D31; D63 \newlineKeywords: local robustness, orthogonal moments,
machine learning, U-statistics, Inequality of Opportunity, AUC, pairwise difference estimators. \newlineR package:
\url{https://joelters.github.io/home/code/}
\end{abstract}
\section{Introduction}
A fast-growing literature in economics and other sciences employs Machine
Learning (ML) methods to estimate high-dimensional first-steps in two-step settings.
The availability of these ML tools brings new
opportunities, but also introduces new challenges. One key challenge is that
regularization and model selection---crucial for variance control in ML---can
induce substantial bias in downstream inferences.
This concern has motivated a recent and expanding literature on Debiased
Machine Learning (DML); see chernozhukov2018double and references therein. DML methods address
second-step bias by constructing Neyman-orthogonal moments---identifying moments whose derivatives with respect to the first-step limit vanish. Combined with cross-fitting---a form of sample splitting---these methods have proven effective in settings where the target parameter is
a population average of a function of the observed data and the first-step ML
limit. Importantly, the existing DML literature focuses almost exclusively
on parameters of this form---such as average treatment effects or other
policy-relevant expectations.
However, many parameters of interest in economics are more complex functionals
of the data than simple averages. This paper takes the next natural step in
complexity by considering target parameters that involve double averages over
the data. These double averages are referred to
as V-statistics, which are closely related to U-statistics---a well-known class of estimators originally motivated by
their favorable bias properties (see halmos1946).\footnote{The bias that originally motivated U-statistics
is a higher-order (second-order) bias, distinct from the first-order regularization and
selection biases that are the main focus of this paper.} We show that
U-statistics involving functions of the data and of an ML first-step estimator provide a flexible and powerful framework for inference in a variety
of relevant settings in economics. These include the measurement of inequality
of opportunity, the evaluation of ML predictions in
classification problems, and inference for flexible sample selection models, among many others.
To explain the main issues, let $W_{i}$ denote a data observation for unit
$i$, and let $\gamma$ denote a high-dimensional or ML population first-step.
For example, in our leading application to inequality, $\gamma(x)$ is the optimal
mean-square prediction (i.e., conditional mean) of income given an individual $i$'s vector of circumstances $x$ (e.g.
parental education, sex, race, etc.). The DML literature has considered target
parameters such as
\[
\theta(\gamma)=\mathbb{E}[m(W_{i},\gamma)],
\]
for a moment function $m$ of the data and the first-step $\gamma$. In contrast, this
paper considers parameters of the form
\[
\theta(\gamma)=\mathbb{E}[g(W_{i},W_{j},\gamma)],
\]
for a symmetric moment function $g$ involving pairs of
independent observations.\footnote{Note that our setting nests the previous setting as a special case since we can always
write $\theta(\gamma)=\mathbb{E}[m(W_{i},\gamma)]=\mathbb{E}[g(W_{i}
,W_{j},\gamma)],$ for $g(W_{i},W_{j},\gamma)=\left[ m(W_{i},\gamma
)+m(W_{j},\gamma)\right] /2,$ provided $W_{i}\ $and $W_{j}$ have the same
distribution.}
In both cases, the main econometrics challenge is that the asymptotic
behaviour of two-step estimators $\hat{\theta}(\hat{\gamma})$ depends
critically on the term
\begin{equation}
\frac{\partial\theta(\gamma_{0})}{\partial\gamma}\left[ \hat{\gamma}
-\gamma_{0}\right] ,
\end{equation}
where $\hat{\gamma}$ is the first-step estimator and $\gamma_{0}$ its
probabilistic limit.\footnote{Heuristically, $\hat{\theta}(\hat{\gamma})\approx\hat{\theta}(\gamma_{0})+\frac{\partial\theta(\gamma_{0})}{\partial\gamma}\left[ \hat{\gamma}
-\gamma_{0}\right]$, see Section 6 in newey1994large.} The precise sense of the derivative in ((ref)) is defined in Section (ref).
A key feature of an ML estimator $\hat{\gamma}$ is that
its estimation error $\hat{\gamma}-\gamma_{0},$ and, in particular, its bias
$\mathbb{E}[\hat{\gamma}-\gamma_{0}],$ is often very large. This is an
unavoidable consequence of the high-dimensionality of the statistical problem.
This means that, unless ((ref)) is zero, the large bias of the first-step will propagate into the second-step, invalidating inference and leading to potentially misleading empirical conclusions.
Moment functions $m$ yielding a zero-derivative property in ((ref)) are called
Neyman-Orthogonal moment functions.\footnote{We use the terms Neyman-orthogonal, orthogonal, locally robust or debiased interchangeably. } chernozhukov2022locally have
shown how to construct such locally robust moment functions for parameters
$\theta(\gamma)$ that solve Generalized Method of Moments (GMM) restrictions. Our paper shows how to achieve this goal when
$\theta(\gamma)$ is a U-statistic functional, i.e., of the form $\theta
(\gamma)=\mathbb{E}[g(W_{i},W_{j},\gamma)].$ Moment functions $g$ with the
zero-derivative property in ((ref)) will be also called debiased,
and their corresponding sample analog (or cross-fitted version) estimators $\hat{\theta}(\hat{\gamma})$ based on an ML first-step
$\hat{\gamma}$ will be referred to as DML U-estimators.
In this paper, we (i) construct debiased moment functions $g$ for
U-statistics; (ii) analyze the asymptotic properties of the resulting DML U-estimators; and (iii) illustrate the framework with applications to
some relevant problems in economics, including inequality of opportunity, ML-based
prediction evaluation metrics for classification, and inferences with ML-sample selection corrections via pairwise difference estimators (see powell1987semiparametric).
Our main application, which motivated the methodology, is to Inequality of Opportunity (IOp)---the share of inequality explained by circumstances that are beyond an individual's
control, such as parental income, education, sex, or social origin.
Equality of opportunity has emerged as an important ideal of distributive justice, both in public policy and academic discourse (cf. gaer1993
, fleurbaey1995equal and roemer1998equality).
The literature on economic inequality distinguishes between two main
approaches to measuring inequality of opportunity: the intergenerational
mobility (IGM) approach and the IOp approach. The IGM approach---particularly
the intergenerational elasticity (IGE) of income---has become more popular in
economics and related fields, largely due to its statistical simplicity and
well-understood properties. However, it rests on a narrow conception of
opportunity, typically focusing only on one single circumstance, namely, parental income. In contrast, the IOp
literature embraces a richer set of circumstances and is grounded in a
strong normative framework. Yet, its empirical methods are often ad
hoc (see, e.g., the identification of types for circumstances), and their statistical properties remain poorly understood and underdeveloped.
This paper provides a step-change to address this shortcoming. It contributes to the IOp literature by introducing, to our knowledge,
the first rigorous econometric approach to IOp measurement. Our approach builds on recent developments in DML to offer a general and statistically sound methodology for
estimation and inference on IOp. To illustrate its practical relevance and enhance widespread adoption, we apply the method to European cross-country data and provide an open-source R package, ineqopp, implementing
our estimators and inference procedures.
Specifically, we apply our methodology to one of the most widely used measures of IOp--the Gini coefficient of predicted income from a set of circumstances.\footnote{The Online Appendix (ref) shows how to debias other popular inequality measures used in the literature, such as the Atkinson index and the Generalized Entropy. These results are new and of independent interest.} Earlier studies estimated these
predictions using simple log-linear regressions or grouped data with a small number of types, see the reviews
in roemer2016equality, ramos2016approaches or
ferreira2016individual, while recent work increasingly relies on ML
methods (e.g., brunori2019inequality, brunori2021evolution, brunori2021roots, carranza2022, hufe2022lower, salas2022inheritances, bernardo2025model). However, ML-based plug-in estimators are subject to large
regularization and model selection biases that distort the estimated share
of inequality explained by circumstances, invalidating the resulting inferences. Additionally, commonly employed standard errors for the Gini do not account for the first-step and are generally inconsistent in our two-step setting.
We propose a novel debiased estimator for Gini-based IOp that is valid across a wide
range of ML methods. It is extremely simple to compute, requires no
additional nuisance parameters, and enables valid inference. Simulations show
that the bias in conventional plug-in estimators can be severe. Our debiased estimator substantially reduces this bias. The new IOp estimator can be seen as a generalization to ML of the Lorenz regression method proposed independently by heuchenne2022inference.
We report the first debiased estimates of IOp across 29 European countries using the
2019 wave of the EU-SILC survey. We find that the most commonly employed ML-based plug-in estimates systematically underestimate IOp, both with simulated and real data. The difference between
debiased and plug-in estimates reaches up to 10 percentage points in some cases. In simulations, our debiased estimator behaves similarly to a well-specified parametric model, while in the application, a common parametric approach yields considerably different results.
Furthermore, we estimate a decomposition of the difference between debiased and plug-in estimates into regularization and overfitting terms. Regularization tends to smooth the fitted-values distribution and biases the plug-in downward, whereas overfitting tends to spread it and biases the plug-in upward. We provide estimates of both sources for all countries. This analysis complements the discussion of over/underfitting biases of the plug-in estimator in brunori2019upward, where only plug-in estimators are discussed. Our estimates show far greater stability across different ML algorithms, in contrast to plug-in estimates, which vary widely. A key advantage of our method is that it enables formal inference and hypothesis testing—such as comparing IOp across countries, time, or circumstance sets—offering rigorous tools for policy evaluation.
Our second application is to the evaluation performance of ML classifiers. One of the most popular metrics is the area under the receiver operating characteristic curve (AUC), see bradley1997use. The AUC represents the likelihood of correctly ranking a random observation from a positive class higher than a random observation from the negative class. This measure is systematically employed across all disciplines to evaluate ML methods, and economics is no exception, see, e.g.,
einav2018predictive, kleinberg2018human, herrera2020political, baron2021banking, jorda2021bank, aydin2022consumption, bazzi2022promise, broner2022fiscal, chan2022selection, dellavigna2024bottlenecks, ludwig2024machine, muller2024credit, to mention but a few.
In applied economic settings, researchers typically report the AUC alongside its standard errors.\footnote{For example, all the aforementioned papers published in the Review of Economic Studies report standard errors for the AUC measures.} The most widely used inference method is the influence function approach of delong1988comparing, often combined with cross-validation. ledell2015computationally derive inference results for a cross-validated AUC parameter, which differs from the population AUC considered in this paper. rotnitzky2006doubly and zhou2025doubly examine AUC inference in the presence of nuisance parameters, but they treat the classifier as fixed and focus on nuisance components arising from identification in missing-data settings. By contrast, our analysis addresses the technical challenge posed by the randomness of the classifier itself, which enters through indicator functions. Despite the widespread use of the AUC measure, to the best of our knowledge the theoretical validity of its asymptotic inference remains an open question, and no clear guidance on best practices is currently available.
Interestingly, we show that the AUC can be written as a known functional of our debiased Gini estimand and the proportion of observations in the positive class. Hence, by our previous results, this representation of the AUC implies that it is locally robust to the ML first-step. Thus, our paper provides the first theoretical justification for the inferences reported in the aforementioned references, and supports the view of the AUC as a robust measure of predictive performance (cf. bradley1997use).
Finally, our third application builds on the literature of pairwise difference estimators to construct a DML U-estimator for the linear model with an ML-sample selection correction. Pairwise estimators were introduced in the seminal work by powell1987semiparametric and have been applied to several relevant econometric models, see the survey in Powell1994 and the more recent contributions in, e.g., kyriazidou1997estimation,blundell2004endogeneity,honore2005identification,chen2010symmetry,hong2010pairwise,aradillas2012pairwise,jochmans2013pairwise,chen2020censored.
In its most general versions, the asymptotic theory for these estimators has been derived for first-steps in a relatively small class (a Donsker class), as in sherman1994u. Unfortunately, Donsker conditions do not hold in high-dimensional settings, see chernozhukov2018double for discussion. To avoid
Donsker conditions, we use cross-fitting as in chernozhukov2018double,
but adapted to a U-statistic setting.
We illustrate with a sample selection model, as in ahn1993semiparametric, but where the selection equation is estimated by an ML method (e.g. Random Forest). A prominent empirical example in econometrics is mroz1987sensitivity, where the researcher specifies a very flexible selection equation with various powers, transformations, and interactions of the original conditioning variables to model female labor force participation. The DML method of chernozhukov2018double is not valid here because it would not account for the estimation of the selection equation, see escanciano2023automatic. Employing a pairwise difference approach leads to a simple plug-in estimator, as suggested in ahn1993semiparametric for a kernel first-step. However, classical pairwise estimators are highly biased when the selection equation is flexibly modeled (e.g., with an ML method). We use our methodology to provide the first debiased pairwise difference estimator that is locally robust to an ML-sample selection correction.
There are many other potential applications of our methodology. These include extensions of our results to: dyadic data (see, e.g., chiang2021dyadic); inequality-aware policy learning (see terschuur2025locally); tests for the equality of conditional distributions (see chen2025biased); semiparametric distance-based estimators (e.g. dominguez2004consistent); Causal U-statistics (see mao2018causal); U-statistics with missing data (e.g. the Gini with missing income data); inferences on measures of fairness in algorithmic decision-making, (see, e.g., rychener2022metrizing) and evaluation metrics for treatment rules such as the QINI (see imai2013estimating, sun2021treatment, yadlowsky2025evaluating); among many others. Some of these applications are further discussed in the Online Appendix (ref).
The rest of the paper is organized as follows. Section (ref) introduces the methodology, using IOp as
a running example. Section (ref) explains how to
construct orthogonal moments in a general U-statistic setting. Section (ref) derives the
asymptotic theory in the general setting and for the IOp example. Section
(ref) shows the performance of the debiased IOp estimator in
Monte Carlo simulations. Section (ref) studies IOp in Europe.
Section (ref) considers additional applications of the methodology. Section (ref) contains two appendices. Appendix (ref) gives a Practitioners' guide to IOp, while Appendix (ref) provides the proofs of the general results. Finally, an Online Appendix gathers the proofs for the examples, and further discussion on simulations and applications.
\section{Methodology}
This section introduces our methodology in a setting slightly more general than the one discussed in the Introduction, using IOp as a running example.
We have independent and identically distributed (i.i.d.) data $W_{i}
=(Y_{i},X_{i})$, for $i=1,...,n$, distributed with unknown distribution
$F_{0}$, an unknown first-step function $\gamma_{0}$ and a finite-dimensional
parameter of interest $\theta_0$ in some set $\Theta\subseteq\mathbb{R}^{p}$. We
can form pairs $(W_{i},W_{j}),$ where $W_{j}$ is an independent copy of
$W_{i}$, with realizations $(w_{i},w_{j})$. Let $\mathbb{E}[\cdot]$ denote the
expectation under $F_{0}$, let $\gamma(F)$ be the probabilistic limit (plim) of a
first-step estimator $\hat{\gamma}$ under $F$, like in
newey1994asymptotic, and $\gamma_{0}=\gamma(F_{0})$. We assume there is
a vector $g(w_{i},w_{j},\gamma,\theta)$ of $p$ known identifying moment
functions such that
\[
\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta)]=0\iff\theta=\theta_{0}\in\Theta.
\]
Without loss of generality, we can assume that $g$ is a symmetric function in
$w_{i}$ and $w_{j}$.\footnote{Otherwise, replace $g(w_{i},w_{j},\gamma
_{0},\theta)$ by $g^{*}(w_{i},w_{j},\gamma_{0}
,\theta)=(1/2)[g(w_{i},w_{j},\gamma_{0},\theta)+g(w_{j},w_{i},\gamma
_{0},\theta)]$.} By independence between $W_{i}$
and $W_{j}$,
\begin{equation}
\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta)]=\int\int g(w_{i},w_{j},\gamma
_{0},\theta)F_{0}(dw_{i})F_{0}(dw_{j}).
\end{equation}
A key feature of this moment is that it is a quadratic functional of $F_{0}$.
This property is what differentiates our analysis from the standard debiasing
literature (see, e.g., chernozhukov2018double), and it leads naturally
to U-statistics (see, e.g., lee2019u for a comprehensive treatment of U-statistics). If we replace $F_{0}$ by the empirical
distribution of $\{W_{i}\}_{i=1}^{n}$ in ((ref)), we obtain a so-called V-statistic for each $\theta\in\Theta$ (or V-process if indexed by $\theta\in\Theta$)
\[
V_{n}g(\cdot,\gamma_{0},\theta)=\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}
^{n}g(W_{i},W_{j},\gamma_{0},\theta),
\]
and its unbiased U-statistic counterpart
\[
U_{n}g(\cdot,\gamma_{0},\theta)=\frac{2}{n(n-1)}\sum_{i<j}g(W_{i},W_{j}
,\gamma_{0},\theta),
\]
where, henceforth, we use the short notation $\sum_{i<j}=\sum_{i=1}^{n-1}
\sum_{j=i+1}^{n}$.\footnote{Note that
$\mathbb{E}[V_{n}g]=\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta)]+b_{n},$ where
$b_{n}=(\mathbb{E}[g(W_{i},W_{i},\gamma_{0},\theta)]-\mathbb{E}[g(W_{i}
,W_{j},\gamma_{0},\theta)])/n$ is a bias term. In contrast, $\mathbb{E}[U_{n}g]=\mathbb{E}[g(W_{i},W_{j},\gamma_{0}
,\theta)].$ Nevertheless, our debiased estimators address a typically larger
bias term than $b_{n}$, related to the use of ML first-steps, so the difference between $V$ and $U$ statistics is
not central to our results.} Plugging in a first-step estimator $\hat{\gamma}$
for $\gamma_{0}$ leads to the plug-in estimator $\hat{\theta}^{P}$ solving
$U_{n}g(\cdot,\hat{\gamma},\hat{\theta}^{P})=0.$ We illustrate with our
running example.
\textsc{Example }(\textbf{Inequality of Opportunity, IOp})\textsc{:}
The leading measure of IOp is given by the Gini coefficient of income's predictions from a set of circumstances,
which can be expressed in the population as\footnote{The primitive definition of the Gini is twice the area between the Lorenz curve and the 45-degree line. This definition works for both discrete and continuous variables. The U-statistic expression we use in this paper is equivalent to the primitive Lorenz-based expression even with discrete variables.}
\begin{equation}
\theta_{0}=\frac{\mathbb{E}[|\gamma_{0}(X_{i})-\gamma_{0}(X_{j})|]}
{\mathbb{E}[\gamma_{0}(X_{i})+\gamma_{0}(X_{j})]},
\end{equation}
where $\gamma_0(x)=\mathbb{E}[Y_{i}|X_{i}=x]$, $x\in\mathcal{X}$, $\mathcal{X}$ is the support of $X_{i}$, $Y_{i}$ is income, and $X_{i}$ is a vector of circumstances individual $i$ did not choose, such
as parental wealth/income, parental education, sex, color of the skin or
social origin, among others. Rearranging ((ref)), we obtain an
identifying moment function
\begin{equation}
g(w_{i},w_{j},\gamma,\theta)=(\gamma(x_{i})+\gamma(x_{j}))\theta-|\gamma
(x_{i})-\gamma(x_{j})|.
\end{equation}
The estimator $\hat{\theta}^{P}$ solving
$U_{n}g(\cdot,\hat{\gamma},\hat{\theta}^{P})=0$ is the plug-in sample Gini
coefficient of fitted values
\[
\hat{\theta}^{P}=\frac{\sum_{i<j}|\hat{\gamma}(X_{i})-\hat{\gamma}(X_{j}
)|}{\sum_{i<j}\left( \hat{\gamma}(X_{i})+\hat{\gamma}(X_{j})\right) },
\]
for a first-step $\hat{\gamma}$ such as the Conditional Inference Forests (CIF),
which is a popular ML method in the IOp literature, see
brunori2021roots, brunori2019inequality or
brunori2021evolution. We show that
$\hat{\theta}^{P}$ is highly biased when ML first-steps, such as the CIF, are used. We develop
inference methods for $\theta_{0}$ and related parameters based on debiased
estimators, which improve upon plug-in estimators.
As an illustration, see Figure (ref), where we simulate the
plug-in and the debiased estimators of IOp which we develop in Section (ref). The Data
Generating Process (DGP) is the same as in the empirical simulation of Section (ref). The
fitted values are estimated with CIF. The histograms approximate the
distributions of the centered estimators and the solid curves are normal p.d.f.s
with the same variance as the estimators but with no bias (i.e. centered at
zero). Even for as large sample sizes as $n=3000$ the plug-in estimator has a
large bias compared to the DML U-estimator proposed in this paper.
$\blacksquare$
\begin{figure}[h]
\caption{Comparison
of plug-in and debiased estimator.}
\end{figure}
More broadly, plug-in estimators are generally biased by model
selection and/or regularization in the first-step. Intuitively, ML algorithms
trade-off bias and variance. Hence, some bias might be allowed in the first-step if it aids prediction. This first-step bias can propagate into the second-step, thereby distorting the estimation of the target parameter.
To explain the source of the bias in the context of our running example, define the sign function $sgn(x)=1(x>0)-1(x<0)$ (equals
$-1,0$ and $1$ for negative, zero and positive numbers, respectively) and
consider an expansion of the numerator of the Gini of fitted values (the
numerator is the most difficult term in an expansion of the bias)
\begin{align*}
\sum_{i<j}|\hat{\gamma}(X_{i})-\hat{\gamma}(X_{j})| & =\sum_{i<j}\left(
\hat{\gamma}(X_{i})-\hat{\gamma}(X_{j})\right) sgn(\hat{\gamma}(X_{i}
)-\hat{\gamma}(X_{j}))\\
& =\sum_{i<j}\left( Y_{i}-Y_{j}\right) sgn(\hat{\gamma}(X_{i})-\hat{\gamma
}(X_{j}))\\
& +\sum_{i<j}\left( \hat{\gamma}(X_{i})-\gamma_{0}(X_{i})-\hat{\gamma}
(X_{j})+\gamma_{0}(X_{j})\right) sgn(\hat{\gamma}(X_{i})-\hat{\gamma}
(X_{j}))\\
& -\sum_{i<j}\left( Y_{i}-\gamma_{0}(X_{i})-Y_{j}+\gamma_{0}(X_{j})\right)
sgn(\hat{\gamma}(X_{i})-\hat{\gamma}(X_{j})).
\end{align*}
The first term in the last equality leads to a
(nonlinear) bias that is bounded and of second order under some conditions. The second
term is the leading first-order bias term. This term is the main source of the
bias problem. With ML estimators the bias $\mathbb{E}[\hat{\gamma}
(X_{i})-\gamma_{0}(X_{i})]$ is so large that makes $\sqrt{n}\left(
\hat{\theta}^{P}-\theta_{0}\right) $ to diverge (see Section (ref) for an illustration of this problem). The third term is another second-order bias term.\footnote{Dealing with and controlling for the nonlinear biases is much more difficult in this application than in other applications commonly studied in the DML literature, see Proposition (ref).}
In the general setting, we can explain the source of the bias problem through a Taylor
expansion argument as in the Introduction. Differentiating $\mathbb{E}[g(W_{i},W_{j},\gamma,\theta(\gamma))]=0$ with respect to $\gamma$
at $\gamma_{0}$ yields
\begin{equation}
\frac{\partial\theta(\gamma_{0})}{\partial\gamma}=-\left[ \frac{\partial
\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta_0)]}{\partial\theta}\right]
^{-1} \frac{\partial}{\partial\gamma}
\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta_0)].
\end{equation}
We show in Proposition (ref) below that this
derivative is not zero for the IOp example, which leads to the bias problem observed
in Figure (ref). The first factor in ((ref)) is a standard derivative, but the second factor is a functional derivative whose formal definition requires some care. We formalize these heuristics and the functional derivative with respect to
$\gamma$ in the next section. The next section also shows how to construct identifying moments with a zero derivative in our U-statistics setting.
\subsection{Construction of debiased moments}
To reduce the first-step bias, we aim to find an adjustment term $\phi$ to the original moment function $g$ such that $\psi=g+\phi$ satisfies $\mathbb{E}[\psi(W_{i},W_{j},\gamma_{0},\theta_0)]=0$ and
\begin{equation}
\frac{\partial}{\partial\gamma}\mathbb{E}[\psi(W_{i},W_{j},\gamma_{0},\theta_0)]=0,
\end{equation}
so the derivative in ((ref)) with $g$ replaced by $\psi$ is zero. This section formalizes this statement and shows how such $\phi$ can be found. To that end, let $F_0$ again
be the cumulative distribution function (cdf) for $W_{i}$ which is unrestricted
except for regularity conditions such as the existence of $\gamma(F_0)$. For
example, for $\gamma(F_0)=\mathbb{E}_{F_0}[Y|X]$ we require $\mathbb{E}
_{F_0}[|Y|]<\infty$, where $\mathbb{E}_{F}$ denotes expectation under $F.$
Henceforth, $d/d\tau$ denotes the derivative from the right with respect to
$\tau$ evaluated at $\tau=0$. Let $F_{\tau}=F_{0}+\tau(H-F_{0})$ be a local deviation from $F_{0}$ along some
alternative distribution $H$, with $\tau\in\lbrack0,1]$. The alternative
distribution $H$ is chosen such that $\gamma_{\tau}=\gamma(F_{\tau})$ exists
for $\tau$ small enough. We aim to find a function $\phi$ such that for all $\theta$ and
all $H$ in $F_{\tau}=F_{0}+\tau(H-F_{0}),$
\begin{equation}
\mathbb{E}_{F_{\tau}}[\phi(W_{i},W_{j},\gamma(F_{\tau}),\alpha(F_{\tau
}),\theta)]=0\text{ for all }\tau\in\lbrack0,\bar{\tau}),\text\bar{\tau}>0,
\end{equation}
and
\begin{equation}
\frac{d}{d\tau}\mathbb{E}[g(W_{i},W_{j},\gamma(F_{\tau}),\theta)]=\int\int
\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)K_{H}\left( dw_{i}
,dw_{j}\right) ,
\end{equation}
with $K_{H}\left( dw_{i},dw_{j}\right) =F_{0}(dw_{i})H(dw_{j})+H(dw_{i}
)F_{0}(dw_{j}).$ Here, $\alpha_{0}$ is an additional nuisance parameter which
satisfies ((ref)). We may assume without loss of generality that $\phi$ is
symmetric in $w_{i},w_{j}$.\footnote{As with $g,$ setting $\phi^{\ast}
(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)=(1/2)[\phi(w_{i},w_{j},\gamma
_{0},\alpha_{0},\theta)+\phi(w_{j},w_{i},\gamma_{0},\alpha_{0},\theta)]$, we
have $\int\int\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)K_{H}
(dw_{i},dw_{j})=\int\int\phi^{\ast}(w_{i},w_{j},\gamma_{0},\alpha_{0}
,\theta)K_{H}(dw_{i},dw_{j})$ and $\phi^{\ast}$ is symmetric.} Equation
((ref)) is a zero mean condition on $\phi.$ The pathwise derivative in
((ref)) formalizes the heuristic derivative $\partial
\mathbb{E}[g(W_{i},W_{j},\gamma_{0},\theta_0)]/\partial\gamma$ of the previous section.
The function $\phi$ can be found by solving the functional equation ((ref)) and it
characterizes the local effect of the first-step $\gamma(F)$ on the functional $\mu(F)=\mathbb{E}[g(W_{i},W_{j},\gamma(F),\theta)]$ as $F$
varies away from $F_{0}$ in any direction $H$. In Section
(ref) we give new results characterizing $\phi$ for
first-steps satisfying orthogonality conditions, including high-dimensional regressions.
\textsc{Example }(\textbf{IOp, cont.})\textsc{:} In this example $\hat{\gamma}$ is a
nonparametric consistent estimator of the conditional mean of income $Y$ given
circumstances $X$, i.e. $\gamma(F)=\mathbb{E}_{F}[Y|X]$. A first technical
challenge we face for computing $\phi$ in this example is the lack of
differentiability of the absolute value. Despite this lack of
differentiability, we are able to compute $\phi$ from an application of a general result in the next section (cf. Lemma (ref)) under the following regularity condition. Again, let $\gamma_{\tau}=\gamma(F_{\tau})$ for $F_{\tau}=F_{0}+\tau(H-F_{0})$.
\begin{assumption}
For each $H$, there exists an $\epsilon>0$ such that $|\gamma_{\tau}(X)-\gamma_{0}(X)|\leq \tau C_{H}(X)$ a.s., for all $\tau\in[0,\epsilon]$. Moreover, $\mathbb{E}[C_{H}(X)]<\infty.$
\end{assumption}
Assumption (ref) is satisfied if $H$ is absolutely continuous with respect to $F_{0}$, and other mild moment conditions, such as $\mathbb{E}
_{H}[|Y|]<\infty$ and $\mathbb{E}
_{H}[|\gamma_{0}(X)|]<\infty$, hold.
\begin{proposition}
Under Assumption (ref), the following function
$\phi$ satisfies ((ref)) and ((ref)) for the Gini of fitted values
\begin{equation}
\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)=\alpha_{01}(y_{i}+y_{j}-\gamma
_{0}(x_{i})-\gamma_{0}(x_{j}))+\alpha_{02}(x_{i},x_{j})(y_{i}-y_{j}-\gamma
_{0}(x_{i})+\gamma_{0}(x_{j})),
\end{equation}
where $\alpha_{0}=(\alpha_{01},\alpha_{02})$, $\alpha_{01}=\theta$ and $\alpha_{02}(x_{i},x_{j})=-sgn(\gamma_{0}(x_{i})-\gamma_{0}(x_{j}))$.
\end{proposition}
The proof of Proposition (ref) is given in Online Appendix (ref). The fact that $\phi\neq0,$ or more precisely that ((ref))
is not zero, explains the bias observed in Figure (ref). Indeed, arguing as in ((ref)), for the Gini of fitted values,
\begin{equation*}
\frac{d\theta(\gamma_{\tau})}{d\tau}=-(\mathbb{E}[2Y])^{-1} \int\int
\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)K_{H}\left( dw_{i}
,dw_{j}\right),
\end{equation*}
where $\phi$ is given in ((ref)). $\blacksquare$
To motivate our definition of the adjustment term $\phi$, note that equation
((ref)) is a zero mean condition and it implies by the chain rule
\begin{align}
\frac{d}{d\tau}\mathbb{E}[\phi(W_{i},W_{j},\gamma(F_{\tau}),\alpha(F_{\tau
}),\theta)] & =-\frac{d}{d\tau}\mathbb{E}_{F_{\tau}}[\phi(W_{i},W_{j}
,\gamma(F_{0}),\alpha(F_{0}),\theta)]\nonumber\\
& =-\int\int\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)K_{H}(dw_{i}
,dw_{j}),
\end{align}
for all $\theta$ and $H$. Equation ((ref)) shows that the effect the first-steps have on $\phi$ \textquotedblleft cancels out\textquotedblright\ with
the effect they have on the original identifying moment in ((ref)). Thus, letting $\psi(w_{i},w_{j},\gamma,\alpha,\theta)=g(w_{i},w_{j}
,\gamma,\theta)+\phi(w_{i},w_{j},\gamma,\alpha,\theta)$, the total first-step effect on $\psi$ is zero, i.e., we have local robustness (by ((ref)) and ((ref)))
\[
\frac{d}{d\tau}\mathbb{E}[\psi(W_{i},W_{j},\gamma(F_{\tau}),\alpha(F_{\tau
}),\theta)]=0.
\]
This local robustness motivates estimators based on the orthogonal moments
$\psi$. This construction led to the expression of the new IOp estimator,
which would not have been possible from the results in the debiased
literature.
\textsc{Example }(\textbf{IOp, cont.}
)\textsc{:} Adding $\phi$ of Proposition (ref) to the
identifying moment $g$ from ((ref)) leads after rearranging to
\begin{align*}
\psi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta) & =(\gamma_0(X_{i}
)+\gamma_0(X_{j}))\theta-|\gamma_{0}(X_{i})-\gamma_{0}(X_{j})|+\theta
(Y_{i}+Y_{j}-\gamma_{0}(X_{i})-\gamma_{0}(X_{j}))\\
& -sgn(\gamma_{0}(X_{i})-\gamma
_{0}(X_{j}))(Y_{i}-Y_{j})+sgn(\gamma_{0}(X_{i})-\gamma
_{0}(X_{j}))(\gamma
_{0}(X_{i})-\gamma_{0}(X_{j}))\\
& =\theta(Y_{i}+Y_{j})-sgn(\gamma_{0}(X_{i})-\gamma
_{0}(X_{j}))(Y_{i}-Y_{j}),
\end{align*}
where we use that for $a,b\in\mathbb{R}$, $sgn(a-b)(a-b)=|a-b|$. In this example, we get a surprisingly simple orthogonal moment function by employing our U-statistic setting. $\blacksquare
$
\subsection{DML U-estimators}
Estimation of nuisance parameters and moment conditions with the same
observations can induce an \textquotedblleft overfitting\textquotedblright
\ bias. Also, machine learning first-steps usually do not satisfy Donsker
conditions (see chernozhukov2018double for discussion). We use cross-fitting to overcome these issues (see
bickel1982adaptive, schick1986asymptotically,
klaassen1987consistent, chernozhukov2018double and
chernozhukov2022locally), but adapted to the U-statistics setting. We
partition the set of pairs $\{(i,j)\in\mathcal\{1,...,n\}^{2}:i<j\}$ into $T$ triangles and $T(T-1)/2$ rectangles. Then, we split each of these rectangles into four smaller rectangles. Hence, we get $L = T + 2T(T-1) = T(2T - 1)$ blocks $\mathcal{I}=\{I_{1},...,I_L\}$. If $n$ is divisible by $2\cdot T$, the blocks all leave $n - n/T$ observations for constructing nuisance parameter estimators. Otherwise, blocks leave approximately (up to a few observations) $n - \lfloor n/T \rfloor$ observations for estimating nuisance parameters. For an illustration, see Figure (ref). For an algorithm for general $(n,T)$ see Online Appendix (ref). The estimators $\hat{\gamma}_{l}$ and $\hat{\alpha}_{l}$ are
constructed using observations not present in the pairs in $I_{l}$, for $l=1,...,L$.
\begin{figure}[h]
\caption{\scriptsize Cross-fitting partition: 1) split $i = 1,...,n$ and $j = 1,...,n$ into $T=2$ equally sized folds, 2) Denote the blocks adjacent to the diagonal $I_1$ and $I_2$, 3) split the remaining rectangle into $4$ equally sized blocks. Balanced splits are always achieved when $n$ is divisible by $2\cdot T$. For other cases see Online Appendix (ref).}
\end{figure}
The debiased sample moment is $\hat{\psi}(\theta)=\hat{g}(\theta)+\hat{\phi}(\theta)$, where\footnote{We could also implement the debiased sample moment as $\hat{\psi}(\theta)=\hat{g}(\theta)+\tilde{\phi}$, where $\tilde{\phi}$ uses a preliminary cross-fitted estimator of $\theta_0$, with no significant changes in the theory.}
\[
\hat{g}(\theta)
=\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}} g(W_{i},W_{j},\hat{\gamma}_{l},\theta),
\quad
\hat{\phi}(\theta)
=\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}} \phi(W_{i},W_{j},\hat{\gamma}_{l},\hat{\alpha}_{l},\theta).
\]
We call an
estimator solving $\hat{\psi}(\hat{\theta})=0$ a DML U-estimator.\footnote{chernozhukov2018double recommend repeatedly estimating with different random cross-fitting splits and taking the median to improve robustness to the partition; this is also possible here but may be computationally expensive.} We give
sufficient conditions in Section (ref) so that the asymptotic variance
of $\sqrt{n}(\hat{\theta}-\theta_{0})$ is $V=B^{-1}\Sigma B^{\prime-1},$ where
\[
B=\frac{\partial\mathbb{E}[\psi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta
)]}{\partial\theta},\quad\Sigma=4\mathbb{V}ar\biggl(\mathbb{E}[\psi
(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta_{0})|W_{i}]\biggr).
\]
The above asymptotic variance can be estimated by $\hat{V}=\hat{B}^{-1}
\hat{\Sigma}\hat{B}^{\prime-1},$ where
\begin{align*}
\hat{B} & =\binom{n}{2}^{-1}\sum_{i<j}\frac{\partial}{\partial\theta}
\psi(W_{i},W_{j},\hat{\gamma},\hat{\alpha},\hat{\theta}),\\
\hat{\Sigma} & =\frac{4}{n(n-1)^{2}}\sum_{i=1}^{n}\biggl[\sum_{j\neq i}
\psi(W_{i},W_{j},\hat{\gamma},\hat{\alpha},\hat{\theta})\biggr]\biggl[\sum
_{j\neq i}\psi(W_{i},W_{j},\hat{\gamma},\hat{\alpha},\hat{\theta
})\biggr]^{\prime}.
\end{align*}
For variance estimation, we can use the whole sample for $\hat{\gamma}$ and $\hat{\alpha}$, since we only need consistency. It is also possible to cross-fit, but it is computationally inconvenient.
\textsc{Example }(\textbf{IOp, cont.})\textsc{:
}Let $\hat{\gamma}_{l}(x)$ denote cross-fitted ML predictions of
income given circumstances $x$. Solving the debiased orthogonal sample moment we get
\[
\hat{\theta}=\frac{\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}sgn(\hat{\gamma}_{l}(X_{i})-\hat{\gamma
}_{l}(X_{j}))(Y_{i}-Y_{j})}{\sum_{i<j}(Y_{i}+Y_{j})}.
\]
The debiased IOp estimator $\hat{\theta}$ resembles the standard Gini coefficient for income,
but rather than weighting by the sign of the differences in income, it weights by the sign of the difference in predictions.
This means that whenever two individuals have the same fitted values their
difference in incomes cannot be attributed to inequality of opportunity. In this example,
\[
h(w) \equiv \mathbb{E}[\psi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta)|W_{i} = w]=\mathbb{E}[\theta
(y+Y_j)- sgn(\gamma_0(x) - \gamma_0(X_j))(y - Y_j)].
\]
The asymptotic variance of the DML U-estimator is given by
$V=\Sigma/B^2,$ where
\[
\Sigma=4\mathbb{V}ar(h(W_i)) \text{ and }B=2\mathbb{E}[Y_{i}],
\]
and it can be consistently estimated by $\hat{V}=\hat{\Sigma}/4\bar{Y}^{2}$, where $\bar{Y}=n^{-1}\sum_{i=1}^{n}Y_{i}$ and
\begin{equation}
\hat{\Sigma} = \frac{4}{n} \sum_{i=1}^n \left( \frac{1}{n-1} \sum_{j\neq i}
\hat{\theta}
(Y_i+Y_j)- sgn(\hat{\gamma}(X_i) - \hat{\gamma}(X_j))(Y_i - Y_j) \right)^2.
\end{equation}
$\blacksquare$
The construction of the orthogonal moment requires deriving an adjustment term $\phi$ for the
given identifying moment $g$ satisfying ((ref)) and ((ref)). We have derived in Proposition
(ref) the adjustment term $\phi$ corresponding to the identifying moment $g$ in
the running IOp example. The next section explains how $\phi$ was obtained for this and other examples, providing a general theory for characterizing $\phi$ for a general class of two-step U-statistics.\footnote{IOp practitioners can safely skip the next section, and jump directly to the Monte Carlo, Section (ref).}
\section{General Construction of Adjustment Terms}
The expression for $\phi$ depends on the original identifying moment function $g$ and how the first-step limit $\gamma(F)$ is identified. Henceforth, $L_{2}(X_{i})$ is the set of measurable functions of $X_{i}$ with finite second moments. We define $L_{2}(X_{i},X_{j})$ analogously, and sometimes we drop the reference to random variables and use simply $L_{2}\ $(the meaning will be clear from the context). We consider first-step functions $\gamma(F)$ in
a linear set $\Gamma\in L_{2}(X_{i})$, such that
\begin{equation}
\mathbb{E}_{F}[\nu(X_{i})\left( Y_{i}-\gamma(F)\right) ]=0\text{ for all
}\nu\in\Gamma.
\end{equation}
This setting is quite general and fits well with the applications we are interested in. It covers linear, nonparametric, and high-dimensional regressions, among others.
Specifically, if $\Gamma=\{\beta^{\prime}X_{i}:\beta\in\mathbb{R}^{d}\}$, with $d$ the dimension of $X_{i}$, then $\gamma(F)(X_{i})=\beta_{F}^{\prime}X_{i}$ is a linear regression fit. If $\Gamma=L_{2}(X_{i})$,
then $\gamma(F)(X_{i})=\mathbb{E}_{F}[Y_{i}|X_{i}]$ is a nonparametric regression fit. When $\Gamma$ is the mean-square limit
of linear combinations $\sum_{k=1}^{K}\beta _{0k}b_{k}(X_{i})$, for some $K$, a sequence of real numbers $(\beta
_{0k})_{k=1}^{\infty }$ and a dictionary $(b_{k})_{k=1}^{\infty }$ of
functions of $X_{i}$ with finite variance, then $\gamma(F)(X_{i})$ is a high-dimensional regression fit. For other examples
included in this setting, see newey1994asymptotic. We require the following assumption on $\Gamma$.
\begin{assumption}
$\Gamma\subseteq L_{2}(X_{i})$ is a closed linear set that contains constant functions.
\end{assumption}
For simplicity of exposition, we take $\gamma(F)$ to be real-valued, though the extension of
our results to multiple first-steps follows straightforwardly from the chain rule.
The two key steps to derive $\phi$ are (i) a linearization step of the original
identifying moment function $g$ with respect to $\gamma$; and (ii) a projection step to provide an integral representation as in ((ref)), see newey1994asymptotic for such a representation in the standard GMM setting. We separate the analysis into these two steps.
\subsection{Linearization}
For this linearization step, we assume there exist functions
$\delta_{m}\in L_{2}$, constants $(c_{1m},c_{2m})$, $m=1,...,M,$ and an integer
$M\geq1$, such that
\begin{equation}
\frac{d}{d\tau}\mathbb{E}[g(W_{i},W_{j},\gamma(F_{\tau}),\theta)]=\sum
_{m=1}^{M}\frac{d}{d\tau}\mathbb{E}[\delta_{m}(X_{i},X_{j},\gamma_{0}
,\theta)(c_{1m}\gamma_{\tau}(X_{i})+c_{2m}\gamma_{\tau}(X_{j}))],
\end{equation}
where $\gamma_{\tau}(x)=\gamma(F_{\tau})(x)$ is the first-step limit corresponding to $F_{\tau}=F_{0}
+\tau(H-F_{0})$. For simplicity of notation, we drop the possible dependence
of $\delta_{m}$, $c_{1m}$ and $c_{2m}$ on $\theta$ and $\gamma_{0}.$ Henceforth, we also use
the short notation $\delta_{ij,m}(\gamma)\equiv\delta_{m}(X_{i},X_{j}
,\gamma)\ $ and $\delta_{ij,m}\equiv\delta_{ij,m}(\gamma_{0})$.
\textsc{Example }(\textbf{IOp, cont.}
)\textsc{:} In the proof of Proposition (ref) we show that
\begin{align}
\frac{d}{d\tau}\mathbb{E}[g(W_{i},W_{j},\gamma(F_{\tau}),\theta)]
&= \frac{d}{d\tau}\mathbb{E}\left[\theta\left(\gamma_{\tau}(X_{i}) + \gamma_{\tau}(X_{j})\right)\right] \nonumber \\
&\quad - \frac{d}{d\tau}\mathbb{E}\left[sgn(\gamma_{0}(X_{i}) - \gamma_{0}(X_{j}))\left(\gamma_{\tau}(X_{i}) - \gamma_{\tau}(X_{j})\right)\right],
\end{align}
so that ((ref)) is satisfied with $M=2$, $\delta_{ij,1}=\theta$, $(c_{11},c_{21})=(1,1),$ $\delta_{ij,2}=-sgn(\gamma_{0}(X_{i})-\gamma_{0}(X_{j}))$, and
$(c_{12},c_{22})=(1,-1)$. $\blacksquare$
\subsection{Projection}
After the linearization, we modify the projection step of newey1994asymptotic to fit our current U-statistic framework. The novel idea is to view $c_{1m}\gamma_{\tau}(X_{i})+c_{2m}\gamma_{\tau}(X_{j})$ as an orthogonal projection of $c_{1m}Y_{i}+c_{2m}Y_{j}$ on a suitable set $\mathcal{S}$ of $L_2(X_{i},X_{j})$ to be defined below.
We consider two cases, depending on the modeling assumptions. In case (i), the researcher follows a fully nonparametric approach and sets $\Gamma=L_{2}(X_{i})$ and $\mathcal{S}=L_2(X_{i},X_{j})$. In the case (ii), $\Gamma\subseteq L_{2}(X_{i})$ is any other set that satisfies Assumption (ref), and the researcher defines $\mathcal{S}=\Gamma+\Gamma \in L_2(X_{i},X_{j})$. This second case is useful when, for example, nonparametric estimation is not feasible (e.g., when $X$ is high-dimensional). Each case defines a different nuisance parameter $\alpha_0 \in \mathcal{S}$. Henceforth $\Pi_{V}(.)$ is the orthogonal projection operator onto $V$, for a closed linear set $V$.
The goal is to express the derivative of the identifying moment function as the integral representation in ((ref)). After the linearization, we show below that $\alpha_{0m}=\Pi_{\mathcal{S}}(\delta_m)$ satisfies
\begin{equation*}
\mathbb{E}[\delta_{m}(X_{i},X_{j},\gamma_{0}
,\theta)(c_{1m}\gamma_{\tau}(X_{i})+c_{2m}\gamma_{\tau}(X_{j}))]=\mathbb{E}[\alpha_{0m}(X_{i},X_{j},\gamma_{0}
,\theta)(c_{1m}\gamma_{\tau}(X_{i})+c_{2m}\gamma_{\tau}(X_{j}))].
\end{equation*}
and
\begin{equation}
\mathbb{E}_{F_{\tau}}[\alpha_{0m}(X_{i},X_{j},\gamma_{0}
,\theta)(c_{1m}Y_{i}+c_{2m}Y_{j}-c_{1m}\gamma_{\tau}(X_{i})-c_{2m}\gamma_{\tau}(X_{j}))]=0.
\end{equation}
We can then apply the chain rule to ((ref)) to express the derivative of interest as an integral representation in ((ref)). The additional nuisance parameters $\alpha_0=(\alpha_{01},...,\alpha_{0M})$ are characterized in the following lemma, which gives the expression for the adjustment term $\phi$. Henceforth, for simplicity of notation we drop the dependence of $\alpha_0$ on $\gamma_{0}$ and $\theta$.
\begin{lemma}
Suppose ((ref)), ((ref)) and Assumption
(ref) hold. Then,
\begin{equation}
\phi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta)=\sum_{m=1}^{M}\alpha_{0m}
(X_{i},X_{j})\left( c_{1m}Y_{i}+c_{2m}Y_{j}-c_{1m}\gamma_{0}(X_{i}
)-c_{2m}\gamma_{0}(X_{j})\right),
\end{equation}
satisfies ((ref))-((ref)) with $\alpha_{0}=(\alpha_{01},...,\alpha_{0M})^{\prime}$, where, for $m=1,...,M$,
\begin{description}
• $\alpha_{0m}=\delta_{m}(\cdot,\gamma_{0})$ for $\Gamma=L_{2}(X_{i})$,
• and for any other $\Gamma,$
\begin{equation}
\alpha_{0m}(X_{i},X_{j})=\Pi_{\Gamma}\mathbb{E}[\delta_{ij,m}|X_{i}]+\Pi_{\Gamma}\mathbb{E}[\delta_{ij,m}|X_{j}
]-\mathbb{E}[\delta_{ij,m}].
\end{equation}
\end{description}
\end{lemma}
\textsc{Example }(\textbf{IOp, cont.}
)\textsc{:} Since we are modelling the prediction fully nonparametricaly (case (i)), then $\alpha_{01}(X_{i},X_{j})=\delta_{ij,1}=\theta$ and $\alpha_{02}(X_{i},X_{j})=\delta
_{ij,2}=-sgn(\gamma_{0}(X_{i})-\gamma_{0}(X_{j}))$. $\blacksquare$
The estimator $\hat{\alpha}_{l}$ for $\alpha_{0}(X_{i},X_{j})$ depends on $\Gamma$. When $\delta$ is known up to the first-step $\gamma_{0},$ as in Lemma (ref)(i), $\alpha_{0}=\delta
(\cdot,\gamma_{0})$ is also known up to $\gamma_{0}$ and we can set
$\hat{\alpha}_{l}=\delta(\cdot,\hat{\gamma}_{l}).$ For other cases, we can
obtain $\hat{\alpha}_{l}$ from any ML method that estimates
the projection of $\delta(\cdot,\hat{\gamma}_{l})$ onto $\Gamma$ (such as Lasso, neural nets, sieves, etc.). For example, using Lemma (ref)(ii) and letting $n_{l}=\sum_{j:(i,j)\notin I_{l}}1$, we can estimate the conditional expectations in ((ref)) by
\[
\tilde{\alpha}_{1l}(X_{i})=\frac{1}{n_{l}}\sum_{j\neq i}\delta_{ij}
(\hat{\gamma}_{l})\text{ and }\tilde{\alpha}_{2l}(X_{j})=\frac{1}{n_{l}}
\sum_{i\neq j}\delta_{ij}(\hat{\gamma}_{l}),
\]
respectively, and then the orthogonal projection of $\tilde{\alpha}
_{rl}(X_{i})$ onto $\Gamma$ in ((ref)), for $r=1,2$, with any ML estimating orthogonal projections and
with observations not in $I_{l}$.
If $\alpha_{0m}=0$ for all $m=1,...,M$, then the original identifying moment is locally robust and no adjustment term is necessary (we set $\phi=0$). We illustrate with debiased IOp.
\textsc{Example }(\textbf{Debiased IOp}
)\textsc{:} If $g(W_{i},W_{j}, \gamma_0,\theta_0)=(Y_i+Y_j)\theta_0-sgn(\gamma_{0}(X_{i})-\gamma_{0}(X_{j}))(Y_i-Y_j)$, then the corresponding adjustment term is $\phi=0$. $\blacksquare$
\section{Asymptotic Theory}
\subsection{The General Case}
The aim of using a debiased moment function and cross-fitting is to be able to perform valid inferences. First, we will show that under some relatively mild conditions,
\begin{equation}
\sqrt{n}\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\psi(W_{i}
,W_{j},\hat{\gamma}_{l},\hat{\alpha}_{l},\theta_{0})=\sqrt{n}\binom{n}{2}
^{-1}\sum_{i<j}\psi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta_{0})+o_{p}(1).
\end{equation}
This means that the first-order asymptotic properties of the debiased
$\hat{\theta}$ are not affected by the estimation of first-steps. In ((ref)), the estimator $\hat{\alpha}_{l}$ is evaluated at $\theta_0$ when $\alpha_{0}(X_{i},X_{j})$ depends on $\theta_0$. This simplifies the verification of the required convergence conditions below.
Let $\left\vert \cdot\right\vert $ and $||\cdot||_q$ be
the Euclidean and $L_{q}$ norms, respectively. Specifically, for a measurable function $m$ of $W_i$, $||m||_q =(\mathbb{E}[|m(W_i)|^q])^{1/q}$, for $q\geq1$. For simplicity of notation, we drop the subindex for $q=2$, and simply write $||m||$. Also define the supremum norm, $||m||_\infty =\sup_{x\in\mathcal{X}}|m(x)|$. The terms $\rightarrow_{p}$ and $\rightarrow_{d}$ denote convergence in probability and distribution, respectively.
\begin{assumption}
$\mathbb{E}[|\psi(W_{i},W_{j},\gamma_{0},\alpha
_{0},\theta_{0})|^{2}]<\infty$ and for each $l=1,...,L,$
\begin{enumerate}
[(i)]
• $\int\int|g(w_{i},w_{j},\hat{\gamma}_{l},\theta_{0})-g(w_{i}
,w_{j},\gamma_{0},\theta_{0})|^{2}F_{0}(dw_{i})F_{0}(dw_{j})\rightarrow_{p}0$;
• $\int\int|\phi(w_{i},w_{j},\hat{\gamma}_{l},\alpha_{0},\theta_{0}
)-\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta_{0})|^{2}F_{0}(dw_{i}
)F_{0}(dw_{j})\rightarrow_{p}0$;
• $\int\int|\phi(w_{i},w_{j},\gamma_{0},\hat{\alpha}_{l},\theta_{0}
)-\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta_{0})|^{2}F_{0}(dw_{i}
)F_{0}(dw_{j})\rightarrow_{p}0$.
\end{enumerate}
\end{assumption}
These are mild mean-square consistency conditions for $\hat{\gamma}_{l}$ and
$\hat{\alpha}_{l}$ separately. A linearization argument like equation (3.2)
often implies that the left hand side of Assumption (ref)(i)-(ii)
are bounded above by a constant times $\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert,$ so $L_{2}$ consistency suffices. Assumption
(ref)(iii) typically follows from $L_{2}$ consistency of
$\hat{\alpha}_{l}$.
There is a large literature checking $L_{q}$-convergence rates for different
machine learners under low level sparsity or smoothness conditions on the
nuisance parameters. The traditional non-parametric literature gives rates for
kernel regression (e.g. tsybakov2009introduction) and sieves/series (e.g.
chen2007large). For $L_{1}$-penalty estimators such as Lasso see, e.g.,
belloni2011 and belloni2013least. Also for low level conditions
for shrinkage and kernel estimators see Appendix B in
sasaki2021estimation. Rates for $L_{2}$-boosting in low dimensions are
found in zhang2005boosting, and more recently in
kueck2023estimation with high-dimensional data. For results on versions
of random forests see wager2015adaptive and athey2019generalized
. Finally, for a modern setting of deep neural networks with rectified linear (ReLU) activation function see, e.g., schmidt2020nonparametric and farrell2021deep. For recent results on $L_q$ rates for machine learners see el_hanchi2023optimal and peng2025adversarial.
Define now the following interaction term
\[
\hat{\xi}_{l}(w_{i},w_{j})=\phi(w_{i},w_{j},\hat{\gamma}_{l},\hat{\alpha}
_{l},\theta_{0})-\phi(w_{i},w_{j},\gamma_{0},\hat{\alpha}_{l},\theta_{0}
)-\phi(w_{i},w_{j},\hat{\gamma}_{l},\alpha_{0},\theta_{0})+\phi(w_{i}
,w_{j},\gamma_{0},\alpha_{0},\theta_{0}).
\]
\begin{assumption}
For each $l=1,...,L$, either i)
\[
\sqrt{n}\int\int\hat{\xi}_{l}(w_{i},w_{j})F_{0}(dw_{i})F_{0}(dw_{j}
)\rightarrow_{p}0,\int\int|\hat{\xi}_{l}(w_{i},w_{j})|^{2}F_{0}(dw_{i}
)F_{0}(dw_{j})\rightarrow_{p}0,
\]
or ii) $\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}|\hat{\xi}_{l}
(W_{i},W_{j})|\rightarrow_{p}0$, or (iii) $\sqrt{n}\binom{n}{2}^{-1}
\sum_{(i,j)\in I_{l}}\hat{\xi}_{l}(W_{i},W_{j})\rightarrow_{p}0$.
\end{assumption}
These are rate conditions on the remainder term $\hat{\xi}_{l}(w_{i},w_{j})$.
For $\phi$ in ((ref)), the interaction term is a sum over $m$ of terms of the
form
\[
\hat{\xi}_{lm}(w_{i},w_{j})=\left( \hat{\alpha}_{lm}(x_{i},x_{j})-\alpha
_{0m}(x_{i},x_{j})\right) \left( c_{1m}\hat{\gamma}_{l}(X_{i})+c_{2m}
\hat{\gamma}(X_{j})-c_{1m}\gamma_{0}(X_{i})-c_{2m}\gamma_{0}(X_{j})\right) .
\]
Therefore, Assumption (ref) follows for first-steps
satisfying orthogonality restrictions if $\sqrt{n}||\hat{\alpha}_{lm}
-\alpha_{0m}||||\hat{\gamma}_{l}-\gamma_{0}||=o_{p}(1).$ This is a product rate
condition which allows for nuisance estimators to converge at slower rates as
long as the product converges faster to zero than the $n^{-1/2}$-rate. Define now $\bar{\psi
}(\gamma,\alpha,\theta) \equiv\mathbb{E}[\psi(W_{i},W_{j},\gamma,\alpha
,\theta)]$, i.e. the expectation is only over $(W_i,W_j)$. Henceforth, $C$ is a generic positive constant that may change
from to expression to expression.
\begin{assumption}
For each $l=1,...,L$ and $\theta$, i) $\int\int
\phi(w_{i},w_{j},\gamma_{0},\hat{\alpha}_{l},\theta)F_{0}(dw_{i})F_{0}
(dw_{j})=0$ with probability approaching one; and ii) $\sqrt{n}\bar{\psi}(\hat{\gamma
}_{l},\alpha_{0},\theta_{0})\rightarrow_{p}0$.
\end{assumption}
Assumption (ref) (i) incorporates the global robustness
property of $\alpha$ and is in most cases easy to check by inspection of
$\phi$. For example, it holds for $\phi$ in ((ref)).
Assumption (ref) (ii) is a small bias condition. A sufficient condition for Assumption (ref) (ii) is that $\bar{\psi}(\gamma
,\alpha_{0},\theta_{0}) \leq C || \gamma - \gamma_0||_q^\rho$ for all $\gamma$ with $||\gamma - \gamma_0||$ small enough and that $|| \hat{\gamma}_l - \gamma_0||_q = o_p(n^{-1/2\rho})$, for $\rho>0$.
\begin{lemma}
If Assumptions (ref)
-(ref) are satisfied then equation ((ref)) holds.
\end{lemma}
For valid inference, we also need convergence of the asymptotic variance
estimators. To simplify the computation, we implemented the variance estimator
without cross-fitting. Define $\hat{g}_{ij}=g(W_{i},W_{j},\hat{\gamma},\hat{\alpha}
,\hat{\theta})$ and $g_{ij}=g(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta_{0}).$
Then define the following leave-one out average $\hat{g}_{-i} = (n-1)^{-1}
\sum_{j\neq i}\hat{g}_{ij}$ and $g_{-i} = (n-1)^{-1}\sum_{j\neq i}g_{ij}$, let
$\hat{\phi}_{-i}$ and $\phi_{-i}$ be defined in the same way.
\begin{lemma}
If $\mathbb{E}[|\psi(W_{i},W_{j},\gamma_{0}
,\alpha_{0},\theta_{0})|^{2}]<\infty$, $n^{-1}\sum_{i=1}^{n}|\hat{g}_{-i} -
g_{-i}|^{2} \to_{p} 0$ and $n^{-1}\sum_{i=1}^{n}|\hat{\phi}_{-i} - \phi
_{-i}|^{2} \to_{p} 0$, then $\hat{\Sigma}\rightarrow_{p}\Sigma$.
\end{lemma}
We also need convergence of the Jacobian of the moment condition
$\hat{B}\rightarrow_{p}B$. Define $\tilde{\psi}_{ij}=\psi(W_{i},W_{j}
,\hat{\gamma},\hat{\alpha},\theta_{0})$.
\begin{assumption}
$B$ exists and there is a neighborhood
$\mathcal{N}$ of $\theta_{0}$ such that i) $||\hat{\gamma}-\gamma
_{0}||\rightarrow_{p}0$; ii) for all $||\gamma-\gamma_{0}||$ and
$||\alpha-\alpha_{0}||$ small enough $\psi(W_{i},W_{j},\gamma,\alpha,\theta)$
is differentiable in $\theta$ on $\mathcal{N}$ with probability approaching
one and there is $C>0$ and $d(W_{i},W_{j},\gamma,\alpha)$ such that
$\mathbb{E}[d(W_{i},W_{j},\gamma,\alpha)] \leq C$ and such that for $\theta
\in\mathcal{N}$ and $||\gamma-\gamma_{0}||$ and $||\alpha-\alpha_{0}||$ small
enough
\[
\left\vert \frac{\partial\psi(W_{i},W_{j},\gamma,\alpha,\theta)}
{\partial\theta}-\frac{\partial\psi(W_{i},W_{j},\gamma,\alpha,\theta_{0}
)}{\partial\theta}\right\vert \leq d(W_{i},W_{j},\gamma,\alpha)|\theta
-\theta_{0}|^{1/C}
\]
iii) For each $k$, $\mathbb{E}[|\partial\tilde{\psi}_{ij}/\partial\theta_{k} -
\partial\psi_{ij}/\partial\theta_{k}|] \to 0$.
\end{assumption}
\begin{lemma}
If Assumption (ref) is
satisfied and if $\bar{\theta} \rightarrow_{p} \theta_{0}$ then $\partial
\hat{\psi}(\bar{\theta})/\partial\theta\to_{p} B$.
\end{lemma}
Now we are ready to give the main asymptotic result of the paper,
\begin{theorem}
If Assumptions (ref)-(ref)
and conditions in Lemma (ref) are satisfied, $\hat{\theta
}\rightarrow_{p}\theta_{0}$ and $V=B^{-1}\Sigma B^{\prime-1}$ is nonsingular,
then
\[
\sqrt{n}(\hat{\theta}-\theta_{0})\rightarrow_{d}\mathcal{N}(0,V).
\]
Also, $\hat{V}\rightarrow_{p}V$.
\end{theorem}
Theorem (ref) assumes consistency of $\hat{\theta}$, which Appendix (ref) establishes under mild conditions.
\textsc{Remark }(\textbf{On degeneracy of orthogonal moments}). When
\[
\mathbb{E}[\psi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta_{0})|W_{i}]=0 \text{ a.s},
\]
then $\Sigma=0$ and $V=0$. In this case, we say that the kernel $\psi$ or the U-statistic is degenerate. In this situation, our results show the estimator converges to $\theta_0$
faster than $n^{-1/2}$. Therefore, our asymptotic distribution theory provides distributional results only for
non-degenerate orthogonal moments. We focus on this case because it is the
most common one in the applications we are interested in. A complete analysis of the
degenerate case is beyond the scope of this paper and is deferred to future
research. $\blacksquare$
\subsection{Asymptotic properties of the debiased IOp}
We give lower-level conditions for the IOp example. To simplify the notation, define the difference of fitted values $\Delta_{\gamma}(X_{i},X_{j})=\gamma(X_{i})-\gamma(X_{j})$, and, in particular, $\Delta_{0}\equiv\Delta_{\gamma_{0}}(X_{i},X_{j})$.
\begin{assumption}
(i) $||\hat{\gamma}_{l} - \gamma_{0}||_1 = o_p(1)$; (ii) $ P(\Delta_0=0| X_i \neq X_j) = 0$; (iii) $P(0 < |\Delta_0| \leq t ) \leq C t^\beta$ \text{ for } $\beta > 0$, and for all $t>0$ sufficiently small.
\end{assumption}
Assumption (ref) (i) requires $\hat{\gamma}_{l}$ to be $L_1$ consistent. Assumption (ref) (ii) requires that, given $X_i \neq X_j$, the mass at zero of $|\Delta_0|$ is zero. This holds if $\Delta_0$ has an absolutely continuous distribution or in the discrete case if $\gamma_0$ is injective. Assumption (ref) (iii) requires that the mass of $|\Delta_0|$ strictly above zero disappears at a given rate. If $\Delta_0$ has a Lebesgue density that is bounded at $0$ then it holds with $\beta = 1$. For $\beta = \infty$, the mass around zero disappears. This is the case when $\Delta_0$ is discrete since there exists $\tilde{t}$ small enough such that $P(0 < |\Delta_0| \leq t) = 0$ for any $t \leq \tilde{t}$. The next result uses (i) and (ii) to show mean-square consistency of the sign difference and (iii) to control the smoothness of the functional $\gamma \rightarrow\mu(\Delta_{\gamma})=\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0]$.
\begin{proposition}
Under Assumption (ref) (i) and (ii),
\[
\int\int|sgn(\Delta_{\hat{\gamma}_l}) - sgn(\Delta_0)|^{2} F_{0}(dw_{i})F_{0}(dw_{j}) \to_{p} 0,
\]
and under Assumption (ref) (iii),
\[
\mu(\Delta_{\gamma}) \leq
\begin{cases}
C ||\gamma - \gamma_0||_q^{\frac{q(1+\beta)}{q+ \beta}} &\text{ if } q \in [1,\infty),\\
C ||\gamma - \gamma_0||_\infty^{1+\beta} &\text{ if } q = \infty.
\end{cases}
\]
\end{proposition}
In $L_2$ norm (i.e., $q=2$), we have that $\mu(\Delta_{\gamma}) \leq C ||\gamma - \gamma_0||^{2 \frac{1+\beta}{2+\beta}}$, so in the continuous case with bounded density at 0 (i.e., $\beta=1$) we achieve a smoothness exponent of $4/3$ while in the discrete case ($\beta = \infty$) we achieve a quadratic bound. This Proposition is an improvement upon the results in clemenccon2011minimax.\footnote{They assume $P(|\gamma_0(X_i) - \gamma_0(x)| \leq t) \leq C t^\beta$ for all $x$. This implies Assumption (ref)(iii) for $\beta \leq 1$. In contrast, we allow for $\beta > 1$, which is key for the rates in Assumption (ref). As a result, we recover Massart’s margin condition from classification (massart2000some; audibert2007fast) in a ranking setting.} The next assumption relates to the trade-off between smoothness properties and required rates of the ML estimators.
\begin{assumption}
$||\hat{\gamma}_{l} - \gamma_{0}||_q = o_p(n^{-\rho_{q\beta}})$ where $\rho_{q\beta} \geq \frac{1}{2q} \cdot \frac{q+\beta}{1+\beta} \text{ if } q \in [1,\infty)$, and $\rho_{q\beta} \geq \frac{1}{2(1+\beta)} \text{ if } q = \infty$.
\end{assumption}
There are several ways in which one can achieve the asymptotic normality results, depending on the margin parameter $\beta$ and the rate of the ML estimator under different norms. Table (ref) shows the most common configurations. In the discrete case ($\beta \to \infty$) and under the supremum norm, we allow for slow convergence rates, even slower than $n^{-1/4}$.
\begin{table}[h!]
\begin{tabular}{@ccc@}
\toprule
\textbf{\(q\)} & \textbf{\(\beta = 1\)} & \textbf{\(\beta = \infty\)} \\
\midrule
\(2\) & \(\displaystyle 3/8 \) & \(\displaystyle 1/4\) \\
\(\infty\) & \(\displaystyle 1/4 \) & \(\displaystyle 0\) \\
\bottomrule
\end{tabular}
\caption{ ML rates required for different margin parameters in $L_2$ and $L_\infty$ norms.}
\end{table}
In our empirical application $\beta = \infty$ so the $L_2$ nonparametric rates of $n^{-1/4}$ typically imposed in the DML literature suffice. For consistent estimation of the variance, we assume the following.
\begin{assumption}
$\mathbb{E}[(\hat{\gamma}(X_{i}) - \gamma
_{0}(X_{i}))^{2}] = o(1)$.
\end{assumption}
This assumption strengthens Assumption (ref) (i) to mean-square convergence (unconditionally). It can be discarded at the cost of cross-fitting the variance estimator.
\begin{proposition}
Let Assumptions (ref), (ref) and
(ref) hold. Assuming further that either $\mathbb{E}[Y_i^2 |X_i] < C < \infty$ a.s. or that $\mathbb{E}[|Y_i - Y_j|^{2+\delta}] < \infty$ for some $\delta > 0$, we have $\sqrt{n}(\hat{\theta}-\theta_{0})\rightarrow_{d}\mathcal{N}(0,V)$ and
$\hat{V}\rightarrow_{p}V$, where $\hat{V}$ is given in ((ref)).
\end{proposition}
\section{Monte Carlo Simulations}
In this section, we evaluate the finite sample performance of the proposed debiased IOp estimator in comparison with the commonly used plug-in estimator. We consider two DGPs: (i) a Gaussian First-Step Error Model that serves to illustrate the main points of the paper; and (ii) an empirically motivated DGP based on the Spanish data used in our empirical application below.
\subsection{A Gaussian First-Step Error Model}
In this DGP, $X \sim \mathcal{N}(0,\sigma_X^2)$ and
\[
Y = \gamma_0(X) + \varepsilon, \text{ where } \gamma_0(X) = 5 + X^2,
\]
where $\varepsilon \sim \mathcal{N}(0, \sigma_\varepsilon^2)$, $\sigma_\varepsilon^2 = \mathbb{V}ar(\gamma_0(X))/StN$, and $StN$ is the Signal to Noise Ratio which is fixed by us. This DGP is called a Gaussian First-Step Error Model because we generate
\[
\hat{\gamma}(x) = \gamma_0(x) + \delta_n(x), \text{ where } \delta_n(x) \sim \mathcal{N}\left(c_1(\mathbb{E}[\gamma_0(X)] - \gamma_0(x))n^{-\rho_1}, \frac{c_2}{StN} n^{-\rho_2}\right),
\]
where $c_1, c_2, \rho_1, \rho_2$ are constants. This DGP helps us check asymptotic theory with a full control over the rates of $\hat{\gamma}$. A nice feature of this DGP is that we can easily check that the key derivative term in ((ref)) is normally distributed. Specifically,
\begin{equation*}
\sqrt{n}\frac{\partial\theta(\gamma_{0})}{\partial\gamma}\left[ \hat{\gamma}
-\gamma_{0}\right]\sim \mathcal{N}\left(C_{1}n^{0.5-\rho_1}, \sigma^2_{n}\right),
\end{equation*}
where $C_{1}=c_{1}\theta_0$, $\sigma^2_{n}=C_{2}n^{1-\rho_2}$ and $C_{2}=c_{2}(\mathbb{E}[2Y])^{-2}\theta^2_{0}/StN$. Computations for the Gaussian First-Step Error Model can be found in Online Appendix (ref). Therefore, the plug-in bias coming from this derivative term diverges when $\rho_1 < 1/2$ and $\rho_2 < 1$. This illustrates the bias problem and the lack of root-$n$ consistency of plug-in estimators based on non-locally robust moment functions, as mentioned in the Introduction.
We can also show that Assumption (ref) (iii) holds with any $\beta < 1$ (see Online Appendix (ref)). Hence, a sufficient condition for valid inference of the debiased estimator requires rates faster than $3/8$ in $L_2$ norm. This implies $\rho_1 > 3/8$ and $\rho_2 > 3/4$. In the simulations, we present slower and faster rates than required. We do so to evaluate the robustness of the finite sample performance to our sufficient conditions. The lower the $\beta$, the less smooth the model and the harder it is to control the nonlinear bias terms. Hence, with this DGP, we are dealing with a difficult case ($\beta<1$).
In Table (ref), we report the bias of the the Plug-in estimator and coverage of the associated confidence intervals in a Monte Carlo simulation for several specifications of the constants $c_1, c_2, \rho_1, \rho_2$, $\sigma_X^2$ and $StN$. We report coverage rates of two different confidence intervals based on the plug-in estimator: (i) a blind application of the plug-in estimator that does not account for the first-step effect on the asymptotic standard error (naive coverage); and (ii) a correction of the standard error for the plug-in estimator that accounts for the first-step error (using the same standard error formula as for the debiased estimator). We also report the confidence interval length based on the debiased estimator.
The results are consistent with our theory. The plug-in method is highly affected by the slow convergence of $\hat{\gamma}$, showing large bias and associated confidence intervals with substantial coverage distortions. In contrast, the debiased IOp estimator has small bias and the confidence interval based on it has an accurate coverage uniformly across different sample sizes and parameter configurations (including deviations from our sufficient conditions).
\begin{table}[!h]
\caption{Gaussian First-step error model simulations}
\fontsize{8}{10}\selectfont
\begin{threeparttable}
\begin{tabular}[t]{cccccccccc}
\toprule
\textbf{StN} & \textbf{$\rho_1$} & \textbf{$\rho_2$} & \textbf{$n$} & \textbf{BP} & \textbf{CP} & \textbf{CPN} & \textbf{BD} & \textbf{CD} & \textbf{CI length}\\
\midrule
0.1 & 0.26 & 0.6 & 1000 & -0.011 & \textcolor{red}{1.000} & \textcolor{red}{0.165} & -0.005 & 0.935 & 0.057\\
& & & 3000 & -0.009 & 0.999 & \textcolor{red}{0.020} & -0.003 & 0.934 & 0.033\\
& & & 6000 & -0.008 & 0.974 & \textcolor{red}{0.001} & -0.002 & 0.935 & 0.023\\
\midrule
& & 0.8 & 1000 & -0.016 & \textcolor{red}{1.000} & \textcolor{red}{0.017} & -0.002 & 0.945 & 0.057\\
& & & 3000 & -0.012 & 0.971 & \textcolor{red}{0.001} & -0.001 & 0.936 & 0.033\\
& & & 6000 & -0.010 & \textcolor{red}{0.751} & \textcolor{red}{0.000} & -0.001 & 0.945 & 0.023\\
\midrule
& 0.5 & 0.6 & 1000 & 0.003 & \textcolor{red}{1.000} & 0.920 & -0.004 & 0.945 & 0.057\\
& & & 3000 & 0.002 & \textcolor{red}{1.000} & \textcolor{red}{0.893} & -0.002 & 0.934 & 0.033\\
& & & 6000 & 0.001 & \textcolor{red}{1.000} & 0.910 & -0.002 & 0.942 & 0.023\\
\midrule
& & 0.8 & 1000 & -0.002 & \textcolor{red}{1.000} & 0.923 & -0.002 & 0.948 & 0.057\\
& & & 3000 & -0.001 & \textcolor{red}{1.000} & 0.915 & -0.001 & 0.936 & 0.033\\
& & & 6000 & -0.001 & \textcolor{red}{1.000} & 0.911 & -0.001 & 0.947 & 0.023\\
\midrule
1 & 0.26 & 0.6 & 1000 & -0.017 & \textcolor{red}{0.115} & \textcolor{red}{0.010} & -0.001 & 0.947 & 0.024\\
& & & 3000 & -0.013 & \textcolor{red}{0.009} & \textcolor{red}{0.001} & 0.000 & 0.935 & 0.014\\
& & & 6000 & -0.011 & \textcolor{red}{0.001} & \textcolor{red}{0.000} & 0.000 & 0.936 & 0.010\\
\midrule
& & 0.8 & 1000 & -0.017 & \textcolor{red}{0.085} & \textcolor{red}{0.006} & 0.000 & 0.949 & 0.024\\
& & & 3000 & -0.013 & \textcolor{red}{0.005} & \textcolor{red}{0.000} & 0.000 & 0.940 & 0.014\\
& & & 6000 & -0.011 & \textcolor{red}{0.000} & \textcolor{red}{0.000} & 0.000 & 0.934 & 0.010\\
\midrule
& 0.5 & 0.6 & 1000 & -0.003 & 0.983 & \textcolor{red}{0.889} & -0.001 & 0.947 & 0.024\\
& & & 3000 & -0.001 & 0.982 & \textcolor{red}{0.893} & 0.000 & 0.938 & 0.014\\
& & & 6000 & -0.001 & 0.982 & \textcolor{red}{0.893} & 0.000 & 0.935 & 0.010\\
\midrule
& & 0.8 & 1000 & -0.003 & 0.979 & \textcolor{red}{0.861} & 0.000 & 0.950 & 0.024\\
& & & 3000 & -0.002 & 0.980 & \textcolor{red}{0.873} & 0.000 & 0.942 & 0.014\\
& & & 6000 & -0.001 & 0.973 & \textcolor{red}{0.865} & 0.000 & 0.933 & 0.010\\
\bottomrule
\end{tabular}
\begin{tablenotes}[para]
• Simulation Results 1000 iterations. sd(X) = 1, $c_1 = 1$, $c_2 = 0.5$. BP: Bias Plug-in, CP: Coverage Plug-in, CPN: Coverage Plug-in Naive, BD: Bias Debiased, CD: Coverage Debiased.
\end{tablenotes}
\end{threeparttable}
\end{table}
\subsection{Empirically-Based DGP}
In this section, we calibrate a DGP to the data used in our empirical application. We take the Spanish data set and run a regression of log income on the education of the mother, the education of the father, and the occupation of the father. In this section, $X_i$ will denote these three variables for individual $i$. We then take the fitted values and the variance of the residuals, which we denote by $\sigma_\varepsilon^2$, to generate new data $Y_i = \exp(\gamma_0(X_{i}) + \varepsilon_i)$, where $\varepsilon_i \sim \mathcal{N}(0,\sigma_\varepsilon^2)$ and $\gamma_0(X_{i})$ is sampled with replacement from the OLS fitted values. The econometrician does not know which circumstances enter the DGP and employs all the circumstances available in our empirical application. She uses grouped Lasso to select relevant variables with the R package \texttt{grpreg} (see breheny2015group).
The DGP employed has $52$ groups with distinct conditional mean and a true IOp of $0.085$. When the econometrician does not know the variables that enter the DGP, the possible number of groups she might entertain is much larger. We feed the MLs $Y$ and not $\ln(Y)$ which in this setting makes it harder for the MLs to detect the signal, this is reflected in a realistically low Signal to Noise Ratio (StN) of 3%.
For the first-step estimator, we use RF and CIF. The latter is popular in the IOp literature (e.g. brunori2021roots) and the former is more popular in general. We use the R packages \texttt{ranger} and \texttt{party}. As a benchmark, we provide results for a parametric approach that runs OLS of $\ln(Y)$ on the correct set of circumstances (see ferreira2011measurement). The fitted values $\hat{\gamma}(x)$ from this log-linear regression consistently estimate $\gamma_0(x)$. As predictions we take $\exp(\hat{\gamma}(X))$. Since the Gini is scale-invariant, the Gini of $\exp(\gamma_0(X))$ is the same as the Gini of $\mathbb{E}[Y|X] = \exp(\gamma_0(X))\sigma^2_\varepsilon/2$. Hence, the parametric model is correctly specified (up to scale) and expected to perform well. For the plug-in confidence intervals, we report two coverages, one based on the correct asymptotic variance, which accounts for the first-step estimation, and a naive one that ignores the first-step estimation.
In Table (ref) we show the results of the simulations. For RF and CIF, the plug-in estimator is generally more biased and the coverages of its associated confidence intervals are far from the nominal level, particularly when using naive standard errors. The debiased estimators show very low bias and lead to accurate inference. As expected, the parametric approach displays very low bias even in its plug-in version. The coverage of the associated confidence intervals is only valid when the first-step estimation effect is taken into account, which we are the first to do. Our debiased ML methods show remarkably similar performance to the well-specified parametric approach. We also show average confidence interval length, which is decreasing in sample size, as expected, and comparable with its parametric counterpart, illustrating the adaptive and efficient features of the proposed debiased inferences.
\begin{table}[H]
\caption{Simulation results with 1000 iterations}
\begin{tabular}[t]{llcccccc}
\toprule
\multicolumn{2}{c} & \multicolumn{3}{c}{Plug-in} & \multicolumn{3}{c}{Debiased} \\
\cmidrule(l{3pt}r{3pt}){3-5} \cmidrule(l{3pt}r{3pt}){6-8}
ML & n & Bias & Coverage & Coverage naive & Bias & Coverage & Avg. CI length\\
\midrule
\addlinespace[0.3em]
\multicolumn{8}{l}{\textbf}\\
\textbf{RF} & 1000 & -0.012 & 0.916 & \textcolor{red}{0.186} & 0.002 & 0.953 & 0.081\\
& 3000 & -0.012 & \textcolor{red}{0.806} & \textcolor{red}{0.140} & -0.001 & 0.937 & 0.047\\
& 6000 & -0.011 & \textcolor{red}{0.741} & \textcolor{red}{0.097} & -0.001 & 0.923 & 0.033\\
\addlinespace[0.3em]
\multicolumn{8}{l}{\textbf}\\
\textbf{CIF} & 1000 & -0.033 & \textcolor{red}{0.621} & \textcolor{red}{0.081} & -0.008 & 0.901 & 0.077\\
& 3000 & -0.022 & \textcolor{red}{0.550} & \textcolor{red}{0.034} & -0.003 & 0.922 & 0.046\\
& 6000 & -0.016 & \textcolor{red}{0.495} & \textcolor{red}{0.037} & -0.002 & 0.934 & 0.032\\
\addlinespace[0.3em]
\multicolumn{8}{l}{\textbf}\\
\textbf{Parametric} & 1000 & 0.022 & \textcolor{red}{0.737} & \textcolor{red}{0.076} & -0.006 & 0.918 & 0.077\\
& 3000 & 0.012 & \textcolor{red}{0.868} & \textcolor{red}{0.106} & -0.003 & 0.926 & 0.046\\
& 6000 & 0.008 & \textcolor{red}{0.887} & \textcolor{red}{0.103} & -0.001 & 0.956 & 0.032\\
\bottomrule
\end{tabular}
\end{table}
\section{Inequality of Opportunity in Europe}
We measure IOp in 29 European countries using the 2019 EUSILC. Income is equivalized household income, the unit is the individual and we restrict the sample to ages 25–59.
Circumstances include parental characteristics and aspects of the respondent’s household and living conditions around age 14: sex, country of birth, whether he/she was living with the father, the number of adults/working adults/kids in the household,
population of the municipality, tenancy of the house, country of birth of the mother, education of the parents, occupational
status of the parents, father's managerial position, father's occupation,
financial situation and whether school materials were accessible. Remember that it refers to when the individual was around
14 years old, so, for instance, financial situation refers to the financial
situation of the household where the individual resided when he/she was around
14 years old. All circumstances are discrete, so there are many different combinations of
the categories. As in the simulations, we first use Grouped Lasso to select among all circumstances.
Figure (ref) shows debiased and plug-in estimates of relative IOp with 95% CIs for the debiased estimator. The wider CIs use the Šidák correction (see vsidak1967rectangular) for simultaneous inference across 29 countries. We also report the (now likely misspecified) parametric approach (see ferreira2011measurement) used in the simulations with bootstrap standard errors. For ML, we use CIF, the most common ML in the IOp literature; results for other MLs appear in Online Appendix (ref). The plug-in estimator systematically underestimates IOp and often falls outside the debiased CI.
\begin{figure}[H]
\caption{IOp with CIF}
\end{figure}
We see a lot of heterogeneity in IOp across countries.
Nordic countries such as Denmark,
Finland, Sweden or Norway are countries with low relative IOp. Netherlands and Germany are also in
the lower range. Eastern
countries are heterogeneous, Czechia or Hungary have lower IOp while Romania or Bulgaria have very high IOp (around 60%). Southern countries have high levels of IOp. We report more detailed results in Table (ref).
\begin{table}[h]
\caption{ Results for CIF. Hyperparameter tuning grids: Number of Trees: \{300, 500, 800\}, CIF Depth: \{1, 3, 6\}, Number of variables to split (mtry): \{5, 6, 7, 8, 9\}}
\fontsize{10}{12}\selectfont
\begin{tabular}[t]{cccccc}
\toprule
Country & Mean & Gini & Plug in/Gini & Debiased/Gini & n\\
\midrule
Austria & 29783 & 0.268 & 22 % & 26 % (22%,30%) & 5150\\
Belgium & 28423 & 0.237 & 34 % & 38 % (35%,41%) & 6368\\
Bulgaria & 6603 & 0.415 & 49 % & 60 % (57%,63%) & 5906\\
Switzerland & 52732 & 0.273 & 19 % & 26 % (21%,31%) & 4715\\
Cyprus & 20281 & 0.299 & 30 % & 37 % (32%,42%) & 4647\\
Czechia & 12200 & 0.232 & 24 % & 29 % (26%,32%) & 6881\\
Germany & 28799 & 0.271 & 16 % & 21 % (17%,26%) & 5981\\
Denmark & 36736 & 0.260 & 5 % & 7 % (1%,13%) & 2076\\
Estonia & 14307 & 0.283 & 22 % & 27 % (24%,30%) & 5669\\
Greece & 9915 & 0.310 & 34 % & 38 % (35%,41%) & 14144\\
Spain & 17707 & 0.327 & 38 % & 42 % (40%,45%) & 16975\\
Finland & 29297 & 0.263 & 10 % & 11 % (5%,17%) & 4403\\
France & 27001 & 0.274 & 27 % & 34 % (31%,37%) & 7924\\
Croatia & 8916 & 0.281 & 31 % & 34 % (31%,37%) & 7120\\
Hungary & 6925 & 0.282 & 26 % & 28 % (22%,33%) & 4568\\
Ireland & 32224 & 0.278 & 31 % & 37 % (32%,42%) & 3660\\
Italy & 20006 & 0.319 & 32 % & 38 % (36%,40%) & 16360\\
Lithuania & 10383 & 0.347 & 29 % & 35 % (30%,40%) & 3643\\
Luxembourg & 45435 & 0.329 & 32 % & 40 % (35%,44%) & 3533\\
Latvia & 10787 & 0.336 & 27 % & 33 % (29%,37%) & 3344\\
Malta & 19263 & 0.266 & 29 % & 35 % (31%,39%) & 3632\\
Netherlands & 30629 & 0.259 & 15 % & 20 % (12%,29%) & 4446\\
Norway & 41869 & 0.245 & 15 % & 24 % (17%,31%) & 2366\\
Poland & 8621 & 0.292 & 34 % & 38 % (35%,40%) & 14293\\
Portugal & 12030 & 0.305 & 39 % & 41 % (38%,43%) & 13783\\
Romania & 4826 & 0.343 & 47 % & 56 % (53%,59%) & 5932\\
Serbia & 3984 & 0.329 & 36 % & 41 % (38%,44%) & 5648\\
Sweden & 29282 & 0.300 & 17 % & 24 % (14%,34%) & 2027\\
Slovakia & 9197 & 0.221 & 28 % & 32 % (28%,35%) & 5727\\
\bottomrule
\end{tabular}
\end{table}
Interpreting the direction and size of the difference between debiased and plug-in estimates is challenging. Although here the plug-in systematically underestimates IOp, Online Appendix (ref) shows that this depends on the ML method employed. A useful heuristic is to decompose the difference between the debiased and plug-in estimators into three terms.
\begin{align*}
\hat{\theta} - \hat{\theta}^P &= \frac{\sum_{l=1}^L \sum_{(i,j) \in I_l} (sgn(\hat{\gamma}_l(X_i) - \hat{\gamma}_l(X_j)) - sgn(\hat{\gamma}(X_i) - \hat{\gamma}(X_j)))(Y_i - Y_j)}{\sum_{i<j} Y_i + Y_j} \\
&+ \frac{\sum_{i<j} sgn(\hat{\gamma}(X_i) - \hat{\gamma}(X_j))(Y_i - \hat{\gamma}(X_i) - Y_j + \hat{\gamma}(X_j))}{\sum_{i<j} Y_i + Y_j} \\
&+ \sum_{i<j} |\hat{\gamma}(X_i) - \hat{\gamma}(X_j)| \left(\frac{1}{\sum_{i<j} Y_i + Y_j} - \frac{1}{\sum_{i<j} \hat{\gamma}(X_i) + \hat{\gamma}(X_j)}\right).
\end{align*}
The first term captures the overfitting that cross-fitting removes. The second corrects the numerator holding fixed the denominator: it can be seen either as the covariance between prediction orderings and prediction errors --- positive when higher/lower predictions are systematically under/overestimated, as in income data --— or as a correction for regularization bias, which tends to compress prediction differences relative to income differences. In both interpretations, the term is expected to be positive, though it may turn negative if prediction orderings are systematically wrong. The third term replaces predicted incomes with actual incomes in the denominator. In Figure (ref) we report these terms divided by the Gini of income, since our interest lies in relative IOp measures. The white circles show the sum of all the terms.
\begin{figure}[H]
\caption{Debiased minus plug-in decomposition for CIF}
\end{figure}
The denominator effect is generally negligible since the mean of predictions is very similar to the mean of the outcome. The regularization term turns out to be always positive (i.e. causes plug-in to underestimate) and the overfitting causes the plug-in to overestimate. For the CIF, the regularization term is larger, causing the plug-in estimator to underestimate in all countries.
In Online Appendix (ref) we show that the signs of the regularization and overfitting terms broadly coincide across MLs. However, the sign of the total difference can differ per ML. In Online Appendix (ref) we also show the relation of the different biases with sample size. For instance, the strong overfitting bias in Scandinavian countries is likely due to small sample sizes. Hence, the resulting difference is a complex interplay between regularization and overfitting, the choice of the ML and the noise in the data.
Finally, in Figure (ref) we illustrate a very attractive feature of our estimator. The debiased estimates are much less sensitive to the choice of the ML. We estimate debiased and plug-in IOp for all countries with 6 different MLs (Lasso, Ridge, Random Forests, CIF, XGBoost and Catboost) and plot all point estimates. We can
see that the debiased estimates are much more concentrated and near to each
other than the plug-in estimates. The plug-in estimates are much more
dispersed and in some cases there are differences of almost 40 percentage points from
one machine learner to the other. Even differences of around 15 percentage points
between plug-in estimates using different machine learners are not uncommon.
This result alone shows the importance of using locally robust debiased
estimators when estimating IOp.
\begin{figure}[h!]
\caption{Sensitivity of Plug-in and Debiased to choice of ML in the first-step.}
\end{figure}
\section{Additional Example Applications}
\subsection{Inference on ML Performance through the AUC}
As mentioned in the Introduction, the Area Under the Curve (AUC), where the curve refers to the ROC (receiver
operating characteristic) curve, is one of the most commonly used measures
to evaluate the performance of an ML\ binary classifier, see bradley1997use. In this section, we show that the AUC falls under our setting. We have a binary outcome $Y_i$ and covariates $X_i$. We restrict our attention to the case where an ML estimator $\hat{\gamma}$ of $\gamma_0(x) = \mathbb{E}[Y_i| X_i = x]$ is used to classify observations. The AUC parameter is defined as
\begin{align*}
\theta_0 &= \mathbb{E}[1(\gamma_0(X_i) > \gamma_0(X_j)) | Y_i = 1, Y_j = 0] + \frac{1}{2}\mathbb{E}[1(\gamma_0(X_i) = \gamma_0(X_j)) | Y_i = 1, Y_j = 0],
\end{align*}
which can be rewritten in the following form to make it a particular case of our theory
\begin{align*}
\theta_0 &= \frac{\mathbb{E}\left[(1(\gamma_0(X_i) > \gamma_0(X_j)) + \frac{1}{2}1(\gamma_0(X_i) = \gamma_0(X_j))) Y_i(1-Y_j) \right]}{2p_{0}(1-p_{0})} \\
&+ \frac{\mathbb{E}\left[(1(\gamma_0(X_i) < \gamma_0(X_j)) + \frac{1}{2}1(\gamma_0(X_i) = \gamma_0(X_j))) Y_j(1-Y_i)\right]}{2p_{0}(1-p_{0})}
\end{align*}
where $p_{0} = P(Y_i = 1)$.
Our next result shows that the AUC is a known functional of our debiased Gini coefficient and the probability $p_{0}$.
\begin{proposition}
The AUC can be written as
\[
\theta_0 = \frac{1}{2} + \frac{\mathbb{E}[sgn(\gamma_0(X_i) - \gamma_0(X_j))(Y_i - Y_j)]}{4p_{0}(p_{0}-1)} = \frac{G(\gamma_0(X_i)) + p_{0} - 1}{2(p_{0}-1)},
\]
\end{proposition}
\begin{proof} The proof follows from the identities
\begin{align*}
1(\gamma_0(x_i) > \gamma_0(x_j)) + \frac{1}{2}1(\gamma_0(x_i) = \gamma_0(x_j)) &= [1+sgn(\gamma_0(x_i) - \gamma_0(x_j))]/2, \\
1(\gamma_0(x_i) < \gamma_0(x_j)) + \frac{1}{2}1(\gamma_0(x_i) = \gamma_0(x_j)) &= [1-sgn(\gamma_0(x_i) - \gamma_0(x_j))]/2,
\end{align*}
\end{proof}
There are several important implications from this result. First, it shows that the AUC is locally robust to the ML first-step. This result provides theoretical support for the view held in applied fields that the AUC is a robust measure of predictive accuracy (see, e.g., bradley1997use). To our knowledge, we are the first to provide theoretical guarantees on the robustness of the AUC (the claim in the literature was based on Monte Carlo experiments). Second, in combination with our previous results, it shows the validity of asymptotic inferences and standard errors based on cross-fitted implementations of the AUC, which we recommend as best practice for applications.
\subsection{A Debiased Difference Estimator for ML-Sample Selection}
In this section, we apply our methodology to propose a debiased pairwise difference estimator for the linear regression model subject to sample selection (see heckman1974shadow, heckman1976common). We consider a two-step pairwise difference estimator, as in ahn1993semiparametric, but with an ML
estimator for the selection equation (e.g., employing modern model selection or regularization tools in the selection equation). We show that the simple
plug-in pairwise estimator of ahn1993semiparametric is biased by the first-step, and provide, to the best of our knowledge, the first debiased pairwise difference estimator for an ML-selection model.
For simplicity of exposition, we focus on a semiparametric partially linear sample selection model, although our results have applications more widely, see, e.g., honore2005identification, jochmans2013pairwise, and other
references cited in the Introduction. The linear sample selection model is
defined by $Y_{i}=Y_{i}^{\ast }D_{i}, $ and the equations
\begin{eqnarray}
Y_{i}^{\ast } &=&S_{i}^{\prime }\theta _{0}+\varepsilon _{i},
\\
D_{i} &=&1(\eta _{0}(X_{i})\geq V_{i}), \notag
\end{eqnarray}
where $\theta _{0}$ is the parameter of interest, $Y_{i}$ is an outcome, $
X_{i}$ is a vector of exogenous covariates that contains $S_{i},$ $D_{i}$ is
a selection indicator, $\eta _{0}$ is an unknown function, and $(\varepsilon
_{i},V_{i})$ are unobservable error terms satisfying an index restriction so that
\begin{equation*}
\mathbb{E}[Y_{i}|X_{i},D_{i}=1]=S_{i}^{\prime }\theta _{0}+\lambda
_{0}(\gamma _{0}(X_{i})),
\end{equation*}
where $\gamma _{0}(X_{i})=\mathbb{E}[D_{i}|X_{i}]$ is the selection
population fitted value and $\lambda _{0}(\gamma _{0}(X_{i}))=\mathbb{E}
[\varepsilon _{i}|X_{i},\eta _{0}(X_{i})\geq V_{i}]$. The data is a random
sample of $W_{i}=(Y_{i},X_{i},D_{i}).$
Let $U_{i}=Y_{i}-S_{i}^{\prime }\theta _{0}-\lambda _{0}(\gamma _{0}(X_{i}))$
denote the prediction error, and notice that
\begin{eqnarray*}
Y_{i}-Y_{j} &=&\left( S_{i}-S_{j}\right) ^{\prime }\theta _{0}+\lambda
_{0}(\gamma _{0}(X_{i}))-\lambda _{0}(\gamma _{0}(X_{j}))+U_{i}-U_{j} \\
&\approx &\left( S_{i}-S_{j}\right) ^{\prime }\theta _{0}+U_{i}-U_{j}\text{
if }\gamma _{0}(X_{i})\approx \gamma _{0}(X_{j}).
\end{eqnarray*}
This motivates a local weighted regression of $Y_{i}-Y_{j}$ on $S_{i}-S_{j}$. Introduce the kernel weights
\begin{equation*}
\varpi _{ij}(\gamma )=\frac{1}{h}K\left( \frac{\gamma (X_{i})-\gamma (X_{j})
}{h}\right) D_{i}D_{j},
\end{equation*}
where $h\downarrow 0$ is a bandwidth parameter, and $K$ is a symmetric (fourth-order) kernel function.
By choosing some \textquotedblleft instruments\textquotedblright\ $
Z_{i}=Z(X_{i})$, ahn1993semiparametric suggested an (approximate)
instrumental variables identifying equation
\begin{equation}
\mathbb{E}[g(W_{i},W_{j},\gamma _{0},\theta _{0})]\approx 0,
\end{equation}
with $g(W_{i},W_{j},\gamma ,\theta )=\varepsilon _{ij}(\theta )\varpi
_{ij}(\gamma )(Z_{i}-Z_{j})\ $and the error term $\varepsilon _{ij}(\theta
)=Y_{i}-Y_{j}-\left( S_{i}-S_{j}\right) ^{\prime }\theta $. For simplicity
of presentation, we consider the just-identified case.
ahn1993semiparametric used the estimating equation ((ref)) as the basis to propose a computationally simple plug-in pairwise difference estimator
\begin{equation*}
\hat{\theta}^{P}=\hat{\Sigma}_{ZS}^{-1}\hat{\Sigma}_{ZY},
\end{equation*}
where for generic (conformable) random vectors $A$ and $B,$ and weights $\hat{\varpi}_{ij} \equiv \varpi _{ij}(\hat{\gamma})$,
\begin{equation*}
\hat{\Sigma}_{AB} =\binom{n}{2}^{-1}\sum_{i<j}(A_{i}-A_{j})(B_{i}-B_{j})^{
\prime }\hat{\varpi}_{ij}.
\end{equation*}
The estimator $\hat{\theta}^{P}$ has several advantages compared to alternative semiparametric estimators, such as having a simple closed form expression and an intuitive matching and instrumental variables interpretation. However, as we show below, $\hat{\theta}^{P}$ is biased by the ML first-step estimation. In this paper, we propose a debiased version of ahn1993semiparametric pairwise difference estimator $\hat{\theta
}^{P}$ when the first-step is an ML estimator $\hat{\gamma}$ (e.g., a Random Forest estimator). The estimator $\hat{\theta}^{P}$ solves a U-statistic estimating equation for each choice of $h,$ and hence falls into
our setting for a fixed $h$. Since we need to make $h\downarrow 0$ to obtain a consistent estimator, we need to generalize our asymptotic expansions above. To simplify the notation, we drop the dependence on $h$ when is
convenient.
Our first result derives the adjustment term for the pairwise difference estimator. The assumptions needed for this section are essentially those of
ahn1993semiparametric and they are all gathered, together with the proofs, in the Online Appendix (ref).
To simplify notation, define the weights
\begin{equation*}
q_{ij}(\gamma )=\frac{1}{h^{2}}\dot{K}\left( \frac{\gamma (X_{i})-\gamma
(X_{j})}{h}\right) \gamma (X_{i})\gamma (X_{j}),
\end{equation*}
where $\dot{K}(u)=\partial K(u)/\partial u$. We consider the nonparametric case where $\hat{\gamma}$ is a
consistent estimator for $\gamma _{0}(X_{i})=\mathbb{E}[D_{i}|X_{i}]$.
\begin{proposition}
Under Assumptions (ref) and (ref) in the Online Appendix (ref), the adjustment term of the pairwise difference estimator
is
\begin{equation*}
\phi (W_{i},W_{j},\gamma _{0},\alpha_{0},\theta )=\alpha
_{0}(X_{i},X_{j})(D_{i}+D_{j}-\gamma _{0}(X_{i})-\gamma _{0}(X_{j})),
\end{equation*}
where $\alpha _{0}(X_{i},X_{j})=\left[ \lambda _{0}(\gamma _{0}(X_{i}))-\lambda
_{0}(\gamma _{0}(X_{j}))\right] q_{ij}(\gamma _{0})(Z_{i}-Z_{j}).$
\end{proposition}
To implement the debiased pairwise estimator, introduce the short notation $
v_{i}:=v_{i}(\gamma _{0})=D_{i}-\gamma _{0}(X_{i})$ and $v_{ij}(\gamma
_{0})=v_{i}(\gamma _{0})-v_{j}(\gamma _{0}).$ The estimated debiased moment function is given by
\begin{equation*}
\psi (W_{i},W_{j},\hat{\gamma}_{l},\hat{\alpha}_{l},\theta )=\varepsilon
_{ij}(\theta )\varpi _{ij}(\hat{\gamma}_{l})(Z_{i}-Z_{j})+\hat{\alpha}
_{l}(X_{i},X_{j})v_{ij}(\hat{\gamma}_{l}),
\end{equation*}
where $\hat{\alpha}_{l}(X_{i},X_{j})=\left[ \hat{\lambda}_{l}(\hat{\gamma}
_{l}(X_{i}))-\hat{\lambda}_{l}(\hat{\gamma}_{l}(X_{j}))\right] q_{ij}(\hat{
\gamma}_{l})(Z_{i}-Z_{j}),$ $\hat{\lambda}_{l}(\cdot )$ is a nonparametric estimator for $\lambda
_{0}(\cdot )$, and $v_{ij}(\hat{\gamma}_{l})=(D_{i}+D_{j}-\hat{\gamma}
_{l}(X_{i})-\hat{\gamma}_{l}(X_{j})).$ For example, $\hat{\lambda}_{l}$
could be obtained from a nonparametric regression of $Y_{i}-S_{i}^{\prime }
\hat{\theta}_{l}\ $on $\hat{\gamma}_{l}(X_{i}),$ for a preliminary estimator
$\hat{\theta}_{l}$ of $\theta _{0}\ $(e.g., a plug-in estimator or
iterations thereof). These nuisance estimators are constructed using
observations not present in the pairs in $I_{l}.$ We allow for generic
estimators $\hat{\gamma}_{l}$ and $\hat{\lambda}_{l}$ satisfying some rate
conditions in Online Appendix (ref).
The debiased pairwise difference U-estimator is $\hat{\theta}=\hat{\Sigma}
_{c,ZS}^{-1}\hat{\Sigma}_{c,ZY},$ where
\begin{eqnarray*}
\hat{\Sigma}_{c,ZY} &=&\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\varpi _{ij}(\hat{
\gamma}_{l})(Z_{i}-Z_{j})(Y_{i}-Y_{j})+\hat{\alpha}_{l}(X_{i},X_{j})v_{ij}(
\hat{\gamma}_{l}) \\
\hat{\Sigma}_{c,ZS} &=&\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\varpi _{ij}(\hat{
\gamma}_{l})(Z_{i}-Z_{j})(S_{i}-S_{j})^{\prime }.
\end{eqnarray*}
We introduce some quantities that appear in the asymptotic variance of $\hat{
\theta}.$ Define
\begin{equation*}
\xi _{i}\equiv U_{i}D_{i}\gamma _{i}\left( Z_{i}-\mu _{z}(\gamma
_{i})\right) f_{i}+\dot{\lambda}_{0}(\gamma _{i})\gamma _{i}^{2}\left(
Z_{i}-\mu _{z}(\gamma _{i})\right) v_{i},
\end{equation*}
where $\gamma _{i}=\gamma (X_{i}),$ $\mu _{z}(\gamma_{i} ):=\mathbb{E}
[D_{i}Z_{i}|\gamma _{i}]/\gamma _{i},$ $f_{i}=f_{\gamma _{0}}(\gamma
_{i})$ is the conditional density of $\gamma _{i}$ given $D_{i}=1,\ $and $
\dot{\lambda}_{0}(u)=\partial \lambda _{0}(u)/\partial u$. Define the matrix
\begin{equation*}
\Sigma _{ZS}:=\mathbb{E}[\gamma _{i}^{2}f_{i}\left( Z_{i}-\mu _{z}(\gamma
_{i})\right) \left( S_{i}-\mu _{s}(\gamma _{i})\right) ],
\end{equation*}
which is assumed to be non-singular, where $\mu _{s}(\gamma_{i} ):=\mathbb{E}
[D_{i}S_{i}|\gamma _{i}]/\gamma _{i}$. Then, we have the following
asymptotic normality result. Define $V=B^{-1}\Sigma B^{\prime -1},$ $B=2\Sigma _{ZS}$ and $\Sigma =4
\mathbb{E}[\xi _{i}\xi _{i}^{\prime }].$
\begin{theorem}
Under Assumptions (ref) to (ref) in the Online Appendix (ref),
\begin{equation*}
\sqrt{n}(\hat{\theta}-\theta _{0})\rightarrow _{d}N(0,V).
\end{equation*}
\end{theorem}
\section{Appendices}
\addcontentsline{toc}{section}{Appendices}
\subsection{Practitioner's Guide to IOp}
We provide an R package so that practitioners can easily apply our IOp estimator. The package is hosted in \href{https://github.com/joelters/ineqopp}{https://github.com/joelters/ineqopp}. At the moment of writing it includes three main functions \texttt{IOp}, \texttt{peffect} and \texttt{IOptest}. \texttt{IOp} provides IOp estimates, \texttt{peffect} provides debiased IOp partial effects where the difference of IOp with all circumstances vs IOp without some circumstances is computed and \texttt{IOptest} tests for equality of IOp between two populations. Examples of usage of these functions for an outcome vector $Y$ and a dataframe with circumstances $X$ are provided within the package documentation.
We stress that we only provide valid inference for \texttt{est_method $=$ TRUE} and \texttt{CFit $=$ TRUE}. Hence, we advise against presenting estimates, standard errors or confidence intervals based on other values of these arguments. Plug-in or debiased without cross-fitting estimators are only implemented for comparing with the proposed cross-fitted debiased estimator. Once we have chosen the debiased estimator with cross-fitting these are the choices left for the practitioner:
\begin{enumerate}
• \textbf{ML choice:} the estimator of the first-step. Our recommendation is to use an ML that non-parametrically estimates the conditional mean. One can choose from lasso, ridge, random forests, conditional inference forests, XGBoosting, Catboosting or neural networks. In practice, one could perform cross-validation to select the ML method, as mentioned in the next point.
• \textbf{Hyperparameter tuning:} the package allows setting several tuning parameters for each method. We recommend performing a cross-validation exercise to choose the parameters. This can be done for the MLs implemented in the \texttt{ineqopp} package with the function \texttt{MLtuning} in the package \texttt{ML} in \href{https://github.com/joelters/ML}{https://github.com/joelters/ML}. For the inference properties to hold, it is essential that the ML estimation is accurate; therefore, we strongly recommend conducting this process carefully. This can be done for several MLs to pick the best one.
• One can decide to get relative IOp as the one reported in the empirical application (the Gini of the predictions over the Gini of $Y$) instead of just the Gini of the predictions. This is done by setting \texttt{IOp_rel $=$ TRUE}.
\end{enumerate}
\subsection{Proofs of General Results}
\subsubsection{Main results}
Without loss of generality, we prove all results for the case in which $M = 1$. The general case follows by applying the same arguments for each $m$ in $\{1,...,M\}$.
\begin{proof}
[Proof of Lemma (ref)]Let $\mathcal{S}$ denote $L_2(X_{i},X_{j})$ in case (i) and $\mathcal{S}=\Gamma+\Gamma$ in case (ii). Set
$\alpha_{0}(X_{i},X_{j})=\Pi_{\mathcal{S}}(\delta(X_{i},X_{j},\gamma_{0}))$. Since $\alpha_{0}(\cdot,x)\in\Gamma$ and
$\alpha_{0}(x,\cdot)\in\Gamma$ for all $x\in\mathcal{X}$, as we show below, it holds by iterated expectations,
\[
\mathbb{E}_{F_{\tau}}[\alpha_{0}(X_{i},X_{j})(c_{1}(Y_{i}-\gamma_{\tau}
(X_{i}))+c_{2}(Y_{j}-\gamma_{\tau}(X_{j})))]=0.
\]
Thus, by the chain rule and because $\alpha_{0}\in\mathcal{S}$,
\begin{align*}
\frac{d}{d\tau}\mathbb{E}[g(W_{i},W_{j},\gamma(F_{\tau}),\theta)] &
=\frac{d}{d\tau}\mathbb{E}[\delta(X_{i},X_{j},\gamma_{0})(c_{1}\gamma_{\tau
}(X_{i})+c_{2}\gamma_{\tau}(X_{j}))]\\
& =\frac{d}{d\tau}\mathbb{E}[\alpha_{0}(X_{i},X_{j})(c_{1}\gamma_{\tau}
(X_{i})+c_{2}\gamma_{\tau}(X_{j}))]\\
& =\frac{d}{d\tau}\mathbb{E}_{F_{\tau}}[\alpha_{0}(X_{i},X_{j})(c_{1}
Y_{i}+c_{2}Y_{j}-c_{1}\gamma_{0}(X_{i})-c_{2}\gamma_{0}(X_{j}))]\\
& =\int\int\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)K_{H}\left(
dw_{i},dw_{j}\right)
\end{align*}
with $\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta)=\alpha_{0}(x_{i}
,x_{j})(c_{1}y_{i}+c_{2}y_{j}-c_{1}\gamma_{0}(x_{i})-c_{2}\gamma_{0}(x_{j}))$.
To complete the proof, note that for (ii), using the short notation $\delta_{ij}(\gamma)\equiv\delta(X_{i}
,X_{j},\gamma),$ and noting that by iterated expectations and independence, we
can write
\begin{align*}
\delta_{ij}(\gamma) & =\mathbb{E}[\delta_{ij}(\gamma)|X_{i}]+\mathbb{E}
[\delta_{ij}(\gamma)|X_{j}]-\mathbb{E}[\delta_{ij}(\gamma)]+U_{ij},\\
& \equiv\tilde{\delta}_{ij}(\gamma)+U_{ij},
\end{align*}
where $\mathbb{E}[U_{ij}|X_{i}]=\mathbb{E}[U_{ij}|X_{j}]=0$. This expansion
implies that only $\tilde{\delta}_{ij}$ matters for the derivative in
((ref)), so $\delta_{ij}$ could be replaced everywhere by
$\tilde{\delta}_{ij}$. Thus, substituting $\delta_{ij}$ from the last display
in the expression for $\alpha_{0}$ we obtain
\begin{align*}
\alpha_{0}(X_{i},X_{j}) & =\Pi_{\mathcal{S}}(\delta(X_{i},X_{j},\gamma
_{0}))\\
& =\Pi_{\mathcal{S}}\mathbb{E}[\delta_{ij}(\gamma_{0})|X_{i}]+\Pi
_{\mathcal{S}}\mathbb{E}[\delta_{ij}(\gamma_{0})|X_{j}]-\mathbb{E}[\delta
_{ij}(\gamma)]\\
& =\Pi_{\Gamma}\mathbb{E}[\delta_{ij}(\gamma_{0})|X_{i}]+\Pi_{\Gamma
}\mathbb{E}[\delta_{ij}(\gamma_{0})|X_{j}]-\mathbb{E}[\delta_{ij}(\gamma)],
\end{align*}
where last expression uses that $\mathcal{S}=\Gamma+\Gamma.$
\end{proof}
\subsubsection{Asymptotic theory proofs}
\begin{proof}
[Proof of Lemma (ref)]Define
\begin{align*}
\hat{R}_{1,ij,l} & =g(W_{i},W_{j},\hat{\gamma}_{l},\theta_{0})-g(W_{i}
,W_{j},\gamma_{0},\theta_{0}),\quad\hat{R}_{2,ij,l}=\phi(W_{i},W_{j}
,\hat{\gamma}_{l},\alpha_{0},\theta_{0})-\phi(W_{i},W_{j},\gamma_{0}
,\alpha_{0},\theta_{0}),\\
\hat{R}_{3,ij,l} & =\phi(W_{i},W_{j},\gamma_{0},\hat{\alpha}_{l},\theta
_{0})-\phi(W_{i},W_{j},\gamma_{0},\alpha_{0},\theta_{0}),\quad(i,j)\in I_{l}.
\end{align*}
Then
\begin{equation}
g(W_{i},W_{j},\hat{\gamma}_{l},\theta_{0})+\phi(W_{i},W_{j},\hat{\gamma}
_{l},\hat{\alpha}_{l},\theta_{0})-\psi(W_{i},W_{j},\gamma_{0},\alpha
_{0},\theta_{0})=\hat{R}_{1,ij,l}+\hat{R}_{2,ij,l}+\hat{R}_{3,ij,l}+\hat{\xi
}_{l}(W_{i},W_{j}).
\end{equation}
Let $\hat{A}_{ij,l} = \hat{R}_{1,ij,l}+\hat{R}_{2,ij,l}+\hat{R}_{3,ij,l}$. Let $N_{l}$ be the observations not in $I_{l}$,\footnote{$N_l$ is a set of observations not of pairs, that is why we do not use the complement of $I_l$, e.g. $n = 4$, $I_l = \{(1,2),(1,3)\}$ and $N_l = \{4\}$.} then we can rewrite
\[
\hat{A}_{ij,l} + \hat{\xi
}_{l}(W_{i},W_{j}) = \left(\hat{A}_{ij,l} - \mathbb{E}[\hat{A}_{ij,l} | N_l]\right) + \mathbb{E}[\hat{A}_{ij,l} | N_l] + \hat{\xi
}_{l}(W_{i},W_{j}).
\]
Since pairs in $I_{l}$ are dependent only when
one or two of the members of the pair coincide (also we omit the fact that we
are dealing with vectors since the convergence of the vector is the
convergence of its elements) we have that for $b = 1,2,3$ and $|A|$ the cardinality of set $A$
\begin{align*}
& \mathbb{E}\biggl[\biggl(\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}
(\hat{R}_{b,ij,l}-\mathbb{E}(\hat{R}_{b,ij,l}|N_{l})\biggr)^{2}
\biggr|N_{l}\biggr]=\\
& n\binom{n}{2}^{-2}\biggl[|I_l|\mathbb{V}ar(\hat{R}_{b,ij,l}
|N_{l})+|I_l|\cdot 2 (p_l - 1) \mathbb{C}ov(\hat{R}_{b,ij,l},\hat{R}_{b,ik,l}
|N_{l})\biggr],
\end{align*}
where $p_l$ is the number of distinct observations that occur as either component of a pair in $I_l$. $n \binom{n}{2}^{-2} |I_l| \to 0$ and by Lemma (ref) in the Online Appendix (ref), $|I_l|\cdot 2 (p_l - 1) \to c < \infty$. Hence
\begin{align*}
& \mathbb{E}\biggl[\biggl(\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}
(\hat{R}_{b,ij,l}-\mathbb{E}(\hat{R}_{b,ij,l}|N_{l})\biggr)^{2}
\biggr|N_{l}\biggr]\\
& \leq c \cdot \mathbb{C}ov(\hat{R}_{b,ij,l},\hat{R}_{b,ik,l}|N_{l})+o_{P}(1)\\
& \leq c \cdot\sqrt{\mathbb{E}(\hat{R}_{b,ij,l}^{2}|N_{l})\mathbb{E}(\hat
{R}_{b,ik,l}^{2}|N_{l})}+o_{P}(1)\rightarrow_{p}0,
\end{align*}
where convergence for $b = 1,2,3$ follows from Assumption
(ref). Then, by the conditional Markov inequality, triangle inequality, and Dominated Convergence Theorem (DCT)
\[
\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}\left(\hat{A}_{ij,l} - \mathbb{E}[\hat{A}_{ij,l} | N_l]\right)\rightarrow_{p}0.
\]
Note that $\mathbb{E}[\hat{R}_{1,ij,l}+\hat
{R}_{2,ij,l}|N_{l}]=\mathbb{E}[\psi(W_{i},W_{j},\hat{\gamma}_{l}
,\alpha_{0},\theta_{0}) | N_l]$. Therefore, by Assumption
(ref)
\[
\biggl|\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}\mathbb{E}(\hat
{R}_{1,ij,l}+\hat{R}_{2,ij,l}|N_{l})\biggr|\leq 2\sqrt{n}|\bar{\psi}(\hat{\gamma}
_{l},\alpha_{0},\theta_{0})|\rightarrow_{p}0.
\]
Hence, by the triangle inequality
\[
\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}(\hat{R}_{1,ij,l}+\hat
{R}_{2,ij,l}+\hat{R}_{3,ij,l}) = \sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}} \left(\hat{A}_{ij,l} - \mathbb{E}[\hat{A}_{ij,l} | N_l]\right) + \mathbb{E}[\hat{A}_{ij,l} | N_l] \rightarrow_{p}0.
\]
Hence, by Assumption (ref) and the triangle inequality
\begin{align*}
& \sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}\psi(W_{i},W_{j},\hat{\gamma
}_{l},\hat{\alpha}_{l},\theta_{0})-\psi(W_{i},W_{j},\gamma_{0},\alpha
_{0},\theta_{0})\\
& =\sqrt{n}\binom{n}{2}^{-1}\sum_{(i,j)\in I_{l}}(\hat{R}_{1,ij,l}+\hat{R}_{2,ij,l}+\hat{R}
_{3,ij,l}+\hat{\xi}_{l}(W_{i},W_{j}))\rightarrow_{p}0.
\end{align*}
The same remains true when we sum across finite folds.
\end{proof}
\begin{proof}
[Proof of Lemma (ref)] Let $\hat{\psi}_{ij} = \psi(W_i,W_j,\hat{\gamma},\hat{\alpha},\hat{\theta})$, $\psi_{ij} = \psi(W_i,W_j,\gamma_0,\alpha_0,\theta_0)$ and define $\hat{\chi}_i = (n-1)^{-1}\sum_{j \neq i} \hat{\psi}_{ij}$, $\tilde{\chi}_i = (n-1)^{-1}\sum_{j \neq i} \psi_{ij}$ and $\chi_i = \mathbb{E}[\psi_{ij}|W_i]$. Without loss of generality, we focus on the scalar case. Note that $\hat{\Sigma} = 4n^{-1}\sum_{i=1}^n \hat{\chi}_i^2$ and we want to show that $\hat{\Sigma} \to_p 4\mathbb{E}[\chi_i^2] = \Sigma$. To do this note that
\[
\frac{1}{n}\sum_{i=1}^n(\hat{\chi}_i - \tilde{\chi}_i + \tilde{\chi´}_i)^2 = \frac{1}{n}\sum_{i=1}^n \tilde{\chi}_i^2 + \frac{1}{n}\sum_{i=1}^n (\hat{\chi}_i - \tilde{\chi}_i)^2 + \frac{2}{n}\sum_{i=1}^n (\hat{\chi}_i - \tilde{\chi}_i)\tilde{\chi}_i.
\]
The first term in the right hand side goes in probability to $\mathbb{E}[\chi_i^2]$ by the strong law of large numbers for U-statistics (see p.190 in serfling1980approximation). By Cauchy-Schwartz, the third term can be bounded as
\[
\frac{2}{n}\sum_{i=1}^n (\hat{\chi}_i - \tilde{\chi}_i)\tilde{\chi}_i \leq 2\biggl(\frac{1}{n}\sum_{i=1}^n (\hat{\chi}_i - \tilde{\chi}_i)^2 \biggr)^{\frac{1}{2}}\biggl(\frac{1}{n}\sum_{i=1}^n \tilde{\chi}_i^2 \biggr)^{\frac{1}{2}}.
\]
Hence, if we show that $n^{-1}\sum_{i=1}^n |\hat{\chi}_i - \tilde{\chi}_i|^2 \to_p 0$ we have that $n^{-1}\sum_{i=1}^n \hat{\chi}_i^2 = \mathbb{E}[\chi_i^2] + o_p(1)$. It follows by the $C_r$ inequality and the assumptions in the lemma that
\begin{align*}
\frac{1}{n}\sum_{i=1}^n |\hat{\chi}_i - \tilde{\chi}_i|^2 &= \frac{1}{n}\sum_{i=1}^n \biggl|\frac{1}{n-1} \sum_{j\neq i}[\hat{g}_{ij} - g_{ij} + \hat{\phi}_{ij} - \phi_{ij}]\biggr|^2 \\
&\leq \frac{2}{n}\sum_{i=1}^n |\hat{g}_{-i} - g_{-i}|^2 + \frac{2}{n}\sum_{i=1}^n |\hat{\phi}_{-i} - \phi_{-i}|^2 \to_p 0.
\end{align*}
\end{proof}
\begin{proof}
[Proof of Lemma (ref)]Define
\[
\hat{B} = \binom{n}{2}^{-1}\sum_{i < j}\partial\psi(W_{i}
,W_{j},\hat{\gamma},\hat{\alpha},\bar{\theta})/\partial\theta, \quad \tilde{B}
=\binom{n}{2}^{-1}\sum_{i<j}\partial\psi(W_{i},W_{j}
,\hat{\gamma},\hat{\alpha},\theta_{0})/\partial\theta.
\]
By ii), with probability approaching 1
\[
\mathbb{E}\biggl[\binom{n}{2}^{-1}\sum_{i<j}d(W_{i},W_{j}
,\hat{\gamma},\hat{\alpha})\biggr]=\mathbb{E}[d(W_{i},W_{j},\hat
{\gamma}, \hat{\alpha})]\leq C,
\]
By Markov inequality , $\binom{n}{2}^{-1}\sum_{i<j}d(W_{i}
,W_{j},\hat{\gamma},\hat{\alpha})=O_{p}(1)$. So by ii) and $\bar{\theta} \to_p \theta_0$
\begin{align*}
\biggl|\hat{B}-\tilde{B}\biggr| & \leq\binom{n}{2}^{-1}
\sum_{i<j}\biggl|\frac{\partial\psi(W_{i},W_{j},\hat{\gamma},\hat{\alpha}
,\bar{\theta})}{\partial\theta}-\frac{\partial\psi(W_{i},W_{j},\hat{\gamma}
,\hat{\alpha},\theta_{0})}{\partial\theta}\biggr|\\
& \leq\binom{n}{2}^{-1}\sum_{i<j}d(W_{i},W_{j},\hat{\gamma},\hat{\alpha}
)|\bar{\theta}-\theta_{0}|^{1/C} \\
&= O_p(1)o_p(1) \to_p 0.
\end{align*}
It follows from Assumption (ref) (iii) and Markov Inequality that
\[
\biggl|\tilde{B}-\binom{n}{2}^{-1}\sum_{i<j
}\partial\psi(W_{i},W_{j},\gamma_{0},\alpha_0,\theta_{0})/\partial\theta \biggr|\rightarrow
_{p}0
\]
and $\binom{n}{2}^{-1}\sum_{i<j}\partial\psi(W_{i},W_{j},\gamma
_{0},\alpha_0,\theta_{0})/\partial\theta\rightarrow_{p}B$ by the strong law of large numbers for U-statistics (see p.190 in serfling1980approximation). Hence, the result follows
from the triangle inequality.
\end{proof}
\begin{proof}
[Proof of Theorem (ref)]By the Mean Value Theorem, ((ref)) and consistency of the Jacobian term, (with $\tilde{\psi}(\theta _{0})$ defined as $\hat{\psi}(\theta _{0})$ but with true nuisance parameters)
\begin{eqnarray*}
0 &=&\sqrt{n}\hat{\psi}(\theta _{0})-\frac{\partial\hat{\psi}(\bar{\theta})}{\partial\theta} \sqrt{n}(\hat{\theta}-\theta _{0})
\\
&=&\sqrt{n}\tilde{\psi}(\theta _{0})-B\sqrt{n}(\hat{\theta}-\theta
_{0})+o_{P}(1).
\end{eqnarray*}
The normality result follows from standard U-statistic theory (see Theorem
12.3 in van2000asymptotic) and Slutsky's lemma. Hence $\sqrt{n}
(\hat{\theta} - \theta_{0}) \rightarrow_{d} \mathcal{N}(0,V)$. Consistency of the variance estimator follows from Lemmas (ref) and (ref) and standard arguments.
\end{proof}
\subsubsection{Consistency}
Re-define the interaction term in Section (ref) as
\[
\hat{\xi}_{l}(w_{i},w_{j},\theta)=\phi(w_{i},w_{j},\hat{\gamma}_{l},\hat{\alpha}_{l},\theta)-\phi(w_{i},w_{j},\gamma_{0},\hat{\alpha}_{l},\theta)-\phi(w_{i},w_{j},\hat{\gamma}_{l},\alpha_{0},\theta_0)+\phi(w_{i},w_{j},\gamma_{0},\alpha_{0},\theta_0).
\]
\begin{theorem}
If (i) $\mathbb{E}[g(W_i,W_j,\gamma_0,\theta)] = 0 \text{ iff } \theta = \theta_0$, (ii) $\Theta$ is compact, (iii) $\int \int |g(w_i,w_j,\hat{\gamma}_l,\theta) - g(w_i,w_j,\gamma_0,\theta)| F_0(dw_i) F_0(dw_j) \to_p 0$ and $\mathbb{E}[|g(W_i,W_j,\gamma_0,\theta)] < \infty$ for all $\theta \in \Theta$, (iv) there is $C > 0$ and $d(W_i,W_j, \gamma)$ such that for $||\gamma - \gamma_0||$ small enough and all $\tilde{\theta}, \theta \in \Theta$
\[
|g(W_i,W_j,\gamma,\tilde{\theta}) - g(W_i,W_j,\gamma,\theta)| \leq d(W_i,W_j,\gamma)|\tilde{\theta} - \theta|^{1/C},
\]
and $\mathbb{E}[d(W_i,W_j,\gamma)] < C$ and (v) Assumption (ref) (ii), (iii), $\int \int |\hat{\xi}_l(w_i,w_j,\theta)| F_0(dw_i) F_0(dw_j) \to_p 0$ and $\mathbb{E}[|\phi(W_i,W_j,\gamma_0,\alpha_0,\theta_0|] < \infty$. Then $\hat{\theta} \to_p \theta_0$.
\end{theorem}
\begin{proof}
Define $\hat{g}(\theta) \equiv \binom{n}{2}^{-1} \sum_{l=1}^L \sum_{(i,j)\in I_l} g(W_i,W_j,\hat{\gamma}_l, \theta)$ and $\bar{g}(\theta) \equiv \mathbb{E}[g(W_i,W_j,\gamma_0,\theta)]$. It follows from the conditional Markov inequality and (iii) that $\hat{g}(\theta) \to_p \bar{g}(\theta)$ for all $\theta \in \Theta$. Let $\tilde{\phi}_{ij} \equiv \binom{n}{2}^{-1} \sum_{i<j} \phi(W_i,W_j,\gamma_0,\alpha_0,\theta_0)$ and $\hat{\phi}_{ij}(\theta) \equiv \binom{n}{2}^{-1} \sum_{l=1}^L \sum_{(i,j) \in I_l} \phi(W_i,W_j,\hat{\gamma}_l,\hat{\alpha}_l,\theta)$. In the notation of Lemma (ref), $\hat{\phi}_{ij}(\theta) - \tilde{\phi}_{ij} = \binom{n}{2}^{-1} \sum_{l=1}^L \sum_{(i,j) \in I_l} \hat{R}_{2,ij,l} + \hat{R}_{3,ij,l} + \hat{\xi}_l(W_i,W_j,\theta)$, so $\hat{\phi}_{ij}(\theta) - \tilde{\phi}_{ij} \to_p 0$ for all $\theta \in \Theta$ by Assumption (v) and the conditional Markov inequality. By consistency of U-statistics we have that $\tilde{\phi}_{ij} \to_p \mathbb{E}[\phi(W_i,W_j,\gamma_0,\alpha_0,\theta_0)] = 0$ so by the triangle inequality we have that $\hat{\phi}_{ij}(\theta) \to_p 0$ for all $\theta \in \Theta$. Therefore, defining $\hat{\psi}_{ij}(\theta) \equiv \binom{n}{2}^{-1} \sum_{l=1}^L \sum_{(i,j) \in I_l}\psi(W_i,W_j,\hat{\gamma}_l,\hat{\alpha}_l,\theta)$, we have that $\hat{\psi}_{ij,l}(\theta) = \hat{g}(\theta) + o_p(1)$.
By triangle inequality and (iv) we know that with probability approaching one
\[
|\hat{g}(\hat{\theta}) - \hat{g}(\theta)| \leq \binom{n}{2}^{-1} \sum_{i<j} |g(W_i,W_j,\hat{\gamma}_l,\hat{\theta}) - g(W_i,W_j,\hat{\gamma}_l,\theta)| \leq \underbrace{\binom{n}{2}^{-1} \sum_{i<j} d(W_i,W_j,\hat{\gamma}_l)}_{\equiv \hat{M}_l}|\hat{\theta} - \theta|^{1/C},
\]
and by (iv) and the conditional Markov inequality $\hat{M}_l = O_p(1)$. By Corollary 2.2 in newey1991uniform we have that $\sup_{\theta \in \Theta} |\hat{\psi}(\theta) - \bar{g}(\theta)| = o_p(1)$. We also know that $\bar{g}(\theta)$ is continuous by (iv). So the conclusion follows from the proof of Theorem 2.6 in newey1994large applied to the Háyek projection of $\hat{g}(\hat{\theta})$.
\end{proof}
\putbib[references]
\pagenumbering{arabic}
\setcounter{page}{1}
center[center omitted — 49 chars of source]
bibunit[ectabib]
\addcontentsline{toc}{section}{Appendices}
Proofs of Inequality of Opportunity
proof[Proof of Proposition (ref)] We use the short notation $\Delta _{\tau }\equiv \gamma _{\tau
}(X_{i})-\gamma _{\tau }(X_{j})$ and, in particular, $\Delta _{0}\equiv
\gamma _{0}(X_{i})-\gamma _{0}(X_{j}).$ By $|\Delta _{\tau }|=sgn(\Delta
_{\tau })\Delta _{\tau }$ and the chain rule,
\begin{eqnarray*}
\frac{d}{d\tau }\mathbb{E}(g(W_{i},W_{j},\gamma (F_{\tau }),\theta )) &=&
\frac{d}{d\tau }\mathbb{E}[\theta (\gamma _{\tau }(X_{i})+\gamma _{\tau
}(X_{j}))]-\frac{d}{d\tau }\mathbb{E}(|\Delta _{\tau }|) \\
&=&\frac{d}{d\tau }\mathbb{E}[\theta (\gamma _{\tau }(X_{i})+\gamma _{\tau
}(X_{j}))]-\frac{d}{d\tau }\mathbb{E}(sgn(\Delta _{0})\Delta _{\tau })-\frac{
d}{d\tau }\mathbb{E}(sgn(\Delta _{\tau })\Delta _{0}).
\end{eqnarray*}
The first two summands are already in the form of ((ref)), so
we can directly apply Lemma (ref) to find that these two terms have
adjustment terms
\begin{eqnarray*}
\phi _{1}(w_{i},w_{j},\gamma _{0},\theta ) &=&\theta (y_{i}+y_{j}-\gamma
_{0}(x_{i})-\gamma _{0}(x_{j})) \\
\phi _{2}(w_{i},w_{j},\gamma _{0},\theta ) &=&-sgn(\Delta
_{0})(y_{i}-y_{j}-\gamma _{0}(x_{i})+\gamma _{0}(x_{j})).
\end{eqnarray*}
Now, we proceed to show that
\[
\frac{d}{d\tau }\mathbb{E}(sgn(\Delta _{\tau })\Delta _{0})=0.
\]
By Assumption (ref),
\[
|sgn(\Delta _{\tau })-sgn(\Delta _{0})|\leq 2\cdot 1(|\Delta _{0}|\leq \tau
V_{H}),
\]
where $V_{H}\equiv V_{H}(X_{i},X_{j})=C_{H}(X_{i})+C_{H}(X_{j}).$ Thus,
\begin{eqnarray*}
\left\vert \mathbb{E}\left[ \Delta _{0}\left\{ \frac{sgn(\Delta _{\tau
})-sgn(\Delta _{0})}{\tau }\right\} \right] \right\vert &\leq &2\mathbb{E}
\left[ \left\vert \Delta _{0}\right\vert \frac{1(|\Delta _{0}|\leq \tau
V_{H})}{\tau }\right] \\
&\leq &2\mathbb{E}\left[ V_{H}1(0<|\Delta _{0}|\leq \tau V_{H})\right] \\
&\rightarrow &0 as \tau \rightarrow 0,
\end{eqnarray*}
where the convergence follows by continuity of the $\sigma $-finite positive
measure
\[
\mu (A)=\int_{A}V_{H}(x_{1},x_{2})F_{0}(dx_{1})F_{0}(dx_{2}),
\]
$\mathbb{E}\left[ V_{H}\right] <\infty $ and the fact that the sets $
A_{n}=\{(x_{1},x_{2})\in \mathcal{X}\times \mathcal{X}:0<|\Delta
_{0}(x_{1},x_{2})|\leq \tau _{n}V_{H}(x_{1},x_{2})\}$ are monotonically
decreasing for a sequence $\tau _{n}\downarrow 0$ as $n\rightarrow \infty ,$
and $\lim_{n\rightarrow \infty }A_{n}=\emptyset .$ Hence, $\phi =\phi
_{1}+\phi _{2}.$
\bigbreak
proof[Proof of Proposition (ref)]
We start with the smoothness part and leave the consistency of the sign for the end of the proof. Since $|sgn(\Delta_{\gamma}) - sgn(\Delta_0)| \leq 2 \cdot 1_A$ a.s. where $A = \{(x_i,x_j): sgn(\gamma(x_i) - \gamma(x_j)) \neq sgn(\gamma_0(x_i) - \gamma_0(x_j)) \}$, we have that
\begin{align*}
|\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0] | \leq 2 \mathbb{E}[|\Delta_0| 1_A].
\end{align*}
For some $t > 0$, we can decompose into a small margin and a large margin term
\[
\mathbb{E}[|\Delta_0| 1_A 1(0 < |\Delta_0| \leq t)] + \mathbb{E}[|\Delta_0| 1_A 1(|\Delta_0| > t)].
\]
Let us now deal with the small margin term. We restrict the analysis to the set $A$ where $\Delta_{\gamma} \cdot \Delta_0 < 0$ or $\Delta_{\gamma} \cdot \Delta_0 = 0$.
If $\Delta_0 > 0$, then $\Delta_{\gamma} < 0$ meaning that
\begin{align*}
\Delta_{\gamma} &= \Delta_0 - (\gamma_0(X_i) - \gamma(X_i)) + \gamma_0(X_j) - \gamma(X_j) < 0\\
\implies &- (\gamma_0(X_i) - \gamma(X_i)) + \gamma_0(X_j) - \gamma(X_j) < - \Delta_0 \\
\implies &|\Delta_0| < |\gamma_0(X_i) - \gamma(X_i)| + |\gamma_0(X_j) - \gamma(X_j)|.
\end{align*}
When $\Delta_0 < 0$, then $\Delta_{\gamma} > 0$ so
\begin{align*}
&- (\gamma_0(X_i) - \gamma(X_i)) + \gamma_0(X_j) - \gamma(X_j) > - \Delta_0 \\
\implies &|-(\gamma_0(X_i) - \gamma(X_i)) + \gamma_0(X_j) - \gamma(X_j)| > - \Delta_0 = |\Delta_0| \\
\implies &|\Delta_0| < |\gamma_0(X_i) - \gamma(X_i)| + |\gamma_0(X_j) - \gamma(X_j)|.
\end{align*}
If $\Delta_0 = 0$ then $\mathbb{E}[|\Delta_0| 1_A] = 0$ and if $\Delta_\gamma = 0$
\[
|\Delta_0| = |\Delta_0 - \Delta_\gamma| = | \gamma_0(X_i) - \gamma(X_i) - (\gamma_0(X_j) - \gamma(X_j)) | \leq | \gamma_0(X_i) - \gamma(X_i)| + | \gamma_0(X_j) - \gamma(X_j) |.
\]
Hence,
\begin{align*}
\mathbb{E}[|\Delta_0| 1_A 1( 0 < |\Delta_0| \leq t)] &\leq 2\mathbb{E}[|\gamma(X_i) - \gamma_0(X_i)| 1(0<|\Delta_0| \leq t)].
\end{align*}
Now, using Hölder's inequality with $q' = q/(q-1)$ we get
\begin{align*}
\mathbb{E}[|\Delta_0| 1_A 1(0 < |\Delta_0| \leq t)] &\leq 2 ||\gamma - \gamma_0||_q P(0 < |\Delta_0| \leq t)^{1/q'},
\end{align*}
so by definition of $q'$ and Assumption (ref) (iii) we get
\[
\mathbb{E}[|\Delta_0| 1_A 1(0 < |\Delta_0| \leq t)] \leq 2C ||\gamma - \gamma_0||_q t^{\frac{\beta(q-1)}{q}}.
\]
For the large margin, let $B_i = |\gamma_0(X_i) - \gamma(X_i)|$ for notational brevity. Using again that on $A$: $|\Delta_0| < |\gamma_0(X_i) - \gamma(X_i)| + |\gamma_0(X_j) - \gamma(X_j)|$ (also inside the indicator) we show that
\begin{align*}
\mathbb{E}[|\Delta_0| 1_A 1(|\Delta_0| > t)] &\leq \mathbb{E}[(B_i + B_j) 1(B_i + B_j > t)] \\
&\leq 2 \mathbb{E}[B_i 1(B_i > t/2)] + 2 \mathbb{E}[B_i 1(B_j > t/2)] \\
&\leq 4C_{\Delta} ||\gamma - \gamma_0||_q P(B_i > t/2)^{\frac{q-1}{q}} \\
&\leq 4 \cdot 2^{q-1}\cdot C_{\Delta} \frac{||\gamma - \gamma_0||_q^q}{t^{q-1}},
\end{align*}
where we have used that $a+b > t$ implies $\{a > t/2\} \cup \{b > t/2\}$ and monotonicity of measures, then Hölder and Markov inequalities. Putting all together (for $q \geq 1$) we have
\begin{align*}
|\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0] | &\leq 4 \cdot 2^{q-1}\cdot C \left( ||\gamma - \gamma_0||_q t^{\frac{\beta (q-1)}{q}} + \frac{||\gamma - \gamma_0||_q^q}{t^{q-1}} \right).
\end{align*}
Optimizing $t$ and plugging in the optimal, we get
\[
|\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0] \leq C(\beta,q) ||\gamma - \gamma_0||_q^{\frac{q(1+\beta)}{q + \beta}}.
\]
Now, for the supremum norm, notice that by the derivations above, we have that $A$ implies
\begin{align*}
|\Delta_0| \leq 2 ||\gamma - \gamma_0||_\infty,
\end{align*}
and hence also $1_A \leq 1(|\Delta_0| \leq 2 ||\gamma - \gamma_0||_\infty)$ a.s. Then
\begin{align*}
|\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0] &= |\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0 1(|\Delta_0| > 0)]| \\
&\leq 2 \mathbb{E}[|\Delta_0| 1_A 1(0<|\Delta_0|)] \\
&\leq 2 ||\gamma - \gamma_0||_\infty P(0 < |\Delta_0| \leq 2 ||\gamma- \gamma_0||_\infty).
\end{align*}
Hence, by Assumption (ref) (iii)
\[
|\mathbb{E}[(sgn(\Delta_{\gamma}) - sgn(\Delta_0))\Delta_0] \leq C ||\gamma - \gamma_0||_\infty^{1 + \beta}.
\]
For the consistency of the sign difference, note that
\begin{align*}
\mathbb{E}[|sgn(\Delta_{\gamma}) - sgn(\Delta_0)|] &= 2\mathbb{E}[1_A 1(X_i \neq X_j)] \\
&= 2 \mathbb{E}[1_A 1(X_i \neq X_j) 1(|\Delta_0| \leq t)] + 2 \mathbb{E}[1_A 1(|\Delta_0| > t)1(X_i \neq X_j)].
\end{align*}
Let now $t \equiv t_n \downarrow0$, then the first term goes to zero by Assumption (ref) (ii). For the second term, we use the same arguments as above to show that
\begin{align*}
\mathbb{E}[1_A 1(|\Delta_0| > t_n)1(X_i \neq X_j)] &\leq \mathbb{E}[1(B_i + B_j > t_n)] \\
&\leq 4 \frac{||\gamma - \gamma_0||_1}{t_n}.
\end{align*}
Hence, using $\hat{\gamma}_l$ instead of $\gamma$ and conditioning on $N_{l}$, the consistency of the sign difference follows by selecting a sequence $t_n \downarrow0$ sufficiently slow.
proof[Proof of Proposition (ref)]
We first want to show that
\begin{align*}
\sqrt{n} \binom{n}{2}^{-1} &\sum_{l=1}^L \sum_{(i,j) \in I_l} \left( \theta_0 (Y_i + Y_j) - sgn(\hat{\gamma}_l(X_i) - \hat{\gamma}_l(X_j))(Y_i-Y_j) \right) = \\
\sqrt{n} \binom{n}{2}^{-1} &\sum_{i<j} \left( \theta_0 (Y_i + Y_j) - sgn(\gamma_0(X_i) - \gamma_0(X_j))(Y_i - Y_j)\right) + o_p(1).
\end{align*}
Define
\[
\hat{A}_{ij,l} = \left[sgn(\gamma_0(X_i) - \gamma_0(X_j)) - sgn(\hat{\gamma}_l(X_i) - \hat{\gamma}_l(X_j))\right](Y_i-Y_j)
\]
Then, adding and subtracting $\mathbb{E}[\hat{A}_{ij,l} | N_{l}]$ and focusing on a single fold $l$ the above boils down to
\begin{align*}
\sqrt{n} \binom{n}{2}^{-1} \sum_{(i,j) \in I_l} \left[ \left( \hat{A}_{ij,l} - \mathbb{E}[\hat{A}_{ij,l} | N_{l}] \right) + \mathbb{E}[\hat{A}_{ij,l} | N_{l}] \right] = o_p(1).
\end{align*}
Following the same arguments as in the general case
\begin{align*}
P\left( \sqrt{n} \binom{n}{2}^{-1} \left|\sum_{(i,j) \in I_l} \left( \hat{A}_{ij,l} - \mathbb{E}[\hat{A}_{ij,l} | N_{l}] \right) \right| > \varepsilon \biggl| N_{l} \right) &\leq 4 \sqrt{\mathbb{E}[\hat{A}_{ij,l}^2 | N_{l}] \mathbb{E}[\hat{A}_{ik,l}^2 | N_{l}]/\varepsilon^2},
\end{align*}
and if conditional second moments of $Y_i$ given $X_i$ are bounded
\begin{align*}
\mathbb{E}[\hat{A}_{ij,l}^2 | N_{l}] &= \mathbb{E}\left( \left[sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l})\right]^2(Y_i-Y_j)^2 | N_{l} \right) \\
&\leq C \mathbb{E}\left( |sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l})| \, | N_{l} \right) \to_p 0,
\end{align*}
or if $2+\delta$ moments of $Y_i - Y_j$ exist
\begin{align*}
\mathbb{E}[\hat{A}_{ij,l}^2 | N_{l}] &= \mathbb{E}\left( \left[sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l})\right]^2(Y_i-Y_j)^2 | N_{l} \right) \\
&\leq C \mathbb{E}\left( |sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l})| \, | N_{l} \right)^{\frac{\delta}{2+\delta}} \mathbb{E}[|Y_i - Y_j|^{2+\delta}] \to_p 0,
\end{align*}
The convergences in probability follow from Assumption (ref) (ii) by Proposition (ref). Now we need to show that
\begin{align*}
\sqrt{n} \binom{n}{2}^{-1} \sum_{(i,j) \in I_l}\mathbb{E}[\hat{A}_{ij,l} | N_{l}] \to_p 0.
\end{align*}
By the same arguments as in the proof of the general case
\begin{align*}
\biggl| \sqrt{n} \binom{n}{2}^{-1} \sum_{(i,j) \in I_l}\mathbb{E}[\hat{A}_{ij,l} | N_{l}] \biggr| &\leq 2\sqrt{n} |\mathbb{E}[(sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l}))\Delta_0 | N_{l}] |.
\end{align*}
By Proposition (ref)
\begin{align*}
\sqrt{n} |\mathbb{E}[(sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l}))\Delta_0 | N_{l}] | &\leq \begin{cases}
\sqrt{n}C(\beta,q) ||\hat{\gamma}_l - \gamma_0||_q^{\frac{q(1+\beta)}{q+ \beta}} & if q \in [1,\infty),\\
\sqrt{n}2C_{\Delta} ||\hat{\gamma}_l - \gamma_0||_\infty^{1+\beta} & if q = \infty.
\end{cases}
\end{align*}
Under Assumption (ref), for any choice of $q$ we get the result that
\begin{align*}
\sqrt{n} |\mathbb{E}[(sgn(\Delta_0) - sgn(\Delta_{\hat{\gamma}_l}))\Delta_0 | N_{l}] | &= o_p(1).
\end{align*}
For the conditions for the consistency of the variance, we first show that
\begin{align*}
\hat{\Sigma} - \Sigma &= \frac{1}{n}\sum_{i=1}^n\left[ \left(\underbrace{\sum_{j\neq i} (Y_i + Y_j)\hat{\theta} - sgn(\Delta_{\hat{\gamma}})(Y_i-Y_j)}_{\equiv \hat{B}_{i}} \right)^2
- \left(\underbrace{\sum_{j\neq i} (Y_i + Y_j)\theta_0 - sgn(\Delta_0)(Y_i-Y_j)}_{\equiv B_{i}} \right)^2 \right] \\
&\leq \frac{1}{n}\sum_{i=1}^n (\hat{B}_{i} - B_{i})^2 + B_{i}(\hat{B}_{i} - B_{i}).
\end{align*}
For the first term, we have that
\begin{align*}
(\hat{B_i} - B_i)^2 \leq (\hat{\theta} - \theta_0)^2\left(\frac{1}{n-1} \sum_{j \neq i} (Y_i + Y_j) \right)^2 +
\left(\frac{1}{n-1} \sum_{j \neq i} (sgn(\hat{\Delta}) - sgn(\Delta_0))(Y_i - Y_j) \right)^2,
\end{align*}
doing some algebraic manipulations
\begin{align*}
\left(\frac{1}{n-1} \sum_{j \neq i} (Y_i + Y_j) \right)^2 &\leq C(Y_i^2 + \bar{Y}_n^2),
\end{align*}
where $\bar{Y}_n = (1/n)\sum_{i=1}^n Y_i$. Hence,
\begin{align*}
\frac{(\hat{\theta} - \theta_0)^2}{n}\sum_{i=1}^n \left(\frac{1}{n-1} \sum_{j \neq i} (Y_i + Y_j) \right)^2 &\leq C(\hat{\theta} - \theta_0)^2\left(\bar{Y}_n + \frac{1}{n} \sum_{i=1}^n Y_i^2 \right) \to_p 0.
\end{align*}
On the other hand
\begin{align*}
&P\left(\frac{1}{n} \left| \sum_{i=1}^n \left( \sum_{j \neq i} (sgn(\hat{\Delta}) - sgn(\Delta_0))(Y_i - Y_j) \right)^2 \right| > \varepsilon \right)\\
&\leq \frac{1}{\varepsilon(n-1)^2}\mathbb{E} \left[ \left| \sum_{j \neq i} (sgn(\hat{\Delta}) - sgn(\Delta_0))(Y_i - Y_j) \right|^2 \right] \\
&= \frac{1}{\varepsilon(n-1)^2}\left((n-1)\mathbb{E} \left[ \hat{A}_{ij}^2 \right] +
(n-1)^2 \mathbb{E} \left[ \hat{A}_{ij} \hat{A}_{ik} \right] \right),
\end{align*}
where $\hat{A}_{ij} \equiv (sgn(\hat{\gamma}(X_i) - \hat{\gamma}(X_j)) - sgn(\gamma_0(X_i) - \gamma_0(X_j)))(Y_i - Y_j)$. The first term goes to zero since its divided by $n-1$ and by the boundedness of the sign and boundedness of the conditional variance of $Y_i$. The second term goes to zero by applying Cauchy Schwarz and noting that by the finiteness of the conditional variance of $Y_i$
\[
\sqrt{\mathbb{E}[\hat{A}_{ij}^2]} \leq C \sqrt{\mathbb{E}\left( |sgn(\hat{\Delta} - \Delta_0)| \right)} \to 0,
\]
by the fact that Proposition (ref) goes through using the exact same arguments unconditionally for any fold under Assumption (ref) and Assumption (ref) (ii). Now,
\begin{align*}
\left| \frac{1}{n} \sum_{i=1}^n B_i(\hat{B}_i - B_i) \right| \leq \sqrt{\frac{1}{n} \sum_{i=1}^n B_i^2} \sqrt{\frac{1}{n} \sum_{i=1}^n (\hat{B}_i - B_i)^2} = O_p(1)o_p(1) = o_p(1).
\end{align*}
Finally, by the continuous mapping theorem $\hat{\Sigma}/\bar{Y}^2_n \to_p \Sigma/\mathbb{E}[Y_i]^2$ so $\hat{V} \to_p V$. Now we need to show the consistency of $\hat{\theta}$. For that note that
\[
\hat{m} \equiv \binom{n}{2}^{-1} \sum_{l = 1}^L \sum_{(i,j) \in I_{l}} sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j) = \binom{n}{2}^{-1} \sum_{l = 1}^L |I_l| \underbrace{\frac{1}{|I_l|} \sum_{(i,j) \in I_{l}} sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j)}_{\equiv \hat{m}_{l}},
\]
where $|I_l|$ is the cardinality of $I_l$. We show first consistency of $\hat{m}_l$ for $l \in \{1,...,L\}$ to $m_0 = \mathbb{E}[sgn(\Delta_0)(Y_i - Y_j)]$. We have that
\[
\hat{m}_l - m_0 = \hat{m}_l - \mathbb{E}[sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j) | N_{l}] + \mathbb{E}[sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j) | N_{l}] - m_0.
\]
The second difference goes to zero in probability by Proposition (ref) with $q = 1$, invoking Assumption (ref) (i). For the first difference, define $p_l$ as the amount of unique observations in $I_l$, for example if $I_l = \{(1,8),(1,9), (2,8), (2,9) \}$ then $p_l = 4$ (i.e. $1,2,8,9$). By conditional Chebyshev and the fact that observations are i.i.d conditional on $N_{l}$ we have that
\begin{align*}
P&\left( \left| |I_l|^{-1} \sum_{(i,j) \in I_{l}} \left(sgn(\Delta_{\hat{\gamma}_l}(Y_i - Y_j) - \mathbb{E}[sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j)] \right) \right| > \varepsilon \biggl| N_{l} \right) \\
&\leq (|I_l| \varepsilon)^{-2} \biggl( |I_l| \mathbb{V}ar(sgn(\Delta_{\hat{\gamma}_l}(Y_i - Y_j) | N_{l}) \\
&+ |I_l|2(p_l - 2) \mathbb{C}ov(sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_j), sgn(\Delta_{\hat{\gamma}_l})(Y_i - Y_k) | N_{l}) \biggr) \to 0,
\end{align*}
where the convergence of the variance term follows directly by $|I_l|/|I_l|^2 \to 0$ and the covariance term goes to zero since $p_l/|I_l| \to 0$ by Lemma (ref). Hence $\hat{m}_l - m_0 \to_p 0$. Then by the Continuous Mapping theorem (CMT) and the fact that $L$ is fixed $\hat{m} \to_p m_0$. Finally since $\theta_0 = m_0/(2\mathbb{E}[Y_i])$ and $\hat{\theta} = \hat{m}/((2/n)\sum_{i=1}^n Y_i)$, $\hat{\theta} \to_p \theta_0$ by the Law of Large numbers and the CMT. The final asymptotic normality result follows exactly as the mean value expansion proof of Theorem (ref) by leveraging the asymptotic equivalence with the infeasible estimator, which has been proven above and by
\begin{align*}
\psi(W_i,W_j,\gamma_0,\theta) &= \theta(Y_i + Y_j) - sgn(\gamma_0(X_i) - \gamma_0(X_j))(Y_i - Y_j) \\
\frac{\partial \psi(W_i,W_j,\gamma_0,\theta)}{\partial \theta} &= (Y_i + Y_j).
\end{align*}
Proofs for the Pairwise Difference Estimator
We shall assume the same conditions as in ahn1993semiparametric, which we repeat here for the sake of completeness.
assumptionThe observed data $W_{i}=(Y_{i},X_{i},D_{i})$ is an
iid sample, with $D_{i}\in \{0,1\}$, and $(Y_{i},X_{i})$ having sixth-order
moments.
assumptionThe data satisfy the restrictions in ((ref)),
with $X_{i}=(S_{i},C_{i})$.
assumptionThe conditional distribution of $\gamma
_{i}:=\gamma _{0}(X_{i})$ given $D_{i}=1$ is absolutely continuous with
respect to the Lebesgue measure, with conditional density function $
f_{\gamma _{0}}(\cdot )$ that is continuous and bounded from above.
assumptionThe matrix $\Sigma _{ZS}:=\mathbb{E}[\gamma _{i}^{2}f_{i}\left( Z_{i}-\mu _{z}(\gamma
_{i})\right) \left( S_{i}-\mu _{s}(\gamma _{i})\right) ] $
is non-singular, where $f_{i}:=f_{\gamma _{0}}(\gamma _{i}),$ $\mu
_{z}(\gamma ):=\mathbb{E}[D_{i}Z_{i}|\gamma _{i}=\gamma ]/\gamma _{i}$ and $
\mu _{s}(\gamma ):=\mathbb{E}[D_{i}S_{i}|\gamma _{i}=\gamma ]/\gamma _{i}$.
assumptionThe kernel function $K$ is symmetric, twice
continuously differentiable, with derivatives $\dot{K}$ and $\ddot{K}$,
respectively; $\int t^{l}K\left( t\right) dt=0$ for $0<l<4$; for some $
l_{0}>0$, $K(t)=0$ for $\left\vert t\right\vert >l_{0}$.
assumptionThe bandwidth $h$ is of the form $h=c_{n}n^{-\delta
},$ where for a positive constant $c_{0},$ $c_{0}<c_{n}<c_{0}^{-1}$ and $
\delta \in (1/8,1/6)$.
assumptionThe conditional density $f_{\gamma _{0}}(\cdot ),$
the functions $\mu _{z}(\cdot )$ and $\mu _{s}(\cdot )$ and the function $
\lambda _{0}(\cdot )\ $and its derivative are all fourth-order continuously
differentiable, with derivatives that are bounded for all $\gamma $ in the
support of $\gamma _{i}$ given $D_{i}=1$. Furthermore, define
\begin{equation*}
\sigma _{ij}^{2}=\mathbb{E}[\varepsilon _{ij}^{2}(\theta _{0})\left\vert
Z_{i}-Z_{j}\right\vert ^{2}D_{i}D_{j}|\gamma _{i},\gamma _{j}],
\end{equation*}
\begin{equation*}
\sigma _{ij,ZS}^{2}=\mathbb{E}[\left\vert Z_{i}-Z_{j}\right\vert
^{2}\left\vert S_{i}-S_{j}\right\vert ^{2}|\gamma _{i},\gamma _{j}],
\end{equation*}
and assume that $\sigma _{ij}^{2},\sigma _{ij,ZS}^{2}<C$ a.s.
These are essentially the same assumptions as in ahn1993semiparametric, and
are extensively discussed there. The only additional assumption we have added is the
bounded conditional variances in Assumption (ref).
We define implicitly the rates $r_{q}$ and $r_{\lambda }$ as $||\hat{\gamma}
_{l}-\gamma _{0}||_{q}=O_{P}(n^{-r_{q}}),$ $q\geq 1,$ and $\left\Vert \hat{
\lambda}_{l}-\lambda _{0}\right\Vert =O_{P}(n^{-r_{\lambda }}),$
respectively. Recall $\delta$ is defined in Assumption (ref).
assumption(i) $r_{2}>\delta ;$ (ii) $r_{\lambda }>3\delta /2;$
(iii) $r_{2q}>[(3+1/q)\delta /4]+1/4$ for some $q>1;$ (iv) $r_{2}+r_{\lambda
}>3\delta /2+1/2;$ (v) $\mathbb{E}[\left\vert
Z_{i}-Z_{j}\right\vert ^{2q/(q-1)}|\gamma _{i},\gamma _{j}]<C$, a.s., for $q$ as in (iii).
These are mild rate conditions on the nuisance parameters. These conditions require that $r_{2}>1/8$, $r_{\lambda }>3/16\ $
and $r_{2}+r_{\lambda }>11/16.$ The latter is a product rate condition$.$
For $q\downarrow 1,$ (iii) requires that $r_{2}>\delta +1/4>3/8$, while at
the other extreme, for $q=\infty ,$ $r_{\infty }>3\delta /4+1/4>11/32,$
which is achieved, for example, by the smoothing kernel estimator of ahn1993semiparametric.
proof[Proof of Proposition (ref)]
To simplify notation, define the weights
\begin{equation*}
\dot{\varpi}_{ij}(\gamma )=\frac{1}{h^{2}}\dot{K}\left( \frac{\gamma
(X_{i})-\gamma (X_{j})}{h}\right) D_{i}D_{j}.
\end{equation*}
By the chain rule and iterated expectations
\begin{eqnarray*}
\frac{d}{d\tau }\mathbb{E}[g(W_{i},W_{j},\gamma _{\tau },\theta )] &=&\frac{d
}{d\tau }\mathbb{E}[\varepsilon _{ij}(\theta )\dot{\varpi}_{ij}(\gamma
_{0})(Z_{i}-Z_{j})\left( \gamma _{\tau }(X_{i})-\gamma _{\tau
}(X_{j})\right) ] \\
&=&\frac{d}{d\tau }\mathbb{E}[\alpha _{0}(X_{i},X_{j},\theta )\left( \gamma _{\tau
}(X_{i})-\gamma _{\tau }(X_{j})\right) ],
\end{eqnarray*}
where $\alpha _{0}(X_{i},X_{j},\theta )=\left[ \lambda_{0}(\gamma_0(X_{i}))-\lambda_{0}(\gamma_0(X_{j}))\right] q_{ij}(\gamma _{0})(Z_{i}-Z_{j}).$ Then, conclude applying Lemma (ref)(i).
We proceed to verify the sufficient conditions for the asymptotic theory, with some adaptations to deal with the local nature of the U-statistics.
lemmaUnder Assumptions (ref) to (ref),
\begin{equation*}
\sqrt{n}\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\psi
(W_{i},W_{j},\hat{\gamma}_{l},\hat{\alpha}_{l},\theta _{0})=\sqrt{n}\binom{n
}{2}^{-1}\sum_{i<j}\psi (W_{i},W_{j},\gamma _{0},\alpha _{0},\theta
_{0})+o_{p}(1).
\end{equation*}
proof[Proof of Lemma (ref)]
Under the conditions on the kernel, it can be shown that
\begin{equation*}
\left\vert K(y)-K(x)\right\vert \leq \left\vert y-x\right\vert K^{\ast }(x),
\end{equation*}
for a bounded and integrable kernel $K^{\ast },$ see A.8 in hansen2008uniform.
Henceforth, we will use the notation $\hat{\Delta}_{ij,l}=\hat{\gamma}
_{l}(X_{i})-\hat{\gamma}_{l}(X_{j})$ and $\Delta _{ij}\equiv \Delta
_{\gamma _{0}}(X_{i},X_{j}).$
To check Assumption 3(i), we use Cauchy-Schwarz and the boundedness of the
kernel $K^{\ast }$ and conditional variances to show, for $
w_{i}=(y_{i},x_{i},d_{i})$ and $w_{j}=(y_{j},x_{j},d_{j}),$
\begin{eqnarray*}
&&\int \int |g(w_{i},w_{j},\hat{\gamma}_{l},\theta
_{0})-g(w_{i},w_{j},\gamma _{0},\theta _{0})|^{2}F_{0}(dw_{i})F_{0}(dw_{j})
\\
&=&\int \int |\varepsilon _{ij}(\theta
_{0})(z_{i}-z_{j})d_{i}d_{j}|^{2}\left\vert \varpi _{ij}(\hat{\gamma}
_{l})-\varpi _{ij}(\gamma _{0})\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&\leq &\int \int |\varepsilon _{ij}(\theta
_{0})(z_{i}-z_{j})d_{i}d_{j}|^{2}\left\vert \hat{\Delta}_{ij,l}-\Delta
_{ij}\right\vert ^{2}\left( h^{-1}K^{\ast }(\Delta _{ij}/h)\right)
^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&\leq &h^{-2}\left\Vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\Vert ^{2} \\
&=&O_{P}\left( n^{2\delta }||\hat{\gamma}_{l}-\gamma _{0}||^{2}\right) \\
&=&o_{P}\left( 1\right) .
\end{eqnarray*}
Similarly, by H\"{o}lder inequality, for any $q>1,$ and with $\alpha_{i,j}=\alpha_0(X_{i},X_{j})$,
\begin{eqnarray}
&&\int \int |\phi (w_{i},w_{j},\hat{\gamma}_{l},\alpha _{0},\theta
_{0})-\phi (w_{i},w_{j},\gamma _{0},\alpha _{0},\theta
_{0})|^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \notag \\
&=&\int \int |\alpha _{ij}|^{2}\left\vert \hat{\Delta}_{ij,l}-\Delta
_{ij}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \notag \\
&\leq &C\int \int \left\vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\vert
^{2}\left( \frac{1}{h^{2}}\dot{K}\left( \frac{\Delta _{ij}}{h}\right)
\right) ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \notag \\
&\leq &Ch^{-4}\left\Vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\Vert
_{2q}^{2}\left( \int \int \left( \dot{K}\left( \frac{\Delta _{ij}}{h}\right)
\right) ^{2q/(q-1)}F_{0}(dw_{i})F_{0}(dw_{j})\right) ^{1-1/q} \notag \\
&\leq &Ch^{-4}\left\Vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\Vert
_{2q}^{2}h^{1-1/q} \notag \\
&=&O_{P}\left( n^{(3+1/q)\delta }\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert _{2q}^{2}\right) \notag \\
&=&o_{P}\left( 1\right) ,
\end{eqnarray}
where we have used that $(3+1/q)\delta /2<[(3+1/q)\delta /4]+1/4$ for any $
q\geq 1,$ and
\begin{equation*}
\mathbb{E}[|\alpha _{ij}|^{2}|X_{i},X_{j}]\leq C\left( \frac{1}{h^{2}}\dot{K}
\left( \frac{\Delta _{ij}}{h}\right) \right) ^{2}.
\end{equation*}
Now, we proceed with Assumption 3(iii). Define
\begin{equation*}
\hat{D}_{ij}(\hat{\gamma}_{l})=\left[ \hat{\lambda}_{l}(\hat{\gamma}
_{l}(X_{i}))-\hat{\lambda}_{l}(\hat{\gamma}_{l}(X_{j}))-\lambda _{0}(\hat{
\gamma}_{l}(X_{i}))-\lambda _{0}(\hat{\gamma}_{l}(X_{j}))\right],
\end{equation*}
and write
\begin{eqnarray*}
\hat{\alpha}_{l}-\alpha _{0} &=&\left[ \lambda _{0}(\hat{\gamma}
_{l}(X_{i}))-\lambda _{0}(\hat{\gamma}_{l}(X_{j}))-\lambda _{0}(\gamma
_{0}(X_{i}))-\lambda _{0}(\gamma _{0}(X_{j}))\right] q_{ij}(\hat{\gamma}
_{l})(Z_{i}-Z_{j}) \\
&&+\left[ \lambda _{0}(\gamma _{0}(X_{i}))-\lambda _{0}(\gamma _{0}(X_{j}))
\right] \left[ q_{ij}(\hat{\gamma}_{l})-q_{ij}(\gamma _{0})\right]
(Z_{i}-Z_{j}) \\
&&+\hat{D}_{ij}(\hat{\gamma}_{l})q_{ij}(\hat{\gamma}_{l})(Z_{i}-Z_{j}),
\end{eqnarray*}
so that, using that $\lambda _{0}$ is Lipschitz,
\begin{eqnarray*}
&&\int \int |\phi (w_{i},w_{j},\gamma _{0},\hat{\alpha}_{l},\theta
_{0})-\phi (w_{i},w_{j},\gamma _{0},\alpha _{0},\theta
_{0})|^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&=&\int \int |\hat{\alpha}_{l}-\alpha _{0}|^{2}\left\vert v_{ij}(\gamma
_{0})\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&\leq &2\int \int |\hat{\alpha}_{l}-\alpha
_{0}|^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&\leq &4\int \int \left\vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\vert
^{2}\left( \frac{1}{h^{2}}\dot{K}\left( \frac{\hat{\Delta}_{ij,l}}{h}\right)
\right) ^{2}\left\vert z_{i}-z_{j}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j})
\\
&&+4\int \int \left[ \lambda _{0}(\gamma _{0}(x_{i}))-\lambda _{0}(\gamma
_{0}(x_{j}))\right] ^{2}\left[ q_{ij}(\hat{\gamma}_{l})-q_{ij}(\gamma _{0})
\right] ^{2}\left\vert z_{i}-z_{j}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j})
\\
&&+4\int \int \hat{D}_{ij}^{2}(\hat{\gamma}_{l})q_{ij}^{2}(\hat{\gamma}
_{l})\left\vert z_{i}-z_{j}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&:&=4(I+II+III).
\end{eqnarray*}
From ((ref)) it follows that $I=o_{P}\left( 1\right) $.
To deal with $II$, write
\begin{equation*}
q_{ij}(\hat{\gamma}_{l})-q_{ij}(\gamma _{0})=\hat{\gamma}_{l}(X_{i})\hat{
\gamma}_{l}(X_{j})(Z_{i}-Z_{j})\frac{1}{h^{3}}\ddot{K}\left( \frac{\Delta
_{ij}}{h}\right) \left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) +R_{ij,q},
\end{equation*}
where $R_{ij,q}$ collects smaller order terms. Then, the leading term in $II$
is bounded by H\"{o}lder inequality, for any $q>1,$ by
\begin{eqnarray}
&&\int \int \left[ \lambda _{0}(\gamma _{0}(x_{i}))-\lambda _{0}(\gamma
_{0}(x_{j}))\right] ^{2}\frac{1}{h^{6}}\ddot{K}^{2}\left( \frac{\Delta _{ij}
}{h}\right) \left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) ^{2}\left\vert
z_{i}-z_{j}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \notag \\
&\leq &Ch^{-6}\left\Vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\Vert
_{2q}^{2}\left( \int \int \left[ \lambda _{0}(\gamma _{0}(x_{i}))-\lambda
_{0}(\gamma _{0}(x_{j}))\right] ^{2p}\ddot{K}^{2p}\left( \frac{\Delta _{ij}}{
h}\right) F_{0}(dw_{i})F_{0}(dw_{j})\right) ^{1/p} \\
&\leq &Ch^{3-1/q}\left\Vert \hat{\Delta}_{ij,l}-\Delta _{ij}\right\Vert
_{2q}^{2} \notag \\
&=&O_{P}\left( n^{-\delta (3-1/q)}\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert _{2q}^{2}\right) \notag \\
&=&o_{P}\left( 1\right) , \notag
\end{eqnarray}
where we have used the Lipschitz property of $\lambda _{0}.$
The term $III$ is bounded by a standard change of variables by
\begin{eqnarray*}
&&\int \int \hat{D}_{ij}^{2}(\hat{\gamma}_{l})\frac{1}{h^{4}}\dot{K}
^{2}\left( \frac{\hat{\gamma}_{l}(X_{i})-\hat{\gamma}_{l}(X_{j})}{h}\right)
\left\vert z_{i}-z_{j}\right\vert ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&\leq &Ch^{-3}\left\Vert \hat{\lambda}_{l}-\lambda _{0}\right\Vert ^{2} \\
&=&O_{P}\left( n^{3\delta }\left\Vert \hat{\lambda}_{l}-\lambda
_{0}\right\Vert ^{2}\right) \\
&=&o_{P}\left( 1\right) .
\end{eqnarray*}
The interaction term has the form
\begin{equation*}
\hat{\xi}_{l}(w_{i},w_{j})=\left( \hat{\alpha}_{l}(x_{i},x_{j})-\alpha
_{0}(x_{i},x_{j})\right) \left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) ,
\end{equation*}
and hence, Assumption 4 follows if $\sqrt{n}||\hat{\alpha}
_{l}-\alpha _{0}||||\hat{\gamma}_{l}-\gamma _{0}||=o_{p}(1).$ From our previous expansions, this product rate condition holds if
\begin{equation*}
\sqrt{n}n^{(3+1/q)\delta /2}\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert _{2q}||\hat{\gamma}_{l}-\gamma _{0}|| \\
=O_{P}\left( n^{(3+1/q)\delta /2+1/2}\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert _{2q}^{2}\right) \\
=o_{p}(1),
\end{equation*}
and
\begin{equation*}
\sqrt{n}n^{3\delta /2}\left\Vert \hat{\lambda}_{l}-\lambda _{0}\right\Vert ||
\hat{\gamma}_{l}-\gamma _{0}||=o_{p}(1).
\end{equation*}
We proceed to verify Assumption 5, $
\sqrt{n}\bar{\psi}(\hat{\gamma}_{l},\alpha _{0},\theta _{0})\rightarrow
_{p}0 $. Define
\begin{eqnarray*}
Q_{ij} &:&=\varpi _{ij}(\hat{\gamma}_{l})-\varpi _{ij}(\gamma _{0}) \\
R_{ij} &:&=Q_{ij}-\frac{1}{h^{2}}\dot{K}\left( \frac{\gamma _{i}-\gamma _{j}
}{h}\right) D_{i}D_{j}\left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) ,
\end{eqnarray*}
and note that by a Taylor's series expansion
\begin{eqnarray*}
Q_{ij} &:&=\frac{1}{h^{2}}\dot{K}\left( \frac{\Delta _{Q}}{h}\right)
D_{i}D_{j}\left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) , \\
R_{ij} &:&=\frac{1}{h^{3}}\ddot{K}\left( \frac{\Delta _{R}}{h}\right)
D_{i}D_{j}\left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right) ^{2},
\end{eqnarray*}
where $\Delta _{Q}$ and $\Delta _{R}$ are intermediate values between $\hat{
\Delta}_{ij,l}$ and $\Delta _{ij}.$ Using this, write,
\begin{eqnarray*}
\bar{\psi}(\hat{\gamma}_{l},\alpha _{0},\theta _{0}) &=&\int \int
\varepsilon _{ij}(\theta _{0})\varpi _{ij}(\hat{\gamma}
_{l})(z_{i}-z_{j})F_{0}(dw_{i})F_{0}(dw_{j}) \\
&&+\int \int \alpha _{ij}v_{ij}(\hat{\gamma}_{l})F_{0}(dw_{i})F_{0}(dw_{j})
\\
&=&\int \int \varepsilon _{ij}(\theta _{0})\varpi _{ij}(\gamma
_{0})(z_{i}-z_{j})F_{0}(dw_{i})F_{0}(dw_{j}) \\
&&+\int \int \varepsilon _{ij}(\theta _{0})(z_{i}-z_{j})\frac{1}{h^{3}}\ddot{
K}\left( \frac{\Delta _{R}}{h}\right) d_{i}d_{j}\left( \hat{\Delta}
_{ij,l}-\Delta _{ij}\right) ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&:&=I+II.
\end{eqnarray*}
To simplify the notation, define the function
\begin{equation*}
m(\gamma _{1},\gamma _{2})=\gamma _{1}\gamma _{2}\left( \lambda _{0}(\gamma
_{1})-\lambda _{0}(\gamma _{2})\right) \left( \mu _{z}(\gamma _{1})-\mu
_{z}(\gamma _{2})\right) .
\end{equation*}
A standard Taylor expansion and Assumptions (ref)-(ref) yield
\begin{eqnarray*}
I &=&\int \int m(\gamma _{1},\gamma _{2})\frac{1}{h}K\left( \frac{\gamma
_{1}-\gamma _{2}}{h}\right) f_{\gamma _{0}}(\gamma _{1})f_{\gamma
_{0}}(\gamma _{2})d\gamma _{1}d\gamma _{2} \\
&=&O(h^{4}).
\end{eqnarray*}
Hence, $\sqrt{n}I=n^{(1-8\delta )/2}\rightarrow 0.$
Additionally, note that
\begin{eqnarray*}
II &=&\int \int m(\gamma _{0}(x_{i}),\gamma _{0}(x_{j}))\frac{1}{h^{3}}\ddot{
K}\left( \frac{\Delta _{R}}{h}\right) \left( \hat{\Delta}_{ij,l}-\Delta
_{ij}\right) ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&=&\int \int m(\gamma _{0}(x_{i}),\gamma _{0}(x_{j}))\frac{1}{h^{3}}\ddot{K}
\left( \frac{\Delta _{ij}}{h}\right) \left( \hat{\Delta}_{ij,l}-\Delta
_{ij}\right) ^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&&+\int \int m(\gamma _{0}(x_{i}),\gamma _{0}(x_{j}))\frac{1}{h^{3}}\left[
\ddot{K}\left( \frac{\Delta _{R}}{h}\right) -\ddot{K}\left( \frac{\Delta
_{ij}}{h}\right) \right] \left( \hat{\Delta}_{ij,l}-\Delta _{ij}\right)
^{2}F_{0}(dw_{i})F_{0}(dw_{j}) \\
&=&:II_{A}+II_{B}.
\end{eqnarray*}
To deal with $II_{A},$ we apply H\"{o}lder inequality, for any $q>1,$ and
argue as in ((ref)) to obtain
\begin{eqnarray*}
\sqrt{n}\left\vert II_{A}\right\vert &=&O_{P}\left( n^{-\delta
(3-1/q)+1/2}\left\Vert \hat{\gamma}_{l}-\gamma _{0}\right\Vert
_{2q}^{2}\right) \\
&=&o_{p}(1),
\end{eqnarray*}
because $r_{2q}>[(3+1/q)\delta /4]+1/4$.
Using again A.8 in hansen2008uniform and H\"{o}lder inequality, for any $q>1,$
\begin{eqnarray*}
\sqrt{n}\left\vert II_{B}\right\vert &\leq &\sqrt{n}h^{-2}\left\Vert \hat{\gamma}_{l}-\gamma _{0}\right\Vert
_{2q}^{3} \\
&\leq &O_{P}\left( n^{2\delta +1/2}\left\Vert \hat{\gamma}_{l}-\gamma
_{0}\right\Vert _{2q}^{3}\right) \\
&=&o_{p}(1).
\end{eqnarray*}
lemmaUnder Assumptions (ref) to (ref),
\begin{equation*}
\sqrt{n}\binom{n}{2}^{-1}\sum_{i<j}\psi (W_{i},W_{j},\gamma _{0},\alpha
_{0},\theta _{0})\rightarrow _{d}N(0,\Sigma ),
\end{equation*}
where $\Sigma =\mathbb{E}[\xi _{i}\xi _{i}^{\prime }],$ and $\xi
_{i}\equiv U_{i}D_{i}\gamma _{i}\left( Z_{i}-\mu _{z}(\gamma _{i})\right)
f_{i}+\dot{\lambda}_{0}(\gamma _{i})\gamma _{i}^{2}\left( Z_{i}-\mu
_{z}(\gamma _{i})\right) v_{i}.$
proof[Proof of Lemma (ref)]
We apply the Hajek projection as in Lemma A.3 in ahn1993semiparametric to
obtain
\begin{equation*}
h_{n}(W_{i}):=\mathbb{E}[\psi (W_{i},W_{j},\gamma _{0},\alpha _{0},\theta
_{0})|W_{i}]=\mathbb{E}[\varepsilon _{ij}(\theta _{0})\varpi _{ij}(\gamma
_{0})(Z_{i}-Z_{j})+\alpha _{0}(X_{i},X_{j})v_{ij}(\gamma _{0})|W_{i}].
\end{equation*}
To that end, we compute, with $K_{ij}=K((\gamma_0(X_{i})-\gamma_0(X_{j}))/h)$,
\begin{eqnarray*}
&&\mathbb{E}[\varepsilon _{ij}(\theta _{0})\varpi _{ij}(\gamma
_{0})(Z_{i}-Z_{j})|W_{i},\gamma _{j}] \\
&=&\left( \lambda _{0}(\gamma _{i})-\lambda _{0}(\gamma _{j})\right)
D_{i}\gamma _{j}\left( Z_{i}-\mu _{z}(\gamma _{j})\right) \frac{1}{h}K_{ij}
\\
&&+U_{i}D_{i}\gamma _{j}\left( Z_{i}-\mu _{z}(\gamma _{j})\right) \frac{1}{h}
K_{ij}.
\end{eqnarray*}
Hence,
\begin{equation*}
\mathbb{E}[\varepsilon _{ij}(\theta _{0})\varpi _{ij}(\gamma
_{0})(Z_{i}-Z_{j})|W_{i}]\equiv U_{i}D_{i}\gamma _{i}\left( Z_{i}-\mu
_{z}(\gamma _{i})\right) f_{i}+\tilde{r}_{1,i}.
\end{equation*}
Similarly,
\begin{eqnarray*}
\mathbb{E}[\alpha _{0}(X_{i},X_{j})v_{ij}(\gamma _{0})|W_{i},\gamma _{j}] &=&
\mathbb{E}[\alpha _{0}(X_{i},X_{j})v_{i}|W_{i},\gamma _{j}] \\
&=&\left( \lambda _{0}(\gamma _{i})-\lambda _{0}(\gamma _{j})\right) \frac{1
}{h^{2}}\dot{K}_{ij}\gamma _{i}\gamma _{j}\left( Z_{i}-\mu _{z}(\gamma
_{j})\right) v_{i},
\end{eqnarray*}
(where we have used that $\mathbb{E}[\gamma _{j}Z_{j}|\gamma _{j}]=\gamma
_{j}\mu _{z}(\gamma _{j})$), and
\begin{equation*}
\mathbb{E}[\alpha _{0}(X_{i},X_{j})v_{ij}(\gamma _{0})|W_{i}]\equiv \dot{
\lambda}_{0}(\gamma _{i})\gamma _{i}^{2}\left( Z_{i}-\mu _{z}(\gamma
_{i})\right) v_{i}+\tilde{r}_{2,i}.
\end{equation*}
Define $\xi _{i}\equiv U_{i}D_{i}\gamma _{i}\left( Z_{i}-\mu _{z}(\gamma
_{i})\right) f_{i}+\dot{\lambda}_{0}(\gamma _{i})\gamma _{i}^{2}\left(
Z_{i}-\mu _{z}(\gamma _{i})\right) v_{i},$ and write
\begin{eqnarray*}
\frac{1}{\sqrt{n}}\sum_{i=1}^{n}h_{n}(W_{i}) &=&\frac{1}{\sqrt{n}}
\sum_{i=1}^{n}\xi _{i}+\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\tilde{r}_{1,i}+\frac{
1}{\sqrt{n}}\sum_{i=1}^{n}\tilde{r}_{2,i} \\
&=&\frac{1}{\sqrt{n}}\sum_{i=1}^{n}\xi _{i}+o_{p}(1) \\
&\rightarrow &_{d}N(0,\Sigma ),
\end{eqnarray*}
where the first equality follows by definition, the second from smoothness,
and Assumptions (ref)-(ref), so that
\begin{eqnarray*}
\left\vert \frac{1}{\sqrt{n}}\sum_{i=1}^{n}\tilde{r}_{l,i}\right\vert
&=&O_{p}(\sqrt{n}h^{4}) \\
&=&o_{p}(1),
\end{eqnarray*}
and the last convergence follows from the standard Central Limit Theorem
(CLT).
lemmaUnder Assumptions (ref) to (ref), if $\bar{\theta}
\rightarrow _{p}\theta _{0}$ then $\partial \hat{\psi}(\bar{\theta}
)/\partial \theta \rightarrow _{p}B$, $B=2\Sigma _{ZS}.$
proof[Proof of Lemma (ref)]
Note $\partial \hat{\psi}(\bar{\theta})/\partial \theta =
\hat{\Sigma}_{c,ZS}$ for all $\bar{\theta}.$ Write
\begin{eqnarray*}
\hat{\Sigma}_{c,ZS} &=&\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in
I_{l}}\varpi _{ij}(\hat{\gamma}_{l})(Z_{i}-Z_{j})(S_{i}-S_{j})^{\prime } \\
&=&\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\varpi _{ij}(\gamma
_{0})(Z_{i}-Z_{j})(S_{i}-S_{j})^{\prime } \\
&&+\binom{n}{2}^{-1}\sum_{l=1}^{L}\sum_{(i,j)\in I_{l}}\left[ \varpi _{ij}(
\hat{\gamma}_{l})-\varpi _{ij}(\gamma _{0})\right]
(Z_{i}-Z_{j})(S_{i}-S_{j})^{\prime } \\
&\equiv &\tilde{\Sigma}_{c,ZS}+R_{c,ZS}
\end{eqnarray*}
Lemma 3.1 in ahn1993semiparametric has shown $\tilde{\Sigma}_{c,ZS}=2\Sigma
_{ZS}+o_{P}(1).$ It remains to show that $R_{c,ZS}=o_{P}(1)$, which follows from the same arguments as those in the proof of Lemma (ref).
proof[Proof of Theorem (ref)]
By the Mean Value Theorem and Lemma (ref), (with $\tilde{\psi}(\theta _{0})$ defined as ($\hat{\psi}(\theta _{0})$) but with true nuisance parameters)
\begin{eqnarray*}
0 &=&\hat{\psi}(\theta _{0})-\tilde{\Sigma}_{c,ZS}(\hat{\theta}-\theta _{0})
\\
&=&\tilde{\psi}(\theta _{0})-\tilde{\Sigma}_{c,ZS}(\hat{\theta}-\theta
_{0})+o_{P}(1)
\end{eqnarray*}
Lemma 3.1 in ahn1993semiparametric has shown $\tilde{\Sigma}_{c,ZS}=2\Sigma
_{ZS}+o_{P}(1).$ Thus,
\begin{eqnarray*}
\sqrt{n}(\hat{\theta}-\theta _{0}) &=&\sqrt{n}2\Sigma _{ZS}^{-1}\hat{\psi}
(\theta _{0})+o_{P}(1) \\
&=&\frac{\Sigma _{ZS}^{-1}}{\sqrt{n}}\sum_{i=1}^{n}\xi _{i}+o_{P}(1) \\
&\rightarrow &_{d}N(0,V)
\end{eqnarray*}
by the CLT, since $\mathbb{E}[|\xi _{i}|^{2}]<\infty \ $.
Computations for the Gaussian First-Step Error Model
In this section, we provide some calculations for the Gaussian First-Step
Model that are used in the main text. Recall that in this DGP $X\sim
\mathcal{N}(0,\sigma _{X}^{2})$ and
equation*[equation* omitted — 84 chars of source]
where $\varepsilon \sim \mathcal{N}(0,\sigma _{\varepsilon }^{2})$, $\sigma
_{\varepsilon }^{2}=\mathbb{V}ar(\gamma _{0}(X))/StN$, and $StN$ is the
Signal to Noise Ratio. This is DGP is called a Gaussian First-Step Error
Model because we generate
equation*[equation* omitted — 217 chars of source]
where $c_{1},c_{2},\rho _{1},\rho _{2}$ are constants. We write $\delta
_{n}(x)=\mu _{n}(x)+\sigma _{\gamma n}U,$ where $U$ is a standard normal,
and
eqnarray*[eqnarray* omitted — 137 chars of source]
We first show that
equation*[equation* omitted — 189 chars of source]
where $C_{1}=\theta _{0}c_{1}$, $\sigma _{n}^{2}=C_{2}n^{1-\rho _{2}}$ and $
C_{2}=c_{2}(\mathbb{E}[2Y])^{-2}\theta _{0}^{2}/StN$. Therefore, the bias
coming from this derivative term diverges when $1/4<\rho _{1}<1/2$. This
illustrates the bias problem and the lack of root-$n$ consistency of plug-in
estimates based on non-locally robust moment functions mentioned in the
introduction.
From our derivative calculations in the proof of Proposition 1,
equation*[equation* omitted — 251 chars of source]
Hence, $\frac{\partial \theta (\gamma _{0})}{\partial \gamma }\left[ \hat{
\gamma}-\gamma _{0}\right] $ follows a normal distribution with mean
equation*[equation* omitted — 156 chars of source]
and variance
equation*[equation* omitted — 89 chars of source]
where we have used that $\mathbb{E}[sgn\left( X_{i}^{2}-X_{j}^{2}\right) ]=0.
$
Next, we verify the conditions for the asymptotic theory (i.e. Assumption
7). Straightforward calculations show that the random variable $\Delta
_{0}=\gamma _{0}(X_{i})-\gamma _{0}(X_{j})=X_{i}^{2}-X_{j}^{2}$ is
absolutely continuous with a density
equation*[equation* omitted — 177 chars of source]
where $K_{0}$ is the modified Bessel function of the second kind of order 0.
It is known that as $z\downarrow 0$ $K_{0}\left( z\right) \sim -\ln z-\xi ,$
where $\xi $ is Euler's constant, hence Assumption 7(ii) and (iii) holds
with any $\beta <1.$ This means that we require $||\hat{\gamma}-\gamma
_{0}||=o_{p}(n^{-\rho _{2\beta }}),$ with
equation*[equation* omitted — 66 chars of source]
Since this holds for any $\beta <1,$ it must be that $\rho _{2\beta } > 3/8$ and $||\hat{\gamma}_l - \gamma_0||^2 = o_p(n^{-3/4})$. Note that
align*[align* omitted — 301 chars of source]
so we need $\min(2\rho_1, \rho_2) > 3/4$ which implies $\rho _{1}>3/8=0.375$
and $\rho _{2}>3/4$. Note that these are sufficient conditions. In the
simulations, we consider values below and above the required to
evaluate the robustness of the finite sample performance to our sufficient
conditions.