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.
71,903 characters
Adaptive Inference in Multivariate Nonparametric Regression Models Under Monotonicity
\maketitle
\begin{abstract}
We consider the problem of adaptive inference on a regression function at a point under a multivariate nonparametric regression setting. The regression function belongs to a H\"older class and is assumed to
be monotone with respect to some or all of the arguments. We derive the
minimax rate of convergence for confidence intervals (CIs) that adapt to the
underlying smoothness, and provide an adaptive inference procedure that
obtains this minimax rate. The procedure differs from that of
\cite{cai2004adaptation}, intended to yield shorter CIs under practically
relevant specifications. The proposed method applies to general linear
functionals of the regression function, and is shown to have favorable
performance compared to existing inference procedures.
\end{abstract}
\newpage
\section{Introduction}
We consider the problem of inference on a regression function at a point under
the nonparametric regression model
\begin{equation*}
y_{i} = f(x_{i}) + u_{i}, \quad u_{i} \overset{\mathrm{i.i.d.}}{\sim} N(0, \sigma^{2}(x_{i})),
\end{equation*}
where $f$ is assumed to lie in a H\"older class with exponent
$\gamma \in (0,1]$. Procedures based on $\gamma$ is conservative (or suboptimal) when the
true regression function in fact lies in a smoother H\"older class with
$\gamma' > \gamma$. Adaptive procedures try to overcome this issue by
automatically adjusting to the (unknown) underlying smoothness class. However,
unlike in the case of estimation, where adaptation to the unknown smoothness
class is in general possible with an additional logarithmic term
(\citealp{lepskii1991ProblemAdaptiveEstimationa}), adaptation is impossible in the case of inference without further
restrictions on the function class (\citealp{low1997nonparametric}) .
Two shape restrictions that can be used to overcome this impossibility have been
discussed in the literature, convexity and monotonicity. In this paper, we
impose monotonicity on the regression function to construct a CI that adapts to
the underlying smoothness of the regression function. The main difference with other papers that consider
adaptation under a monotonicity condition (\citealp{cai2013AdaptiveConfidenceIntervals}; \cite{armstrong2015AdaptiveTestingRegression}) is our general treatment of the
dimension of $x_{i}$. To our knowledge, this is the first
paper to construct adaptive CIs, under
a multivariate nonparametric regression setting.
We consider coordinate-wise monotonicity with respect to all or some of the
coordinates. A function $f$ is \textit{coordinate-wise monotone} with respect to
$\mathcal{V} \subseteq \{1, \dots, k\}$ if $x_{j}\geq z_{j}$ for all
$j \in\mathcal{V}$ and $x_{j}= z_{j}$ for all $j\notin\mathcal{V}$ imply
$f(z) \geq f(z)$. The minimax expected length of a CI over the H\"older class
with exponent $\gamma$ converges to $0$ at the well-known rate of
$n^{-1/(2 + k/\gamma)}$. When the regression function is monotone in all
variables, i.e., $\mathcal{V} = \{1, \dots, k\}$, we can construct a CI that
achieves this minimax rate over all $\gamma \in (0,1]$ just as in the univariate
case. Also, again as in the univariate case, if the regression is not monotone
to any of the variables so that $\mathcal{V} = \emptyset$, there is no scope for
adaptation.
An interesting case is when the function is monotone with respect to only some
of the variables so that $k_{+} := \lvert \mathcal{V} \rvert < k$, which can arise due to the multivariate nature of the problem. In this case, we
show that for a CI that maintains coverage over the H\"older class with exponent
$\gamma$, the minimax expected length over a smoother class $\gamma' > \gamma$
converges to $0$ at the rate $n^{-1/(2 + k_{+}/\gamma' +
(k-k_{+})/\gamma)}$. The denominator of the exponent can be written as
$2 + k/\gamma - k_{+}( 1/\gamma- 1/\gamma')$. This is the sum of a term that comes from the minimax rate over $\gamma$, $2 + k/\gamma$, and $- k_{+}( 1/\gamma- 1/\gamma')$. In this sense,
$k_{+}( 1/\gamma- 1/\gamma')$ exactly quantifies the possible gain from
monotonicity, indicating larger gains if the regression function is monotone in
more variables and/or smoother.
We propose a CI that obtains this minimax rate (of adaptation) for a sequence of
H\"older exponents $\{\gamma_{j}\}_{j=1}^{J} \subset (0, 1]$. While the method
provided by \cite{cai2004adaptation} can be used to construct such a CI, we
provide an alternative method that builds upon the one-sided CI proposed by
\cite{armstrong2018optimal}. Their one-sided CI ``directs power'' to a smoother class while maintaining
coverage over a larger class of functions. Our CI is constructed by combining
the lower and upper versions of their one-sided CI to create a two-sided CI, and
then taking the intersection of a sequence of such two-sided CIs that direct
power to each $\gamma_{j}$. An appropriate Bonferroni correction is used to
obtain correct coverage. This CI can be used in more general nonparametric
regression settings, as long as the parameter of interest is a linear functional
of the regression function and the regression functions lies in a convex function
class.
While the proposed CI obtains the minimax length over $\gamma_{j}$ for each $j$
up to a constant factor that does not depend on the sample size, this constant
does depend on the number of parameter spaces $J$ the CI adapts to. This is in
contrast with the CI of \cite{cai2004adaptation}, which gives a multiplicative
constant that does not depend on $J$. However, the multiplicative constant of
our CI grows slowly with $J$ at a $(\log J)^{1/2}$ rate, and is smaller than the
constant given by \cite{cai2004adaptation} for any reasonable specification of
$J$. Even if one wishes to adapt to $J= 10^{3}$ parameter spaces, our CI obtains
the minimax expected length of each parameter space within a multiplicative
constant of $4.14$, whereas this constant is $16$ for the CI by
\cite{cai2004adaptation}. A simulation study confirms that our CI can be
significantly shorter in practice as well. Nonetheless, the uniform constant
that \cite{cai2004adaptation} obtain is theoretically attractive and allows one
to adapt to the continuum of H\"older exponents $(0,1]$ in this context.\\
\noindent \textbf{Related literature. }An adaptation theory for CIs in a nonparametric regression setting was developed
by \cite{cai2004adaptation}. \cite{cai2013AdaptiveConfidenceIntervals} provide a
procedure for constructing adaptive CIs that adapt to each individual function
under monotonicity and convexity. \cite{armstrong2015AdaptiveTestingRegression}
provides an inference method for the regression function at a point, possibly on
the boundary of the support, that adapts to the underlying H\"older classes
under a monotonicity assumption. As noted earlier, the main difference of our
paper is that we consider a multivariate regression setting where there is no
restriction on the dimension of the independent variable as long as it is fixed
and finite. The adaptation theory for CIs builds upon the more classical minimax
theory for CIs, which has been developed in \cite{Donoho1994} and
\cite{low1997nonparametric}. \cite{cai2012MinimaxAdaptiveInference} provides an
excellent review on the theory of minimax and adaptive CIs, along with the
minimax and adaptive estimation problems.
While the focus of this paper is on adaptive CIs, there are other forms of
confidence sets that are of interest in the context of nonparametric regression
setting. Adaptive confidence balls have been considered in
\cite{genovese2005ConfidenceSetsNonparametric},
\cite{cai2006AdaptiveConfidenceBalls} and
\cite{robins2006AdaptiveNonparametricConfidencea}. An adaptation theory for
confidence bands has been considered in, for example,
\cite{dumbgen1998NewGoodnessoffitTests},
\cite{genovese2008AdaptiveConfidenceBands}, and
\cite{cai2014AdaptiveConfidenceBands}. In the context of density estimation,
adaptive confidence bands have also been considered in
\cite{hengartner1995FiniteSampleConfidenceEnvelopes},
\cite{gine2010ConfidenceBandsDensitya}, and
\cite{hoffmann2011AdaptiveInferenceConfidencea}.
Recently, there has been interest in isotonic regression in general
dimensions. The monotonicity condition imposed in such models is the same as the
one we impose here with $\mathcal{V} = \{1, \dots, k\}$.
\cite{han2019IsotonicRegressionGeneral} derive minimax rates for the least
squares estimation problem. \cite{deng2020ConfidenceIntervalsMultiple} provide a
method for constructing CIs at a point based on block max-min and min-max
estimators.\\
\noindent \textbf{Outline.} Section \ref{sec:Setup} describes the nonparametric
regression model and the function class we consider. Section
\ref{sec:adapt-conf-interv} introduces the notion of adaptivity in more detail
and describes our procedure for constructing adaptive CIs. Section
\ref{sec:adapt-infer-f0} presents the main result of the paper, the minimax rate of adaptation, and an adaptive CI that obtains this rate by solving the corresponding
modulus problem. Section \ref{sec:Simulation} provides a simulation study, and
Section \ref{sec:empir-illustr} illustrates our method in the
context of production function estimation.
Any proof omitted in the main text can be found in the appendix. Appendix
\ref{appdx-sec:Proofs_lemmas} collects the proofs for lemmas and
corollaries. Appendix \ref{sec: proof main thm} contains the proof for our main
theoretical result, Theorem \ref{thm:rate_result}.
\section{Nonparametric Regression Under Monotonicity\label{sec:Setup}}
We observe $\left\{ \left(y_{i},x_{i}\right)\right\} _{i=1}^{n}$ and consider a
nonparametric regression model,
\begin{equation}
y_{i}=f(x_{i})+u_{i},\label{eq:npr_model}
\end{equation}
where $x_{i}\mathcal{\in X}\subset\mathbb{R}^{k}$ is a (fixed) regressor,
$f:\mathbb{R}^{k}\to\mathbb{R}$ is the unknown regression function that lies in
some function class $\mathcal{F}$, and $u_{i}$'s are independent with
$u_{i}\sim N(0,\sigma^{2}(x_{i}))$ and $\sigma^{2}(\cdot)$ known. The parameter
of interest is $f(x_{0})$. For the rate results provided in Section \ref{sec:Asymptotic-rate-results}, we require that $x_{0} \in \mathrm{Int}\,
\mathcal{X}$. However, we note that the solution to the modulus problem given in Section \ref{sec:solut-modul-probl} does not depend on whether $x_0$ is on the boundary or not. Without loss of generality, we normalize $x_{0}$ to be $0$.
We take the $\mathcal{F}$ to be the class of functions that are H\"older
continuous and nondecreasing in all or some of the variables. Let
$\Lambda(\gamma,C)$ denote the set of functions from $\mathbb{R}^{k}$ to
$\mathbb{R}$ that are H\"older continuous with H\"older constants $(\gamma,C)$,
\begin{equation*}
\Lambda(\gamma,C):=\left\{ f\in\mathcal{\mathcal{F}}(\mathbb{R}^{k},
\mathbb{R}):\left|f(x)-f(z)\right\rvert \leq C\left\Vert x-z \right\Vert
^{\gamma}\text{ for all }x,z\in\mathcal{X}\right\} ,
\end{equation*}
where $\mathcal{F}\left(\mathbb{R}^{k}, \mathbb{R}\right)$ is the set of
functions from $\mathbb{R}^{k}$ to $\mathbb{R}$, $\gamma\in[0,1],$ $C \geq 0$
and $\left\Vert \cdot\right\Vert $ is a norm on $\mathbb{R}^{k}$. For notational
simplicity, we omit the dependence of the function class on the choice of the
norm $\lVert \cdot \rVert$. We impose the following restriction that
$\lVert \cdot \rVert$ is monotone in the magnitude of each element, which is
satisfied by most norms used in practice. such as the $\ell_{p}$ norm or a
weighted version of it. We discuss the relationship between this assumption and
the monotonicity of the regression function in Remark \ref{rem:norm}.
\begin{assumption}
\label{assu:norm}$\left\Vert \cdot\right\Vert $ is a norm on $\mathbb{R}^{k}$
such that $\left\Vert z \right\Vert $ is nondecreasing in
$\left\lvert z_{j}\right\rvert $ for each $j\in\{1,...,k\}.$
\end{assumption}
We now define the (coordinate-wise) monotone H\"older class. For a subset of the
covariate indices $\mathcal{V}\subset\left\{ 1,\dots,k\right\} $, write
\begin{equation*}
\Lambda_{+,\mathcal{V}}(\gamma,C):=\left\{ f\in\Lambda(\gamma,C):f(x)\geq
f(z)\,\,\text{if }x_{j}\geq z_{j}\text{ }\forall
j \in\mathcal{V} \text{ and }x_{j}= z_{j}\,\,\forall
j\notin\mathcal{V}\right\}.
\end{equation*}
This is the set of H\"older continuous functions that are nondecreasing,
coordinate-wise, with respect to the $j$th element for $j \in
\mathcal{V}$. Define $k_{+}:=\left\lvert \mathcal{V}\right\rvert .$ By a
relabeling argument, it is without loss of generality to write
$\mathcal{V}:=\left\{ 1,\dots,k_{+}\right\} $. If $k_{+}=k$, then
$\Lambda_{+,\mathcal{V}}(\gamma,C)$ is the set of nondecreasing and H\"older
continuous functions where the monotonicity is with respect to the
coordinate-wise partial ordering on $\mathbb{R}^{k}$.
\section{Adaptive Confidence Intervals}
\label{sec:adapt-conf-interv}
\subsection{Notion of Adaptivity}
In this section, we discuss the problem of inference for a general linear
functional of the regression function, $Lf$. Consider a sequence of convex
parameter spaces $\mathcal{F}_{1}$, \dots, $\mathcal{F}_{J}$, with the
requirement that $\mathcal{F}_{j} \subset \mathcal{F}_{J}$ for all $j \leq
J$. Note that the parameter spaces are not necessarily nested, but there is a
largest convex parameter space that nests all the other parameter spaces. Here,
$\mathcal{F}_{J}$ reflects a conservative choice of the parameter space where
the researcher believes the true regression function to lie in. Hence, the CI we
construct will be required to maintain correct coverage over this space. An
adaptive CI maintains this correct coverage over the largest parameter space
$\mathcal{F}_{J}$ while having good performance (e.g. shorter expected length)
when the true function happens to lie in the smaller parameter space
$\mathcal{F}_{j}$, simultaneously for all $j \leq J$.
Then, a natural question is how well a CI that maintains coverage over
$\mathcal{F}_{J}$ can perform over $\mathcal{F}_{j}$, which is one of the main
questions that \cite{cai2004adaptation} raise and address in detail in the
context of two-sided CIs. The case of one-sided CIs has been considered by
\cite{armstrong2018optimal}, along with other questions.
\subsubsection{Two-sided Adaptive CIs}
\label{sec:two-sided-adaptive}
Let $\mathcal{I}_{\alpha,2}^{J}$ denote the set of all two-sided CIs that have
coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$ Following
\citet{cai2004adaptation}, the performance criterion we consider for two-sided
CIs is the worst-case expected length. That is, the performance of a CI, $CI$,
over the parameter space $\mathcal{F}_{j}$ is measured by
$\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}}_{f}\mu(CI)$ with smaller values of this quantity
meaning better performance. Here, $\operatorname{\mathbf{E}}_{f}$ denotes the expectation when the true
regression function is $f$ and $\mu$ is the Lebesgue measure on the real
line. Then, the shortest possible worst-case expected length a CI can achieve
over $\mathcal{F}_{j}$ (while maintaining correct coverage over
$\mathcal{F}_{J}$) is characterized by the quantity
\begin{equation*}
L_{j,J}^{\ast}:=\inf_{CI\in\mathcal{I}_{\alpha,2}^{J}}\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}}_{f}\mu(CI).
\end{equation*}
Following \citet{cai2004adaptation}, we say a CI is \textit{adaptive} if it
achieves $L_{j,J}^{\ast}$ for all $j \leq J$ up to a multiplicative constant
that does not depend on the sample size. Let $z_{q}$ denote the $q$--quantile of
the standard normal distribution. \cite{cai2004adaptation} show that
$L_{j,J}^{\ast} \asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J}),$
with $\asymp$ denoting asymptotic equivalence\footnote{We write
$a_{n}\asymp b_{n}$ if
$
0<\underset{n\to\infty}{\lim\inf}({a_{n}/}{b_{n}})\leq\underset{n\to\infty}{\lim\sup}({a_{n}}/{b_{n}})<\infty.$}
and $\omega_{+}(\delta,\mathcal{F}_{j},\mathcal{F}_{J})$ is the \textit{between
class modulus of continuity} defined as
\begin{align*}
&\omega_{+}(\delta,\mathcal{F}_{j},\mathcal{F}_{J}) \\
:=&\sup\left\{
\left\lvert Lf_{J}-Lf_{j}\right\rvert :\textstyle\sum_{i=1}^{n}\left((f_{J}(x_{i})-f_{j}(x_{i}))/\sigma(x_{i})\right)^{2}\leq\delta^{2},f_{j}\in\mathcal{F}_{j},f_{J}\in\mathcal{F}_{J}\right\},
\end{align*}
for $\delta \geq 0.$\footnote{Note that the definition
is slightly different with \cite{cai2004adaptation} due to the $\sigma(x_{i})$
term that appears in the denominator of the summand. This is because we divide
both sides of \eqref{eq:npr_model} by the (known) $\sigma(x_{i})$ to convert
the model into the same form as that of \cite{cai2004adaptation}.} In general,
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$ is more tractable than
$L_{j,J}^{\ast}$, and thus the strategy is to construct a CI that has worst case
length over $\mathcal{F}_{j}$ bounded by
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$, up to a multiplicative
constant. We refer to the rate at which
$\omega_{+}(z_{1-\alpha}, \mathcal{F}_{j}, \mathcal{F}_{J}) $ converges to $0$ as the
\textit{minimax rate of adaptation} (of $\mathcal{F}_{j}$ over
$\mathcal{F}_{J}$). If $\mathcal{F}_{j}=\mathcal{F}_{J}$, this is the minimax
rate over $\mathcal{F}_{J}$, which is the fastest rate at which the worst-case
expected length over $\mathcal{F}_{J}$ of a CI that maintains correct coverage
over the same space $\mathcal{F}_{J}$ can achieve.
\subsubsection{One-sided CIs}
\label{sec:one-sided-cis}
While our main focus is on adaptive two-sided CIs, the construction of our
adaptive CI relies heavily on the one-sided CI proposed by
\cite{armstrong2018optimal}. Hence, we briefly describe the notion of adaptivity
in the context of one-sided CIs. For one-sided CIs,
we follow \citet{armstrong2018optimal} and consider the $\beta$th quantile of
excess length as the performance criterion. More specifically, for a one sided
lower CI, $[\hat{c},\infty)$, we denote the $\beta$th quantile of the excess
length at $f$ as $q_{\beta,f}(Lf-\hat{c})$, where $q_{\beta,f}(\cdot)$ denotes
the $\beta$th quantile function when the true regression function is $f$. Under
this criterion, the best possible performance over $\mathcal{F}_{j}$ is
quantified by
\begin{equation*}
\ell_{j,J}^{\ast}:=\inf_{\hat{c}:[\hat{c},\infty)\in\mathcal{I}_{\alpha,\ell}^{J}}\sup_{f\in\mathcal{F}_{j}}q_{\beta,f}(Lf-\hat{c}),
\end{equation*}
where $\mathcal{I}_{\alpha, \ell}^{J}$ denotes the set of all one-sided lower
CIs that have coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$
\citet{armstrong2018optimal} showed that
$\ell_{j,J}^{\ast}=\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{J},\mathcal{\mathcal{F}}_{j})$,
where
$\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{J},\mathcal{\mathcal{F}}_{j})$
is the \textit{ordered class modulus of continuity} defined as
\begin{align*}
&\omega(\delta,\mathcal{F}_{j},\mathcal{F}_{k}) \\
:=&\sup\left\{
Lf_{k}-Lf_{j}:\textstyle\sum_{i=1}^{n}\left((f_{k}(x_{i})-f_{j}(x_{i}))/\sigma(x_{i})\right)^{2}
\leq \delta^{2},f\in\mathcal{F}_{j},f_{J}\in\mathcal{F}_{J}\right\},
\end{align*}
for any $\delta \geq 0$ and $j,k \leq J$. We refer to the optimization problem
in the definition as the ordered modulus problem. Naturally, an analogous result
holds for upper one-sided CIs so that
$u_{j,J}^{\ast}=\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{j},\mathcal{\mathcal{F}}_{J})$,
where
\begin{equation*}
u_{j,J}^{\ast}:=\inf_{\hat{c}:(-\infty,\hat{c}]\in\mathcal{I}_{\alpha,u}^{J}}\sup_{f\in\mathcal{F}_{j}}q_{\beta,f}(\hat{c}-Lf),
\end{equation*}
with $\mathcal{I}_{\alpha, u}^{J}$ denoting the set of all one-sided upper CIs
that have coverage at least $1-\alpha$ over $\mathcal{F}_{J}.$
We say a one-sided lower CI, $[\hat{c}^{\ast},\infty)$, is adaptive if there
exists some $c>0$ that does not depend on $n$ such that
\begin{equation*}
\sup_{f\in\mathcal{F}_{j}}q_{\beta,f}(Lf-\hat{c}^{\ast})\leq
c\,\omega(z_{1-\alpha}+z_{\beta},\mathcal{\mathcal{F}}_{j},\mathcal{\mathcal{F}}_{J})
\end{equation*}
for all $j \leq J$, and similarly for one-sided upper CIs.
\subsubsection{Modes of Adaptation}
\label{sec:scope-adaptation}
Note that it must be the case that
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j}) \leq
\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$ (and similarly for the
ordered moduli) because
$ \omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})$ takes the supremum over a
smaller set. However, if it happens to be the case that
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})$,
an adaptive CI, $CI^{\ast}$, satisfies
$\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}}_{f}\mu(CI^{\ast})\leq\overline{c}\,L_{j,j}^{\ast}$
for all $j\leq J$. \cite{cai2004adaptation} define such CI to be \textit{strongly
adaptive.} This is an ideal case because we obtain $L_{j,j}^{\ast}$, up to a
multiplicative constant, which is the minimax length we could have achieved if
we ``knew'' that our true regression function lied in the smaller class
$\mathcal{F}_{j}$ (i.e., if we made a stronger assumption that the true
regression function lies in this smaller class). While adaptive CIs exist in
general, strong adaptation is possible only when
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})$
for all $j \leq J$. This is not a property of a given procedure, but of the
given statistical model.
The least desirable case is when
$\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})\asymp\omega_{+}(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{J})$,
because this leaves no scope of adaptation. An intermediate case is when
\begin{equation*}
\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{j})\prec\omega_{+}(z_{1-\alpha},\mathcal{F}_{j},\mathcal{F}_{J})\prec\omega_{+}(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{J}),
\end{equation*}
so that the minimax rate of adaptation is better than the worst-case minimax rate
over $\mathcal{F}_{J}$ but not as good as the minimax rate over
$\mathcal{F}_{j}$.\footnote{For positive sequences $\left\{ a_{n}\right\} $ and
$\left\{ b_{n}\right\} $, we write $a_{n} \prec b_{n}$ if
$\underset{n\to\infty}{\lim\inf}({b_{n}/}{a_{n}}) = \infty$.} That is, one can
do better than simply taking the most conservative parameter space as the true
space but not quite as good as knowing that the true function actually lies in
the smaller parameter space. Hence, the minimax adaptation rate plays an
important role in determining whether sharp adaptation is possible. In Section
\ref{sec:Asymptotic-rate-results}, we derive the minimax rates of adaptation under
the model given in Section \ref{sec:Setup}.
\subsection{Construction of Adaptive CIs}
\label{sec:constr-adapt-cis}
\citet{cai2004adaptation} provide a general method of constructing adaptive CIs
of $Lf$ under the general model \eqref{eq:npr_model}. Here, we provide an
alternative method that is intuitive and gives smaller constants in the case of
non-nested parameter spaces.\footnote{For a given adaptive CI, $CI^{\ast}$, we
refer to the positive number $c$ (that does not depend on $n$) such that
$ \sup_{f\in\mathcal{F}_{j}}E\mu(CI^{\ast})\leq
c\,\omega_{+}\left(z_{\alpha},\mathcal{F}_{j},\mathcal{F}_{J}\right), $ as the
``constant'' of $CI^{\ast}$.} For the nested case, the CI of
\cite{cai2004adaptation} has a bounded constant even as $J \to \infty$, which is
an attractive theoretical property. For the CI we propose, the constant will
grow with $J$ in general. In practice, however, one can only adapt to finitely
many parameter spaces due to computational constraints. The proposed procedure
gives a smaller constant than that of \cite{cai2004adaptation} even for unrealistically large values of $J$ (e.g.,
$J = 10^{10}$).
The main building block for our adaptive CI is the minimax one-sided CI proposed
by \citet{armstrong2018optimal}, which relies on the ordered modulus. We say
that $(f_{j}, f_{k}) \in \mathcal{F}_{j} \times \mathcal{F}_{k} $ is a solution
to $\omega(\delta, \mathcal{F}_{j}, \mathcal{F}_{k})$ if $(f_{j}, f_{k})$ solves
the optimization problem corresponding to
$\omega(\delta, \mathcal{F}_{j}, \mathcal{F}_{k})$. Let
$(f_{J,{\delta}}^{*,Jj},
g_{j,{\delta}}^{*,Jj})\in\mathcal{F}_{J}\times\mathcal{F}_{j}$ be a solution to
the ordered modulus
$\omega\left({\delta},\mathcal{F}_{J},\mathcal{F}_{j}\right),$ and define the
estimator
\begin{align}
\begin{aligned}
\hat{L}_{\delta}^{\ell,j} = &\frac{1}{2}\, L({g_{j,\delta}^{*,Jj}
+f_{J,\delta}^{*,Jj}}) \\
&+
\frac{\omega'\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)}{\delta}\sum_{i=1}^{n}(g_{j,\delta}^{*,Jj}(x_{i})-f_{J,\delta}^{*,Jj}(x_{i}))\Bigg(\frac{y_{i}}{\sigma(x_{i})}-\frac{g_{j,\delta}^{*,Jj}(x_{i})+f_{J,\delta}^{*,Jj}(x_{i})}{2}\Bigg),
\end{aligned}
\label{eq:lhat_lower}
\end{align}
where $\omega'(\cdot, \mathcal{F}_{J}, \mathcal{F}_{j})$ is the derivative of $\omega(\cdot, \mathcal{F}_{J}, \mathcal{F}_{j})$. Based on this estimator, define a lower one-sided CI by subtracting the maximum
bias and an appropriately scaled normal quantile:
\begin{equation}\label{eq:opt lower CI}
\hat{c}_{\alpha, \delta}^{\ell,j}:= \hat{L}_{\delta}^{\ell,j}-\frac{1}{2}\omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)+\frac{1}{2}\delta\omega'\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)-z_{1-\alpha}\omega'\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right).
\end{equation}
The following theorem from \citet{armstrong2018optimal} shows that for a
specific choice of $\delta$, this CI is optimal in the sense that it achieves
$\ell^{\ast}_{j,J}.$
\begin{lemma}[Theorem 3.1 of \cite{armstrong2018optimal}]
\label{thm:ArmKol}Let $\underline{\delta}=z_{\beta}+z_{1-\alpha}$. Then,
\begin{equation*}
\underset{f\in\mathcal{F}_{j}}{\sup}q_{f,\beta}(Lf-\hat{c}_{\alpha,
\underline{\delta}}^{\ell,j})= \ell^{\ast}_{j,J}= \omega(\underline{\delta},\mathcal{F}_{J},\mathcal{F}_{j}).
\end{equation*}
\end{lemma}
The excess length $Lf-\hat{c}_{\alpha, \underline{\delta}}^{\ell,j}$ follows a
Gaussian distribution because it is a affine transformation of the data, which
follows a Gaussian distribution by assumption. Hence, the median and mean of the
excess length are the same. Taking $\beta=1/2$, we can replace $q_{f,\beta}$
with the expectation under $f$, which gives
\begin{equation}
\underset{f\in\mathcal{F}_{j}}{\sup}\operatorname{\mathbf{E}}_{f}\left(Lf-\hat{c}_{\alpha}^{\ell,j}\right)
= \omega\left(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{j}\right),\label{eq:exp_ex_len}
\end{equation}
where we define
$\hat{c}_{\alpha}^{\ell,j}:= \hat{c}_{\alpha, z_{1-\alpha}}^{\ell,j}.$ Likewise,
we can define an optimal upper one-sided CI
$(-\infty, \hat{c}_{\alpha, \underline{\delta}}^{\ell,j}]$ such that
\begin{equation}
\sup q_{f,\beta}(\hat{c}_{\alpha, \underline{\delta}}^{u,j}-Lf) =
u^{\ast}_{j,J} = \omega(\underline{\delta},\mathcal{F}_{j},\mathcal{F}_{J})\label{eq:one-sd-optim-up},
\end{equation}
where the precise definition of $\hat{c}_{\alpha, \underline{\delta}}^{\ell,j}$
is given in Appendix \ref{sec:defint-optim-upper}. Similarly, let
$\hat{c}_{\alpha}^{u,j}$ denote the upper counterpart of
$\hat{c}_{\alpha}^{\ell,j}.$
Using the optimal one-sided CIs, we first show how a naive Bonferroni procedure
leads to a two-sided adaptive CI. We then provide a method that improves upon
this naive Bonferroni CI by taking into account the correlation among the
CIs. The naive Bonferroni CI is defined as
\begin{equation}\label{eq:naive_Bon_CI}
CI_{\alpha}^{Bon,J}:=\cap_{j=1}^{J}[\hat{c}_{\alpha/2J}^{\ell,j},\hat{c}_{\alpha/2J}^{u,j}].
\end{equation}
This has coverage at least $1-\alpha$ over $\mathcal{F}_{J}$ because each
$[\hat{c}_{\alpha/2J}^{\ell,j},\hat{c}_{\alpha/2J}^{u,j}]$ has coverage
$1 - \alpha/J$ over $\mathcal{F}_{J}$ and $CI_{\alpha}^{Bon,J}$ is simply the
intersection of such CIs. The following theorem shows that this CI is indeed
adaptive.
\begin{theorem}
\label{thm:BonfCI} For any $j=1,\dots,J$, we have
\begin{equation}
\sup_{f\in\mathcal{F}_{j}}\operatorname{\mathbf{E}} \mu(CI_{\alpha}^{Bon,\mathcal{J}})\leq\frac{2z_{1-\frac{\alpha}{2J}}}{z_{1-\frac{\alpha}{2}}}\,\omega_{+}(z_{1-\frac{\alpha}{2}},\mathcal{F}_{j},\mathcal{F}_{J}).\label{eq:BonBound}
\end{equation}
\end{theorem}
The constant $2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$ increases with
the number of parameter spaces $J$.\footnote{The constant,
$z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$, grows with $J$ at the rate
$(\log J)^{1/2}$. This is the same rate that \cite{cai2004adaptation} find in
their analysis of the case with non-nested parameter spaces. Their constant is
at least eight times greater than what we provide here, but does not require
that the largest space in consideration is convex.} On the other hand, the
constant given in \citet{cai2004adaptation} is 16 and thus does not depend on
the number of parameter spaces. However, we note that
$2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}$ is not too large, in fact
smaller than $16$, for reasonable specifications of $J.$ For example, when
$\alpha=0.05$ and $J=50$, we get
$2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}\approx3.36$, which is
considerably smaller than the constant given in \citet{cai2004adaptation}. Even
for unrealistically large $J$ such as $J=10^{10}$, we have
$2z_{1-\frac{\alpha}{2J}}/z_{1-\frac{\alpha}{2}}<8$, which is still less than
half of the constant given by \citet{cai2004adaptation}. Simulation results
given in Section \ref{sec:Simulation} confirm that not only the upper bound, but
also the actual length itself is often much shorter for our CI.
\begin{remark}
Suppose one is interested in constructing the one-sided CI in an adaptive
way. Note that Lemma \ref{thm:ArmKol} implies that any one-sided CI
$[\hat{c}_{\alpha}^{\ell,J},\infty)$ with coverage probability $1-\alpha$
should satisfy
\[
\underset{f\in\mathcal{F}_{j}}{\sup}\operatorname{\mathbf{E}} (Lf-\hat{c}_{\alpha}^{\ell,J})\geq\omega\left(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{j}\right).
\]
Define $\hat{c}_{\alpha}^{\ell,J}=\max _{j}\hat{c}_{\alpha/J}^{\ell,j}.$ Then,
by an analogous argument to Theorem \ref{thm:BonfCI}, we have
\[
\underset{f\in\mathcal{F}_{j}}{\sup}\operatorname{\mathbf{E}} (Lf-\hat{c}_{\alpha}^{\ell,J})\leq\frac{z_{1-\frac{\alpha}{J}}}{z_{1-\alpha}}\omega\left(z_{1-\alpha},\mathcal{F}_{J},\mathcal{F}_{j}\right).
\]
Therefore, $\left[\hat{c}_{\alpha}^{\ell,J},\infty\right)$ is an adaptive
one-sided CI in a similar sense with the two-sided case.
\end{remark}
The naive CI given in \eqref{eq:naive_Bon_CI} does not take into account the
possible correlation among the CIs that we take the intersection of. However, if
parameter spaces are ``close'' to each other, the corresponding CIs will be
correlated, implying that there is room for improvement over the Bonferroni
procedure. Consider the CIs of the form
$CI^{\tau,\mathcal{J}}=\cap_{j=1}^{J}[\hat{c}_{\tau}^{\ell,j},\hat{c}_{\tau}^{u,j}]$. If
we take $\tau = \alpha/(2J)$, this is precisely the CI given in
\eqref{eq:naive_Bon_CI}. The CI that gives the smallest constant among CIs of
such forms is $CI^{\tau^{\ast},\mathcal{J}},$ where $\tau^{\ast}$ is the largest
possible $\tau$ such that $CI^{\tau,\mathcal{J}}$ has correct coverage over
$\mathcal{F}_{J}$:
\begin{equation*}
\tau^{\ast}:= \sup_{\tau} \tau \,\, \text{ s.t. } \inf_{f \in
\mathcal{F}_{J}}\operatorname{\mathbf{P}} _{f}(Lf \in CI^{\tau,\mathcal{J}}) \geq 1- \alpha.
\end{equation*}
We know that $\tau = \alpha/(2J)$ satisfies the constraint, and also that any
$ \tau > \alpha$ does not because then $[\hat{c}_{\tau}^{\ell,j}, \infty)$ will
have coverage probability $1-\tau < 1 - \alpha$. Hence, we can restrict $\tau$
to lie in $[\alpha/(2J), \alpha]$.
However, the coverage probability
$ \inf_{f \in \mathcal{F}_{J}}\operatorname{\mathbf{P}} _{f}(Lf \in CI^{\tau,\mathcal{J}})$ is unknown in
general, rendering $CI^{\tau^{\ast},\mathcal{J}}$ infeasible. Instead, we
replace this coverage probability with a lower bound that we can calculate
either analytically or via simulation. Then, we take $\tau^{\ast}$ as the largest value that makes this lower bound at least $1-\alpha$. As we show later, using $\tau^{\ast}$ rather than $\alpha/(2J)$ can only make the resulting CI shorter.
Let $(V(\tau)', W(\tau)')'$ be a centered Gaussian random vector with unit
variance. The covariance terms for
$V(\tau)=\left(V_{1}(\tau),...,V_{J}(\tau)\right)' $ is given by
\[
\text{Cov}\left(V_{j}(\tau),V_{\ell}(\tau)\right)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*,
Jj}(x_{i})-f_{J,z_{1-\tau}}^{*,
J j}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*,
J\ell}(x_{i})-f_{J,z_{1-\tau}}^{*, J \ell}(x_{i})\big).
\]
Likewise, the covariance terms for
$ W(\tau)=\big(W_{1}(\tau),...,W_{J}(\tau)\big)'$ is given by
\[
\text{Cov}\big(W_{j}(\tau),W_{\ell}(\tau)\big)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*,
jJ}(x_{i})-f_{J,z_{1-\tau}}^{*,
jJ}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*, \ell
J}(x_{i})-f_{J,z_{1-\tau}}^{*, \ell J}(x_{i})\big).
\]
Finally, the covariance terms across $V(\tau)$ are $W(\tau)$ given as
\[
\text{Cov}\big(V_{j}(\tau),W_{\ell}(\tau)\big)=\frac{1}{z_{1-\tau}^{2}}\sum_{i=1}^{n}\big(g_{j,z_{1-\tau}}^{*,
Jj}(x_{i})-f_{J,z_{1-\tau}}^{*,
Jj}(x_{i})\big)\big(g_{\ell,z_{1-\tau}}^{*,
\ell J}(x_{i})-f_{J,z_{1-\tau}}^{*, \ell J}(x_{i})\big).
\]
This Gaussian random vector can be used to tune the critical value, as the
following lemma implies.
\begin{lemma}
\label{lem:Bonf-improv}
Let $\tau^{*}\in\left[\frac{\alpha}{2J},\alpha\right]$ to be the largest value
of $\tau$ such that
\begin{equation}
\operatorname{\mathbf{P}} \left(\max\left\{ V\left(\tau\right)',W\left(\tau\right)'\right\}
>z_{1-\tau}\right) \leq \alpha\label{eq:VWmaxCond}.
\end{equation}
Then, we have
$\sup_{f\in\mathcal{F}_{J}}\operatorname{\mathbf{P}} \left(f(0)\notin
CI_{\tau^{*}}^{\mathcal{J}}\right)\leq\alpha$.
\end{lemma}
Such a $\tau^{\ast}$ always exists because the inequality \eqref{eq:VWmaxCond}
holds with $\tau = \alpha/(2J)$ due to the union bound. A solution $\tau^{*}$
can be found via numerical simulation. By construction, its length will be also
bounded by (\ref{eq:BonBound}). In Section \ref{sec:constr-adapt-ci}, we show
that as $n \to \infty$ the distribution of $(V(\tau)', W(\tau)')'$ does not
depend on $\tau$, under our setting of $Lf = f(0)$ with $f$ belonging to a
H\"older class. Hence, finding $\tau^{*}$ boils down to simply finding the
$1-\alpha$ quantile of the maximum of a Gaussian vector in this case.
\section{Adaptive Inference for $f(0)$}
\label{sec:adapt-infer-f0}
In this section, we provide an adaptive inference procedure for $f(0)$. To
construct the adaptive CI introduced in Section \ref{sec:adapt-conf-interv}, we
first solve the corresponding modulus problem. By using this solution to the
modulus problem, we derive the minimax rate of adaptation. Finally, we provide
a CI that obtain this rate, using the method described in Section
\ref{sec:adapt-conf-interv}.
\subsection{Solution to the Modulus Problem}
\label{sec:solut-modul-probl}
Let
$\Lambda_{+,\mathcal{V}}(\gamma_{j},C_{j})\subset\Lambda_{+,\mathcal{V}}(\gamma_{J},C_{J})$
with $\gamma_{j}\geq\gamma_{J}$ and $C_{j}\leq
C_{J}$.
To construct the adaptive CI, we first calculate the ordered moduli,
$\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{J},C_{J}\right)\right)$
and
$\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{J},C_{J}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right)\right),$
for each $j=1,\dots,J.$ For notational simplicity, we consider the case with
$J=2$ and solve
$\omega_{+}\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)$, from which the general solution follows immediately.
Recall the definition of the ordered modulus of continuity
\begin{align*}
&\sup
f_{2}(0)-f_{1}(0) \\
\text{s.t. }
&\sum_{i=1}^{n}\left(\left(f_{2}(x_{i})-f_{1}(x_{i})\right)/\sigma(x_{i})\right)^{2}\leq\delta^{2}, \,\,
f_{j}\in\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right) \text{ for }j = 1,2,
\end{align*}
with the maximized value denoted by
$\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right).$
It is convenient to solve the inverse modulus problem instead, which is defined
as
\begin{align}
\begin{aligned}
\inf \,
&\sum_{i=1}^{n}\left(\left(f_{2}(x_{i})-f_{1}(x_{i})\right)/\sigma(x_{i})\right)^{2} \\
\text{s.t. }& f_{2}(0)-f_{1}(0)=b,\,\,
f_{j}\in\Lambda_{+,\mathcal{V}}\left(\gamma_{j},C_{j}\right) \text{ for }j =
1,2,
\end{aligned}\label{eq:InvMod}
\end{align}
for $b>0$, with the square root of the maximized value denoted by the inverse
(ordered) modulus
$
\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right).$
We provide a closed form solution for the this inverse problem, from which we can recover the solution to the
original problem by finding $b$ such that
$\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)=\delta$.
Note that this is simply a search problem on the positive real line.
To characterize the solution to \eqref{eq:InvMod}, we show two simple lemmas
about the properties of the class
$\Lambda_{+,\mathcal{V}}\left(\gamma,C\right)$. For
$z=(z_{1},\dots,z_{k})\in\mathbb{R}^{k}$, define
\[
\left(z\right)_{\mathcal{V}+}=\begin{cases}
\max\left\{ z_{i},0\right\} & i\in\mathcal{V}\\
z_{i} & i\notin\mathcal{V}
\end{cases}
\]
and $\left(z\right)_{\mathcal{V}-}=\left(-z\right)_{\mathcal{V}+}.$
\begin{lemma}
\label{lem:h_plus} Suppose Assumption \ref{assu:norm} holds, and let
$\gamma\in[0,1]$ and $C>0$. Define
\begin{equation*}
h_{+}(x) = C\left\lVert (x)_{\mathcal{V}+}\right\rVert ^{\gamma}\quad\text{and}\quad
h_{-}(x) = -C\left\lVert (x)_{\mathcal{V}-}\right\rVert ^{\gamma}.
\end{equation*}
Then, $h_{+},\ h_{-}\in\Lambda_{+,\mathcal{V}}\left(\gamma,C\right)$.
\end{lemma}
The following lemma asserts that the class of functions we consider is closed
under the maximum operator.
\begin{lemma}
\label{lem:holder_max_holder}Suppose
$h_{1},h_{2}\in\Lambda_{+,\mathcal{V}}\left(\gamma,C\right)$. Then,
$\max \left\{ h_{1},h_{2}\right\}
\in\Lambda_{+,\mathcal{V}}\left(\gamma,C\right)$.
\end{lemma}
The next lemma can be used to establish the solutions to the problem
(\ref{eq:InvMod}). This is a generalization of Proposition 4.1 of
\citet{beliakov2005monotonicity}, which gives the same result for the special
case of $\gamma=1$.
\begin{lemma}
\label{lem:InvModMax}Given $f_{0}\in\mathbb{R}$ and $0<\gamma\leq1$, define
\[
\Lambda_{+,\mathcal{V}}^{f_{0}}\left(\gamma,C\right)=\left\{
f\in\Lambda_{+,\mathcal{V}}\left(\gamma,C\right):\ f(0)=f_{0}\right\} .
\]
Then, for any $x\in\mathbb{R}^{k}$, we have
\begin{eqnarray*}
\underset{f\in\Lambda_{+,\mathcal{V}}^{f_{0}}\left(\gamma,C\right)}{\max }f(x) & = & f_{0}+C\left\lVert (x)_{\mathcal{V}+}\right\rVert ^{\gamma}\\
\underset{f\in\Lambda_{+,\mathcal{V}}^{f_{0}}\left(\gamma,C\right)}{\min }f(x) & = & f_{0}-C\left\lVert (x)_{\mathcal{V}-}\right\rVert ^{\gamma}.
\end{eqnarray*}
\end{lemma}
We are now ready to characterize the solution to the inverse modulus problem
(\ref{eq:InvMod}). For $r\in\mathbb{R}$, define
$\left(r\right)_{+}:= \max\left\{ r,0\right\} $.
\begin{prop}
\label{prop:InvModRes} Suppose Assumption \ref{assu:norm} holds, and define
\begin{eqnarray*}
f_{1}^{*}(x) & = & C_{1}\left\lVert \left(x\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{1}}\\
f_{2}^{*}(x) & = & \max \left\{ b-C_{2}\left\lVert \left(x\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{2}},\ C_{1}\left\lVert \left(x\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{1}}\right\} .
\end{eqnarray*}
Then, $\left(f_{1}^{*},f_{2}^{*}\right)$ solves the inverse modulus problem
(\ref{eq:InvMod}), and the inverse modulus is given by
\begin{align}
\begin{aligned}
&\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right) \\
= &\Big( \sum\nolimits_{i=1}^{n}\left(\left(b-C_{1}\left\lVert
\left(x_{i}\right)_{\mathcal{V}+}\right\rVert
^{\gamma_{1}}-C_{2}\left\lVert
\left(x_{i}\right)_{\mathcal{V}-}\right\rVert
^{\gamma_{2}}\right)/\sigma\left(x_{i}\right)\right)_{+}^{2} \Big)^{1/2}.
\end{aligned}
\end{align}
\end{prop}
\begin{proof}
To solve (\ref{eq:InvMod}), note that it is without loss of generality to
restrict attention to the functions with $f_{1}(0)=0$ and $f_{2}(0)=b$, which
is satisfied by $f_{1}^{*}$ and $f_{2}^{*}$. To simplify notation, write
$\mathcal{F}_{1}^{0}=\Lambda_{+,\mathcal{V}}^{0}\left(\gamma_{1},C_{1}\right)$
and
$\mathcal{F}_{2}^{b}=\Lambda_{+,\mathcal{V}}^{b}\left(\gamma_{2},C_{2}\right)$.
Since $f_{2}(0)>f_{1}(0)$, we want
$ f_{1}(x) = \max_{f\in\mathcal{F}_{1}^{0}}f(x)$ and
$ f_{2}(x) = \min_{f\in\mathcal{F}_{2}^{b}}f(x)$ as long as $x\in\mathcal{X}$
satisfies $\min_{f\in\mathcal{F}_{2}} f(x)\geq \max_{f\in\mathcal{F}_{1}}f(x)$, and
$f_{1}(x)=f_{2}(x)$ otherwise. Note that $f_{1}^{*}$ and $f_{2}^{*}$ are designed
exactly to achieve this goal, which follows by Lemma \ref{lem:InvModMax}.
It remains to check whether $f_{1}^{*}\in\mathcal{F}_{1}$ and
$f_{2}^{*}\in\mathcal{F}_{2}$. The former case is trivial. For the latter
case, note that
$f_{1}^{\ast}\in\Lambda_{+,\mathcal{V}}(\gamma_{1},C_{1})\subseteq\Lambda_{+,\mathcal{V}}(\gamma_{2},C_{2})$.
Now, by Lemma \ref{lem:holder_max_holder}, we have
$f_{2}^{*}\in\mathcal{F}_{2}.$
\end{proof}
The following corollary states an analogous result regarding the inverse modulus
$\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right)$.
\begin{corollary}
\label{cor:InvModRes(21)} Define
\begin{eqnarray*}
g_{1}^{*}(x) & = & b-C_{1}\left\lVert \left(x\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{1}}\\
g_{2}^{*}(x) & = & \min \left\{ b-C_{1}\left\lVert \left(x\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{1}},\ C_{2}\left\lVert \left(x\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{2}}\right\} .
\end{eqnarray*}
Then, $\left(g_{1}^{*},g_{2}^{*}\right)$ solves the inverse modulus
$\omega^{-1}\left(b,\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right)$.
\end{corollary}
\begin{remark}[Role of Assumption \ref{assu:norm}]\label{rem:norm} Proposition
\ref{prop:InvModRes} requires Assumption \ref{assu:norm} due to the specific
form of monotonicity we consider. By considering coordinate-wise monotonicity, we
must take a norm that is ``aligned'' with this direction of monotonicity. The
assumption precisely imposes this. This is a unique feature that arises in the
multivariate setting. To allow for more general norms, let $\mathcal{B}$ be an
orthonormal basis of $\mathbb{R}^{k}$, and denote by $z^{\mathcal{B}}$ the
coordinate vector of $z \in \mathbb{R}^{k}$ with respect to $\mathcal{B}$ and
$z^{\mathcal{B}}_{j}$ its $j$th component. Suppose the regression function is
monotone in the coefficients with respect to this basis $\mathcal{B}$, so that
the monotone H\"older class is given as
\begin{equation*}
\Lambda_{+,\mathcal{V}}(\gamma,C):=\left\{ f\in\Lambda(\gamma,C):f(x)\geq
f(z)\,\,\text{if }x^{\mathcal{B}}_{j}\geq z^{\mathcal{B}}_{j}\text{ }\forall
j \in\mathcal{V} \text{ and }x^{\mathcal{B}}_{j}\geq z^{\mathcal{B}}_{j}\,\,\forall
j\notin\mathcal{V}\right\}.
\end{equation*}
Then, the condition we want to impose on the norm $\lVert \cdot \rVert$ is
monotonicity with respect to the magnitude of $z^{\mathcal{B}}_{j}$. A special
case is the Mahalanobis distance.
\end{remark}
\subsection{Minimax Rate of Adaptation}
\label{sec:Asymptotic-rate-results}
Using this solution to the inverse modulus, we derive the rate of convergence of
the between class of modulus, which characterizes how fast the worst-case
expected length of the adaptive CIs can go to 0 as $n\rightarrow\infty$. We
derive the rates under the assumption that the sequence of design points
$\{x_{i}\}_{i=1}^{\infty}$ is a realization of a sequence of independent and
identically distributed random vectors $\{ X_{i} \}_{i=1}^{\infty}$ drawn from a
distribution that satisfies some mild regularity conditions. This gives an
intuitive restriction on the design points, and also shows that the result
applies under random design points as well.\footnote{Consider the model
$y_{i}=f(X_{i})+\varepsilon_{i}$, for $i=1,\dots,n$, with the
$X_{i}\overset{\mathrm{i.i.d.}}{\sim}p_{X}$ with
$\varepsilon_{i}|X_{i}\sim N(0,\sigma^{2}(X_{i}))$. Then, conditional on
$\{ X_{i} \}_{i=1}^{n}=\{x_{i}\}_{i=1}^{n}$, this model is equivalent with our
model.} Define
$r(\gamma_{1}, \gamma_{2}) =
({2+k_{+}/\gamma_{1}+(k-k_{+})/\gamma_{2}})^{-1}$. The following theorem fully
characterizes the minimax rate of adaptation.
\begin{theorem}\label{thm:rate_result}
Let $\left\{ X_{i}\right\} _{i=1}^{\infty}$ be an i.i.d. sequence of random
vectors with support $\mathcal{X}$. Suppose $X_{i}$ admits a probability
density function $p_{X}(\cdot)$ that is continuous at $0$ with $p_{X}(0)>0,$
and assume $\sigma(\cdot) = 1$. Then, for almost all realizations
$\left\{ x_{i}\right\} _{i=1}^{\infty}$ of
$\left\{ X_{i}\right\} _{i=1}^{\infty}$ and for all $\delta > 0$, we have
\begin{align*}
\lim_{n\to\infty}n^{r(\gamma_{1},\gamma_{2})}\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)
& =\left(\delta^{2}/c_{1,2}^{\ast}\right)^{r(\gamma_{1}, \gamma_{2})},\text{ and}\\
\lim_{n\to\infty}n^{r(\gamma_{1},\gamma_{2})}\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right) & =\left(\delta^{2}/c_{2,1}^{\ast}\right)^{r(\gamma_{1}, \gamma_{2})},
\end{align*}
where $c_{1,1}^{\ast}$ and $c_{2,1}^{\ast}$ are constants that depend only on
the function spaces.
\end{theorem}
\begin{remark}
The result immediately implies the rate of convergence for the between class
modulus
\begin{align*}
& \lim_{n\to\infty}n^{r(\gamma_{1},\gamma_{2})}\omega_{+}\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)\right)\\
= & \max\left\{ \delta^{2}/c_{1,2}^{\ast},\delta^{2}/c_{2,1}^{\ast}\right\} ^{r(\gamma_{1},\gamma_{2})}.
\end{align*}
Hence, if a CI maintains coverage over
$\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)$, the best possible
worst-case length of this CI over
$\Lambda_{+,\mathcal{V}}\left(\gamma_{2},C_{2}\right)$ goes to $0$ at the same
rate as $n^{-r(\gamma_{1},\gamma_{2})}$.
\end{remark}
\begin{remark}[Heteroskedasticity]
For simplicity, the theorem imposes a homoskedasticity condition (i.e.,
$\sigma(\cdot) = 1$). However, allowing for general $\sigma(\cdot)$ is
straightforward and requires only weak regularity conditions on
$\sigma(\cdot)$. See Appendix \ref{apdx:Heteroskedasiticy} for details.
\end{remark}
Theorem \ref{thm:rate_result} shows how the monotonicity restriction plays a
role in determining the minimax rates of adaptation to H\"older coefficients
under the multivariate nonparametric regression setting. When $k_{+}=k$, the
minimax rate of adaptation is $n^{-\frac{1}{2+k/\gamma_{1}}}$, which equals the
minimax convergence rate over
$\omega\left(\delta,\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right),\Lambda_{+,\mathcal{V}}\left(\gamma_{1},C_{1}\right)\right)$.
This shows that strong adaptation is possible if the regression function is
monotone with respect to all the variables, just like in the univariate case. On
the other hand, when $k_{+}=0$, the rate becomes
$n^{-\frac{1}{2+k/\gamma_{2}}}$, consistent with the previous findings that
there is no scope of adaptation for general H\"older classes without any shape
constraint. Importantly, Theorem \ref{thm:rate_result} characterizes the
convergence rate for the case where $0<k_{+}<k$, where it gives an intuitive
intermediate rate between the two extreme.
\subsection{Construction of the Adaptive CI}
\label{sec:constr-adapt-ci}
Here, we give the explicit formula of the CIs for our parameters of interest,
now that we have derived the form of the moduli of continuity and the solutions
to the modulus problems in the previous section. We first consider
$L_{0}f$. Before stating the result, it is convenient to define the following
functions
\begin{align*}
D_{Jj,\delta}(x_{i}) := &
(\omega(\delta,\mathcal{F}_{J},\mathcal{F}_{j})-C_{j}\lVert
(x_{i})_{\mathcal{V}-}\rVert
^{\gamma_{j}}-C_{J}\lVert
(x_{i})_{\mathcal{V}+}\rVert ^{\gamma_{J}})_{+},
\text{ and}\\
D_{jJ, \delta}(x_{i}) := &( \omega(\delta,\mathcal{F}_{j},\mathcal{F}_{J})-C_{J}\lVert (x_{i})_{\mathcal{V}-}\rVert ^{\gamma_{J}}-C_{j}\lVert (x_{i})_{\mathcal{V}+}\rVert ^{\gamma_{j}})_{+}.
\end{align*}
\begin{corollary}
\label{cor:chat formula}For $Lf=f(0)$ and $\beta=1/2$, the lower CI defined in \eqref{eq:opt lower CI} is given by
\[
\hat{c}_{\delta}^{\ell,j}=\hat{L}_{\delta}^{\ell,j}-\frac{1}{2}\left( \omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)+\frac{\delta^{2}}{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)} \right),
\]
where
\begin{align*}
\hat{L}_{\delta}^{\ell,j} = &\, \frac{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)y_{i}}{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)}+\frac{\omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)}{2}\\
& -\frac{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)(\omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right)-C_{j}\left\lVert \left(x_{i}\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{j}}+C_{J}\left\lVert \left(x_{i}\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{J}})}{2\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)}.
\end{align*}
Likewise, the upper bound of the CI is given by
\[
\hat{c}_{\delta}^{u,j}=\hat{L}_{\delta}^{u,j}+\frac{1}{2}\left(\omega\left(\delta,\mathcal{F}_{j},\mathcal{F}_{J}\right)+\frac{\delta^{2}}{\sum_{i=1}^{n}D_{jJ, \delta}\left(x_{i}\right)}\right),
\]
where
\begin{align*}
\hat{L}_{\delta}^{u,j} = & \, \frac{\sum_{i=1}^{n}D_{jJ, \delta}\left(x_{i}\right)y_{i}}{\sum_{i=1}^{n}D_{jJ, \delta}\left(x_{i}\right)}+\frac{\omega\left(\delta,\mathcal{F}_{j},\mathcal{F}_{J}\right)}{2}\\
& -\frac{\sum_{i=1}^{n}D_{jJ, \delta}\left(x_{i}\right)(\omega\left(\delta,\mathcal{F}_{j},\mathcal{F}_{J}\right)-C_{J}\left\lVert \left(x_{i}\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{J}}+C_{j}\left\lVert \left(x_{i}\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{j}})}{2\sum_{i=1}^{n}D_{jJ, \delta}\left(x_{i}\right)}.
\end{align*}
\end{corollary}
The first terms in the formula of $\hat{L}_{\delta}^{\ell,j}$ and
$\hat{L}_{\delta}^{u,j}$ are the random terms linear in $y_{i}$ while the
remaining terms are non-random fixed terms. If $\mathcal{V}=\{1,...,k\}$ (so
the function is monotone in every coordinate), the random terms can be viewed
as a kernel estimator with a data-dependent bandwidth. Too see this, if we
define
\[
k(x)=\left[1-C_{j}\left\lVert \left(x\right)_{\mathcal{V}-}\right\rVert ^{\gamma_{j}}-C_{J}\left\lVert \left(x\right)_{\mathcal{V}+}\right\rVert ^{\gamma_{J}}\right]_{+},
\]
and
\[
h_{mn}\left(x\right)=\begin{cases}
\omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right){}^{1/\gamma_{J}} & \text{if the }m\text{th coordinate of }x\geq0\\
\omega\left(\delta,\mathcal{F}_{J},\mathcal{F}_{j}\right){}^{1/\gamma_{j}}
& \text{otherwise},
\end{cases}
\]
we have
\[
\frac{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)y_{i}}{\sum_{i=1}^{n}D_{Jj,\delta}\left(x_{i}\right)}=\frac{\sum_{i=1}^{n}k\left(x_{1i}/h_{1n}\left(x_{i}\right),...,x_{ki}/h_{kn}\left(x_{i}\right)\right)y_{i}}{\sum_{i=1}^{n}k\left(x_{1i}/h_{1n}\left(x_{i}\right),...,x_{ki}/h_{kn}\left(x_{i}\right)\right)}.
\]
Hence, the CI can be considered to be based on a Nadaraya-Watson type
estimator, correcting for the bias.
As described in Section \ref{sec:constr-adapt-cis}, the proposed CI is given by
$ \cap_{j=1}^{J} [\hat{c}_{z_{1-\tau^{\ast}}}^{\ell,j},
\hat{c}_{z_{1-\tau^{\ast}}}^{\ell,j}]$, where $\tau^{\ast}$ is defined in Lemma
\ref{lem:Bonf-improv}. Here, we show that the distribution of
$(V(\tau)',W(\tau)')'$ does not depend on $\tau$ as $n \to \infty$. The
implication of this invariance with respect to $\tau$, is that calculating
$\tau^{\ast}$ boils down to calculating the quantile of the maximum of Gaussian
vectors. The variance matrix of this limiting Gaussian random vector is known,
and thus the said quantile can be easily simulated. Moreover, when
$\gamma_{1} = \cdots = \gamma_{J}$ so that the parameters spaces differs only in
$C_{j}$, $\tau^{\ast}$ can be shown to be bounded away from zero by a constant
that does not depend on $J$, for large $n$. Hence, the constant of the CI does
not grow to infinity as $J \to \infty$ in this case.\footnote{This is especially
useful when one wishes to adapt to $C$ while keeping $\gamma$ fixed. For
example, \cite{kwon2020InferenceRegressionDiscontinuity} take $\gamma_{j} = 1$
and consider the problem of adapting to the Lipschitz constant in a regression
discontinuity setting.}
\begin{lemma}\label{lem:tau}
Under the same set of conditions of Theorem \ref{thm:rate_result},
$(V\left(\tau\right)',W\left(\tau\right)' )' \overset{d}{\to}
(V_{\infty}',W_{\infty}')'$ as $n \to \infty$, where
$(V_{\infty}',W_{\infty}')'$ is a Gaussian random vector that does not depend
on $\tau.$ Moreover, if $\gamma_{1}= \cdots = \gamma_{J}$, then, for large
$n$, $\tau^{\ast} > \eta$ for some $\eta > 0$ that does not depend on $J$.
\end{lemma}
\begin{remark}[Dependence on $J$]
The proof reveals that when all $J$ parameter spaces correspond to different
H\"older exponents (i.e., $\gamma_{1} > \cdots > \gamma_{J}$), the dependence
of $\tau^{\ast}$ on $J$ does not vanish and in fact results in CIs whose
constants grow at the same rate as the naive Bonferroni CI, $(\log
J)^{1/2}$. However, some finite sample improvement in terms of the length of
the resulting CI compared to the naive Bonferroni CI is shown in the empirical
exercise. When some of the parameter spaces have the same H\"older exponent,
the improvement can be significant. As an extreme case, when
$\gamma_{1}= \cdots = \gamma_{J}$, $\tau^{\ast}$ can is bounded away from $0$
by a constant that does not depend on $J$, which is exactly what the second
part of the lemma asserts.
\end{remark}
\section{Simulation Results}
\label{sec:Simulation}
In this section, we compare the performances of the adaptive CI of
\citet{cai2004adaptation} and the adaptive CI constructed using the naive Bonferroni procedure described in Section \ref{sec:constr-adapt-cis}. As a benchmark, we also provide the lengths
of the shortest fixed length confidence intervals of \citet{Donoho1994}, referred to as minimax CIs. We consider inference for $f(0)$, given some regression function $f$. We consider the case where the researcher is uncertain about the value of the H\"older exponent $\gamma$, and thus tries to adapt to its value.
First, we construct adaptive CIs with respect to two smoothness parameters
$\left(\gamma_{1},\gamma_{2}\right)=\left(1,10^{-3}\right)$ while fixing $C = 1$, which gives $J = 2$. We vary $n$ over $\{10^2,\ 5 \times 10^2,\ 10^3,\ 5 \times 10^3,\ 10^4 \}$
to investigate the rate of adaptation as the sample size grows. The true
regression function is over $\mathbb{R}^2$ and given by either $f_{1}$ or $f_{2}$, defined as
\begin{align*}
f_{1}(x_{1},x_{2}) = 0, \quad
f_{2}(x_{1},x_{2}) = \left\lVert \left(x_{1},x_{2}\right)_{\mathcal{V}+}\right\rVert _{2}^{\gamma_{2}},
\quad
\mathcal{V}=\{1,2\}.
\end{align*}
By construction, we have $f_{j}\in\Lambda_{+,\mathcal{V}}(\gamma_{j},1)$. The covariates are drawn from a uniform distribution over $[-1/(2\sqrt{2}),
1/(2\sqrt{2})]^{2}$, and the noise terms, $\{u_{i}\}_{i=1}^{n}$, are drawn
from a standard normal distribution. The outcome variable is given as
$y_{i} = f(x_{i}) + u_{i}$, for $f \in \{f_1, f_2\}$. We fix the draw of
$\left\{ x_{i} \right\}_{i=1}^{n}$ within each simulation iteration. We run 500 iterations to calculate the average lengths and coverage probabilities of CIs. The nominal coverage probability is $.95$ for
all CIs.
\begin{table}
\centering{}
\caption{Lengths of CIs when $f=f_{1}$ with $J = 2$
\label{tab:Lengths-of-CIs_f1}}
\begin{tabular}{@{\extracolsep{5pt}} lcccc} \\[-1.8ex]\hline
\hline \\[-1.8ex] & AdaptBonf & CL & Minimax ($\gamma = \gamma_2$) & Minimax ($\gamma = \gamma_1$) \\
\hline \\[-1.8ex] n = 100 & $0.925$ & $5.092$ & $1.459$ & $0.617$
\\ n = 500 & $0.478$ & $1.857$ & $1.212$ & $0.391$ \\
n = 1,000 & $0.382$ & $1.530$ & $1.150$ & $0.329$ \\ n = 5,000 & $0.256$ & $1.037$ & $1.064$ & $0.222$ \\ n =
10,000 & $0.212$ & $0.868$ & $1.045$ & $0.187$ \\
\hline \\[-1.8ex]
\end{tabular}
\end{table}
\begin{table}
\caption{Lengths of CIs when $f=f_{2}$ with $J = 2$
\label{tab:Lengths-of-CIs_f2}}
\centering{}
\begin{tabular}{@{\extracolsep{5pt}} lccc} \\[-1.8ex]\hline \hline \\[-1.8ex]
& AdaptBonf & CL & Minimax ($\gamma = \gamma_2$) \\ \hline \\[-1.8ex] n = 100 &
$1.599$ & $5.092$ & $1.459$ \\ n = 500 &
$1.279$ & $3.504$ & $1.212$ \\ n = 1,000 & $1.197$ & $3.222$ & $1.150$ \\ n =
5,000 & $1.086$ & $2.860$ & $1.064$ \\ n = 1,0000 &
$1.060$ & $2.703$ & $1.045$ \\ \hline
\\[-1.8ex]
\end{tabular}
\end{table}
\begin{table}
\caption{Coverage probabilities of adaptive CIs ($J = 2$)\label{tab:CovPRob}}
\centering{}
\begin{tabular}{@{\extracolsep{5pt}} lcccccc} \\[-1.8ex]\hline \hline
\\[-1.8ex] & \multicolumn{2}{c}{$f = f_1$} & \multicolumn{2}{c}{$f = f_2$} \\
\hline \\[-1.8ex] & AdaptBonf & CL & AdaptBonf & CL \\
\hline \\[-1.8ex] n = 100 & $0.988$ & $1.000$ & $0.992$ & $1.000$
\\ n = 500 & $0.972$ & $1.000$ & $0.976$ & $0.998$ \\ n =
1,000 & $0.970$ & $1.000$ & $0.982$ & $1.000$ \\ n = 5,000 & $0.972$ & $1.000$ & $0.970$ & $1.000$ \\ n = 10,000 & $0.974$ & $1.000$ & $0.976$ & $1.000$ \\
\hline \\[-1.8ex]
\end{tabular}
\end{table}
Table
\ref{tab:Lengths-of-CIs_f1} shows the results for the case where $f=f_{1}$. Each column corresponds to 1) our proposed (naive) Bonferroni adaptive procedure (AdaptBonf), 2) the adaptive CI of \cite{cai2004adaptation} (CL, henceforth), 3) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_2, 1)$, and 4) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_1, 1)$. Regarding the last two minimax procedures, we refer to them as the ``conservative minimax CI'' and the ``oracle minimax CI'', respectively. Note that the oracle minimax CI is an optimal benchmark, which is only feasible when we actually know the true regression function is in the smaller parameter space $\Lambda_{+,\mathcal{V}}(\gamma_1, 1)$. In Table \ref{tab:Lengths-of-CIs_f1}, the average lengths of both adaptive confidence intervals decrease considerably as $n$ increases from 100 to 10,000. In comparison, the length
of the conservative minimax CI (column 3) decreases only about 28\% for the same change in the
sample size. This shows the lengths of the adaptive confidence intervals decrease more sharply when the true function is smooth, as predicted by the theory.
To compare the performances of different adaptive inference procedures, note that the average lengths of the CI of CL adapting to the H\"older exponents (column 2) are often wider than the conservative minimax CI (column 3). When $n = 100$, the former is more than three times
wider than the latter, and the adaptive procedure starts to dominate the minimax procedure only when $n$ is greater than 5,000. In comparison, our proposed
Bonferroni adaptive procedure (column 1) yields shorter CIs than those by CL, as
predicted in Section \ref{sec:constr-adapt-cis}.
To compare the Bonferroni adaptive CI with the conservative minimax CI, the lengths of the former are always
exceeded by those of the minimax CI, even for the relatively small sample size of $n=100$. Moreover, the length of the adaptive CI becomes only 20\% of the length of the conservative minimax CI for the sample size of $n=10^4$.
The Bonferroni procedure also performs well even when compared to the infeasible oracle minimax CI (column 4), with the length of the former only 13\% wider than the latter when $n = 10^4$. This demonstrates the strong adaptivity property of the adaptive procedure when the regression function is monotone with respect to all variables, as shown in Section \ref{sec:Asymptotic-rate-results}.
Table \ref{tab:Lengths-of-CIs_f2} demonstrates the analogous simulation results when $f=f_{2}$. In this case, the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_2, 1)$ (column 3) is referred to as the oracle minimax CI. While
the lengths of the oracle minimax procedure are considerably shorter than the CIs of CL for
various values of $n$, the performance of the Bonferroni CIs almost matches that of
the oracle minimax procedure. Especially, the performance of the Bonferroni adaptive
procedure becomes extremely close to the oracle minimax procedure when $n$ is greater than 500.
Table \ref{tab:CovPRob} shows the coverage probabilities of
adaptive CIs for both of the cases when $f = f_1$ and $f = f_2$. While all the CIs achieve the correct coverage, none of those CIs exactly achieves the nominal coverage of $.95$, reflecting the conservative nature of the adaptive CIs. We can see that the adaptive procedure of CL is particularly conservative, almost always yielding 100\% coverage probabilities.
So far we considered adapting to the smoothness parameters at two extremes,
$\gamma\in(0.001,1)$. Since the multiplicative constant for the Bonferroni procedure increases with $J$, a concern is that the performance of the Bonferroni procedure relative to the CL procedure might get worse when $J$ is larger. To investigate the possibility, we consider adapting to a wider set of parameters, $\{\gamma_j\}_{j = 1}^6$, where $\gamma_j = 1 - (j - 1)/5$ for $j = 1,...,5$ and $\gamma_6 = 10^{-3}$. Moreover, rather than taking the extreme value of $\gamma$ as the true parameter, we consider the case where $\gamma$ takes an intermediate value, $\gamma=1/2$. The true regression function is given by
\[
f_3(x_{1},x_{2})=\left\lVert \left(x_{1},x_{2}\right)_{\mathcal{V}+}\right\rVert _{2}^{1/2},
\quad
\mathcal{V} = \{1,2\},
\]
so that $f_3 \in\Lambda_{+,\mathcal{V}}(1/2,1)$.
Table \ref{tab:Lengths-of-CIs-MultiGam_Full} displays the simulation results corresponding to this specification. Each column corresponds to 1) our proposed Bonferroni adaptive procedure, 2) the adaptive CI of CL, 3) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(\gamma_6, 1)$, and 4) the minimax CI with respect to $\Lambda_{+,\mathcal{V}}(1/2, 1)$. As before, we refer to the last two CIs as the conservative minimax CI and the oracle minimax CI, respectively. We observe the same pattern as in the case of adapting to two parameters---adaptive CIs shrink faster than the conservative minimax CI as the sample size increases, and the Bonferroni adaptive CIs are shorter than the ones of CL. While the ratio of the length of the Bonferroni CI to that of the CI of CL is larger in this case compared to the case where $J = 2$, especially when $n$ is large, the Bonferroni CI is still more than 50 \% narrower than the CI of CL, and not much wider than the oracle minimax CI.
\begin{table}
\caption{Lengths of CIs when $f=f_{3}$ with $J = 6$
\label{tab:Lengths-of-CIs-MultiGam_Full}}
\centering{}
\begin{tabular}{@{\extracolsep{5pt}} lcccc} \\[-1.8ex]\hline \hline
\\[-1.8ex] & AdaptBonf & CL & Minimax ($\gamma = \gamma_6$) & Minimax ($\gamma = 0.5$) \\
\hline \\[-1.8ex] n = 100 & $1.495$ & $5.092$ & $1.459$ & $0.908$ \\
n =
500 & $0.972$ & $2.417$ & $1.212$ & $0.649$ \\
n = 1,000 & $0.833$ & $2.441$ & $1.150$ & $0.580$ \\ n =
5,000 & $0.615$ & $1.816$ & $1.064$ & $0.447$ \\ n =
10,000 & $0.543$ & $1.145$ & $1.045$ & $0.400$ \\
\hline \\[-1.8ex]
\end{tabular}
\end{table}
\section{Empirical Illustration}
\label{sec:empir-illustr}
In this section, we apply our procedure to the production function estimation
problem for the Chinese chemical industry. Specifically, we use the firm-level data of
\citet{jacho2010identification} for the year
2001, which was also used by \citet{horowitz2017nonparametric} to illustrate
their method of constructing the uniform confidence band for the production
function under shape restrictions.
In the dataset, the dependent variable is the logarithm of value-added real
output ($y$), and the explanatory variables are the logarithms of the net value
of the real fixed asset ($k$) and the number of employees ($\ell$). After
removing the outliers for $y,k$ and $\ell$, the remaining sample size was
$n=1,636$.\footnote{We used the conventional way of outlier detection, removing
the observations that are greater than the third quantile plus IQR times 1.5,
or less than the first quantile minus IQR times 1.5. Our resulting sample size
is close to \citet{horowitz2017nonparametric}, who have $n=1,638$. } Table
\ref{tab:Summary-statistics:-Chinese} shows the brief summary of the variables
used in our analysis. We are interested in construction of the confidence
interval for
$f(k_{0},\ell_{0}) := \operatorname{\mathbf{E}} \left[y|k=k_{0},\ell=\ell_{0}\right]$. We
take $\left(k_{0},\ell_{0}\right)$ to be medians of each variable.
\begin{table}[t]
\caption{Summary statistics - Chinese chemical industry dataset for the year
2001}
\label{tab:Summary-statistics:-Chinese}
\centering
\begin{tabular}{@{\extracolsep{5pt}} lcccc}
\\[-1.8ex]\hline \hline \\[-1.8ex] & Min & Mean & Median & Max \\
\hline
\\[-1.8ex]
Log output & $6.472$ & $9.952$ & $9.937$ & $13.233$
\\
Log fixed asset & $7.463$ & $10.818$ & $10.759$ &
$14.226$ \\ Log labor & $3.664$ & $6.352$ & $6.386$ & $9.142$ \\
\hline \\[-1.8ex]
\end{tabular}
\end{table}
The first step is to estimate the variance of the error term. We assume homoskedastic errors for simplicity. The variance estimator is defined as
\[
\hat{\sigma}^{2}=\frac{\sum_{i=1}^{n}\left(y_{i}-\hat{r}(k_{i},\ell_{i})\right)}{n-2\nu_{1}+\nu_{2}},
\]
where $\hat{r}(k_{i},\ell_{i})$ is the estimator for the conditional mean using
kernel regression, $\nu_{1}=\text{tr}(L)$, $\nu_{2}=\text{tr}(L'L)$, where $L$ is
the weight matrix for the kernel estimator. Refer to \citet{wassermann2006all}
for a justification for this variance estimator. We used the Gaussian kernel with the bandwidth chosen by
expected Kullback-Leibler cross validation as in \citet{hurvich1998smoothing}.
For the function space, we consider adapting to a sequence of parameter
spaces $\{\Lambda_{+,\mathcal{V}}(\gamma_{j},C) \}_{j = 1}^6$ with $\gamma_j = 1 - (j - 1)/5$ for $j = 1,...,5$ and $\gamma_6 = 10^{-3}$. We take $\mathcal{V}=\{1,2\}$, assuming that the
production function is nondecreasing in both fixed assets and labor, which is
consistent with economic theory. To make $\Lambda_{+,\mathcal{V}}(\gamma_{j},C) \subset \Lambda_{+,\mathcal{V}}(\gamma_{6},C)$ hold for all $j = 1,...,5$, we only use observations in a restricted support, and the effective sample size is given by $n_{\text{eff}} = 272$.
For the norm, we use the Euclidean norm weighted by the inverse of the standard
deviation of each input,
$\left\lVert (k,\ell)\right\rVert
=(({k}/{s_{k}})^{2}+({\ell}/{s_{\ell}})^{2})^{1/2}$ where $s_{k}$ and
$s_{\ell}$ are standard deviations of $k$ and $\ell$, respectively. We take
conservative values of $C$ by setting
\[
C=2\times\underset{(i,j)\in\left\{ 1,...,n_{\text{eff}}\right\}
^{2}}{\max }\frac{\left\lvert y_{j}-y_{i}\right\rvert }{\left\lVert (k_{j},\ell_{j})-(k_{i},\ell_{i})\right\rVert ^{\gamma_{6}}}.
\]
\begin{table}[t!]
\caption{95\% confidence intervals for $f(k_{0},\ell_{0})$}
\label{tab:CI_empirical}
\centering
\begin{tabular}{@{\extracolsep{5pt}} lcc}
\hline
\hline
\\[-1.8ex]
& CI & Length \\
\hline \\[-1.8ex]
Minimax ($\gamma = \gamma_6$) & [7.501, 13.005] & $5.504$
\\
Minimax ($\gamma = \gamma_1$)
& [9.922, 10.484] & $0.562$ \\
AdaptBonf (Naive) & \textbf{[9.766, 11.264]}
& $\mathbf{1.498}$ \\
AdaptBonf (Calibrated) &
\textbf{[9.864, 11.188]} & $\mathbf{1.324}$ \\
Cai and Low & [7.134, 12.049] & $4.915$ \\
\hline
\\[-1.8ex]
\end{tabular}
\end{table}
We compare different procedures to construct CIs. The methods in
comparison are the minimax CI with respect to the largest space $\Lambda_{+,\mathcal{V}}(\gamma_{6},C)$ (row 1), the restricted minimax CI with respect to the smallest space $\Lambda_{+,\mathcal{V}}(\gamma_{1},C)$ (row 2), the adaptive Bonferroni
CI adapting to $\{\gamma_j\}_{j = 1}^6$ (row 3), the same adaptive CI, but
taking into account the correlations between different CIs (fourth row), and the adaptive CI of \citet{cai2004adaptation} (henceforth CL) adapting to $\{\gamma_j\}_{j = 1}^6$ (fifth row). Note that all the CIs maintain correct coverage over the largest space $\Lambda_{+,\mathcal{V}}(\gamma_{6},C)$, except for the second one, which is valid only over the smallest space $\Lambda_{+,\mathcal{V}}(\gamma_{1},C)$. We refer to the first minimax CI as the conservative minimax CI.
Table \ref{tab:CI_empirical} demonstrates the 95\% confidence intervals for
$f(k_{0},\ell_{0})$ produced by different inference methods.
First of all, the lengths of the
adaptive Bonferroni CIs are much shorter than the conservative minimax CI, while the
procedure of CL yields a wider CI, almost as long as the conservative minimax
CI. We can also observe that the adaptive Bonferroni CI using the calibrated
value of $\tau^*$ (fourth row) is relatively narrower than its naive version
taking $\tau = 0.05/2J$ (third row). Lastly, while the length of the second
minimax procedure (second row) is the shortest, it is only valid when we are confident that the true regression function is in the smallest function space we consider, $ \Lambda_{+,\mathcal{V}}(\gamma_{1},C)$. Together with the simulation results in the previous section, our empirical analysis demonstrates the advantage of using an adaptive procedure when the monotonicity restriction is plausible as well as good finite sample performance of our proposed Bonferroni adaptive procedure.
\newpage
\bibliography{KK_rfap_Nov262020}
\newpage