The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
104,232 characters
Debiased Machine Learning of Set-Identified Linear Models
\maketitle
\begin{abstract}
This paper provides estimation and inference methods for an identified set's boundary (i.e., support function) where the selection among a very large number of covariates is based on modern regularized tools. I characterize the boundary using a semiparametric moment equation. Combining Neyman-orthogonality and sample splitting ideas, I construct a root-$N$ consistent, uniformly asymptotically Gaussian estimator of the boundary and propose a multiplier bootstrap procedure to conduct inference. I apply this result to the Partially Linear Model, the Partially Linear IV Model and the Average Partial Derivative with an interval-valued outcome.
\end{abstract}
\section{Introduction and Motivation.}
Interval-valued outcomes are ubiquitous in economic research. Examples of such outcomes include bidders' valuation in English auctions (\cite{HaileTamer}), income and wages (\cite{Trostel}, \cite{Gafarov}), house prices (\cite{GRT}, \cite{BerSasaki}), and county-level employment rates (\cite{Dorn}). An outcome is interval-valued if the actual outcome $Y$ is missing, but there exist an observable lower bound $Y_L$ and an upper bound $Y_U$ so that
\begin{align}
\label{eq:ylyu}
Y_L \leq Y \leq Y_U \text{ a.s. }
\end{align}
When an outcome is interval-valued, the parameter of interest is a set, where each point corresponds to a possible random variable $Y$ in the band \eqref{eq:ylyu}.
The main contribution of this paper is to provide an estimator of the identified set's boundary, where the selection among high-dimensional controls is based on modern machine learning/regularized methods. The paper focuses on identified sets whose boundary can be represented by a moment equation as in \cite{BM}, \cite{BMM}. In this paper, the equation depends on an identified functional nuisance parameter, for example, a conditional mean function. A naive approach would be to plug-in a machine learning estimate of the nuisance parameter into the moment equation and solve for the boundary. However, modern regularized methods (machine learning techniques) have bias converging slower than the parametric rate, which cannot be made small by classic techniques (e.g., undersmoothing). As a result, plugging such estimates into the moment equation produces a biased, suboptimal estimate of the boundary itself.
To overcome the transmission of the bias into the second stage, I adjust the moment equation to make it insensitive or, formally, Neyman-orthogonal, to the biased estimation of the nuisance parameter. While orthogonality has been extensively studied in the point-identified case, set identification presents several challenges. The first one is the non-smoothness of $x \rightarrow \min (x,0)$ function at $x=0$, often occurred in censored LAD (e.g., \cite{Powell}). The second one is to establish uniformly valid inference over the boundary in addition the pointwise one.
When the identified set is multi-dimensional, its boundary consists of continuum points. As a result, economists are interested in uniform inference in addition to the pointwise inference. Establishing uniform inference is not trivial. To control the speed at which an empirical sample average concentrates around the population mean, I invoke maximal inequalities of \cite{CCKAS} instead of Markov inequality that is typically sufficient in the point-identified case. I propose multiplier bootstrap algorithm to conduct inference. By virtue of orthogonality, only the moment function (the second stage), not the nuisance parameter estimate (the first stage), needs to be resampled in simulation. As a result, this multiplier bootstrap is faster to compute than the weighted bootstrap of \cite{CCMS}, which is based on a non-orthogonal moment and involves resampling of both stages. I demonstrate the proposed approach in a simulation exercise and give a brief empirical illustration.
\paragraph{Literature review. Set identification} This paper bridges the gap between three literatures: set-identified models, debiased/orthogonal machine learning, and non-smooth models. Set identification is a vast area of research, encompassing a wide variety of approaches: linear and quadratic programming, random set theory, support function, and moment inequalities
(\cite{Manski90}, \cite{ManskiPepper}, \cite{Manski:2002}, \cite{HaileTamer}, \cite{CHT}, \cite{BM}, \cite{Molinari2008}, \cite{CilibertoTamer}, \cite{LeeBound}, \cite{Stoye}, \cite{AndrewsShiECMA}, \cite{BMM2}, \cite{CCMS}, \cite{BMM3}, \cite{CherRigStoker}, \cite{BMM}, \cite{CLR}, \cite{FanPark}, \cite{KaidoWhite}, \cite{KaidoSantos}, \cite{KaidoWhite2}, \cite{Pakesetal}, \cite{ShiShum}, \cite{Kaido:2016}, \cite{Kasy2016}, \cite{KlineTamer}, \cite{AndrewsShi}, \cite{CanayBugniShi}, \cite{Kaido}, \cite{ChenTamerChristensen}, \cite{GafarovMeierOlea}, \cite{Shi}, \cite{Gafarov}, \cite{KaidoMolinariStoye}, \cite{SyrgkanisTamer}, \cite{Torgovitsky},\cite{MolinariStoye}, \cite{BerSasaki}, \cite{Honore}, \cite{AndrewsRothPakes}, \cite{kallus2020localized}, \cite{FanTao}, \cite{FanTao2}, \cite{MolinariMolchanovPeng}, \cite{HsiehShiShum}, \cite{DongHsiehShum}), see e.g. \cite{Tamer:2010} or \cite{Molinari:2018} for a review. This paper generalizes the \cite{BMM}'s model by allowing its components to depend on a functional nuisance parameter, covering e.g., \cite{CCMS} and \cite{Kaido} as special cases.
\paragraph {Orthogonality. }Next, this paper contributes to a large body of work on debiased inference for parameters following regularization or model selection (\cite{Neyman:1959}, \cite{Neyman:1979}, \cite{HardleStoker1989}, \cite{NeweyStoker}, \cite{Newey1994}, \cite{Robins}, \cite{robinson:88}, \cite{ZhangZhang}, \cite{JM}, \cite{chernozhukov2016double}, \cite{LRSP}, \cite{Program}, \cite{sasaki2018estimation}, \cite{sasaki2020unconditional}, \cite{Sasaki}, \cite{chiang2019multiway}, \cite{ning2020doubly}, \cite{chernozhukov2021debiased}, \cite{chernozhukov2021automatic}, \cite{CherSem}, \cite{NSS}, \cite{singh2020debiased}, \cite{Colangelo}, \cite{Lieli}, \cite{ZimLech}). A basic idea is to make the moment condition insensitive, or, formally, Neyman-orthogonal, to the biased estimation of the nuisance parameter. For a semiparametric GMM setting, the work by \cite{AckerbergChen} derives an orthogonal moment condition whose nuisance functions are identified by conditional moment restriction, such as conditional mean and conditional quantiles. Combining Neyman-orthogonality and sample splitting, \cite{LRSP} and \cite{chernozhukov2016double} derive a root-$N$ consistent and asymptotically normal estimator for a single target parameter. This has idea has been extended for many functional parameters in $Z$-estimation framework, in the context of distribution regression (\cite{BelCherWei}) and quantile regression (\cite{sasaki2020unconditional}). Next, the paper is related to literature on the non-smooth estimating equations (\cite{Powell}, \cite{Powell2}, \cite{PowellStockStoker}, \cite{kaplan_sun_2017}, \cite{franguridi2021conditional}). Finally, the paper contributed to a small, but growing literature on machine learning for bounds and partially identified models (\cite{kallus2019assessing}, \cite{jeong2020robust}, \cite{SemSupp2}, \cite{Bonvini_2021}).
\paragraph{Structure of the paper.} The paper is organized as follows. Section \ref{sec:setup1} demonstrates main points for the partially linear model of \cite{robinson:88}. Section \ref{sec:theory} states theoretical results. Section \ref{sec:application} applies the results to models with an interval-valued outcome. Section \ref{sec:montecarlo} presents finite-sample evidence. Section \ref{sec:empirical} contains an empirical illustration. Section \ref{sec:proofs} contains the proofs of main results.
\section{Set-Up.}
\label{sec:setup1}
\subsection{ General Framework }
\label{sec:setup}
I focus on parameters that are linear in an unobserved scalar outcome $Y$. The identified set takes the form
\begin{align}
\label{eq:idset}
\mathcal{B} = \{ \beta = \Sigma^{-1} \mathbb{E} V (\eta_0) Y, \quad Y_L \leq Y \leq Y_U\},
\end{align}
where the random vector $V(\eta_0) \in \mathrm{R}^{d}$ depends on a nuisance function $\eta_0$. Examples of the nuisance functions
$$
\eta_0 = \eta_0(X)
$$
include the propensity score, the conditional density, and the regression function, among others. The matrix $\Sigma \in \mathrm{R}^{d \times d}$ is identified by a moment equation
\begin{align}
\label{eq:sigma}
\Sigma = \mathbb{E} A(W,\eta_0)
\end{align}
and is assumed invertible. The key innovation of this framework is to allow the vector $V(\eta)$, the matrix function $A(W, \eta)$, and the bounds $Y_L, Y_U$ to depend on a functional nuisance parameter $\eta$, covering the models in \cite{BMM}, \cite{CCMS}, \cite{Kaido}, and many others as special cases.
\subsection{Examples}
\begin{example}[Partially Linear Model]
\label{ex:plp}
Consider the partially linear model of \cite{robinson:88}
\begin{align}
\label{eq:plm}
Y &= D' \beta_0 + f_0(X) + U, \quad \mathbb{E} [U\mid D,X] = 0,
\end{align}
where $D \in \mathrm{R}^d$ is a treatment (policy) variable, $\beta_0$ is a causal (structural) parameter, $$X = (1,X_1, X_2, \dots, X_{p_X}) $$ is a vector of covariates whose dimension may be large relative to the sample size (e.g., $p_X \gg N$), and $f_0(\cdot)$ is an integrable function. The parameter $\beta_0$ can be represented as the minimizer of the least squares criterion function
\begin{align}
\label{eq:plppointlong}
\beta_0& =\arg \min_{b \in \mathrm{R}^d, f \in L_2(P)} \mathbb{E} (Y - D' b - f(X))^2.
\end{align}
As in Lemma \ref{lem:fromlongtoshort}, $\beta_0$ coincides with the minimizer of a shorter criterion function
\begin{align}
\label{eq:plppoint1}
\beta_0& =\arg \min_{b \in \mathrm{R}^d} \mathbb{E} (Y - (D - \eta_0(X))' b)^2.
\end{align}
Thus, the identified set $\mathcal{B}$ is a special case of model \eqref{eq:idset}-\eqref{eq:sigma} with $V(\eta)$ and $A(W,\eta)$ defined as follows. The treatment regression function is
\begin{align}
\label{eq:cexp}
\eta_0(X) = \mathbb{E}[ D \mid X],
\end{align}
the treatment residual is
\begin{align}
\label{eq:resid}
V (\eta) = D-\eta(X),
\end{align}
the matrix function is
\begin{align}
\label{eq:aweta}
A(W,\eta) = (D - \eta(X))(D - \eta(X))'.
\end{align}
\end{example}
\begin{example}[Partially Linear IV Model]
\label{ex:plpiv}
Consider the following partially linear IV model
\begin{align}
\label{eq:plpiv}
Y &= D' \beta_0 + f_0(X) + U, \quad \mathbb{E}[ U \mid Z, X]=0 \\
D &= m_0(X) + E, \quad \mathbb{E} [ E \mid X ] =0, \label{eq:m0} \\
Z &= \eta_0(X) + V, \quad \mathbb{E}[ V \mid X ] =0, \label{eq:eta0}
\end{align}
where $Z \in \mathrm{R}^d$ is the instrument for $D$. The parameter $\beta_0$ can be characterized by a moment equation
\begin{align*}
\mathbb{E} (Z - \eta_0(X)) ( Y - (D - m_0(X))' \beta_0) = 0,
\end{align*}
which gives a closed-form expression for $\beta_0$
\begin{align*}
\beta_0 = (\mathbb{E} (Z - \eta_0(X)) (D - m_0(X))' )^{-1} \mathbb{E} (Z - \eta_0(X)) Y.
\end{align*}
The model is a special case of \eqref{eq:idset}-\eqref{eq:sigma} with
\begin{align}
V (\eta) &= Z - \eta(X) \label{eq:plpivspec} \\
A (W, \eta,m) &= (Z - \eta(X)) (D - m(X))', \label{eq:plpivspec2}
\end{align}
where $V = V(\eta_0)$ in \eqref{eq:eta0}. If $Z = D$, the model \eqref{eq:plpiv}-\eqref{eq:eta0} coincides with \eqref{eq:plm}-\eqref{eq:cexp}.
\end{example}
\begin{example}[Average Partial Derivative]
\label{ex:apd} An important parameter in economics is the average partial derivative. This parameter shows the average effect of a small change in a variable of interest $D$ on the outcome $Y$ conditional on the covariates $X$. To describe this change, define the conditional expectation function of an outcome $Y$ given the variable $D$ and exogenous variable $X$ as $$\mu(D,X):= \mathbb{E} [Y|D,X]$$ and its partial derivative with respect to $D$ as $\nabla_D\mu(D,X):=\nabla_D\mu(d,X)|_{d=D}$. Then, the average partial derivative is defined as
\begin{align}
\label{eq:apd}
\beta &= \mathbb{E} \nabla_D\mu(D,X).
\end{align}
For example, when $Y$ is the logarithm of consumption, $D$ is the logarithm of price, and $X$ is the vector of other demand attributes, the average partial derivative stands for the average price elasticity.
Assume that the variable $D$ has bounded support $\mathcal{D} \subset \mathcal{R}^d$. Furthermore, the conditional density $f ( D \mid X)$ has positive density on this support a.s. in $X$. \cite{HardleStoker1989} have shown that the average partial derivative can be represented as $\beta = \mathbb{E} V Y,$
where $V = - \nabla_D \log f(D|X) = - \frac{\nabla_D f(D|X)}{f(D|X)}$ is the negative partial derivative of the logarithm of the density $f(D|X)$. When $Y$ is interval-valued, the identified set $\mathcal{B}$ for $\beta$ is a special case of \eqref{eq:idset}-\eqref{eq:sigma} with $A(W, \eta) = \Sigma = I_d$, the nuisance function $\eta_0(D, X) = \frac{\nabla_D f(D|X)}{f(D|X)}$ and the vector $V ( \eta) = -\eta$. \cite{Kaido} studies a special case of this problem without covariates.
\end{example}
\subsection{Single treatment}
\label{sec:singletreat}
In this section, I derive an orthogonal moment equation for the upper bound $\beta_U$ on the causal parameter $\beta_0$ in Example \ref{ex:plp}. Because the bias due to sign mistake in $Y^{\text{best}}(\eta) \neq Y^{\text{best}}(\eta_0)$ proves to be second-order, I derive the orthogonal moment treating $Y^{\text{best}}(\eta_0)$ as observed. Next, I formally control the bias due to the sign mistake.
The following two subsections elaborate on these points.
\paragraph{Moment equation for $\beta_U$. }Consider Example \ref{ex:plp} with a single treatment (i.e., $d=1$). The identified set $\mathcal{B}$ becomes a closed interval $[\beta_L, \beta_U]$. The upper bound $\beta_U$ is
\begin{align}
\label{eq:betu1d}
\beta_U = \max_{\{ Y: Y_L \leq Y \leq Y_U \}} \bigg\{ \dfrac{\mathbb{E} (D - \eta_0(X)) \cdot Y }{\mathbb{E} (D - \eta_0(X))^2 } \bigg\}.
\end{align}
To maximize the numerator of \eqref{eq:betu1d}, take $Y=Y_U$ for positive values of $D - \eta_0(X)$ and $Y=Y_L$ otherwise. Define the best-case outcome
\begin{align}
\label{eq:ubg}
Y^{\text{best}} (\eta) &= \begin{cases} Y_L, \quad D - \eta(X) \leq 0, \\
Y_U, \quad D - \eta(X) > 0.
\end{cases}
\end{align}
Plugging \eqref{eq:ubg} into \eqref{eq:betu1d} gives the moment function for $\beta_U$
\begin{align}
\label{eq:naive:1d}
m (W, \beta_U, \eta) := (Y^{\text{best}}(\eta) - (D- \eta(X)) \beta_U) (D -\eta(X)).
\end{align}
\paragraph{ Establishing orthogonality. } In what follows, $\eta$ corresponds to an instance of the nuisance parameter whose true value is $\eta_0$. Consider an infeasible moment function
\begin{align}
\label{eq:naive:1d0}
m_0(W, \beta_U, \pmb{\eta} ) := (Y^{\text{best}}(\eta_0) - (D- \pmb{\eta} (X)) \beta_U) (D - \pmb{\eta} (X)),
\end{align}
where the best-case outcome $Y^{\text{best}}(\eta_0) $ is treated as if it was observed, and $ \pmb{\eta}$ indicates the instance of the nuisance parameter treated as unknown. The moment equation \eqref{eq:naive:1d0} is not orthogonal to the perturbations of $\widehat{\eta}(X) - \eta_0(X)$
\begin{align*}
\partial_{r} \mathbb{E} [ m_0 (W, \beta_U, r(\eta - \eta_0) + \eta_0) ] |_{r=0} = -\mathbb{E}[ Y^{\text{best}}(\eta_0) ( \eta(X) - \eta_0(X)) ] \neq 0.
\end{align*}
Therefore, the bias of the estimation error $\widehat{\eta}(X)-\eta_0(X)$ translates into the moment \eqref{eq:naive:1d0}. To overcome the transmission of this bias, \cite{robinson:88} proposes an orthogonal moment equation
\begin{align*}
g_0 (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \})= (Y^{\text{best}}(\eta_0)- \pmb{\gamma}_U(X) - (D- \pmb{\eta} (X)) \beta_U) (D - \pmb{\eta} (X)),
\end{align*}
where the true value of $\gamma_U(x)$ is $$\gamma_{U,0}(X) = \mathbb{E}[ Y^{\text{best}}(\eta_0) \mid X]. $$ The orthogonality condition is
\begin{align}
\label{eq:derivative}
&\partial_{r} \mathbb{E} [ g_0 (W, \beta_U, r(\eta - \eta_0) + \eta_0, \gamma_{U,0} ) ] |_{r=0} \\
&=- \mathbb{E}[ (Y^{\text{best}}(\eta_0) - \mathbb{E}[ Y^{\text{best}}(\eta_0)\mid X]) ( \eta(X) - \eta_0(X)) ] =0. \nonumber
\end{align}
Thus, the bias of the estimation error, $\widehat{\eta}(X) - \eta_0(X)$, does not translate into the moment \eqref{eq:naive:1d0}. The orthogonality condition with respect to $\gamma_{U,0}$ can be verified in a similar way.
\paragraph{ Bias due to sign mistake. } I now discuss whether the mistake in the best-case outcome $Y^{\text{best}} (\eta) \neq Y^{\text{best}} (\eta_0)$, which has been so far ignored, has any effect on the support function estimate. Consider the difference between the feasible moment $g (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \})$ and its infeasible analog
\begin{align*}
& g (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \}) - g_0 (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \}) = (D - \eta(X)) ( Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) ) \\
&= (D - \eta_0(X)) ( Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) ) \\
&+ ( \eta_0(X) - \eta(X)) ( Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) ).
\end{align*}
Define the first-order bias $B_1( \eta, \eta_0)$ as
\begin{align*}
B_1( \eta, \eta_0):= \mathbb{E} [ (D - \eta_0(X)) ( Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) )]
\end{align*}
and the second-order one
\begin{align*}
B_2( \eta, \eta_0):= \mathbb{E} [ ( \eta_0(X) - \eta(X)) ( Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) )].
\end{align*}
Below, I describe the conditions under which $B_1( \eta, \eta_0)$ and $ B_2( \eta, \eta_0)$ are negligible, (i.e., $o(N^{-1/2})$). If they hold, the feasible moment equation $g (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \})$ is insensitive to the biased estimation of $\eta$. As a result, the support function estimator based on $g (W, \beta_U, \{ \pmb{\eta}, \pmb{\gamma}_U \})$ is asymptotically unbiased under plausible conditions.
Consider the best-case outcome $Y^{\text{best}} (\eta)$
$$
Y^{\text{best}} (\eta) = Y_L + (Y_U - Y_L) 1 \{ D- \eta(X) \geq 0 \}.
$$
The sign mistake $Y^{\text{best}} (\eta) \neq Y^{\text{best}} (\eta_0)$ occurs on the events
\begin{align}
\label{eq:e-}
\mathcal{E}_{-} :&= \bigg\{ D - \eta_0(X) < 0 < D - \eta(X) \bigg\}, \\
\label{eq:e+} \mathcal{E}_{+} :&= \bigg\{ D - \eta(X) < 0 < D - \eta_0(X) \bigg\}.
\end{align}
On these events, the residual cannot exceed estimation error in absolute value
\begin{align}
\label{eq:small}
\bigg\{ \mathcal{E}_{-} \cup \mathcal{E}_{+} \bigg\} &\Rightarrow \bigg\{ 0 < | D -\eta_0(X) | < | \eta(X) - \eta_0(X) | \bigg\}.
\end{align}
If the width $Y_U - Y_L$ is bounded by $M_{UL}$, the estimation error of $Y^{\text{best}}(\eta) $ is bounded as
\begin{align}
\label{eq:bestcaseerror}
&|Y^{\text{best}}(\eta) - Y^{\text{best}}(\eta_0) |=(Y_U - Y_L) 1 \{ \mathcal{E}_{-} \cup \mathcal{E}_{+} \} \\
&\leq M_{UL} 1 \{ 0 < | D -\eta_0(X) | < | \eta(X) - \eta_0(X) | \}. \nonumber
\end{align}
Suppose the conditional density $ h_{V (\eta_0) \mid X} (t,X) $ of $V(\eta_0)$ is bounded by $M_h$ a.s.. Invoking \eqref{eq:small} gives
\begin{align}
\mathbb{E} | \eta (X) - \eta_0(X) | 1{\{\mathcal{E}_{+ } \cup \mathcal{E}_{-}\}} &\leq \mathbb{E}_{X} \int_{- | \eta (X) - \eta_0(X) | }^{| \eta (X) - \eta_0(X) |} t h_{ V (\eta_0) | X} (t) dt \nonumber \\
&\leq 2 M_h \mathbb{E}_{X} ( \eta (X) - \eta_0(X))^2. \label{eq:mainbound2}
\end{align}
As a result, the bias terms $B_1( \eta, \eta_0)$ and $B_2( \eta, \eta_0)$ shrink at the quadratic rate. Combining \eqref{eq:mainbound2} and \eqref{eq:derivative} gives a feasible moment function
\begin{align}
\label{eq:feasible2}
g (W, \{ \pmb{\eta}, \pmb{\gamma}_U \})&= (Y^{\text{best}}(\pmb{\eta}) - \pmb{\gamma}_U (X) - (D- \pmb{\eta}(X)) \beta_U) (D -\pmb{\eta}(X)).
\end{align}
\subsection{Multi-dimensional case}
\label{sec:multid}
In this section, I derive an orthogonal moment for the support function, starting from a non-orthogonal one due to \cite{BMM}, \cite{BM}.
\paragraph{Moment Equation for Support Function. } As shown in \cite{BM}, the identified set $\mathcal{B}$ in \eqref{eq:idset} is a compact and convex set. Thus, it can be described by its projections onto a unit sphere
\begin{align}
\label{eq:unitsphere}
\mathcal{S}^{d-1} := \{q \in \mathrm{R}^{d}, \quad \| q\| = 1\}.
\end{align}
For any direction $q \in \mathcal{S}^{d-1}$, define the support function as the upper bound on $q' \beta_0$
\begin{align}
\label{eq:suppfun}
\sigma(q):= \sup_{b \in \mathcal{B}} q' b.
\end{align}
As proposed in \cite{BM} and \cite{BMM}, define the projected weighting vector
\begin{align}
\label{eq:zp}
z(p,\eta) = p' V (\eta),
\end{align}
the best-case outcome $Y(p,\eta)$
\begin{align}
\label{eq:yq}
Y (p,\eta) &= Y_L + (Y_U - Y_L) 1\{ z(p,\eta)>0\},
\end{align}
and the projection parameter $p(q)$
\begin{align}
\label{eq:pq}
p(q) = \Sigma^{-1} q.
\end{align}
Then, the moment equation for $\sigma(q)$ is
\begin{align}
\label{eq:zqwq}
\sigma(q) &= \mathbb{E} [z (p, \eta_0) Y (p,\eta_0)] \big|_{p = p(q)}.
\end{align}
\paragraph{Orthogonal Moment for Support Function. } The moment equation \eqref{eq:zqwq} is sensitive to the biased estimation of the nuisance parameter $\eta$. To avoid the transmission of this bias into the second stage, I construct another moment function $g(W, p, \xi),$ where $\xi(p)$ is the nuisance parameter, in the following steps.
\begin{enumerate}
\item Starting from an infeasible, smooth moment
\begin{align}
\label{eq:smooth:d0}
m_0(W, p, \pmb{\eta}) := z (p,\pmb{\eta}) Y (p,\eta_0)
\end{align}
derive an infeasible orthogonal moment $g_0(W, p, \xi(p)) $ obeying \eqref{eq:derivative} for each $p$.
\item Invoke Lemma \ref{lem:powell} to bound the bias
\begin{align}
\label{eq:ubias}
\sup_{p \in \mathcal{P}} \mathbb{E} | z(p,\eta) \left(Y (p,\eta)- Y (p,\eta_0) \right) | = O( \mathbb{E} \| \eta(X) - \eta_0(X) \|^2),
\end{align}
where $p$ belongs to a compact bounded set $\mathcal{P}$ defined below. When $\Sigma = I_d$, $\mathcal{P} = \mathcal{S}^{d-1}$.
\item Combine \eqref{eq:ubias} and \eqref{eq:smooth:d0} to obtain the feasible orthogonal moment
\begin{align}
\label{eq:feasible}
g(W, p, \xi(p)) = g_0(W, p, \xi(p)) + z(p,\eta) (Y (p,\eta)- Y (p,\eta_0)).
\end{align}
\end{enumerate}
\begin{example*}[Example \ref{ex:plpiv}, cont.]
Consider Example \ref{ex:plpiv}. The projected weighting vector \eqref{eq:zp} is
\begin{align*}
z(p, \eta) = p' (Z - \eta(X)).
\end{align*}
The infeasible orthogonal moment $g_0(W, p, \xi(p)) $ is
\begin{align}
\label{eq:orthomom1}
g_0(W, p, \xi(p)) = z(p,\eta) (Y (p,\eta_0) - \gamma (p, X)),
\end{align}
where
\begin{align}
\label{eq:rieszgamma}
\gamma_0(p,x) &= \mathbb{E} [ Y (p,\eta_0) \mid X] \\
&= \mathbb{E} [ Y_L \mid X=x] + \mathbb{E} [ (Y_U - Y_L) 1\{ p' V(\eta_0) > 0 \} \mid X=x] \nonumber
\end{align}
The nuisance function $\xi_0(p) = (\eta(\cdot), \gamma(p,\cdot))$. Invoking \eqref{eq:feasible} gives a feasible orthogonal moment
\begin{align}
\label{eq:orthomom}
g(W, p, \xi(p)) = z(p,\eta) (Y (p,\eta) - \gamma (p, X)).
\end{align}
Corollary \ref{cor:plpiv} establishes the asymptotic theory for the support function estimator based on \eqref{eq:orthomom}.
\end{example*}
\begin{example*}[Example \ref{ex:apd}, cont.]
Consider Example \ref{ex:apd}. The projected weighting vector is
\begin{align*}
z(q, \eta) = - q' \partial_D \log f(D \mid X).
\end{align*}
The infeasible orthogonal moment $g_0(W, p, \xi(p)) $ is
\begin{align}
\label{eq:rho:apd}
g_0(W,q,\xi(q)) = z(q, \eta) Y(q, \eta_0) + q' \partial_D \log f(D \mid X) \mu(q,D,X) + q' \nabla_D \mu(q,D,X),
\end{align}
where
\begin{align*}
\mu_0(q,D,X) &= \mathbb{E} [ Y_L \mid D, X] + \mathbb{E} [ (Y_U - Y_L) \mid D,X ] 1\{ -q' \partial_D \log f(D \mid X) > 0 \} \\
&= \gamma_{L,0}(D, X) + \gamma_{UL,0}(D, X) 1\{ -q' \partial_D \log f(D \mid X) > 0 \}.
\end{align*}
Invoking \eqref{eq:feasible} gives a feasible moment equation
\begin{align}
\label{eq:orthomom:apd}
g(W, p, \xi(p)) &= g_0(W,q,\xi(q)) + z(q,\eta) (Y (q,\eta) -Y (q,\eta_0)) \\
&= z(q,\eta) Y(q, \eta) +q' \partial_D \log f(D|X) \mu(q,D,X) + q' \nabla_D \mu(D,X) \nonumber.
\end{align}
In absence of the conditioning covariates, \eqref{eq:orthomom:apd} coincides with the efficient score in \cite{Kaido}. Corollary \ref{cor:apd} establishes the asymptotic theory for the support function estimator based on \eqref{eq:rho:apd}.
\end{example*}
\subsection{ Overview of Main Results}
\label{sec:overview}
The Support Function Estimator $\widehat{\sigma}(q)$ has two stages. In the first stage, I construct an estimate $\widehat{\xi}$ of the nuisance parameter $\xi_0$ using some regularized estimator. In the second stage, I compute the estimated values $(\widehat{\xi}_i)_{i=1}^N$ and the support function estimate. I use different samples in the first and the second stage in the form of cross-fitting. The number $K$ of cross-fit partitions is assumed to be fixed/finite relative to $N$.
\begin{definition}[Cross-Fitting]
\label{sampling}
\mbox{}
\begin{compactenum}
\item For a random sample of size $N$, denote a $K$-fold random partition of the sample indices $[N]=\{1,2,...,N\}$ by $(J_k)_{k=1}^K$, where $K$ is the number of partitions and the sample size of each fold is $n = N/K$. For each $k \in [K] = \{1,2,...,K\}$ define $J_k^c = \{1,2,...,N\} \setminus J_k$.
\item For each $k \in [K]$, construct estimates $\widehat{\xi}_k = \widehat{\xi}( W_{i \in J_k^c})$ and $\widehat{\eta}_k = \widehat{\eta}( W_{i \in J_k^c})$ of the nuisance parameters $\xi_0$ and $\eta_0$ using only the data $\{ W_{j}: j \in J_k^c \}$. For any observation $i \in J_k$, define $\widehat{\xi}_i = \widehat{\xi}_k (W_i)$ and $\widehat{\eta}_i = \widehat{\eta}_k (W_{i})$.
\end{compactenum}
\end{definition}
Definition \ref{sampling} introduces cross-fitting. Cross-fitting plays an essential role in modern debiased inference in semi-parametric models; see, e.g., \cite{bch:2010,zheng:laan,chernozhukov2016double} for recent examples and \cite{hasminskii:debiased} and \cite{schick1986asymptotically} for early, classical uses of simpler sample-splitting methods for debiased inference. Cross-fitting is proposed for the cases where the nuisance parameter $\xi_0(\cdot)$ does not depend on $p \in \mathcal{P}$, such as Example \ref{ex:apd} and Examples \ref{ex:plp}--\ref{ex:plpiv} under Assumption \ref{ass:suffcond}.
\begin{definition}[Support Function Estimator]
\label{def:estimate:psi}
Let $\widehat{\xi}$ and $\widehat{\eta}$ be the estimates of $\xi_0$ and $\eta_0$. Define
\begin{align}
\widehat{\Sigma} &= \dfrac{1}{N} \sum_{i=1}^N A(W_i, \widehat{\eta}_i), \quad \widehat{p}(q) = (\widehat{\Sigma}^{-1})' q \label{sigma:est} \\
\widehat{\sigma}(q)&= \dfrac{1}{N} \sum_{i=1}^N g(W_i, \widehat{p}(q) , \widehat{\xi}_i( \widehat{p}(q) )). \label{eq:psi:est}
\end{align}
\end{definition}
\begin{definition}[Multiplier Bootstrap]
\label{def:bb}
Let $(e_i)_{i=1}^N: e_i $ are i.i.d. truncated exponential random variables $\text{Tr Exp}(1)$ on $[0, \bar{M}]$ independent of the data. Define the bootstrap analog of $ \widehat{\sigma}(q)$ as
\begin{align}
\widetilde{\Sigma} &= \dfrac{1}{N} \sum_{i=1}^N \dfrac{e_i}{\bar{e}} A(W_i, \widehat{\eta}_i), \quad \widetilde{p}(q) = (\widetilde{\Sigma}^{-1})' q \label{sigma:est:boot} \\
\widetilde{\sigma}(q)&= \dfrac{1}{N} \sum_{i=1}^N \dfrac{e_i}{\bar{e}} g(W_i, \widetilde{p}(q) , \widehat{\xi}_i( \widetilde{p}(q) )). \label{eq:psi:boot}
\end{align}
\end{definition}
Under mild conditions on $\xi$, the Support Function Estimator delivers a high-quality estimate $\widehat{\sigma}(q)$ of the support function $\sigma(q)$ with the following properties
\begin{enumerate}
\item With probability (w.p.) $\rightarrow 1$, the estimator converges uniformly over the unit sphere $\mathcal{S}^{d-1}$
\begin{align}
\label{eq:urate}
\sup_{ q \in \mathcal{S}^{d-1} } | \widehat{\sigma} (q) - \sigma(q) | = O_P(1/\sqrt{N}) = o_P(1).
\end{align}
\item The estimator $\widehat{\sigma}(q)$ is asymptotically Gaussian
\begin{align}
\label{eq:limit}
S_N(q) :=\sqrt{N} (\widehat{\sigma}(q) - \sigma(q)) = \mathbb{G}_N(q)+ o_P(1) \quad \text{ uniformly in } \mathcal{S}^{d-1},
\end{align}
where the empirical process $\mathbb{G}_N(q)$ is approximated by a Gaussian process $\mathbb{G}(q)$, which is a tight $P$-Brownian bridge in $ \ell^{\infty}(\mathcal{S}^{d-1})$.
\end{enumerate}
Define the bootstrap statistic
\begin{align*}
\widetilde{S}_N(q):=\sqrt{N} (\widetilde{\sigma}(q) - \widehat{\sigma}(q))
\end{align*}
\paragraph{Pointwise asymptotics. } The sharp identified set for $q' \beta_0$ is $[ - \sigma(-q), \sigma(q)]$. Its $(1-\tau)$-pointwise confidence region (CR) is
\begin{align*}
[\underline{i}(q), \bar{i}(q)]:=[- \widehat{\sigma}(-q) + N^{-1/2} \widehat{C}_{\tau/2}(q), \quad \widehat{\sigma}(q) +N^{-1/2} \widehat{C}_{1-\tau/2}(q) ],
\end{align*}
where the critical values $\widehat{C}_{\tau/2} (q)$ and $\widehat{C}_{1-\tau/2} (q) $ are the $\tau/2$ and $1-\tau/2$ quantiles of the bootstrapped statistic $|\widetilde{S}_N(q)|$. Plugging $$q= \mathbf{e}_k = (0,0,\dots, 0, \underbrace{1}_{k}, 0, \dots, 0) \in \mathrm{R}^d, \quad k = 1,2, \dots, d$$ gives the CR for the projection of the identified set $\mathcal{B}$.
\paragraph{Uniform asymptotics. } The $(1-\tau)$-uniform confidence region (CR) for $[ - \sigma(-q), \sigma(q)]$ is
\begin{align}
\label{eq:uniform}
[\underline{i}_u(q), \bar{i}_u(q)]:=[- \widehat{\sigma}(-q) + N^{-1/2} \widehat{C}^{*}_{\tau/2}, \quad \widehat{\sigma}(q) +N^{-1/2} \widehat{C}^{*}_{1-\tau/2} ],
\end{align}
where $\widehat{C}^{*}_{\tau/2}$ and $\widehat{C}^{*}_{1-\tau/2}$ are the quantiles of the bootstrapped statistic $\sup_{q \in \mathcal{S}^{d-1}} |\widetilde{S}_N(q)|$. Likewise, for any function $f(\cdot)$ and a critical value $\widehat{c}_N = c_N + o_P(1)$ and $c_N = O_P(1)$,
\begin{align*}
{\mathrm{P}} ( f(S_N) \leq \widehat{c}_N) - {\mathrm{P}}^e ( f(\widetilde{S}_N) \leq \widehat{c}_N) \rightarrow_P 0,
\end{align*}
where ${\mathrm{P}}^{e} (\cdot)$ is the probability conditional on the data.
\section{ Theoretical Results.}
\label{sec:theory}
\setcounter{example}{0}
\paragraph{Notation.} I use the empirical process notation. For a generic function $f$ and a generic sample $(W_i)_{i=1}^N$, denote the empirical sample average by $$ {\mathbb{E}_{N}} f(W_i) := \dfrac{1}{N} \sum_{i=1}^N f(W_i) $$ and the
scaled, demeaned sample average by $$\mathbb{G}_N f(W_i) := 1/\sqrt{N} \sum_{i=1}^N [f(W_i) - \int f(w) d P(w)].$$
For two sequences of random variables $\{ a_N, b_N, N \geq 1\}: a_N \lesssim_{P} b_N$ means $a_N = O_{P} (b_N)$. For two sequences of numbers $\{a_N, b_N, N \geq 1\}$, $a_N \lesssim b_N$ means $a_N = O (b_N)$. Let $a \wedge b = \min \{ a, b\}, a \vee b = \max \{ a, b\} $. The $\ell_2$ norm of a vector is denoted by $\| \cdot \|$, the $\ell_1$ norm is denoted by $\| \cdot \|_1$, the $\ell_{\infty}$ norm is denoted by $\| \cdot \|_{\infty}$, and $\ell_0$ norm is denoted by $\| \cdot \|_{0}$. For a matrix $Q$, let $\|Q\|$ be the maximal eigenvalue of $Q$ and $\| Q \|_F$ be the Frobenius norm of $Q$. For a random vector $W$, let $\| W \|_{P,c}:= (\int |W|^c d P)^{1/c}$. The random sample $(W_i)_{i=1}^N$ is a sequence of independent copies of a random element $W$ taking values in a measurable space $(\mathcal{W}, \mathcal{A}_{\mathcal{W}})$ according to a probability law $P$. The $\| f \|_{P_N, 2}$ is the empirical $\ell_2$-norm, denoted as $\| f \|_{P_N, 2}:=(N^{-1} \sum_{i=1}^N f^2(W_i))^{1/2}$. Define the projection set
\begin{align}
\label{eq:pee}
\mathcal{P} = \bigg\{ p \in \mathrm{R}^d: \quad 1/2 \min \operatorname{eig} (\Sigma^{-1}) \leq \| p \| \leq 2 \max \operatorname{eig} (\Sigma^{-1}) \bigg\}
\end{align}
and let $C_P:= 2 \max \operatorname{eig} (\Sigma^{-1}) $. Let $\mathcal{F}_c$ be the space of continuous functions obeying two conditions: (1) $f(Z)$ has a continuous distribution when $Z$ is a tight Gaussian process with non-degenerate covariance function and (b) $f(\xi_N + c) -f(\xi_N) = o(1)$ for any $c=o(1)$ and any $\| \xi_N \|= O_P(1)$. Note that $[p_X]:=\{1,2,\dots, p_X\}$.
\subsection{Assumptions}
Assumption \ref{ass:sigma} is a standard identification condition. It ensures that the eigenvalues of $\Sigma$ are bounded from above and below.
\begin{assumption}[Identification]
\label{ass:sigma}
There exist constants $\lambda_{\min} >0$ and $\lambda_{\max} < \infty $ that bound the eigenvalues of $\Sigma$ in \eqref{eq:sigma} from above and below $0 < \lambda_{\min} \leq \min \operatorname{eig} (\Sigma) \leq \max \operatorname{eig}(\Sigma) \leq \lambda_{\max}$.
\end{assumption}
Assumption \ref{ass:jacobian} ensures that the support function is differentiable on the unit sphere. It requires the distribution of the weighting vector $V(\eta_0)$ to be sufficiently smooth. For example, if the vector $V(\eta_0)$ has a symmetric (i.e., spherical) distribution around zero, the normalized vector $\| V(\eta_0) \|^{-1} V(\eta_0)$ is uniformly distributed on the surface of the unit sphere $\mathcal{S}^{d-1}$, and Assumption \ref{ass:jacobian} holds. Assumption \ref{ass:jacobian} is a common regularity condition in set-identified models (e.g., Condition C.1 in \cite{CCMS}) and censored median regression (e.g., Assumption R.2 in \cite{Powell}).
\begin{assumption}[Smooth boundary]
\label{ass:jacobian}
Let $d \geq 2$. There exists a finite constant $\mathcal{C}_V$ such that
\begin{align}
\label{eq:smooth}
\sup_{q \in \mathrm{S}^{d-1}} {\mathrm{P}} \left( \dfrac{ |q' \Sigma^{-1/2} V (\eta_0) | }{ \| \Sigma^{-1/2} V (\eta_0) \| } \leq \delta \right) \leq \mathcal{C}_V \delta.
\end{align}
\end{assumption}
\begin{example}[Gaussian Noise]
\label{ex:gaussian}
Consider Example \ref{ex:plpiv} with $$V(\eta_0) = Z - \eta_0(X) \sim N(0, \Sigma)$$ independent of $X$. Then, $\Sigma^{-1/2} V(\eta_0) \sim N(0, I_d)$ is standard Gaussian vector and
$\Sigma^{-1/2} V(\eta_0)/ \| \Sigma^{-1/2} V (\eta_0) \|$ is uniformly distributed over the unit sphere $\mathcal{S}^{d-1}$. For any $q \in \mathcal{S}^{d-1}$, $|q'\Sigma^{-1/2} V(\eta_0)|/ \| \Sigma^{-1/2} V (\eta_0) \|$ is uniformly distributed on $[0,1]$ (\cite{Pitman}). As a result, \eqref{eq:smooth} holds with $\mathcal{C}_V = 1$ conditional on $X$ uniformly in $X$.
\end{example}
Assumption \ref{ass:jacobian} is a sufficient condition for the smoothness of the boundary. If it holds, the moment equation \eqref{eq:orthomom} is differentiable in $p$. The gradient
\begin{align}
\label{eq:gp}
G(p):= \mathbb{E} V (\eta_0) Y(p, \eta_0)
\end{align}
is a uniformly continuous function of $p$ (see Lemma \ref{lem:uniderivative} in Online Appendix). As a result, there exists a uniform Gaussian approximation for the support function estimator.
\begin{remark}
\label{rm:discrete}
Consider Example \ref{ex:plp}. If $D$ and $X$ consist of discrete variables only, the distribution of $V (\eta_0)$ cannot be continuous, and Assumption \ref{ass:jacobian} fails. Discrete distributions imply flat surfaces on the identified set, which may not be compatible with uniform Gaussian approximation. In this case,
\cite{CCMS} suggests adding a small amount of continuously distributed noise and work with a slightly expanded set with smooth boundary, while \cite{Gafarov} provides an alternative approach.
\end{remark}
The moment functions $A(W, \eta)$ in \eqref{eq:sigma} and $g(W, p, \xi(p))$ depend on the nuisance parameters $\eta_0$ and $\xi_0$, respectively. Definition \ref{def:nearorthog} introduces a sequence of nuisance realization sets $\Xi_N \subseteq \Xi$ and $\mathcal{T}_N \subseteq \mathcal{T}$ that contain $\xi_0$ and $\eta_0$, as well as their estimators $\widehat{\xi}$ and $\widehat{\eta}$, with probability $1-\phi_N$. As the sample size increases, the sets $\Xi_N$ and $\mathcal{T}_N$ shrink. The shrinkage speed is measured by the rates below.
\begin{definition}[Moment Rates]
\label{def:nearorthog}
Let $\{ \Xi_N, N \geq 1\}$ and $\{ \mathcal{T}_N, N \geq 1\}$ be sequences of subsets of $\Xi$ and $\mathcal{T}$, respectively, obeying the following conditions. (1) The true values $\xi_0$ and $\eta_0$ belong to $\Xi_N$ and $\mathcal{T}_N$ for all $N \geq 1$. There exists a sequence of numbers $\phi_N = o(1)$ such that the first-stage estimators $\widehat{\xi}$ of $\xi_0$ and $\widehat{\eta}$ of $\eta_0$ belong to $\Xi_N$ and $\mathcal{T}_N$ with probability at least $1-\phi_N$. Define the sequences $\mu_N, r_N'', r_N', A_N, \delta_N$ as
\begin{align}
\sup_{\xi \in \Xi_N} \sup_{p \in \mathcal{P}} | \mathbb{E} [ g(W, p, \xi(p)) - g (W,p, \xi_0(p)) ] | = \mu_N \label{eq:mun} \\
\sup_{\eta \in \mathcal{T}_N} \| \mathbb{E} [ A (W, \eta) - A(W, \eta_0) ] \| = A_N \label{eq:an} \\
\sup_{\xi \in \Xi_N} \sup_{p \in \mathcal{P}} (\mathbb{E} (g(W, p, \xi(p)) - g (W,p, \xi_0(p)) )^2 )^{1/2} = r_N'' \label{eq:rnprimeprime} \\
\sup_{p, p_0 \in \mathcal{P} \| p - p_0 \| \lesssim \tau_N} (\mathbb{E} ( g(W, p, \xi_0(p)) - g(W,p_0, \xi_0(p_0)) )^2 )^{1/2} = r_N' \label{eq:rnprime} \\
\sup_{\eta \in \mathcal{T}_N} (\mathbb{E} \| A(W, \eta) - A(W, \eta_0) \|^2)^{1/2} = \delta_N, \label{eq:delta}
\end{align}
where $\tau_N:= N^{-1/2} \log N$ for the case when $A(W,\eta_0)$ is being estimated and $\tau_N:= r_N':=0$ when $A(W, \eta_0) = \Sigma$ is known.
\end{definition}
\begin{assumption}[Regularity Conditions]
\label{ass:ratebound2}
(1) There exist absolute constants $c'>2$ and $\bar B_A< \infty$ such that the following matrix norms are bounded
$$\sup_{\eta \in \mathcal{T}_N} \mathbb{E} \| A(W, \eta) \|^{c'} \leq \bar B_A, \quad \mathbb{E} \| A(W, \eta_0) \|_F^{2} \leq \bar B_A$$ (2) There exists a sequence $v_N = o(N^{-1/4})$ such that $$ \| \mathbb{E}_N A(W_i, \eta_0) - \Sigma \| = O_P (v_N) = o_P(1). $$
\end{assumption}
Assumption \ref{ass:ratebound} requires the rates of Definition \ref{def:nearorthog} to decay sufficiently fast. The matrix moment function $A(W, \eta)$ is assumed to be already orthogonal with respect to $\eta$. This is the case for covariance matrices in Examples \ref{ex:plp} and \ref{ex:plpiv}.
\begin{assumption}[Convergence Rates]
\label{ass:ratebound}
(1) For $c_3>2$ in Assumption \ref{ass:ratebound2}, suppose the numerator moment rates $\mu_N, r_N'', r_N'$ obey the following bounds: $\mu_N = o(N^{-1/2})$ and $$(r_N'' + r_N' ) \log^{1/2} (1/(r_N'' + r_N')) + N^{-1/2+1/c_3} \log N = o(1).$$
(2) The matrix rates $A_N$ and $\delta_N$ obey the following bounds: $A_N = o(N^{-1/2})$ and $\delta_N = o(1)$.
\end{assumption}
Assumption \ref{ass:concentration:chap1} bounds the complexity of the function class $$\mathcal{G}_{\xi} = \{ g (W,p, \xi(p)), p \in \mathcal{P}\}.$$
\begin{assumption}[Complexity Conditions]
\label{ass:concentration:chap1}
(1) There exists a measurable envelope function $G_{\xi} = G_{\xi}(W)$ that almost surely bounds all elements in the class $$ \sup_{p \in \mathcal{P}} | g(W,p,\xi(p)) | \leq G_{\xi}(W)
\quad \text{a.s.}$$
There exists $c>2$ such that $\|G_{\xi}\|_{P,c} := \left(\int_{w \in \mathcal{W}} (G_{\xi}(w))^c \right)^{1/c} < \infty$. (2) There exist constants $a,v$ that do not depend on $N$ such that the uniform covering entropy of the function class $\mathcal{G}_{\xi}$ is bounded
\begin{align}
\label{eq:complexityq}
\log \sup_{Q} N(\epsilon \| G_{\xi} \|_{Q,2}, \mathcal{G}_{\xi} , \| \cdot \|_{Q,2}) \leq v \log (a/\epsilon), \quad \text{ for all } 0 < \epsilon \leq 1.
\end{align}
\end{assumption}
\subsection{Results}
Define the influence function $h_g(W,q)$ as
\begin{align}
\label{eq:hwgq}
h_g(W, q) &= g(W, p(q), \xi_0(p(q))) - \sigma(q)
\end{align}
and the influence function for the matrix estimation
\begin{align}
\label{eq:hwaq}
h_A(W, q) &= -G(p(q))' \Sigma^{-1} ( A(W, \eta_0) - \Sigma) \Sigma^{-1} q ,
\end{align}
where $G(p)$ is the gradient defined in \eqref{eq:gp}. Finally, define
\begin{align}
\label{eq:hwq}
h(W,q) = h_g(W, q)+h_A(W, q).
\end{align}
\begin{theorem}[Limit Theory for the Support Function Process]
\label{thm:limit}
Suppose Assumptions \ref{ass:sigma}-\ref{ass:concentration:chap1} hold. Then, the support function process $S_N(q) = \sqrt{N} (\widehat{\sigma} (q) - \sigma(q))$ is asymptotically linear uniformly on $ \mathcal{S}^{d-1}$
$$S_N(q) = \mathbb{G}_N [h(W,q)] + o_{P} (1)\text{ uniformly on } \mathcal{S}^{d-1},$$
where $h(W, q)$ is as in \eqref{eq:hwq}. Furthermore, the process $S_N(q)$ admits the following approximation
$$
S_N(q) =_d \mathbb{G}[h(q) ] + o_P(1) \quad \text{ in } \ell^{\infty} (\mathcal{S}^{d-1}),$$
where the process $ \mathbb{G}[h(q)] $ is a tight $P$-Brownian bridge in $ \ell^{\infty} (\mathcal{S}^{d-1})$ with a non-degenerate covariance function
\begin{align*}
\Omega(q_1,q_2 ) = \mathbb{E} [h(W,q_1)h(W,q_2 )] - \mathbb{E}[h(W,q_1)] \mathbb{E}[h(W,q_2 )], \quad q_1, q_2 \in \mathcal{S}^{d-1}.
\end{align*}
\end{theorem}
Theorem \ref{thm:limit} is my first main result. It says that the Support Function Estimator is asymptotically equivalent to a tight Gaussian process with a non-degenerate covariance function. Previous work (e.g., \cite{CCMS}, \cite{Kaido}) has derived similar Gaussian approximations for the estimators based on classic nonparametric methods. Combing Neyman-orthogonality and sample splitting, Theorem \ref{thm:limit} allows to accommodate both classic nonparametric and modern regularized/machine learning estimators. Corollary \ref{cor:limit} states that the inference properties of support function estimator.
\begin{corollary}[Limit Inference on Support Function Process]
\label{cor:limit}
Suppose Assumptions \ref{ass:sigma}--\ref{ass:concentration:chap1} hold. For any $f \in \mathcal{F}_c$, $\widehat{c}_N = c_N + o_P(1)$ and $c_N = O_P(1)$,
\begin{align*}
{\mathrm{P}} ( f(S_N) \leq \widehat{c}_N) - {\mathrm{P}} ( f( \mathbb{G}[h(q)]) \leq \widehat{c}_N) \rightarrow 0.
\end{align*}
If $c_N (1-\tau)$ is the $(1-\tau)$-quantile of $ f( \mathbb{G}[h(q)]) $ and $\widehat{c}_N(1-\tau) = c_N(1-\tau) + o_P(1)$ is any consistent estimate of this quantile, then
\begin{align*}
{\mathrm{P}} ( f(S_N) \leq \widehat{c}_N(1-\tau)) \rightarrow 1-\tau.
\end{align*}
\end{corollary}
\begin{theorem}[Limit Theory for the Bootstrap Support Function Process]
\label{thm:bb}
Under conditions of Theorem \ref{thm:limit}, the bootstrap support function process $\widetilde{S}_N(q) = \sqrt{N} (\widetilde{\sigma} (q) - \widehat{\sigma}(q))$ is asymptotically linear uniformly on $ \mathcal{S}^{d-1}$
\begin{align*}
\widetilde{S}_N(q) &= \mathbb{G}_N [ (e-1) h (W,q)]+ o_{P}(1) \text{ uniformly on } \mathcal{S}^{d-1}.
\end{align*}
Furthermore, the bootstrap support function process admits an approximation conditional on the data:
$$
\widetilde{S}_N(q) = \widetilde{\mathbb{G}[h(q)]} + o_{P^e}(1) \quad \text{ in } \ell^{\infty} (\mathcal{S}^{d-1}), \text{ in probability } P,$$
where $\widetilde{\mathbb{G}[h(q)]}$ is a tight $P$-Brownian bridge in $\ell^{\infty} (\mathcal{S}^{d-1})$ with the same distribution as the process $\mathbb{G}[h(q)]$ defined in Theorem \ref{thm:limit}, and independent of $ \mathbb{G}[h(q)] $.
\end{theorem}
Theorem \ref{thm:bb} is my second main result. It establishes the validity of multiplier bootstrap for uniform inference on the support function. In contrast to the weighted bootstrap of \cite{CCMS}, the multiplier bootstrap does not require re-estimating the first-stage nuisance parameter in each bootstrap repetition. Instead, the nuisance parameter is estimated on an auxiliary sample once and plugged into the bootstrap sampling procedure.
\begin{corollary}[Limit Inference on Bootstrap Support Function Process]
\label{cor:bb}
For any $\widehat{c}_N = c_N + o_P(1)$ and $c_N = O_P(1)$,
\begin{align*}
{\mathrm{P}} ( f(S_N) \leq \widehat{c}_N) - {\mathrm{P}}^e ( f(\widetilde{S}_N) \leq \widehat{c}_N) \rightarrow_P 0.
\end{align*}
If $\widehat{c}_N(1-\tau)$ is the $(1-\tau)$-quantile of $ f(\widetilde{S}_N) $ under $P^e$, then
\begin{align*}
{\mathrm{P}} (f(\widetilde{S}_N) \leq \widehat{c}_N(1-\tau)) \rightarrow_P 1-\tau.
\end{align*}
\end{corollary}
\section{Partially Linear IV Model}
\label{sec:application}
\setcounter{assumption}{0}
\setcounter{lemma}{0}
\setcounter{corollary}{0}
\subsection{The First Stage }
In this section, I give examples of the first-stage estimators as well as the nuisance realization sets. I focus on the partially linear IV model of Example \ref{ex:plpiv}. This discussion automatically covers Example \ref{ex:plp} that is a special case of Example \ref{ex:plpiv} with $D=Z$.
Definition \ref{def:fsrate} introduces sequences of nuisance realization sets $\{ \mathcal{T}_N, N \geq 1\}$ and $\{ \mathcal{M}_N, N \geq 1\}$ that contain the true value of $\eta_0$ and $m_0$ and their estimates $\widehat \eta$ and $\widehat m$ with probability $1-o(1)$. As the sample size $N$ increases, the sets $\mathcal{T}_N$ and $\mathcal{M}_N$
shrink. The shrinkage speed is measured by mean square rates $\eta_N$ and $m_N$, respectively.
\begin{definition}[Mean Square Rates]
\label{def:fsrate} Define the mean square rate for the expectation functions $\eta_0(x)$ and $m_0(x)$
\begin{align*}
\sup_{\eta \in \mathcal{T}_N} \left(\mathbb{E} \|\eta (X)- \eta_0(X)\|^2\right)^{1/2} = \eta_N \\
\sup_{m \in \mathcal{M}_N} \left(\mathbb{E} \| m (X)-m_0(X)\|^2\right)^{1/2} = m_N \end{align*}
\end{definition}
Definition \ref{def:fsrate} introduces mean square convergence rates for the expectation functions $\eta_0$ and $\gamma_0$. The bounds on $\eta_N$ are established for
a wide variety of regularized estimators, including partition estimators (\cite{CattaneoFarrell}, \cite{CattaneoFarrellFeng}), $\ell_2$-boosting (\cite{Luo}), deep neural networks (\cite{Schmidt}, \cite{Farrell}), linear and nonlinear sieve estimators \cite{Chen2007}, penalized sieve estimators \cite{Chen2011}, random forest in small (\cite{WagerWalther}) and high (\cite{syrgkanis2020estimation}) dimensions with sparsity structure. For the sake of brevity, the examples of nuisance realization sets $ \mathcal{M}_N$ and $\mathcal{T}_N$ are omitted in this text, but they can be found in \cite{SemCher}, Appendix B or \cite{CGST}, Section 5.
\paragraph{First-Stage Fitted Values: Basic Case.} In this paragraph, I consider the case when the nuisance parameter $\xi_0$ does not depend on $p$.
\begin{assumption}[Independent and Symmetric Residual]
\label{ass:suffcond}
The following conditions hold. (1) The interval width $Y_U - Y_L$ is independent of $V(\eta_0)$ conditional on $X$:
\begin{align}
\label{eq:indep}
(Y_U - Y_L) \quad \rotatebox[origin=c]{90}{$\models$} \quad V(\eta_0) \mid X.
\end{align}
(2) The vector $V(\eta_0)$ is independent of $X$, for any $p \in \mathcal{P}$. (3) $V(\eta_0)$ has a spherical distribution, which implies ${\mathrm{P}} (p' V(\eta_0) >0) = 1/2$ for any $p \in \mathcal{P}$.
\end{assumption}
Suppose Assumption \ref{ass:suffcond} holds. Define
$$
\gamma_{L,0}(X):= \mathbb{E}[ Y_L \mid X], \quad \gamma_{UL,0} (X) = \mathbb{E}[ (Y_U - Y_L) \mid X].
$$
Then, the Riesz representer function takes the form
\begin{align}
\gamma_0(p,X) &= \gamma_{L,0} (X)+1/2 \gamma_{UL,0} (X) \label{eq:midpoint}
\end{align}
and the first-stage fitted values take the form
\begin{align}
\label{eq:crossfitgammal}
\widehat {\gamma} (X_i):= \widehat{\gamma}_L (X_i) + 1/2 \widehat{\gamma}_{UL} (X_i), \quad i=1,2,\dots, N,
\end{align}
where $(\widehat{\gamma}_L (X_i), \widehat{\gamma}_{UL}(X_i))_{i=1}^N$ are the cross-fit first-stage estimates of Definition \ref{sampling}.
\paragraph{First-Stage Fitted Values: High-Dimensional Sparse Case.} In this paragraph, I sketch a possible estimator of the Riesz representer function without relying on Assumption \ref{ass:suffcond}. Define
$$
\rho_0(p,X) = \mathbb{E}[ \mathcal{Y} (p,\eta_0) \mid X] = \mathbb{E}[ (Y_U - Y_L) 1\{ p' V(\eta_0) > 0 \} \mid X].
$$
Assume that the expectation function is approximated as
\begin{align}
\label{eq:maindecomp}
\rho_0(p,X) = \Lambda(Z(X)' \nu_0(p)) + R_p(\eta_0,X),
\end{align}
where $Z: \mathcal{X} \rightarrow \mathrm{R}^{p_X}$ is a set of $p_X$ measurable basis functions of the covariates $X$, $\nu_0(p)$ is a $p_X$-dimensional vector, and $R_p(\eta_0,X)$ is an approximation error, and $\Lambda: \mathrm{R} \rightarrow \mathrm{R}$ is the linear link function (Example \ref{ex:linearlasso}) and logistic link function (Example \ref{ex:lassologistic}).
\begin{example}[Linear Lasso Estimator]
\label{ex:linearlasso}
Let $$\mathcal{Y}_i(p, \widehat \eta):= (Y_{U,i} - Y_{L,i}) 1 \{ p' V_i(\widehat \eta) >0 \}, \quad i=1,2,\dots, N$$ be an estimate of $\mathcal{Y}_i(p, \eta_0)$. Given the penalty level
$$
\lambda_{Z} = \dfrac{1.1}{N^{1/2}} \Phi^{-1} \left(1 - \dfrac{\bar{\gamma}}{ 2 p_X N^d} \right), \quad \bar{\gamma} = .1/ \log N
$$
and the diagonal matrix of penalty loadings $\widehat \Psi_{p}$, define
$$
\widehat \nu (p) \in \arg \min_{\nu \in \mathrm{R}^{p_X}} N^{-1} \sum_{i=1}^N ( \mathcal{Y}_i(p, \widehat \eta)- Z(X_i)' \nu )^2 + \lambda_{Z} \| \widehat \Psi_{p} \nu \|_1
$$
and the fitted value as
$$
\widehat \rho(p,X_i):= Z(X_i)' \widehat \nu (p).
$$
Given the cross-fit estimate $ \widehat \gamma_L (\cdot)$ of the $\gamma_{L,0}(\cdot)$, define
\begin{align}
\label{eq:finalfitted}
\widehat \gamma(p, X_i) = \widehat \gamma_L (X_i) +Z(X_i)' \widehat \nu (p).
\end{align}
\end{example}
An important special case occurs when the brackets have constant width
$$
Y_U - Y_L = \Delta \quad \text{a.s.},
$$
and the nuisance parameter reduces to the conditional probability
$$
\mathbb{E} [ \mathcal{M} (p, \eta_0) \mid X]:={\mathrm{P}} (p' V(\eta_0) > 0 \mid X).
$$
Suppose this probability can be approximated as
\begin{align}
\label{eq:maindecomplogistic}
{\mathrm{P}} (p' V(\eta_0) > 0 \mid X) = \Lambda (Z(X)' \nu_0(p)) + R_p(\eta_0,X),
\end{align}
where $ \Lambda (t) = \exp t/ (\exp t+1)$ is the logistic link function.
\begin{example}[Logistic Lasso Estimator]
\label{ex:lassologistic}
Let $$L(y,t):= -1 (1\{ y=1 \} \log \Lambda(t) + 1\{ y=0 \} \log (1-\Lambda(t)) )$$ be the logistic loss function and let $\mathcal{M}_i(p, \widehat \eta) = 1\{ p' V_i (\widehat \eta)>0\}$ be the estimated outcome.
Given the penalty level $\lambda_Z$ and the diagonal matrix of penalty loadings $\widehat \Psi_{p}$, define
$$
\widehat \nu (p) \in \arg \min_{\nu \in \mathrm{R}^{p_X}} N^{-1} \sum_{i=1}^N L (\mathcal{M}_i(p, \widehat \eta) , Z(X_i)' \nu ) + \lambda_{Z} \| \widehat \Psi_{p} \nu \|_1
$$
The final estimator of fitted values is
\begin{align}
\label{eq:finalfittedlogis} \widehat \gamma(p, X_i) &= \widehat \gamma_L (X_i) + \widehat \rho(p,X_i)=\widehat \gamma_L (X_i) + \Delta \Lambda( Z(X_i)' \widehat \nu(p)).
\end{align}
\end{example}
Lemma \ref{lem:riesz} provides the first-stage convergence rates for linear Lasso estimator. In contrast to regular linear Lasso model, the outcome $\mathcal{Y}(p, \eta_0)$ depends on the nuisance parameter $\eta_0$ and therefore has to be estimated. I show that the linear Lasso estimator remains valid as long as the first-stage mean square rate decays sufficiently fast.
\begin{lemma}[Validity of Linear Lasso with Estimated Outcome]
\label{lem:riesz}
The following conditions hold for $N$ large enough and a sequence $\zeta_N=o(1)$ and $\Lambda(t)=t$. (i) The model \eqref{eq:maindecomp} is approximately sparse with $s_{\nu} =s_{\nu}(N)$
\begin{align*}
\sup_{ p \in \mathcal{P}} \| \nu_0(p) \|_0 \leq s_N
\end{align*}
and $\log (p_X \vee N) \leq \zeta_N N^{1/3}$. (ii) Heteroscedasticity. There exists a constant $c_{\zeta} >0$ so that $0<c_{\zeta} \leq \mathbb{E}[ \zeta^2(p, \eta_0) \mid X ] \text{ a.s. }$. (iii) Lipschitz property of $\nu_0(p)$. For some finite constant $C_L, \sup_{p_1, p_2 \in \mathcal{P}} \| \nu_0(p_1) - \nu_0(p_2) \|_1 \leq C_L \| p_1 - p_2 \|$. (iv) The approximation error decays fast:
$$ \sup_{p \in \mathcal{P}} (\mathbb{E} R^2_0(p,X))^{1/2} = o (\bar{\sigma}_N), \quad \bar{\sigma}_N:=\sqrt{s_N \log (p_X \vee N)/N}.$$
(vi) First Stage. (a) The first-stage mean square rate is fast enough $\eta_N = o (\bar{\sigma}_N)$ and (b) There exists $\eta^{\infty}_N = o(1)$ so that
$\sup_{\eta \in \mathcal{T}_N} \sup_{x \in \mathcal{X}} \| \eta(x) - \eta_0(x) \| \leq \eta^{\infty}_N = o(1)$. Then, under additional conditions on $Z(X)$ in Assumption \ref{lem:rieszappendix1},
the estimate $\widehat \nu (p)$ of Example \ref{ex:linearlasso} is uniformly sparse, that is $\sup_{p \in \mathcal{P}} \| \widehat \eta(p) \|_0 \leq C_X s_N$,
and the following performance bounds hold:
\begin{align}
\sup_{ p \in \mathcal{P}} \| Z(X)' (\widehat \nu (p) - \eta_0(p)) \|_{P_N, 2} &\leq C_X \sqrt{\dfrac{s_N \log p_X }{N}} \label{eq:fsstage} \\
\sup_{ p \in \mathcal{P}} \| \widehat \nu (p) - \eta_0(p) \|_1 &\leq C_X \sqrt{\dfrac{s^2_N \log p_X }{N}}.\label{eq:fsstage2}
\end{align}
\end{lemma}
\begin{lemma}[Validity of Logistic Lasso with Estimated Outcome]
\label{lem:rieszlogistic}
Suppose the conditions of Lemma \ref{lem:riesz} hold for \eqref{eq:maindecomp} with logistic link function $\Lambda(t)$. Then, under Assumption \ref{lem:rieszappendix2},
the estimate $\widehat \nu (p)$ of Example \ref{ex:lassologistic} is uniformly sparse, that is $\sup_{p \in \mathcal{P}} \| \widehat \eta(p) \|_0 \leq \widetilde{C} s_N$, and the bounds \eqref{eq:fsstage}--\eqref{eq:fsstage2} hold.
\end{lemma}
\subsection{The Second Stage }
In this section, I describe the Support Function Estimator for the Partially Linear IV Model. The estimator for Example \ref{ex:plp} is obtained by replacing Steps 1 and 2 and 4 of Algorithm \ref{alg:plpiv} by their analogs in Algorithm \ref{alg:plp}. As an input, Algorithm \ref{alg:plp} takes a direction $q \in \mathcal{S}^{d-1}$ and the first-stage fitted values $(\widehat{\eta}(X_i), \widehat{m}(X_i), \widehat{\gamma}(p,X_i))_{i=1}^N$.
\begin{algorithm}[H]
\begin{algorithmic}[1]
\STATE The treatment residual $$\widehat{V}_i:= Z_i - \widehat{\eta}(X_i), \quad \widehat{E}_i:= D_i - \widehat{m}(X_i) \qquad i=1,2,\dots, N.$$
\STATE The sample covariance matrix: $$\widehat{\Sigma}:= \dfrac{1}{N} \sum_{i=1}^N \widehat{V}_i \widehat{E}_i' .$$
\STATE The best-case outcome as a function of $p$ $$Y_{i}(p, \widehat{\eta}) :=Y_{L,i} + (Y_{U,i}-Y_{L,i})1\{ p' \widehat{V}_i>0 \}, \quad i=1,2,\dots, N. $$
\STATE The IV coefficient of the second-stage residual $Y_{i}(\widehat{p}(q), \widehat{\eta}) - \widehat{\gamma}(\widehat{p}(q), X_i)$ on the instrument residual $\widehat{V}_i$
\begin{align}
\label{eq:betaq}
\widehat{\beta}_q = \widehat{\Sigma}^{-1} \dfrac{1}{N} \sum_{i=1}^N\widehat{V}_i [ Y_{i}(\widehat{p}(q) , \widehat{\eta}) - \widehat{\gamma} (\widehat{p}(q) ,X_i)], \quad \widehat{p}(q) =\widehat{\Sigma}^{-1}q
\end{align}
\STATE Report: the projection of $\widehat{\beta}_q$ onto the direction $q$: $$\widehat{\sigma}(q) = q' \widehat{\beta}_q.$$
\end{algorithmic}
\caption{Support Function Estimator for Partially Linear IV Model.}
\label{alg:plpiv}
\end{algorithm}
\begin{algorithm}[H]
Input: a direction $q \in \mathcal{S}^{d-1}$, estimated values $(\widehat{\eta}(X_i), \widehat{\gamma}(p,X_i))_{i=1}^N$. Estimate the following quantities:
\begin{itemize}[1]
\item[$1'$, $2'$] The treatment residual and the sample covariance matrix $$\widehat{V}_i:= D_i - \widehat{\eta}(X_i), \qquad i=1,2,\dots, N, \quad \widehat{\Sigma}:= \dfrac{1}{N} \sum_{i=1}^N \widehat{V}_i V_i' .$$
\item[$4'$] The OLS coefficient of the second-stage residual $Y_{i}(\widehat{p}(q), \widehat{\eta}) - \widehat{\gamma}(\widehat{p}(q), X_i)$ on the treatment residual $\widehat{V}_i$ as in \eqref{eq:betaq}.
\end{itemize}
\caption{Support Function Estimator for Partially Linear Model.}
\label{alg:plp}
\end{algorithm}
\subsection{Results }
\label{sec:ver}
In this section, I verify Assumptions \ref{ass:ratebound}-- \ref{ass:concentration:chap1} for Example \ref{ex:plpiv}. Then, I establish asymptotic Gaussian approximation for the Support Function Estimator.
\begin{assumption}[Bounded Width $Y_U - Y_L$]
\label{ass:boundedwidth}
The width $Y_U - Y_L$ is bounded by a finite constant $M_{UL}$ a.s., that is
$$
Y_U - Y_L \leq M_{UL} \text{ a.s. }
$$
\end{assumption}
\begin{assumption}[Regularity Conditions]
\label{ass:regularity}
(1) For every $p \in \mathcal{P}$, the residual vector $V(\eta_0) $ has a conditional density $h_{\{p' V(\eta_0) \mid X=x\}}(\cdot,x)$ that is bounded uniformly over $X$ by $M_h$. (2) The residual vector norms $\| V(\eta_0) \| =\| Z-\eta_0(X)\|$ and $ \| E(m_0) \|= \|D-m_0(X)\|$ are subGaussian random variables conditional on $X$. (3) The random variable $Y_L$ has finite conditional second moment $\sup_{x \in \mathcal{X}} \mathbb{E}[ Y^2_L \mid X=x] \leq \bar C_L$ and $\| Y_L \|_{P,4} \leq \bar C_L$ for some finite $ \bar C_L< \infty$. Finally, for some $c'>2$, $\| \| V(\eta_0) \| (Y_L - \gamma_{L,0} (X))\|_{P,c'}$ and $\| V (\eta_0) \|_{P,c'}$ are finite.
\end{assumption}
Lemma \ref{lem:powell} shows that the moment equation \eqref{eq:zqwq} incurs only a second-order bias due to the sign mistake of $V(\eta_0)$ as long as $V(\eta_0)$ is continuously distributed.
\begin{lemma}[First-Order Bias]
\label{lem:powell}
Under Assumptions \ref{ass:boundedwidth} and \ref{ass:regularity} (1), the first-order bias shrinks at the quadratic speed
\begin{align*}
\sup_{p \in \mathcal{P}} | \mathbb{E} [ z(p, \eta_0) ( Y(p, \eta) - Y(p, \eta_0) )] | \leq 2 C_P^2 M_{UL} M_h \mathbb{E} \| \eta(X) - \eta_0(X) \|^2.
\end{align*}
Likewise, the second-order bias shrinks at the quadratic speed
\begin{align*}
\sup_{p \in \mathcal{P}} | \mathbb{E} [ (z(p, \eta) - z(p, \eta_0)) ( Y(p, \eta) - Y(p, \eta_0) )] | \leq 2 C_P^2 M_{UL} M_h \mathbb{E} \| \eta(X) - \eta_0(X) \|^2.
\end{align*}
\end{lemma}
\begin{lemma}[Verification of Assumption \ref{ass:ratebound2}]
\label{cor:powell2}
Suppose Assumption \ref{ass:regularity}(2) holds. For any sequence $\ell_N \rightarrow \infty$, Assumption \ref{ass:ratebound2} is satisfied with $v_N = \sqrt{ \ell_N /N}$, in particular, one can take $\ell_N = \log N$ to ensure that $v_N = o(N^{-1/2} \log N)$.
\end{lemma}
Lemma \ref{cor:powell} verifies the Assumption \ref{ass:ratebound} for Example \ref{ex:plpiv}. Define the mean square rates for the expectation functions
\begin{align*}
\sup_{\gamma_L \in \Gamma_{L,N}} \left(\mathbb{E} (\gamma_{L} (D,X) - \gamma_{L,0}(D,X))^2 \right)^{1/2} =: \gamma_{L,N} \\
\sup_{\gamma_{UL} \in \Gamma_{UL,N}} \left(\mathbb{E} (\gamma_{UL} (D,X) - \gamma_{UL,0}(D,X))^2 \right)^{1/2} =: \gamma_{UL,N}
\end{align*}
\begin{lemma}[Verification of Assumptions \ref{ass:ratebound}--\ref{ass:ratebound2}]
\label{cor:powell}
Let $\eta_0(X)$ be as in \eqref{eq:eta0}, $V(\eta)$ be as in \eqref{eq:plpivspec}, the matrix function $A(W, \eta, m)$ be as in \eqref{eq:plpivspec2} and the orthogonal moment function be as in \eqref{eq:orthomom}. Suppose Assumptions \ref{ass:boundedwidth} and \ref{ass:regularity} and \ref{ass:suffcond} hold. Furthermore, suppose (1) the elements of $\Gamma_{L,N}$ are bounded by finite constant $M_{\gamma}$ in the sup-norm: $\sup_{ \gamma_L \in \Gamma_{L,N} } \sup_{x \in \mathcal{X}} | \gamma_L(x) | \leq M_{\gamma}$ and $\sup_{ \gamma_{UL} \in \Gamma_{UL,N} } \sup_{x \in \mathcal{X}} | \gamma_{UL}(x) | \leq M_{UL\gamma}$ (2) the elements of $\mathcal{M}_N$ and $\mathcal{T}_N$ are bounded by finite constant $M_{\eta}< \infty$ in the sup-norm: $$\sup_{ m \in \mathcal{M}_N } \sup_{x \in \mathcal{X}} \| m(x) \| \leq M_{\eta} , \quad \sup_{ \eta \in \mathcal{T}_N } \sup_{x \in \mathcal{X}} \| \eta(x) \| \leq M_{\eta}. $$ Then, Assumption \ref{ass:ratebound2}(1) holds. Furthermore, the bias rates in Definition \ref{def:nearorthog} can be bounded as follows with $\gamma^2_N := 2(\gamma^2_{L,N} + \gamma^2_{UL, N})$
\begin{align*}
\mu_N &\leq 4 C^2_P M_{UL} M_h \eta^2_N + C_P \eta_N \cdot \gamma_N \\
A_N &\leq \eta_N \cdot m_N.
\end{align*}
In addition, for $\tau_N := N^{-1/2} \log N$, the sequences $r_N'', r_N', \delta_N$ obey
$$
r_N'' = O (\eta_N + \gamma_N), \quad r_N' = O ( N^{-1/4} \log^{1/2} N), \quad \delta_N = O (\eta_N + m_N).
$$
Thus, if $\eta_N = o(N^{-1/4})$ and $ (m_N + \gamma_N) \cdot \eta_N = o(N^{-1/2})$ and $m_N = o(1)$ and $\gamma_N = o(1)$, Assumption \ref{ass:ratebound} holds.
\end{lemma}
Combining the statements in Lemmas \ref{cor:powell2}--\ref{cor:powell}, I obtain the following corollary.
\begin{corollary}[Asymptotic Theory for Partially Linear IV Model with Interval-Valued Outcome]
\label{cor:plpiv}
Suppose Assumptions \ref{ass:sigma}, \ref{ass:jacobian} and \ref{ass:suffcond} and the conditions of Lemma \ref{cor:powell} hold with the fast enough first-stage rates $\eta_N, m_N$ and $\gamma_N:= 2(\gamma_{L,N} + \gamma_{UL, N})$ such that Assumption \ref{ass:ratebound} holds. Then, Theorems \ref{thm:limit} and \ref{thm:bb} and Corollaries \ref{cor:limit} and \ref{cor:bb} hold for the Support Function Estimator of the Algorithm \ref{alg:plpiv} with the first-stage fitted values \eqref{eq:crossfitgammal} and with the influence function $h(W,q)$ equal to \eqref{eq:hwq}.
\end{corollary}
\begin{corollary}[Asymptotic Theory for Partially Linear IV Model with Interval-Valued Outcome]
\label{cor:plpiv1}
Suppose Assumptions \ref{ass:sigma}, \ref{ass:jacobian} and the conditions of Lemma \ref{lem:riesz} and Lemma \ref{cor:powell} hold with the fast enough first-stage rates $\eta_N, m_N$ and $\gamma_N:= \gamma_{L,N}$ such that Assumption \ref{ass:ratebound} holds. Then, Theorems \ref{thm:limit} and \ref{thm:bb} and Corollaries \ref{cor:limit} and \ref{cor:bb} hold for the Support Function Estimator of the Algorithm \ref{alg:plpiv} with the influence function $h(W,q)$ equal to \eqref{eq:hwq}, where the first-stage fitted values are as in \eqref{eq:finalfitted} (assuming Lemma \ref{lem:riesz} holds with $\bar{\sigma}_N = o(N^{-1/4})$ ) or as in \eqref{eq:finalfittedlogis} (assuming Lemma \ref{lem:rieszlogistic} holds with $\bar{\sigma}_N = o(N^{-1/4})$).
\end{corollary}
\begin{remark}
\label{rm:sharp}
Consider Example \ref{ex:plp}. Note that the short least squares regression \eqref{eq:plppoint1} uses only $d$ out of (infinitely) many restrictions implied by the exogeneity restriction \eqref{ex:plp}. Therefore, the identified set
$$
\mathcal{B}_1:=\bigg\{ b \in \mathrm{R}^d: \exists f_b \in \mathcal{L} \text{ and } Y \in [Y_L, Y_U]: \quad \mathbb{E}[ Y - D b - f_b(X) \mid D, X] =0 \bigg\}
$$
is a convex subset of $\mathcal{B}$, but may not coincide with $\mathcal{B}$.
\end{remark}
\section{Average Partial Derivative}
\label{sec:apd}
In this section, I present the identification, estimation and inference results for Example \ref{ex:apd}.
\setcounter{assumption}{0}
\setcounter{lemma}{0}
\setcounter{corollary}{0}
\setcounter{example}{0}
\setcounter{definition}{0}
Assumption \ref{ass:regcond:apd} states the sufficient conditions for compactness and convexity of the identified set $\mathcal{B}$.
\begin{assumption}[Regularity Conditions for Average Partial Derivative]
\label{ass:regcond:apd}
The following conditions hold. (1) There exists a compact and convex set $\mathcal{D}$ with nonempty interior containing the support of $D$, such that $f_0 (d \mid X) = 0$ on the boundary of $\mathcal{D}$ a.s. in $X$. (2) The random variables $D$ and $f_0(D \mid X)$ and $\nabla f_0 (D \mid X)$ have the density conditional on $X$. (3) The random vector $$\eta_0(D,X) = \nabla f_0 (D \mid X)/ f_0 (D \mid X)$$ is $L_{P,2}$-integrable, that is, there exists a finite constant $C_{\text{APD}}< \infty$ such that
$$
\| \| V(\eta_0) \| \|_{P,2} := \| \| \nabla f_0 (D \mid X)/ f_0 (D \mid X) \| \|_{P,2} \leq C_{\text{APD}}.
$$
(4) The functions $\gamma_{L,0}(D,X)$ and $\gamma_{UL,0}(D, X)$ are bounded by some finite constant $C_{\gamma}< \infty$ on the support of $D$ and $X$. (5) Furthermore, $\gamma_{L,0}(D,X)$ and $\gamma_{UL,0}(D, X)$ are continuously differentiable w.r.t to $D$ with bounded derivatives a.s. in $D$ and $X$. (5) The random variables $\| V(\eta_0) \|$ and $
\| V(\eta_0) \| (Y_L - \gamma_{L,0}(D,X))
$ are $L_{P,2}$-integrable and $L_{P,c'}$-integrable for $c'>2$.
\end{assumption}
\begin{lemma}
Suppose Assumption \ref{ass:regcond:apd} (1)-(4) hold. Then, (a) the identified set $\mathcal{B}$ for $\beta_0$ is compact and convex and (b) the support function of $\mathcal{B}$ is given by \eqref{eq:zqwq}.
In addition, if Assumption \ref{ass:jacobian} holds, which implies that $\mathcal{B}$ is strictly convex.
\label{lem:apd}
\end{lemma}
Lemma \ref{lem:apd} extends the Theorem 2.1 of \cite{Kaido} to allow the density $D$ to be conditioned on $X$.
\begin{assumption}[Margin Condition]
\label{ass:margin}
(1) There exists an absolute constant $\bar C_f <\infty$ so that, in some neighborhood $(0, \bar t)$ of zero,
\begin{align*}
\sup_{q \in \mathcal{S}^{d-1}} {\mathrm{P}} \left( |q' \nabla_D f_0(D \mid X)/ f_0 (D \mid X) | \leq \delta \right) \leq \bar C_f \delta , \quad \delta \in (0, \bar{t}).
\end{align*}
\end{assumption}
The margin condition is commonly used in classification analysis (\cite{MammenTsybakov}, \cite{Tsybakov}) and empirical welfare maximization (\cite{KitagawaTetenov}, \cite{MbakopTabord}) and bounds \cite{SemSupp2}. This paper is the first one to introduce it for the study of average partial derivative with an interval-valued variable.
\begin{definition}[Worst-Case Rates]
Let $f_0(D \mid X)$ and $\nabla_D f_0(D \mid X)$ be the conditional density and its derivative. Let $\{ F_N, \quad N \geq 1 \}$ and $\{ F^1_N, \quad N \geq 1 \}$ be sequences of realization sets of the estimates of $f_0 (D \mid X)$ and $\nabla_D f_0 (D \mid X)$, respectively. Assume that these sequences shrink at the following worst-case rates $f^{\infty}_N$ and $f^{\infty}_{1,N}$
\begin{align*}
\sup_{ f \in F_N } \sup_{d, x} | f (d \mid x) - f_0( d \mid x) | \leq f^{\infty}_N \\
\sup_{ f \in F^1_N } \sup_{d, x} \| \nabla_{d} f (d \mid x) - \nabla_d f_0( d \mid x) \| \leq f^{\infty}_{1,N}
\end{align*}
Define
$$
\eta^{\infty}_N:= f^{\infty}_{1,N} + f^{\infty}_N.
$$
Furthermore, assume that the elements of $F_N$ are bounded by some $\bar B_f< \infty$:
$$
\sup_{ f \in F_N } \sup_{d, x} \{ | f^{-1} (d \mid x ) |, | f (d \mid x ) |\} \leq \bar B_f.
$$
and
$$
\sup_{ f \in F^1_N } \sup_{d, x} \| \nabla_{d} f (d \mid x) \| \leq \bar B_f.
$$
\end{definition}
\subsection{The Algorithm}
\begin{algorithm}[H]
\begin{algorithmic}[1]
\STATE The best-case outcome as a function of $q$ $$Y_{i}(q, \widehat{\eta}) :=Y_{L,i} + (Y_{U,i}-Y_{L,i})1\{ -q' \widehat{\eta}(D_i, X_i)>0 \}, \quad i=1,2,\dots, N. $$
\STATE The expectation function of $Y(q, \eta)$ and its derivative
\begin{align*}
\widehat \mu (q, D_i, X_i) = \widehat{ \gamma}_L (D_i, X_i) + \widehat{ \gamma}_{UL} (D_i, X_i) 1\{ -q' \widehat{\eta}(D_i, X_i)>0 \} \\
\nabla_D \widehat \mu (q, D_i, X_i) = \nabla_D \widehat{ \gamma}_L (D_i, X_i) + \nabla_D \widehat{ \gamma}_{UL} (D_i, X_i) 1\{ -q' \widehat{\eta}(D_i, X_i)>0 \}
\end{align*}
\STATE The sample moment
\begin{align}
\label{eq:betaqapd}
\widehat{\beta}_q = \dfrac{1}{N} \sum_{i=1}^N -\widehat{\eta}(D_i, X_i) [ Y_{i}(q, \widehat{\eta}) - \widehat{\mu} (q,D_i, X_i)] + \nabla_D \widehat{\mu} (q,D_i, X_i) \end{align}
\STATE Report: the projection of $\widehat{\beta}_q$ onto the direction $q$: $$\widehat{\sigma}(q) = q' \widehat{\beta}_q.$$
\end{algorithmic}
\caption{Support Function Estimator for Average Partial Derivative.}
\label{alg:apd}
\end{algorithm}
\subsection{Results}
Lemma \ref{lem:bias2} shows that the moment equation \eqref{eq:zqwq} incurs only a second-order bias due to the sign mistake of $V(\eta_0)$ as long as $V(\eta_0)$ is continuously distributed. In contrast to the setup of Lemma \ref{lem:powell}, the estimation error is not orthogonal to the space of $V(\eta_0)$. As a result, the first-order bias is bounded by invoking the margin assumption and the $\ell_{\infty}$-rate.
\begin{lemma}[First-Order Bias]
\label{lem:bias2}
Suppose Assumptions \ref{ass:boundedwidth} and \ref{ass:margin} holds. Then, the first-order bias shrinks at quadratic speed, that is, for some constant $\bar C$ large enough,
\begin{align*}
\sup_{ \eta \in \mathcal{T}_N } | \mathbb{E} [z(q, \eta_0) (Y(q, \eta) - Y(q, \eta_0)) ] | = O ((\eta^{\infty}_N)^2).
\end{align*}
Likewise, the second-order bias shrinks at quadratic speed
\begin{align*}
\sup_{ \eta \in \mathcal{T}_N } |\mathbb{E} [ (z(q, \eta) - z(q, \eta_0)) (Y(q, \eta) - Y(q, \eta_0)) ] | = O ((\eta^{\infty}_N)^2).
\end{align*}
\end{lemma}
\begin{definition}[Mean Square Rates]
\label{def:sqrate}
Define the mean square rates for the expectation functions
\begin{align*}
\sup_{\gamma_L \in \Gamma_{L,N}} \left(\mathbb{E} (\gamma_{L} (D,X) - \gamma_{L,0}(D,X))^2 \right)^{1/2} =: \gamma_{L,N} \\
\sup_{\gamma_{UL} \in \Gamma_{UL,N}} \left(\mathbb{E} (\gamma_{UL} (D,X) - \gamma_{UL,0}(D,X))^2 \right)^{1/2} =: \gamma_{UL,N}
\end{align*}
and for their derivatives
\begin{align*}
\sup_{\gamma_L \in \Gamma_{L,N}} \left(\mathbb{E} \| \nabla_D \gamma_{L} (X) -\nabla_D \gamma_{L,0}(X) \|^2 \right)^{1/2} =: \gamma^1_{L,N} \\
\sup_{\gamma_{UL} \in \Gamma_{UL,N}} \left(\mathbb{E} \| \nabla_D \gamma_{UL} (X) -\nabla_D \gamma_{UL,0}(X) \|^2 \right)^{1/2} =: \gamma^1_{UL,N}
\end{align*}
and let $$
\gamma_N := \gamma_{L,N} + \gamma_{UL,N} + \gamma^1_{L,N} +\gamma^1_{UL,N}.
$$
\end{definition}
\begin{corollary}[Asymptotic Theory for Average Partial Derivative with an Interval-Valued Outcome]
\label{cor:apd}
Suppose Assumptions \ref{ass:regcond:apd} and \ref{ass:boundedwidth} and \ref{ass:margin} hold. In addition, suppose $| \gamma_{L} (D,X) | + \| \nabla_D \gamma_{L} (D,X) \| + | \gamma_{UL} (D,X) | + \| \nabla_D \gamma_{UL} (D,X) \|< M_{\gamma} \text{ a.s.}$. Then, the sequences $\mu_N$ and $r_N'$ can be bounded as follows
$$
\mu_N = O ( (\eta^{\infty}_N)^{3/2}+ \eta^{\infty}_N \cdot \gamma_N).
$$
and $r_N' = O (\gamma_N + (\eta^{\infty}_N)^{1/2} ) $. In particular, if $\eta^{\infty}_N = o (N^{-1/3})$ and $\eta^{\infty}_N \cdot \gamma_N = o (N^{-1/2})$ and $\gamma_N = o(1)$, Assumption \ref{ass:ratebound} holds. Assumption \ref{ass:ratebound2} holds automatically since $A(W, \eta_0) = \Sigma = I_d$ is a known matrix. Finally, Assumption \ref{ass:concentration:chap1} holds. Then, Theorems \ref{thm:limit} and \ref{thm:bb} hold for the Support Function Estimator of Algorithm \ref{alg:apd} and the influence function equal to
\begin{align*}
h(W,q) = g(W,q,\xi(q)) - \mathbb{E} [ g(W,q,\xi(q)) ],
\end{align*}
where $g(W,q,\xi(q))$ is given in (\ref{eq:rho:apd}).
\end{corollary}
\section{Simulation Study}
\label{sec:montecarlo}
In this section, I compare the performance of the classic (series-based non-orthogonal) and the proposed (lasso-based orthogonal) approaches in a moderate-dimensional sparse design. The first approach is to plug a least squares series first-stage estimator into the non-orthogonal moment equation \eqref{eq:zqwq}. The second one is to plug an $\ell_1$-regularized least squares series into the orthogonal moment equation \eqref{eq:orthomom}. Regularization helps to leverage the sparsity assumption and to reduce risk.
Consider the partially linear model of Example \ref{ex:plp} with $d=2$ treatments. The outcome equation is generated as
\begin{align}
Y &=D_1 \beta_1 + D_2 \beta_2 + X \cdot (c_{\theta} \theta_0) + U, \label{eq:theta0}
\end{align}
where $D=(D_1, D_2)$ is the treatment vector, $X \in \mathrm{R}^{p_X}$ is the $p_X$-covariate vector with $p_X=50$, and $U \sim N(0, \sigma_U^2)$ is a normal shock independent of $D$ and $X$ with $\sigma_U =1$. For each treatment $m=1,2$, its reduced form is
\begin{align}
D_m &= X\cdot (c_{D_m} \alpha_m) + V_m, \label{eq:alpha1} \quad m=1,2,
\end{align}
where the coefficients $$\theta_0 = \alpha_1 = \alpha_2 = (1, 1/2^2, \dots, 1/j^2, \dots, 1/p_X^2). $$ The parameters in \eqref{eq:theta0} and \eqref{eq:alpha1} are chosen as
\begin{align*}
c_D& = (c_{D_1}, c_{D_2}) = (2,1), c_{\theta}=1, \beta_0 = (1,1).
\end{align*}
The covariates $X$ are generated from $N(0, \Omega)$, where $\Omega$ is a Toeplitz matrix with correlation coefficient $\rho = 0.5$. That is, for every $(i,j) \in \{1,2,\dots, p_X\}^2$, $\Omega_{ij} = \rho^{|i-j|}$. The vector $V$ is independent of $(X,U)$ and is drawn from the bivariate normal distribution
\begin{align}
\label{eq:vee}
V = (V_1, V_2) \sim N(0, \Sigma), \quad \Sigma = \begin{pmatrix} 1 & 0.5 \\ 0.5 & 1 \end{pmatrix}.
\end{align}
The outcome $Y$ is not included into the data. Instead, the support of $Y$ is partitioned into the bins $ \cup_{s=1}^{S} [b_s, b_{s+1})$ of width $\Delta$:
$$
b_{s+1} = b_s + \Delta, \quad s=1, 2, \dots, S-1,
$$
where $b_0 = -\infty$ and $b_{S+1} = \infty$. The observed bounds $Y_L$ and $Y_U$ are taken to be
\begin{align*}
[Y_L, Y_U]:&= \sum_{s=1}^{S} [b_s, b_{s+1}) \cdot 1 \{ Y \in [b_s, b_{s+1}) \}.
\end{align*}
Thus, the observed data vector $W=(X, D, Y_L, Y_U)$ but does not contain $Y$.
I now derive the true (population) support function. Plugging $Y_U - Y_L = \Delta$ into \eqref{eq:zqwq} gives
$$
\sigma(q) = q' \Sigma^{-1} \mathbb{E} V(\eta_0) Y(q, \eta_0) = q' \Sigma^{-1} \mathbb{E} V Y_L + \Delta \mathbb{E} \max (q' \Sigma^{-1} V,0).
$$
Invoking \eqref{eq:vee} gives $q' \Sigma^{-1} V \sim N(0, q' \Sigma^{-1} \Sigma \Sigma^{-1} q) = \sqrt{q' \Sigma^{-1} q} N(0, 1)$, which implies
$$
\mathbb{E} \max (q' V, 0) = \sqrt{q' \Sigma^{-1} q/2\pi}.
$$
Thus, the support function is
\begin{align}
\label{eq:sigmatrue}
\sigma(q) = q'\kappa_0 + \Delta \sqrt{q' \Sigma^{-1} q/2 \pi} , \quad \kappa_0 := \Sigma^{-1} \mathbb{E} V Y_L.
\end{align}
The classic and the proposed approaches are implemented via Algorithm \ref{alg:plpiv} with different sets of the first-stage fitted values. For each treatment $m \in \{1,2\}$, the first-stage regression parameter $\alpha$ is estimated as
\begin{align}
\label{eq:lasso}
\widehat{\alpha}_{m} :&= \arg \min_{b \in \mathrm{R}^{p_X}} \dfrac{1}{N} \sum_{i =1}^N (D_{im} - X_i' b)^2 + \lambda_{D_m} \| b \|_1,
\end{align}
where $\lambda_{D_m} = 0$ in the series-based case and $\lambda_{D_m} > 0$ in the lasso-based case. For the lasso estimator, the penalty parameter $\lambda_{D_m}$ is chosen according to Algorithm 1 in \cite{Program} (i.e., the default value of \url{rlasso} package). In both cases, the treatment fitted values are
\begin{align*}
\widehat{\eta}(X)&= X' (\widehat{\alpha}_1, \widehat{\alpha}_2).
\end{align*}
In the orthogonal case (the proposed approach), the best-case outcome fitted values are
\begin{align*}
\widehat{\gamma}_U(X)&= X' \widehat{\gamma}_L + \dfrac{1}{2} \Delta,
\end{align*}
where $\gamma_L$ is estimated by Lasso regression of $Y_L$ on $X$ similarly to \eqref{eq:lasso}. Since the classic case does not require partialling out,
the fitted values $\widehat{\gamma}_U(X)$ are set to zero.
The estimator's performance is summarized in terms of its risk and coverage. The total risk of the estimator is defined as
\begin{align}
\label{eq:total}
R_H := \sup_{q \in \mathcal{S}^{d-1}} | \widehat{\sigma}(q) - \sigma(q) |.
\end{align}
It coincides with the Hausdorff distance $d(\widehat{\mathcal{B}},\mathcal{B})$ between the estimated ($\widehat{\mathcal{B}}$) and the true ($\mathcal{B}$) sets. Furthermore, the outer and the inner risks are defined as
\begin{align}
\label{eq:outerinner}
R_O:= \sup_{q \in \mathcal{S}^{d-1}} \max (\widehat{\sigma}(q) - \sigma(q), 0), \quad R_I:= \sup_{q \in \mathcal{S}^{d-1}} \max ( \sigma(q) - \widehat{\sigma}(q), 0).
\end{align}
The rejection frequency is the share of rejected simulation draws
\begin{align}
\label{eq:boot}
\dfrac{1}{N_S} \sum_{s=1}^{N_S} 1\{ R_H^s > c^{*}_{1- \alpha} \},
\end{align}
where $c^{*}_{1- \alpha} $ is the $(1-\alpha)$-quantile of $R_H^b$ of the bootstrap process
\begin{align}
\label{eq:totalboot}
R^b_H := \sup_{q \in \mathcal{S}^{d-1}} | \widehat{\sigma}^b(q) - \widehat{\sigma}(q) |.
\end{align}
Table \ref{tab:sims} compares the classic and the proposed estimators in terms of risk and coverage. Across the board, the proposed estimator has a smaller total risk than the classic one, by a factor of ranging from $2.0$ for $\Delta=1$ to $3.0$ for $\Delta=3$. The classic estimator has higher risk due to the excessive noise of series estimators $\widehat{\alpha}_1$ and $\widehat{\alpha}_2$ in a regime with $p_X =50$ covariates.
I investigate the importance of the first-stage regularization and second-stage orthogonalization, applying one at a time. Table \ref{tab:simsapp} compares the non-orthogonal (Columns (1)-(4)) and the orthogonal (Columns (5)-(8)) estimators based on the true treatment first stage. Across the board, the orthogonal estimator has smaller total risk than the non-orthogonal one, by a factor of $2.5$ on average. The variance reduction occurs due to partialling out the relevant controls from the outcome $Y^{\text{best}}(\eta_0)$. In contrast, the risks of series-based non-orthogonal (Table \ref{tab:sims}, Column (3)) and orthogonal (Table \ref{tab:simsapp}, Column (7)) estimators are close to each other. In particular, the high risk of the series-based estimator cannot be improved by orthogonalization.
Next, I compare the ortho (Table \ref{tab:sims}, Columns 5--8) and non-ortho (Table \ref{tab:simsapp2}, Columns 1--4) estimators based on the lasso-based first stage.
Across the board, the risk of the ortho version is substantially smaller than the non-ortho one, by a factor ranging from $2$ to $5.5$. The non-ortho version has a higher risk, because the estimates of non-zero coefficients in $\alpha_m$ are shrunk to zero. While the shrinkage bias could be reduced by invoking the post-lasso instead of the lasso estimator, post-single-selection inference may not be robust to moderate deviations from zero.
\section{Empirical illustration}
\label{sec:empirical}
This section demonstrates the proposed approach by estimating the gender wage gap with a bracketed wage variable. First, I show that a frequent empirical practice -- midpoint regression -- gives biased results. Instead of this approach, I propose reporting identified set (i.e., the lower and the upper bound) for the parameter of interest and demonstrate how to estimate the set.
The sample for the analysis comes from the U.S. March Supplement of the Current Population Survey (CPS) in 2015, as studied in \cite{MulliganRubinstein} and \cite{CCMS}. The selected sample consists of white non-hispanic individuals, aged 25 to 64 years, working more than 35 hours per week during at least 50 weeks of the year. The resulting sample comprises $32,523$ workers, including $18,137$ men and $14,386$ women. The object of interest is the gender wage gap -- the average difference in log wages between men and women after controlling for the observed characteristics. The data set is augmented by the lower bound $Y_L$ and the upper bound $Y_U$ defined as
\begin{align*}
[Y_L, Y_U]:&= \sum_{s=1}^{S} [b_s, b_{s+1}) \cdot 1 \{ Y \in [b_s, b_{s+1}) \},
\end{align*}
where $b_0=1$ and $$ b_{s+1} - b_s = \Delta \quad s=1,2,\dots, S.$$
I assume that the log wage variable follows the partially linear regression of Example \ref{ex:plp}
$$
Y = D' \beta_0 + X' \gamma_0+ U, \quad \mathbb{E}[ U \mid X, D] =0.
$$
Here, the outcome variable $Y$ is the logarithm of the hourly wage rate, the treatment/policy variable $D$ is an indicator for female gender, and the vector of $p_X=260$ controls includes demographic indicators, region, and experience indicators, as well as their interactions. The coefficients $\beta_0$ and $\gamma_0$ are the target and the nuisance parameters, respectively. The parameter $\beta_0$ is identified as a minimizer of the least squares loss function
$$
\beta_0:=\arg \min_{b \in \mathrm{R}} \mathbb{E} (Y - (D- \eta_0(X))b)^2.
$$
When $Y$ is unobserved, a frequent approach is to replace $Y$ by the bracket midpoint
$$Y_M:=(Y_L+Y_U)/2.$$
The ``mid-point'' regression parameter is taken to be
$$
\beta_{M}:=\arg \min_{b \in \mathrm{R}} \mathbb{E} (Y_M - (D- \eta_0(X))b)^2,
$$
which can differ from $\beta_0$. I report the estimate and the $95 \%$ CI of $\beta_0$ and $\beta_M$. The first-stage treatment and the outcome expectation functions are estimated via logistic and linear Lasso regression of \cite{Program} implemented in the hdm $R$ package, respectively. In addition, I also consider the random forest estimator as implemented in the ranger $R$ package. The final estimator is taken to be the Double Machine Learning estimator of \cite{chernozhukov2016double} with $K=2$-fold cross-fitting. The Lasso-based first-stage fitted values (Columns (1)-(2)) and the random-forest-based (Columns (3)-(4)).
Instead of reporting $\beta_M$ -- which is a biased measure of $\beta_0$ -- I propose reporting the lower and the upper bound on $\beta_0$. The moment equation for $\beta_U$ is given in \eqref{eq:betu1d} and its estimate is defined in the Algorithm \ref{alg:plp} with the fitted values described below. The symmetry-based specification (SYM) imposes Assumption \ref{ass:suffcond}, and the first-stage fitted values are taken to be \eqref{eq:crossfitgammal}. The sparsity-based specification (SPRS) imposes the model \eqref{eq:maindecomp}, and the first-stage are taken to be as in Example \ref{ex:lassologistic}, \eqref{eq:finalfittedlogis}. The treatment expectation function $\eta_0(\cdot)$ is estimated the same as in the point-identified specifications.
Table \ref{tab:empapp} summarizes the findings for the bracket width $\Delta \in \{1,2,3\}$. The ``ground-truth'' gender wage gap ranges between $18 \%$ (Columns (3)-(4)) and $20 \%$ (Columns (1)-(2)). I call this estimate ``ground-truth'' because this is the estimate to be reported if the wage $Y$ was observed. When the bracket width $\Delta=1$ is small, the midpoint estimate $\beta_M$ is close to the estimate of $\beta_0$. When $\Delta=2$, the midpoint estimate $\beta_M$ ranges between $25\%$ (RF) and $28 \%$ (Lasso). Furthermore, the 95 $\%$ CI for $\beta_M$ does not contain the ``ground-truth`` estimate of $\beta_0$ for either RF or Lasso first-stage method. When $\Delta=3$, the magnitude of the bias remains substantial, which speaks against midpoint regression for large values of bracket width.
To estimate $[\beta_L, \beta_U]$, I consider four specifications: SYM-Lasso, SYM-RF, SPRS-Lasso and SPRS-RF. Across the board, the bounds contain $\beta_0$. Furthermore, the specifications yield close results despite being based on very different starting assumptions. That said, the bounds $[\beta_L, \beta_U]$ come out wide. As discussed in Remark \ref{rm:sharp}, the bounds $[\beta_L, \beta_U]$ may not be sharp for $\beta_0$, since they are utilize only $d$ (out of infinitely many) moment restriction implied by \eqref{ex:plp}. The derivation of sharp bounds for $\beta_0$ is left for the future work.
\begin{table}
\centering
\caption{Finite-sample performance of the classic (series-based) and the proposed (Lasso-based) methods}
\begin{tabular}{c|cccccccc}
\toprule
& \multicolumn{4}{c}{Series-based}& \multicolumn{4}{c}{Lasso-based} \\
\hline \\
$N$ & Total & Outer & Inner & Rej.freq & Total & Outer & Inner & Rej.freq \\
\\
& \multicolumn{8}{c}{Panel A: Bracket width $\Delta = 1$}\\
\\
$250$ & 0.24 & 0.18 & 0.23 & 0.03 & 0.13 & 0.11 & 0.11 & 0.09 \\
$500$ & 0.19 & 0.09 & 0.19 & 0.25 & 0.09 & 0.08 & 0.08 & 0.08 \\
$700$ & 0.19 & 0.07 & 0.19 & 0.43 & 0.08 & 0.07 & 0.07 & 0.04 \\
$1,000$ & 0.18 & 0.06 & 0.18 & 0.99 & 0.07 & 0.06 & 0.06 & 0.10 \\
\\
& \multicolumn{8}{c}{Panel B: Bracket width $\Delta = 2$}\\
\\
$250$ & 0.35 & 0.24 & 0.34 & 0.05 & 0.17 & 0.14 & 0.13 & 0.09 \\
$500$ & 0.33 & 0.13 & 0.33 & 0.85 & 0.12 & 0.10 & 0.09 & 0.09 \\
$700$ & 0.32 & 0.09 & 0.32 & 0.99 & 0.10 & 0.08 & 0.08 & 0.03 \\
$1,000$ & 0.32 & 0.07 & 0.32 & 1.00 & 0.09 & 0.07 & 0.07 & 0.10 \\
\\
& \multicolumn{8}{c}{Panel C: Bracket width $\Delta = 3$}\\
\\
$250$ & 0.48 & 0.32 & 0.47 & 0.10 & 0.21 & 0.17 & 0.16 & 0.07 \\
$500$ & 0.46 & 0.16 & 0.46 & 0.98 & 0.15 & 0.12 & 0.12 & 0.07 \\
$700$ & 0.46 & 0.12 & 0.46 & 1.00 & 0.13 & 0.11 & 0.10 & 0.04 \\
$1,000$ & 0.46 & 0.09 & 0.46 & 1.00 & 0.11 & 0.09 & 0.09 & 0.06 \\
\bottomrule
\end{tabular}
\label{tab:sims}
\caption*{Notes. Results are based on 10, 000 simulation runs. Panels A, B and C correspond to the bin width $\Delta = 1, 2, 3$. Table shows the total risk \eqref{eq:total}, the outer and inner risks \eqref{eq:outerinner}, and the rejection frequency \eqref{eq:boot} for the nominal size $\alpha = 0.05$. The supremum over $\mathcal{S}^1$ is approximated by the maximum over the grid consisting of $ 50$ evenly spaced points on unit circumference $\mathcal{S}^1$. Columns (1--4) and (5--8) correspond to the classic and the proposed approach. The number of bootstrap repetitions $B=2, 000$. The true support function $\sigma(q)$ is in \eqref{eq:sigmatrue}. For the description of estimators, see text. } \end{table}
\newpage
\begin{table}[H]
\centering
\small
\caption{Bounds on gender wage gap with bracketed log wage}
\begin{tabular}{c|cc|cc|}
\toprule
& \multicolumn{2}{c}{Lasso}& \multicolumn{2}{c}{RF} \\
& Estimated Set & 95 $\%$ CI & Estimated Set & 95 $\%$ CI \\
\\
& \multicolumn{4}{c}{Panel A: Bracket width $\Delta = 1$}\\
\\
True $\beta_0$ & -0.200 & (-0.213, -0.186) & -0.180 & (-0.194, -0.167) \\
$\beta_M$ & -0.199 & (-0.215, -0.184) & -0.181 & (-0.196, -0.165) \\
$[\beta_L, \beta_U]$ (SYM) & [ -1.195, 0.796] & (-1.211, 0.813) & [-1.133, 0.771] & (-1.150, 0.788) \\
$[\beta_L, \beta_U]$ (SPRS) & [-1.193, 0.794] & (-1.209, 0.810) & [-1.115, 0.754] & (-1.131, 0.770) \\
\\
& \multicolumn{4}{c}{Panel B: Bracket width $\Delta = 2$}\\
\\
True $\beta_0$ & -0.200 & (-0.213, -0.186) & -0.180 & (-0.194, -0.167) \\
$\beta_M$ &-0.279 & (-0.302, -0.256) & -0.251 & (-0.273, -0.228) \\
$[\beta_L, \beta_U]$ (SYM) &[-2.270, { }1.712] & (-2.295, 1.737) & [-2.155, 1.654] & (-2.180, 1.679) \\
$[\beta_L, \beta_U]$ (SPRS) &[-2.266 { }1.708] & (-2.290, 1.731) & [-2.120, 1.618] & (-2.144, 1.643) \\
\\
& \multicolumn{4}{c}{Panel C: Bracket width $\Delta = 3$}\\
\\
True $\beta_0$ & -0.200 & (-0.213, -0.186) & -0.180 & (-0.194, -0.167) \\
$\beta_M$ & -0.158 & (-0.179, -0.137) & -0.131 & (-0.152, -0.110) \\
$[\beta_L, \beta_U]$ (SYM) & [-3.145, {} 2.828] & (-3.172, 2.854) & [-2.988, 2.725] & (-3.015, 2.752) \\
$[\beta_L, \beta_U]$ (SPRS) & [-3.139, {}2.822] & (-3.162, 2.844) & [-2.935, 2.672] & (-2.959, 2.697) \\
\bottomrule
\end{tabular}
\label{tab:empapp}
\caption*{Notes. Estimated parameter (square brackets) and the $95\%$ confidence bands (parentheses) for the parameter. The true parameter $\beta_0$ is based on the observed outcome $Y$. The midpoint parameter $\beta_M$ is based on the midpoint outcome $Y_M$. The upper bound $\beta_U$ is as defined in \eqref{eq:betu1d}, and $\beta_L$ is its analog. The first-stage treatment expectation function is estimated by linear Lasso (Columns (1)-(2)) and random forest (Columns (3)-(4)). The symmetry-based specification (SYM, Row 3) is based on Assumption \ref{ass:suffcond}, and the first-stage fitted values are given in
\eqref{eq:crossfitgammal}. The sparsity-based specification (SPRS, Row 4) is based on Example \ref{ex:lassologistic}, and the first-stage fitted values are given in \eqref{eq:finalfittedlogis}. For more details, see text. }
\end{table}
\newpage
\section{Proofs}
\label{sec:proofs}
\paragraph{Empirical process notation. } Let $\widehat{\eta}_k$ and $\widehat{\xi}_k, \quad k=1,2,\dots, K$ be as in Definition \ref{sampling}. Define an event
\begin{align*}
\mathcal{E}_N &:= \{ \widehat{\eta}_k, (\widehat{\xi}_k(p))_{p \in \mathcal{P}} \in \Xi_N \quad \forall k=1,2,\dots,K \}.
\end{align*} By union bound, this event holds with probability approaching one $${\mathrm{P}} (\mathcal{E}_N ) \geq 1- K \epsilon_N = 1-o(1).$$
For a given partition $k$ in $\{1,2, \dots, K\}$, define the partition-specific averages
\begin{align*}
{\mathbb{E}_{n,k}} f(W_i) :&= \dfrac{1}{n} \sum_{i \in J_k} f(W_i), \\
{\mathbb{G}_{n,k}} f(W_i) :&= \dfrac{1}{\sqrt{n}} \sum_{i \in J_k} [f(W_i) - \int f(w) dP (w)].
\end{align*}
Define the function $\psi_0(p)$
\begin{align}
\psi_0(p) &= \psi(p, \xi_0) = \mathbb{E} [ g(W, p, \xi_0(p))] \label{eq:psi0}
\end{align}
and observe that plugging $p_0(q)$ into $\psi_0(p)$ gives the support function $ \sigma(q) $ at $q$:
$$ \psi_0(p_0(q))= \psi_0(\Sigma^{-1}q)=\sigma(q).$$
For $i \in J_k$ and $k=1,2,\dots, K$, define the partition-specific conditional expectation
\begin{align}
\label{eq:psi}
\psi(p, \widehat{\xi}_k)&:= \mathbb{E} [ g (W_i, p, \widehat{\xi}_k(p)) \mid (W_i)_{i \in J_k^c} ], \quad i \in J_k
\end{align}
and its weighted sample analog
\begin{align}
\label{eq:psihat}
\widehat{\psi}_k^v (p, \widehat{\xi}_k) := {\mathbb{E}_{n,k}} v_i g (W_i, p, \widehat{\xi}_k(p)).
\end{align}
Finally, define the matrix error terms
\begin{align*}
\widehat{\Sigma}_k(\widehat{\eta}_k)&:= {\mathbb{E}_{n,k}} A(W_i, \widehat{\eta}_k), \quad \widehat{\Sigma}(\widehat{\eta}) := \dfrac{1}{K} \sum_{k=1}^K \widehat{\Sigma}_k(\widehat{\eta}_k)
\end{align*}
and the weighted matrix error
\begin{align*}
\widehat{\Sigma}^v_k(\widehat{\eta}_k)&:= {\mathbb{E}_{n,k}} v_i A(W_i, \widehat{\eta}_k), \quad \widehat{\Sigma}^v(\widehat{\eta}) := \dfrac{1}{K} \sum_{k=1}^K \widehat{\Sigma}^v_k(\widehat{\eta}_k).
\end{align*}
\paragraph{Empirical process remainder terms.} Define the remainder term
\begin{align*}
R_{1,k} (p)&= \sqrt{N}( (\widehat{\psi}_k(p, \widehat{\xi}_k) - \widehat{\psi}_k(p_0, \xi_0)) - (\psi(p, \widehat{\xi}_k)- \psi_0(p_0)) ), \quad k=1,2,\dots, K,
\end{align*}
the bias term
\begin{align*}
R_{2,k} (p)&= \sqrt{N}( \psi(p, \widehat{\xi}_k)- \psi_0(p) ), \quad k=1,2,\dots, K,
\end{align*}
the second-order remainder term
\begin{align*}
R(p, p_0)&= \sqrt{N}( \psi_0(p) - \psi_0(p_0) -G(p_0)' (p- p_0))
\end{align*}
and the bootstrap term
\begin{align*}
R_{1,k}^v (p)&= \sqrt{N}( (\widehat{\psi}^v_k(p, \widehat{\xi}_k) - \widehat{\psi}^v_k(p_0, \xi_0)) - (\psi(p, \widehat{\xi}_k)- \psi_0(p_0)) ), \quad k=1,2,\dots, K.
\end{align*}
Thus, the support function process $S_N(q)$ of Theorem \ref{thm:limit} can be decomposed as
\begin{align*}
S_N(q):&= \sqrt{N} (\widehat{\sigma}(q) - \sigma(q)) \\
&= \sqrt{N} \left(\dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k (\widehat{p}(q), \widehat{\xi}_k) - \sigma(q) \right) \\
&=\sqrt{N} \left(\dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k(p_0(q), \xi_0) - \sigma(q) + G(p_0)' (\widehat{p}(q) - p_0(q)) \right) \\
&+\dfrac{1}{K} \sum_{k=1}^K [R_{1,k} (\widehat{p}(q))+ R_{2,k} (\widehat{p}(q)) ] + R(\widehat{p}(q), p_0(q)).
\end{align*}
\paragraph{Misclassification events.} For $p, p_0 \in \mathcal{P}$, define the events $\mathcal{E}_{+}(p), \mathcal{E}_{-}(p), \mathcal{E}_{-}(p,p_0), \mathcal{E}_{+}(p,p_0) $
\begin{align}
\label{eq:e+p}
\mathcal{E}_{+}(p) &= \bigg\{ p' V ( \eta) < 0 < p' V ( \eta_0) \bigg\}, \\
\mathcal{E}_{-}(p) &= \bigg\{ p' V ( \eta_0) < 0 < p' V ( \eta) \bigg\} \label{eq:e-p}
\end{align}
and
\begin{align}
\mathcal{E}_{-}(p,p_0) &= \{ p_0' V(\eta_0) < 0 < p' V(\eta_0) \} \label{eq:e-pp0} \\
\mathcal{E}_{+}(p,p_0) &= \{ p' V(\eta_0) < 0 < p_0' V(\eta_0) \}. \label{eq:e+pp0}
\end{align}
\subsection{Proofs of Main Results.}
\begin{proof} [Proof of Lemma \ref{lem:powell}]
Define
\begin{align*}
z(p, \eta):&= p' V(\eta) \\
B_1 (W, \eta, p) :&= p' V ( \eta_0) ( Y(p, \eta) - Y(p, \eta_0) ) \\
B_2 (W, \eta, p) :&= p' ( V (\eta) - V ( \eta_0)) ( Y(p, \eta) - Y(p, \eta_0) ).
\end{align*}
Observe that
\begin{align*}
z(p, \eta) ( Y(p, \eta) - Y(p, \eta_0)) = B_1 (W, \eta, p) + B_2 (W, \eta, p).
\end{align*}
The mistake in $Y(p, \eta)$ can only occur if $p' V ( \eta_0)$ is small enough
\begin{align}
\bigg\{ Y(p, \eta) \neq Y(p, \eta_0) \bigg\} &\Leftrightarrow \bigg\{ \mathcal{E}_{+}(p) \text{ or } \mathcal{E}_{-}(p)\bigg\} \nonumber \\
&\Rightarrow \bigg\{ 0 < | p' V ( \eta_0) | < | p' (V ( \eta) - V (\eta_0))| \bigg\} \nonumber \\
&\Rightarrow \bigg\{ 0 < | p' V ( \eta_0) | < C_P \| V ( \eta) - V (\eta_0) \| \bigg\}=: \mathcal{E}_{+-}(p) . \label{eq:peta0bound}
\end{align}
Recall that
\begin{align}
\label{eq:pp0}
Y(p, \eta) -Y(p, \eta_0)=(Y_U - Y_L) 1\{ \mathcal{E}_{+}(p) \cup \mathcal{E}_{-}(p)\}
\end{align}
Invoking $| \mathbb{E} X | \leq \mathbb{E} | X| $ and \eqref{eq:pp0} gives
\begin{align}
| \mathbb{E} B_1 (W, \eta, p) |&= | \mathbb{E} p' V ( \eta_0) (Y(p, \eta) -Y(p, \eta_0) ) | \nonumber \\
&\leq \mathbb{E} | p' V ( \eta_0) | (Y_U - Y_L) 1\{ \mathcal{E}_{+}(p) \cup \mathcal{E}_{-}(p)\} \nonumber \\
&\leq \mathbb{E} | p' V ( \eta_0) | (Y_U - Y_L) 1\{ \mathcal{E}_{+-}(p) \} \label{eq:peta0bound2}
\end{align}
Invoking definition of $\mathcal{E}_{+-}(p)$ in \eqref{eq:peta0bound} and $Y_U - Y_L \leq M_{UL} \text{ a.s. }$ gives
\begin{align}
&\mathbb{E} | p' V ( \eta_0) | (Y_U - Y_L) 1\{ \mathcal{E}_{+-}(p) \} \nonumber \\
&\leq C_P M_{UL} \mathbb{E} \| V ( \eta) - V (\eta_0) \| 1\{ \mathcal{E}_{+-}(p) \}. \label{eq:peta0bound3}
\end{align}
The second-order bias term is bounded as
\begin{align}
| \mathbb{E} B_2 (W, \eta, p) | \leq C_P M_{UL} \mathbb{E} \| V ( \eta) - V (\eta_0)\| 1\{ \mathcal{E}_{+}(p) \cup \mathcal{E}_{-}(p)\}. \label{eq:peta0bound4}
\end{align}
Invoking Assumption \ref{ass:regularity} gives
\begin{align}
&\mathbb{E} \| \eta_0(X) - \eta(X) \| 1\{ \mathcal{E}_{+}(p) \cup \mathcal{E}_{-}(p)\} \nonumber \\
&=\mathbb{E}_{X} \| \eta(X) - \eta_0(X) \| \int_{- C_P \| \eta(X) - \eta_0(X) \| }^{C_P \| \eta(X) - \eta_0(X) \| }h_{p' V(\eta_0) \mid X} (t, X) dt \nonumber \\
&\leq 2 C_P M_h \mathbb{E}_{X} \| \eta(X) - \eta_0(X) \|^2 \label{eq:lipbound}.
\end{align}
Combining the bounds gives
\begin{align}
| \mathbb{E} [B_1 (W, \eta, p) + B_2 (W, \eta, p)] | &\leq |\mathbb{E} B_1 (W, \eta, p) | + | \mathbb{E} B_2 (W, \eta, p) | \nonumber \\
&\leq 4 C_P^2 M_{UL} M_{h} \mathbb{E} \| \eta(X) - \eta_0(X) \|^2 \label{eq:powell}.
\end{align}
\end{proof}
\begin{proof}[Proof of Theorem \ref{thm:limit}]
\textbf{Step 1.} This step is required only if $\Sigma$ is unknown, such as in Example \ref{ex:plpiv}. As shown in Lemma \ref{lem:matrixlinarization}, for $v=1$ (regular case) and $v=e$ (bootstrap case),
\begin{align}
\label{eq:sigma4main}
(\widehat{\Sigma}^v (\widehat{\eta}))^{-1} - \Sigma^{-1} = -\Sigma^{-1} ( \widehat{\Sigma}^v (\eta_0) - \Sigma ) \Sigma^{-1} + M^v,
\end{align}
where the remainder matrix $M^v$ obeys $\| M^v \| = o_P (N^{-1/2})$ for both cases. Take $v=1$. Post-multiplying the LHS above by $q$ gives
\begin{align}
\sqrt{N} (\widehat{p}(q) - p_0(q)) &= - \Sigma^{-1} ( \widehat{\Sigma} (\eta_0) - \Sigma ) \Sigma^{-1}q + o_P(1) \nonumber
\end{align}
Likewise, taking $v=e$ gives
\begin{align*}
\sqrt{N} (\widetilde{p}(q) - p_0(q)) &= ((\widehat{\Sigma}^e (\widehat{\eta}))^{-1} - \Sigma^{-1})'q \\
&= - \Sigma^{-1} ( \widetilde{\Sigma}^e (\eta_0) - \Sigma ) \Sigma^{-1}q + o_P(1).
\end{align*}
For some $N$ large enough, $ \widehat{p}(q) \in \mathcal{P} \quad \forall q \in \mathcal{S}^{d-1}$ and $ \widetilde{p}(q) \in \mathcal{P} \quad \forall q \in \mathcal{S}^{d-1}$ with probability $1-o(1)$.
\textbf{ Step 2. } Let $k=1,2,\dots, K$ denote the partition index. We bound $R_{1,k} (\widehat{p}(q))$ and $R^e_{1,k} (\widehat{p}(q))$. Define the function class
\begin{align*}
\mathcal{F}_{2k}^v := \{ v\cdot(g(\cdot,p,\widehat{\xi}_k(p)) - g(\cdot,p_0,\xi_0(p_0))), \quad p , p_0 \in \mathcal{P}, \quad \| p - p_0 \| \leq \tau_N \},
\end{align*}
where $v=1$ (regular case) and $v=e$ (bootstrap case). The class $\mathcal{F}_{2k}$ is obtained as $$\mathcal{F}_{2k} \subset \mathcal{G}_{\widehat{\xi}_k} - \mathcal{G}_{\xi_0},$$ where $\mathcal{G}_{\xi_0}$ and $\mathcal{G}_{\widehat{\xi}_k} $ are defined in Assumption \ref{ass:concentration:chap1}. On the event $\mathcal{E}_N$,
\begin{align*}
&\sup_{p \in \mathcal{P} } |g(W_i,p,\widehat{\xi}_k(p)) - g(W_i,p_0,\xi_0(p_0)) | \\
&\leq \sup_{p \in \mathcal{P}} | g(W_i,p,\widehat{\xi}_k(p)) | + \sup_{p \in \mathcal{P}} | g(W_i,p,\xi_0(p)) | \\
&\leq G_{\widehat{\xi}_k} + G_{\xi_0},
\end{align*}
and $G_{\widehat{\xi}_k \xi_0} := G_{\widehat{\xi}_k} + G_{\xi_0}$ is a measurable envelope for the class $\mathcal{F}_{2k}$. Note that $\| G_{\widehat{\xi}_k \xi_0} \|_{P,c} \leq \| G_{\xi_0} \|_{P,c} + \| G_{\widehat{\xi}_k } \|_{P,c} \leq 2 C_1$. The uniform covering entropy of the function class $ \mathcal{G}_{\widehat{\xi}_k} - \mathcal{G}_{\xi_0}$ is bounded as
\begin{align*}
& \log \sup_{Q} N(\epsilon \| G_{\xi} + G_{\xi_0} \|_{Q,2}, \mathcal{G}_{\widehat \xi} - \mathcal{G}_{\xi_0} , \| \cdot \|_{Q,2}) \\
&\leq \log \sup_{Q} N(\epsilon/2 \| G_{\xi} \|_{Q,2}, \mathcal{G}_{\xi} , \| \cdot \|_{Q,2}) + \log \sup_{Q} N(\epsilon/2 \| G_{\xi_0} \|_{Q,2}, \mathcal{G}_{\xi_0} , \| \cdot \|_{Q,2}) \\
&\leq 2v \log (2a/\epsilon)
\end{align*}
by the proof of Theorem 3 in \cite{andrews:1994b} and Assumption \ref{ass:concentration:chap1}. The class $\mathcal{F}_{2k}^e $ is obtained by multiplication of $\mathcal{F}_{2k}$ by an integrable random variable independent of the data, and therefore retains $P$-Donsker and uniform covering properties of the class $\mathcal{F}_{2}$. In particular,
$G^v_{\widehat{\xi}_k \xi_0} := |v| (G_{\widehat{\xi}_k} + G_{\xi_0})$ is a valid envelope for $ \mathcal{F}_{2k}^v $. Next, for $v=1$ and $v=e$,
\begin{align*}
&\sup_{\xi \in \Xi_N} \mathbb{E} v^2 (g(W,p,\xi(p)) - g(W,p_0,\xi_0(p_0)) )^2 \\
&\leq 2 \bigg( \sup_{\xi \in \Xi_N} \sup_{p \in \mathcal{P}} \mathbb{E} v^2 ( g(W,p,\xi(p)) - g(W,p,\xi_0(p)))^2 \\
&+ \sup_{p_0, p \in \mathcal{P}, \quad \| p - p_0 \| \leq \tau_N} \mathbb{E} v^2 (g(W,p,\xi_0(p)) - g(W,p_0,\xi_0(p_0)) )^2 \bigg) \\
&\leq 2((r_N'')^2 + (r_N')^2).
\end{align*}
Invoking Lemma \ref{lem:maxineq} conditional on $(W_i)_{i \in J_k}$ and taking $\sigma = 2 (r_N' + r_N'')$ gives
\begin{align*}
\sup_{f \in \mathcal{F}_2^v} | {\mathbb{G}_{n,k}} [f] | &= \sup_{ p_0, p \in \mathcal{P}} | \widehat{\psi}_k (p, \widehat{\xi}_k) - \widehat{\psi}_k(p_0, \xi_0) - (\psi(p, \widehat{\xi}_k) - \psi(p_0, \xi_0)) | \\
&\lesssim_P (r_N''+r_N' )\log^{1/2} (1/(r_N''+r_N' )) + N^{-1/2+1/c} \log N \\
&= o_P(1),
\end{align*}
where the last equality follows from Assumption \ref{ass:concentration:chap1}. By Lemma \ref{lem:cond}, $\sup_{f \in \mathcal{F}_2^v} | {\mathbb{G}_{n,k}} [f] | = o_P(1)$ holds unconditionally. as shown in Step 1, wp $1-o(1)$, $\widehat{p}(q) \in \mathcal{P} \quad \forall q$, and
\begin{align*}
\sup_{q \in \mathcal{S}^{d-1}}| R_{1,k} (\widehat{p}(q)) | \lesssim_P \sup_{p \in \mathcal{P}}| R_{1,k} (p) | \lesssim_P \sup_{f \in \mathcal{F}_2} | {\mathbb{G}_{n,k}} [f] | = o_P(1).
\end{align*}
\textbf{ Step 3. Bound on $R_{2,k} (\widehat{p}(q))$. } On the event $\mathcal{E}_N$, for any $k=1,2,\dots, K$
\begin{align*}
\sup_{q \in \mathcal{S}^{d-1}}| R_{2,k} (\widehat{p}(q)) | &\leq \sup_{p \in \mathcal{P}} \sqrt{N} | \psi(p, \widehat{\xi}_k) - \psi_0(p) | \\
&\leq \sup_{ \xi \in \Xi_N} \sup_{p \in \mathcal{P}}\sqrt{N} | \psi(p, \xi) - \psi_0(p) | \lesssim \sqrt{N} \mu_N.
\end{align*}
Combining the bounds over a finite set $k=1,2,\dots, K$ gives
\begin{align*}
\dfrac{1}{K} \sum_{k=1}^K \bigg[R_{1,k} (\widehat{p}(q)) + R_{2,k} (\widehat{p}(q)) \bigg] = o_P(1).
\end{align*}
\textbf{ Step 4. Bound on $R(\widehat{p}(q), p_0)$. }
By Assumption \ref{ass:jacobian}, Lemma \ref{lem:uniderivative} and Step 1, on the event $\mathcal{E}_N$,
\begin{align*}
\sup_{q \in \mathcal{S}^{d-1}} | R(\widehat{p}(q), p_0(q)) | = o (\sqrt{N} \| (\widehat{\Sigma} (\widehat{\eta}))^{-1} - \Sigma^{-1} \|) = o_P(1).
\end{align*}
\textbf{ Step 5. Conclusion. } For the influence function $h(W,q)$ in \eqref{eq:hwq}, the function class
\begin{align*}
\mathcal{H} &= \bigg\{ h(\cdot, q), \quad q \in \mathcal{S}^{d-1} \bigg\} \subseteq \mathcal{G}_{\xi_0} +\mathcal{H}_A
\end{align*}
is included into the sum $\mathcal{G}_{\xi_0} + \mathcal{H}_A $. The function classes $\mathcal{G}_{\xi_0}$ and $\mathcal{H}_A $ are Donsker classes with square integrable envelopes (by Assumption \ref{ass:concentration:chap1} and Lemma \ref{lem:maxineq:suppfun}, respectively). The statement of the lemma follows from the Skorohod-Dudley-Wichura construction, as in \cite{Skorohod}, \cite{Dudley} and \cite{Wichura}.
\end{proof}
The bootstrap support function process can be decomposed as
\begin{align*}
\widetilde{S}_N(q) :&=\sqrt{N} \left( (\widetilde{\sigma}(q) - \sigma(q)) - (\widehat{\sigma}(q) - \sigma(q)) \right).
\end{align*}
The first summand can be decomposed as
\begin{align*}
\widetilde{\sigma}(q) - \sigma(q) &=\dfrac{1}{K} \sum_{k=1}^K \bar{e}^{-1} \widehat{\psi}_k^v(\widetilde{p}(q), \widehat{\xi}_k) - \sigma(q) \\
&=\dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k^v(\widetilde{p}(q), \widehat{\xi}_k) / (1+ o_P(1)) - \sigma(q) \\
&=\dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k^v(\widetilde{p}(q), \widehat{\xi}_k) - \sigma(q) + o_P(1).
\end{align*}
The weighted moment can be decomposed as
\begin{align*}
\sqrt{N} \dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k^v(\widetilde{p}(q), \widehat{\xi}_k) &= \sqrt{N} \dfrac{1}{K} \sum_{k=1}^K \widehat{\psi}_k^v(p_0(q), \xi_0) + \sqrt{N} G(p_0(q))' (\widetilde{p}(q) - p_0(q)) \\
&+ \dfrac{1}{K} \sum_{k=1}^K \bigg[R_{1,k}^e (\widetilde{p}(q)) + R_{2,k} (\widetilde{p}(q)) \bigg] + R(\widetilde{p}(q), p_0(q)).
\end{align*}
\begin{proof}[Proof of Theorem \ref{thm:bb}]
By Comment B.1 in \cite{CCMS}, if the bootstrap random element converges in probability $P$ unconditionally (i.e, $Z_N=o_P(1)$), then $Z_N = o_{P^e} (1)$
in $L^1(P)$ sense and hence in probability $P$, where $P^e$ denotes the probability measure conditional on the data.
\textbf{ Step 1. } As shown in Step 1 of the proof of Theorem \ref{thm:limit}, $\widetilde{p}(q) \in \mathcal{P} \quad \forall q$. The bound on $R_{1,k}^e (\widetilde{p}(q))$ is established in the proof of Theorem \ref{thm:limit}, Step 2. By construction,
\begin{align*}
\sqrt{N}(\psi(p, \widehat{\xi}_k) - \psi(p, \xi_0)) = R_{2,k} (p) = \mathbb{E} [ v \cdot (g(W, p, \widehat{\xi}_k) - g(W, p, \xi_0)) \mid (W_i)_{i \in J_k} ].
\end{align*}
Thus, the bounds on $R_{2,k}(p)$ and $R(p, p_0)$ are established in Steps 3 and 4 of Theorem \ref{thm:limit}.
\textbf{ Step 2. } As shown in Step 5 of Theorem \ref{thm:limit}, the function class $\mathcal{H} \subseteq \mathcal{G}_{\xi_0} + \mathcal{H}_A$ is a Donsker class with
square-integrable envelopes. Then by the Donsker theorem for exchangeable bootstraps, weak convergence holds conditional on the
data,
\begin{align*}
\mathbb{G}_N [ (e-1) h(W,q) ] /\bar{e} \Rightarrow \widetilde{\mathbb{G}[h(q)]} \text{ under} P^e \text{ in probability } P,
\end{align*}
where $ \widetilde{\mathbb{G}[h(q)]} $ is a P-Brownian bridge independent of $G[h(q)]$ with the same distribution as $G[h(q)]$.
\end{proof}
\newpage
\begin{abstract}
This appendix contains supplementary statements for the paper ``Debiased Machine Learning of Set-Identified Linear Models'' by Vira Semenova. Appendix A contains useful technical statements. Appendix B contains additional proofs. Appendix C contains supplementary tables.
\end{abstract}