# Working title: “Optimal prediction bands”

Jianqing Fan, Jiawei Ge and Debarghya Mukherjee  Jianqing Fan†, Jiawei Ge† and
Debarghya Mukherjee†

†Department of Operations Research and Financial Engineering,
Princeton University  This work is partially supported by ONR grant
N00014-19-1-2120 and the NSF grants DMS-2052926, DMS-2053832, and DMS-2210833

# UTOPIA: Universally Trainable Optimal Prediction Intervals Aggregation

Jianqing Fan, Jiawei Ge and Debarghya Mukherjee  Jianqing Fan†, Jiawei Ge† and
Debarghya Mukherjee†

†Department of Operations Research and Financial Engineering,
Princeton University  This work is partially supported by ONR grant
N00014-19-1-2120 and the NSF grants DMS-2052926, DMS-2053832, and DMS-2210833

###### Abstract

Uncertainty quantification for prediction is an intriguing problem with
significant applications in various fields, such as biomedical science,
economic studies, and weather forecasts. Numerous methods are available for
constructing prediction intervals, such as quantile regression and conformal
predictions, among others. Nevertheless, model misspecification (especially in
high-dimension) or sub-optimal constructions can frequently result in biased
or unnecessarily-wide prediction intervals. In this paper, we propose a novel
and widely applicable technique for aggregating multiple prediction intervals
to minimize the average width of the prediction band along with coverage
guarantee, called Universally Trainable Optimal Predictive Intervals
Aggregation (UTOPIA). The method also allows us to directly construct
predictive bands based on elementary basis functions. Our approach is based on
linear or convex programming which is easy to implement. All of our proposed
methodologies are supported by theoretical guarantees on the coverage
probability and optimal average length, which are detailed in this paper. The
effectiveness of our approach is convincingly demonstrated by applying it to
synthetic data and two real datasets on finance and macroeconomics.

Keywords: average predictive widths, predictive interval constructions,
coverage probability, basis expansions, kernel tricks, neural networks.

##  1 Introduction

Uncertainty quantification plays a crucial role in a wide range of scenarios
where decision-makers require not only accurate point prediction but also the
uncertainty quantification associated with the prediction. For instance, in
finance, a stock price prediction interval can help investors make informed
decisions and manage market risks, while in health care, a disease progression
prediction interval is important for making treatment decisions and improving
patient outcomes. Furthermore, reliable prediction intervals with provable
coverage guarantees that quantify uncertainty in machine learning-based
predictive models are now in high demand. In addition to providing a high-
level coverage guarantee, a prediction interval with a small width is also
important because it conveys more precise information about the uncertainties
in prediction, enabling better decision-making.

Several attempts have been made to construct prediction intervals owing to the
importance of uncertainty quantification. One of the earliest and perhaps most
popular methods is quantile regression. Given $\tau>0$, quantile regression
estimates $Q_{\tau}(Y\mid X)$, i.e. the conditional quantile of $Y$ given $X$
at quantile $\tau$, by minimizing a check loss. There is abundance of research
on the estimation of the quantile function, i.e. Koenker and Bassett Jr,
(1978); Koenker and Hallock, (2001); Belloni et al., (2019) to name a few. As
most of statistical modeling, model mis-specifications are of primary concern,
which can easily lead to invalid predictive intervals. These concerns amplify
with high dimensionality. In addition, quantile regression often suffers from
quantile crossing problem (Chernozhukov et al.,, 2010). Though the quantile
crossing problem can be alleviated by estimating several quantile functions
together along with a monotonicity constraint (e.g. see Fakoor et al., (2021);
Moon et al., (2021) and references therein), proper theoretical justification
along with computational efficiency is still lacking.

To address the above concerns, a relatively new and popular approach to
constructing a prediction band that has gained significant attention is the
conformal prediction method. This method was originally proposed in Vovk et
al., (2005) and later analyzed in Vovk et al., (2009); Lei et al., (2013); Lei
and Wasserman, (2014) and many other papers. A key advantage of the conformal-
based method is that it does not require any strong assumption on the
structure of the underlying distribution: The only requirement is for the data
to be exchangeable, which is later modified to incorporate time series
data(Barber et al.,, 2022; Stankeviciute et al.,, 2021; Xu and Xie,, 2020) and
under distribution shift(Tibshirani et al.,, 2019; Hu and Lei,, 2020). In the
vanilla conformal prediction, one needs to re-estimate the mean function to
create the conformity score for each plausible value of the response variable
$y$, which is quite computationally intensive. To alleviate this issue, Lei et
al., (2013) used a split conformal method, where the first mean function is
estimated using a part of the data and then the estimated mean function is
used in the remaining part of the data to create a prediction band. However,
the probability statement in their coverage guarantee is taken over both the
randomness of the past data and the new observation, which is weaker than the
conditional coverage guarantee (as in our Theorems 3.5, 3.7 and 3.9).
Furthermore, the conformal prediction approach can be conservative, i.e., the
resulting prediction intervals may be wider than necessary, leading to
imprecise quantification of uncertainty. Very recently, Liang, (2021) proposes
a kernel-based method to construct prediction bands for uncertainty
quantification. This method is easy to implement but can be computationally
expensive for large sample sizes. In addition, it requires certain conditions
on the distribution and the kernel function (Assumptions [S2] and [S3] in
Liang, (2021)). Furthermore, the objective of this method is to minimize the
complexity of the underlying function class and it is not immediate how this
is related to the width of the prediction interval.

From the above discussion, it is clear that we have several methods available
in the literature for constructing prediction intervals, each of which has its
own pros and cons. This begs an immediate question, if we have several methods
available, can we efficiently aggregate them to create a prediction interval
that has both coverage guarantee and small width? To address this issue, we
propose a method, called UTOPIA, for aggregating prediction bands via a simple
linear or convex programming problem that directly aims to minimize the width
of the aggregated confidence interval along with coverage guarantee. Let us
illustrate here with a simple example (details can be found in Setup 1 of
Section 5). Suppose that $(X,\xi)\overset{i.i.d.}{\sim}\mathrm{Unif}(-1,1)$
and $Y=\xi\sqrt{1+25X^{4}}$. The following picture compares the prediction
performance of UTOPIA and split conformal prediction (SplitCF):

(a) UTOPIA

(b) SplitCF

Figure 1: The blue dots and light blue band represent the testing data and the
prediction band, respectively.

As evident from the above picture, UTOPIA provides a coverage guarantee with a
smaller average width. Going beyond the aggregation of predictive intervals,
our UTOPIA method can be directly used to construct predictions band from a
general function class $\mathcal{F}$ as well as estimate $\mathbb{E}[Y\mid
X]$. One distinguished feature of our method is that we will have the coverage
probability guarantees, whether the function class $\mathcal{F}$ is
misspecified or not. In addition, we attempt to keep the average widths of
bands as small as we can via optimization through linear programming. To
summarize, in this work, we propose an approach that:

  1. 1.

Learns a prediction band with a small average width and achieves a provable
coverage guarantee simultaneously.

  2. 2.

Is easy to implement, only solving a linear/quadratic optimization problem.

  3. 3.

Is demonstrated convincingly through various applications both in synthetic
and real data.

The rest of the paper is organized as follows: Section 2 presents a general
methodology for constructing the prediction intervals. Section 3 provides
theoretical guarantees of our method. Section 4 elaborates on the theoretical
properties of the various important special cases of our general method. In
Section 5 we provide extensive simulations to support our theoretical findings
and to compare and contrast with the existing methods. Last, but not the
least, we further demonstrate the effectiveness of our approach by applying it
to two real datasets (Fama-French data and FRED-MD data) in Section 6. Details
of the proofs are deferred to the Appendix in the interest of space.

##  2 UTOPIA: Methods of Aggregations and Constructions of Predictive
Intervals

We start by assuming that we have a pair of random variables $(X,Y)\sim P$,
where $X\in\mathcal{X}$ and $Y\in\mathcal{Y}$. In this paper we typically
assume $\mathcal{X}={\mathbb{R}}^{p}$ and $\mathcal{Y}={\mathbb{R}}$. For any
such pair, we have the following decomposition:

$Y=m_{0}(X)+\sqrt{v_{0}(X)}\xi\,,$

where $m_{0}(x)=\mathbb{E}[Y\mid X=x]$ is the conditional mean function and
$v_{0}(x)={\sf var}(Y\mid X=x)$ is the conditional variance function. By
definition, we have $\mathbb{E}[\xi\mid X]=0$ and ${\sf var}(\xi\mid X)=1$,
albeit $\xi$ needs not be independent of $X$. A special case is when
$\xi\perp\\!\\!\\!\perp X$, which is known as the _location-scale_ family. As
mentioned in the introduction, our goal is to construct a prediction interval
of the form $[l(X),u(X)]$ of $Y$ given $X$ based on $n$ i.i.d random
observations $(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$ from $P$. Here we aim for
_expected guarantee_ instead of _conditional guarantee_ , i.e.
$\mathbb{P}(\ell(X)\leq Y\leq u(X))\geq 1-\alpha$ where the probability is
taken over the joint distribution of $(X,Y)$. Constructing any such prediction
interval is easy; given any positive function $f$, choose $c\equiv c(f)$ such
that:

$\mathbb{P}(-cf(X)\leq Y\leq cf(X))\geq 1-\alpha\,.$

Such a $c$ always exists when the distribution of $(X,Y)$ is continuous since
the probability on the left-hand side goes to $0$ as $c$ goes to $0$, goes to
$1$ as $c$ goes to $\infty$. In other words, any positive function can be
converted to a prediction interval with $(1-\alpha)$ coverage guarantee just
by a simple re-scaling.

###  2.1 Minimum Average Length Predictive Intervals

Another desirable property of a prediction interval is the small width. For
prediction intervals with average guarantee, this leads to the following
optimization problem:

$\displaystyle\min_{l,u}$ $\displaystyle\mathbb{E}[u(X)-l(X)]$ (2.1)
$\displaystyle{\sf subject\ to}$ $\displaystyle\mathbb{P}(Y\in[l(X),u(X)])\geq
1-\alpha\,.$

The solution to the above optimization problem is the _optimal_ prediction
interval in the sense that it has the smallest average width among all the
prediction intervals with $(1-\alpha)$ average guarantee. However, the above
optimization problem is hard to solve without further restriction for a
general $\alpha>0$. To understand why, recall that it is not possible to find
the optimal $100(1-\alpha)\%$ confidence interval of a general univariate
distribution $F$ unless we assume symmetry and unimodality. However, if
$\alpha=0$, i.e. we want $100\%$ confidence interval, then the optimal range
is $[F^{-1}(0),F^{-1}(1)]$ no matter what the distribution $F$ is. This
property can also be extrapolated to the general problem (2.1), which we
present in the following lemma:

######  Lemma 2.1.

For $\alpha=0$, the optimal solution of (2.1) is

$\displaystyle l^{\star}(x)$ $\displaystyle=F_{Y\mid
X=x}^{-1}(0)=m_{0}(x)+\sqrt{v_{0}(x)}F_{\xi\mid X=x}^{-1}(0)\,,$
$\displaystyle u^{\star}(x)$ $\displaystyle=F_{Y\mid
X=x}^{-1}(1)=m_{0}(x)+\sqrt{v_{0}(x)}F_{\xi\mid X=x}^{-1}(1)\,.$

Furthermore, when the conditional distribution of $\xi$ given $X$ is symmetric
around $0$, then $F_{\xi\mid X=x}^{-1}(0)=-F_{\xi\mid X=x}^{-1}(1)$ and
consequently the optimal prediction interval for $Y$ is:

$\mathrm{PI}^{\star}_{100\%}(x)=\left[m_{0}(x)-\sqrt{v_{0}(x)}F_{\xi\mid
X=x}^{-1}(1),\ m_{0}(x)+\sqrt{v_{0}(x)}F_{\xi\mid X=x}^{-1}(1)\right]\,.$

with optimal length being $2\sqrt{v_{0}(x)}F_{\xi\mid X=x}^{-1}(1)$.

Although, the solution obtained in Lemma 2.1 is non-trivial only if $Y$ is
compactly supported (i.e., $F_{\xi\mid X}^{-1}(0)>-\infty$ and $F_{\xi\mid
X}^{-1}(1)<\infty$), one can restrict themselves in a suitable compact for
unbounded responses. For simplicity of exposition, we work with a compactly
supported response. The above observation leads to the following conclusion:
solving (2.1) for $\alpha=0$ is easier than $\alpha>0$ as long as $(X,Y)$ is
compactly supported. Therefore, given $n$ i.i.d. observations
$(X_{1},Y_{1}),\dots,(X_{n},Y_{n})$, we begin with empirical $100\%$ coverage
guarantee at the observed data; see constraint in (2.2) under the following
assumption:

######  Assumption 2.2.

The random variable $(X,Y)$ is compactly supported. The distribution of $Y$
given $X$ is symmetric around $m_{0}(X)$, i.e. the conditional distribution of
$\xi$ given $X$ is symmetric around the origin almost surely on $X$.

The compact assumption is quite standard in the literature of non-parametric
analysis and it is primarily assumed to simplify some technicalities without
hurting the main message. The assumption of the symmetric distribution of
$\xi$ given $X$ is to ensure that we only need to estimate $F_{\xi\mid
X}^{-1}(1)$ instead of both $F_{\xi\mid X}^{-1}(0)$ and $F_{\xi\mid
X}^{-1}(1)$ as $F_{\xi\mid X}^{-1}(0)=-F_{\xi\mid X}^{-1}(1)$. However,
symmetry assumption is not necessary; our general results continue to hold
without symmetry assumption, see Remark 2.3 for more details.

As evident from Lemma 2.1, we need to estimate $m_{0}$ (the conditional mean
function) and $v_{0}(x)(F_{\xi\mid X}^{-1}(1))^{2}$ to create a prediction
interval under Assumption 2.2. For notational simplicity, henceforth we define
$f_{0}(x)=v_{0}(x)(F_{\xi\mid X}^{-1}(1))^{2}$.

There are several ways to estimate each of them as $m_{0}$ is the conditional
mean of $Y$ given $X$ and $f_{0}(x)$ is the essential supremum of
$(Y-m_{0}(X))^{2}$ given $X$. Therefore, one may apply various methods (both
parametric and non-parametric) to estimate each of the functions, e.g. kernel-
based methods (Nadaraya-Watson estimator or local polynomial regression),
basis expansion-based methods using spline basis or wavelet basis (Fan and
Gijbels,, 1996) or employing neural networks for $m_{0}$ and several quantile
regression related methods for estimating $f_{0}$. In fact for location-scale
model where $\xi\perp\\!\\!\\!\perp X$, then $F_{\xi\mid
X}^{-1}=F_{\xi}^{-1}(1)$ does not depend on $X$. In that case, one can simply
estimate $v_{0}$, which is the conditional mean of $(Y-m_{0}(X))^{2}$ and then
adjust the constant $F_{\xi\mid X}^{-1}(1)$. However, we also need to make
sure that the average width of the prediction interval is as small as possible
and the concerns of model mis-specification for both $m_{0}(\cdot)$ and
$f_{0}(\cdot)$ can not be understated.

###  2.2 Optimal predictive intervals aggregation and construction

Given a collection of candidate functions $\mathcal{F}$ (resp. $\mathcal{M}$)
for $f_{0}$ (resp. $m_{0}$) and $n$ i.i.d. observations
$\\{(X_{1},Y_{1}),\dots,(X_{n},Y_{n})\\}$ we propose to solve the following
optimization problem:

$\displaystyle\min_{m\in\mathcal{M},f\in\mathcal{F}}$ $\displaystyle\ \ \ \
\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$ (2.2) $\displaystyle{\sf subject\ to}$
$\displaystyle\ \ \ \ f(x_{i})\geq(y_{i}-m(x_{i}))^{2}\ \ \ \forall\ 1\leq
i\leq n\,.$

The above optimization problem has the following implication: we want to
choose estimates $(\hat{m},\hat{f})$ of $(m_{0},f_{0})$ such that i)
$Y_{i}\in[\hat{m}(X_{i})-\sqrt{\hat{f}(X_{i})},\hat{m}(X_{i})+\sqrt{\hat{f}(X_{i})}]$
for all the observed samples (i.e. $100\%$ empirical coverage) and ii)
$\hat{f}$ has minimum empirical average bandwidth among all such estimates.
Given that we’re optimizing on both $(m,f)$, one may wonder why we do not have
any explicit loss function for estimating $m_{0}$ in the objective of (2.2)
and how does the estimation of $m_{0}$ affect the width of the estimated
prediction interval. We elaborate on this in Lemma 2.4.

The general formation (2.2) suggests a universal approach for constructing
high-coverage confidence bands with small width by choosing $\mathcal{M}$ and
$\mathcal{F}$ as desired, such as linear combination of some basis functions
(model aggregation, see Section 2.3) and Reproducing Kernel Hilbert Space (see
Section 2.4). Upon obtaining $(\hat{m},\hat{f})$ by solving (2.2), our
estimated prediction interval is:

$\displaystyle\widehat{\mathrm{PI}}_{100\%}(x)=\left[\hat{m}(x)-\sqrt{\hat{f}(x)},\
\hat{m}(x)+\sqrt{\hat{f}(x)}\right].$ (2.3)

If the underlying function is complex (e.g. non-VC class, precisely defined
later), then often sample level guarantee is not sufficient for the population
level guarantee. In these cases, we recommend slightly extending the
prediction interval:

$\displaystyle\widehat{\mathrm{PI}}_{100\%}(x,\delta)=\left[\hat{m}(x)-\sqrt{\hat{f}(x)+\delta},\
\hat{m}(x)+\sqrt{\hat{f}(x)+\delta}\right],$ (2.4)

for some small $\delta>0$ that depends on $n$ and typically decreases to $0$
as $n\to\infty$. For more details, see Remark 3.6.

Although so far we have studied the case $\alpha=0$, in most of the practical
scenarios, we only need to cover the variable $Y$ with a probability of at
least $1-\alpha$ for some $\alpha>0$. In that case, the above-obtained
interval could be much wider than needed. One simple way out is to _scale
down_ $\hat{f}$ (i.e. multiplying $\hat{f}$ by some $\lambda\in(0,1)$ so that
$\lambda\hat{f}$ has desired empirically coverage on a testing sample), as
summarized in Algorithm 1. The key steps are i) first estimate
$\hat{f},\hat{m}$ to guarantee $100\%$ empirical coverage and then ii) scale
it down appropriately for desired $100(1-\alpha)\%$ coverage.

Input : $\alpha$, sample $\mathcal{D}=\\{(X_{i},Y_{i})\\}_{i=1}^{n}$, function
classes $\mathcal{M},\mathcal{F}$, and $\delta$, where $\delta=0$ when
$(\mathcal{F},\mathcal{M})$ are VC class of functions.

Output : $\hat{m},\hat{f},\hat{\lambda}_{\alpha}$

1 Use the data to solve (2.2) and obtain $\hat{m},\hat{f}$.

2 Set $\hat{\lambda}_{\alpha}$ as the $(1-\alpha)^{th}$ empirical quantile of
$(Y_{i}-\hat{m}(X_{i}))^{2}/(\hat{f}(X_{i})+\delta)$, i.e.
$\displaystyle\hat{\lambda}_{\alpha}:=\inf_{\lambda}\left\\{\lambda:\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\\{(Y_{i}-\hat{m}(X_{i}))^{2}>\lambda(\hat{f}(X_{i})+\delta)\\}}\leq\alpha\right\\}\,.$

Return $(\hat{\lambda}_{\alpha},\hat{m},\hat{f})$ with the estimated
prediction interval:
$\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(x,\delta)=\left[\hat{m}(x)-\sqrt{\hat{\lambda}_{\alpha}(\hat{f}(x)+\delta)},\
\hat{m}(x)+\sqrt{\hat{\lambda}_{\alpha}(\hat{f}(x)+\delta)}\right]\,.$

[-0.15in] Aggregation with coverage probability $1-\alpha$

Algorithm 1 UTOPIA: Universally Trainable Optimal Prediction Intervals

[-0.15in] Aggregation with coverage probability $1-\alpha$

It is tempting to think that the scale-down procedure is sub-optimal, however,
our method has two significant advantages applicable to many real-life
scenarios:

  1. 1.

Our proposal is universally applicable, independent of model mis-
specification. It is easy to implement, i.e. once we obtain a good $100\%$
prediction interval, all we need to do is to scale it down properly.
Furthermore, we have theoretical guarantees for the coverage of our estimated
prediction interval (See Section 3.4), regardless of the correct or incorrect
specifications of the model.

  2. 2.

Our procedure yields the optimal quantile function (at $(1-\alpha)$ level) for
the location-scale family, i.e. when $\xi$ independent of $X$. This follows
from the observation that for the location-scale family:

$Q_{1-\alpha}((Y-m_{0}(x))^{2}\mid X=x)=\lambda^{\star}(\alpha)f_{0}(x)$

where
$\lambda^{\star}(\alpha)=F_{\xi^{2}}^{-1}(1-\alpha)/F_{\xi^{2}}^{-1}(1)$.

When $\mathcal{M}$ is not a singleton, then we need to estimate both $f_{0}$
and $m_{0}$ simultaneously and consequently, it may be computationally
challenging. To alleviate this, one may follow a two-step procedure, which
bears the essence of split-conformal prediction band (Lei et al.,, 2013): i)
First estimate $m_{0}$ from $\mathcal{M}$ via parametric or non-parametric
regression and then ii) plugin $\hat{m}$ in (2.2) to estimate $\hat{f}$. We
summarize this procedure in the Algorithm 2. Algorithm 2 primarily differs
from Algorithm 1 in Step 2, where we first estimate $m_{0}$ and then $f_{0}$
separately. This method has a computational advantage as regression is easy to
implement (e.g. using the kernel-based methods, spline-based approaches, or
via deep neural networks).

Input : $\alpha$, sample $\mathcal{D}=\\{(X_{i},Y_{i})\\}_{i=1}^{n}$, function
classes $\mathcal{M},\mathcal{F}$, and $\delta$, where $\delta=0$ when
$(\mathcal{F},\mathcal{M})$ are VC class of functions

Output : $\hat{m},\hat{f},\hat{\lambda}_{\alpha}$

1 Use the data to first estimate $m_{0}$ via non-parametric regression, i.e.
$\hat{m}={\arg\min}_{m\in\mathcal{M}}\left[\sum_{i\in\mathcal{D}}(Y_{i}-m(X_{i}))^{2}+\rho_{n}(m)\right]\,,$
where $\rho_{n}(\cdot)$ is a penalty function.

2 Use the data to solve (2.2) with $\hat{m}$ fixed and obtain $\hat{f}$.

3 Compute the $(1-\alpha)^{th}$ empirical quantile of
$(Y_{i}-\hat{m}(X_{i}))^{2}/(\hat{f}(X_{i})+\delta)$:
$\displaystyle\hat{\lambda}_{\alpha}:=\inf_{\lambda}\left\\{\lambda:\frac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{\\{(Y_{i}-\hat{m}(X_{i}))^{2}>\lambda(\hat{f}(X_{i})+\delta)\\}}\leq\alpha\right\\}\,.$

Return $(\hat{\lambda}_{\alpha},\hat{m},\hat{f})$ with the estimated
prediction interval:
$\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(x,\delta)=\left[\hat{m}(x)-\sqrt{\hat{\lambda}_{\alpha}(\hat{f}(x)+\delta)},\
\hat{m}(x)+\sqrt{\hat{\lambda}_{\alpha}(\hat{f}(x)+\delta)}\right]\,.$

Algorithm 2 Two-step UTOPIA

######  Remark 2.3 (Extension to non-symmetric distribution).

Although we have assumed the distribution of $Y\mid X$ to be symmetric around
its mean $m_{0}$, our method can also be extended to asymmetric distribution.
In that case, if we are only interested in estimating a prediction interval of
$Y$ given $X$, not its mean, then we can modify (2.2) as follows:

$\displaystyle\min_{f_{1}\in\mathcal{F}_{1},f_{2}\in\mathcal{F}_{2}}$
$\displaystyle\ \ \frac{1}{n}\sum_{i}(f_{2}(X_{i})-f_{1}(X_{i}))$ (2.5)
$\displaystyle{\sf subject\ to}$ $\displaystyle\ \ f_{1}(X_{i})\leq Y_{i}\leq
f_{2}(X_{i}),\ \ \forall\ 1\leq i\leq n\,.$

Our theories can be extended easily to incorporate this situation. If $m_{0}$
estimation is necessary, then we can resort to a one-step or two-step method
by replacing $Y_{i}$ by $Y_{i}-m(X_{i}))$ and optimizing over
$(m,f_{1},f_{2})$.

We now address the question that how the estimation of $m_{0}$ affects the
width of the obtained prediction interval through the following lemma.

######  Lemma 2.4.

Let $m(x)$ be an arbitrary mean estimator and $\mathcal{F}$ is the set of all
non-negative functions. Then, the following problem

$\displaystyle\min_{f\in\mathcal{F}}$ $\displaystyle\ \ \ \ \mathbb{E}[f(X)]$
(2.6) $\displaystyle{\sf subject\ to}$ $\displaystyle\ \ \ \
\mathbb{P}\left(f(X)\geq(Y-m(X))^{2}\right)=1.$

has a unique explicit solution

$\displaystyle f_{m}(x)=\left(\sqrt{f_{0}(x)}+|m_{0}(x)-m(x)|\right)^{2}.$
(2.7)

Lemma 2.4 ascertains that $f_{m}$ is strictly larger that $f_{0}(x)$ for all
$x$ where $m_{0}\neq m$ and consequently
$\mathbb{E}[f_{m}(X)]\geq\mathbb{E}[f_{0}(X)]+\|m-m_{0}\|_{2}^{2}$. Thus, a
_bad replacement_ of $m_{0}$ (in the sense $\|m-m_{0}\|$ is large) affects the
width of the optimal prediction interval (though coverage probability is
correct), stipulating the necessity of estimating $m_{0}$ well. Since
estimating $m_{0}$ is important, one may wonder why we do not include any term
corresponding to the estimating error of $m$ into the objective function in
(2.2). The basic intuition is that, we here indirectly minimize
$\sum_{i}(Y_{i}-m(X_{i}))^{2}$ by minimizing an upper bound
$\sum_{i}f(X_{i})$. In other words, if $\sum_{i}f(X_{i})$ is small and
$f(X_{i})\geq(Y_{i}-m(X_{i}))^{2}$, then we automatically ensure
$\sum_{i}(Y_{i}-m(X_{i}))^{2}$ to be small.

###  2.3 Aggregation of Confidence Bands

As elaborated previously, one may use several different techniques for
constructing prediction intervals, i.e. various parametric or non-parametric
methods for estimating $(m_{0},f_{0})$. Furthermore, for a location-scale
family, it is sufficient to estimate $(m_{0},v_{0})$. Additionally, modern
machine learning techniques like neural networks can be employed to create
prediction intervals by utilizing their ability to capture complex patterns in
data. A special case of (2.2) is model aggregation which aims to aggregate
several candidates for prediction interval to guarantee adequate coverage as
well as small width.

Suppose we have $\mathcal{M}_{0}=\\{m_{1},\dots,m_{L}\\}$ and
$\mathcal{F}_{0}=\\{f_{1},\dots,f_{K}\\}$, i.e. we have some candidate
estimates for $(m_{0},f_{0})$, and we take their linear spans for
$\mathcal{M}$ and $\mathcal{F}$ respectively as linear aggregation. Then (2.2)
reduces to the following quadratic optimization problem:

$\displaystyle\min_{\begin{subarray}{c}\alpha_{1},\dots,\alpha_{K}\\\
\beta_{1},\dots,\beta_{L}\end{subarray}}$ $\displaystyle\ \ \ \
\frac{1}{n}\sum^{n}_{i=1}\Bigl{\\{}\sum_{j=1}^{K}\alpha_{j}f_{j}(x_{i})\Bigr{\\}}$
(2.8) $\displaystyle{\sf subject\ to}$ $\displaystyle\ \ \ \
\sum_{j=1}^{K}\alpha_{j}f_{j}(x_{i})\geq\Bigl{(}y_{i}-\Bigl{\\{}\sum_{j=1}^{L}\beta_{j}m_{j}(x_{i})\Bigr{\\}}\Bigr{)}^{2}\
\ \ \forall\ 1\leq i\leq n\,,$ $\displaystyle\ \ \ \ \alpha_{j}\geq 0,\ \ \
\forall\ 1\leq j\leq K\,.$

The first curly bracket indicates the “squared width” at point $X_{i}$ and the
optimization object is the average “squared width”, subject to constraint of
empirical 100% coverage. The second curly bracket is adjusting the center of
the predictive intervals to make the width as small as possible. In other
words, (2.8) aims to find a suitable linear combination of the available
estimators that covers the observed responses as well as has the minimum width
of the prediction interval. In matrix notation, if we define
$\mathbf{M}_{n}\in{\mathbb{R}}^{n\times L}$ and
$\mathbf{F}_{n}\in{\mathbb{R}}^{n\times K}$ as
$\mathbf{M}_{n,ij}=m_{j}(x_{i})$ and $\mathbf{F}_{n,ij}=f_{j}(x_{i})$, Then we
can rewrite (2.8) as

$\displaystyle\min_{{\bm{\alpha}}\in{\mathbb{R}}^{K}_{+},{\bm{\beta}}\in{\mathbb{R}}^{L}}$
$\displaystyle\ \ \ \mathbf{1}^{\top}\mathbf{F}_{n}{\bm{\alpha}}$
$\displaystyle{\sf subject\ to}$ $\displaystyle\ \ \ \
e_{i}^{\top}\left(\mathbf{Y}-\mathbf{M}_{n}{\bm{\beta}}\right)\left(\mathbf{Y}-\mathbf{M}_{n}{\bm{\beta}}\right)^{\top}e_{i}-e_{i}^{\top}\mathbf{F}_{n}{\bm{\alpha}}\leq
0\ \ \ \forall\ 1\leq i\leq n\,.$

As the constraint function is quadratic in ${\bm{\beta}}$ and linear in
${\bm{\alpha}}$, it is jointly convex in $({\bm{\alpha}},{\bm{\beta}})$ and
consequently can be solved efficiently via standard convex optimization
techniques. Let $\widehat{{\bm{\alpha}}}$ and $\widehat{{\bm{\beta}}}$ be the
solution to the above optimization problem. Then the proposed prediction
interval is

$\displaystyle\widehat{\mathrm{PI}}_{100\%}(x)$
$\displaystyle=\left[\sum_{j=1}^{L}\hat{\beta}_{j}m_{j}(x)-\sqrt{\sum_{j=1}^{K}\hat{\alpha}_{j}f_{j}(x)},\
\
\sum_{j=1}^{L}\hat{\beta}_{j}m_{j}(x)+\sqrt{\sum_{j=1}^{K}\hat{\alpha}_{j}f_{j}(x)}\right]\,.$
(2.9)

As the linear span of finitely many functions is a VC class, we do not need to
extend the interval by $\delta$ here, namely taking $\delta=0$. For
$100(1-\alpha)\%$ coverage with some $\alpha>0$, we can scale down
$\hat{f}_{0}$ as per Step 2 of Algorithm 1. For theoretical guarantees see
Section 4.1.

###  2.4 Construction of Prediction Bands by Kernel Basis

In this subsection we estimate $(m_{0},f_{0})$ via kernel basis assuming they
lie on a reproducing kernel Hilbert space. For simplicity, we assume $m_{0}$
is known and aim to approximate $f_{0}$ by some non-negative function over a
RKHS $\mathcal{H}$ with kernel $K$ and feature map $\Phi$, where typically we
take $\Phi(x)=K(x,\cdot)$. As $f_{0}$ is non-negative, it is tempted to
approximate $f_{0}$ by
$\mathcal{F}=\\{\sum_{j}\alpha_{j}K(x_{j},\cdot):\alpha_{j}\geq 0\\}$.
However, as elaborated in Marteau-Ferey et al., (2020), they may not be
universal approximators in the sense that they may not approximate all non-
negative functions and suggested to use the quadratic form
$\\{\langle\Phi(x),\mathcal{A}(\Phi)(x)\rangle_{\mathcal{H}}\\}$ where
$\mathcal{A}$ is bounded linear Hermitian operator from
$\mathcal{H}\mapsto\mathcal{H}$ and $\langle\cdot,\cdot\rangle_{\mathcal{H}}$
is the inner product associate with the Hilbert space $\mathcal{H}$.

Liang, (2021) proposed the following method to estimate $f_{0}$:

$\displaystyle\min_{\mathcal{A}:\mathcal{H}\rightarrow\mathcal{H}}$
$\displaystyle\|\mathcal{A}\|_{\star}$ (2.10) $\displaystyle{\sf subject\ to}$
$\displaystyle\langle\Phi(x_{i}),\mathcal{A}(\Phi)(x_{i})\rangle_{\mathcal{H}}\geq(y_{i}-m_{0}(x_{i}))^{2}$
$\displaystyle\mathcal{A}\succeq 0\,,$

where $\|\cdot\|_{\star}$ denotes the nuclear norm of an operator. The aim of
this objective is to minimize the rank of $\mathcal{A}$. However as our
objective is to minimize the expected empirical width, we start with a simpler
function class
$\mathcal{F}=\\{f(x)=\langle\Phi(x),\mathcal{A}(\Phi)(x)\rangle_{\mathcal{H}},\|\mathcal{A}\|_{\star}\leq
r,\mathcal{A}\succeq 0\\}$ for some given rank $r$ and consequently (2.2)
boils down to the following:

$\displaystyle\min_{\mathcal{A}:\mathcal{H}\rightarrow\mathcal{H}}$
$\displaystyle\frac{1}{n}\sum_{i}\langle\Phi(x_{i}),\mathcal{A}(\Phi)(x_{i})\rangle_{\mathcal{H}}$
(2.11) $\displaystyle{\sf subject\ to}$
$\displaystyle\langle\Phi(x_{i}),\mathcal{A}(\Phi)(x_{i})\rangle_{\mathcal{H}}\geq(y_{i}-m_{0}(x_{i}))^{2}$
$\displaystyle\|\mathcal{A}\|_{\star}\leq r,\mathcal{A}\succeq 0\,.$

By representer theorem in Marteau-Ferey et al., (2020), the solution to the
above infinite-dimensional problem can be written as
$\mathcal{A}=\sum^{n}_{i,j=1}B_{ij}\Phi(x_{i})\otimes\Phi(x_{j})$ for some
positive semi-definite matrix $B$, which implies an equivalent finite-
dimensional problem:

$\displaystyle\min_{B\in\mathbb{S}^{n\times n}}$ $\displaystyle{\sf tr}(KBK)$
(2.12) $\displaystyle{\sf subject\ to}$ $\displaystyle\left\langle
K_{i},BK_{i}\right\rangle\geq(y_{i}-m_{0}(x_{i}))^{2}\ \forall\ 1\leq i\leq n$
$\displaystyle{\sf tr}(KB)\leq r,\ \ B\succeq 0\,.$

Here $K\in\mathbb{S}^{n\times n}$ is the corresponding kernel matrix with
$K_{ij}=K(x_{i},x_{j})$ and $K_{i}\in\mathbb{R}^{n}$ represents the $i$-th
column of the matrix $K$. Notice that (2.12) is semi-definite programming as
we have the cone constraint ($B\succeq 0$) and both the objectives and
constraints are linear in $B$. Therefore, it can be easily implemented via
standard convex programming package (e.g. cvx in MATLAB). Let $\hat{B}$ be the
solution to the above optimization problem. Then the estimated prediction
interval is:

$\displaystyle\widehat{\mathrm{PI}}_{100\%}(x,\delta)$
$\displaystyle=\left[m_{0}(x)-\sqrt{\langle
K_{x},\hat{B}K_{x}\rangle+\delta},\ \ m_{0}(x)+\sqrt{\langle
K_{x},\hat{B}K_{x}\rangle+\delta}\right]\,,$

where $K_{x}\in\mathbb{R}^{n}$ with $K_{x,i}=K(x_{i},x)$. For
$100(1-\alpha)\%$ coverage with some $\alpha>0$, we can scale down $\langle
K_{x},\hat{B}K_{x}\rangle+\delta$ as per Step 2 of Algorithm 1. Furthermore,
when $m_{0}$ is unknown, one can replace $m_{0}$ either by some given estimate
$\hat{m}$ or simultaneously estimate it via (2.2) and a suitable function
class $\mathcal{M}$ (e.g., one may also take $\mathcal{M}$ to be a RKHS). For
more details and the theoretical guarantees see Section 4.3.

##  3 Statistical Properties

In this section, we present our theoretical analysis of UTOPIA. Recall that we
obtain our prediction interval (2.3) by solving (2.2) where
$\hat{f}\in\mathcal{F}$ and $\hat{m}\in\mathcal{M}$. Therefore the convergence
properties of this estimated prediction interval rely on the complexity of the
underlying function classes, especially a trade-off between the number of
available samples and the complexity of the function classes. Before delving
deep, first, we present notations and preliminaries on the complexity of the
underlying function classes that are going to be used throughout this section.

###  3.1 Notations and preliminaries

##### Notations

We use $m_{0}(x)$ and $v_{0}(x)$ to denote the conditional mean function and
conditional variance function of $Y\mid X=x$, respectively. Given any guess
$m$ of $m_{0}$, we denote by $f_{m}(x)$ the supremum of $(Y-m(X))^{2}\mid
X=x$, i.e.

$f_{m}(x):=\inf_{t}\left\\{t\,\big{|}\,\mathbb{P}\left((Y-m(X))^{2}\leq
t\big{|}X=x\right)=1\right\\}$

and $f_{0}(x):=f_{m_{0}}(x)$. We use $\mathcal{D}$ to denote $n$ observed
samples. We assume there exists $\mathcal{B}_{\mathcal{M}}>0$ and
$\mathcal{B}_{\mathcal{F}}>0$ such that $\|m\|_{\infty}\leq B_{\mathcal{M}}$
for all $m\in\mathcal{M}$ and $\|f\|_{\infty}\leq B_{\mathcal{F}}$ for all
$f\in\mathcal{F}$.

##### Complexity of function class

There are various notions available in the literature to quantify the
complexity such as VC-dimension (Vapnik and Chervonenkis,, 2015), metric
entropy (Kolmogorov and Tikhomirov,, 1959), pseudo-dimension (Pollard,, 1990),
fat-shattering dimension (Kearns and Schapire,, 1994) etc. To cover a wide
array of scenarios, we here quantify the complexity of $\mathcal{F}$ through
its Rademacher complexity, defined as follows.

######  Definition 3.1.

Let $\mathcal{F}$ be a function class and $\\{X_{i}\\}^{n}_{i=1}$ be a set of
samples drawn i.i.d. from a distribution $\mathcal{D}$. The Rademacher
complexity of $\mathcal{F}$ is defined as

$\displaystyle\mathcal{R}_{n}(\mathcal{F})=\mathbb{E}_{{\epsilon},\mathcal{D}^{\otimes
n}}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum^{n}_{i=1}{\epsilon}_{i}f(X_{i})\right],\
\ \
{\epsilon}_{1},\dots,{\epsilon}_{n}\overset{\mathrm{i.i.d}}{\sim}\mathrm{Rad}(1/2)\,.$
(3.1)

Rademacher complexity of a function class $\mathcal{F}$ describes how closely
a function class can be spuriously correlated with independent Rademacher
random variables. Typically, a bound on this complexity can be obtained via
covering a number of $\mathcal{F}$ via Dudley’s chaining argument, which we
state here for the convenience of the readers.

######  Definition 3.2.

Given a function class $\mathcal{F}$ with some (semi)-metric $d$ on
$\mathcal{F}$, the ${\epsilon}$-covering number of $\mathcal{F}$ is defined
as:

$\mathcal{N}({\epsilon},d,\mathcal{F})=\inf\left\\{m\geq 0:\exists
f_{1},\dots,f_{m}\in\mathcal{F}\ \ \mathrm{such\ that}\ \
\sup_{f\in\mathcal{F}}\min_{1\leq j\leq
m}d(f_{j},f)\leq{\epsilon}\right\\}\,.$

The ${\epsilon}$-metric entropy of $\mathcal{F}$ is defined as
$\log{\mathcal{N}({\epsilon},d,\mathcal{F})}$.

Clearly $\mathcal{N}({\epsilon},d,\mathcal{F})$ is a non-increasing function
of ${\epsilon}$ and as ${\epsilon}\downarrow 0$,
$\mathcal{N}({\epsilon},d,\mathcal{F})\uparrow\infty$. Depending on the growth
rate of $N({\epsilon},\mathcal{F},d)$ near origin, we have two broad
categories:

  1. 1.

Polynomial growth:  When $\mathcal{N}({\epsilon},d,\mathcal{F})$ grows
polynomially in ${\epsilon}^{-1}$ as ${\epsilon}\downarrow 0$, i.e.

$\mathcal{N}({\epsilon},d,\mathcal{F})\lesssim\left(\frac{1}{{\epsilon}}\right)^{\alpha}\,.$
(3.2)

  2. 2.

Exponential growth: When $\mathcal{N}({\epsilon},d,\mathcal{F})$ grows
exponentially in ${\epsilon}^{-1}$ as ${\epsilon}\downarrow 0$, i.e.

$\log{\mathcal{N}({\epsilon},d,\mathcal{F})}\lesssim\left(\frac{1}{{\epsilon}}\right)^{\alpha}\,.$
(3.3)

As it is evident a function class with polynomial growth is _less complex_
than a function class with exponential growth. Consequently, we expect a
smaller Rademacher complexity for the former class. The relation between the
Rademacher complexity and the covering number is captured through the
following proposition:

######  Proposition 3.3.

Let $\mathcal{F}$ be a function class and $B>0$. Suppose that
$\|f\|_{\infty}\leq B$ for any $f\in\mathcal{F}$. Consider the covering number
of $\mathcal{F}$ in terms of $L_{\infty}$ norm. Then we have the following
bound on the Rademacher complexity in terms of the growth of the covering
number:

$\mathcal{R}_{n}(\mathcal{F})\leq\begin{cases}C_{1}n^{-\frac{1}{2}},&\text{ if
the growth is polynomial or exponential with }0<\alpha\leq 2\,,\\\
C_{2}n^{-\frac{1}{\alpha}},&\text{ if the growth is exponential with
}\alpha>2\,.\end{cases}$

Here $C_{1},C_{2}$ depend only on $\alpha$ and $B$.

Thus, if the covering number of $\mathcal{F}$ (w.r.t. $L_{\infty}$ norm) grows
polynomially or exponentially with $0<\alpha<2$ then the Rademacher complexity
is of the order $n^{-1/2}$ (this is also known as Donsker class of functions).
The Rademacher complexity increases with $\alpha$ when the growth of the
covering number is exponential with $\alpha>2$. We will use this in our
subsequent theorems.

One special sub-class of functions with polynomial growth is Vapnik-
Chervonenkis (VC) class of functions. We present a rigorous definition below:

######  Definition 3.4 (Vapnik-Cervonenkis(VC) class).

Suppose $\mathcal{C}$ is a collection of sets. The collection $\mathcal{C}$ is
said to shatter $\\{x_{1},\ldots,x_{n}\\}$ if
$|\\{\\{x_{1},\ldots,x_{n}\\}\cap C\mid C\in\mathcal{C}\\}|=2^{n}$. The VC
dimension of $\mathcal{C}$ denoted by $\mathrm{VC}(\mathcal{C})$ is the
smallest $n$ for which no set of size $n$ is shattered by $\mathcal{C}$. If
$\mathrm{VC}(\mathcal{C})$ is finite, then $\mathcal{C}$ is called VC class. A
collection of measurable functions
$\mathcal{F}:\mathcal{X}\mapsto{\mathbb{R}}$ on a sample space $\mathcal{X}$
is called a VC class of functions if the collection of all subgraphs of the
functions in $\mathcal{F}$ forms a VC class of sets in
$\mathcal{X}\times\mathbb{R}$. And the VC dimension of $\mathcal{F}$ denoted
by $V_{\mathcal{F}}$ is defined as the VC dimension of the collection of all
subgraphs.

If a function class $\mathcal{F}$ is uniformly bounded by a constant $B$ and
has VC dimension $V_{\mathcal{F}}<\infty$, then we can bound its covering
number by Haussler’s bound (see Theorem 2.6.7 of Van Der Vaart and Wellner,
(1996)):

$\displaystyle\sup_{Q}N({\epsilon}F,\mathcal{F},L_{2}(Q))$ $\displaystyle\leq
CV_{\mathcal{F}}(16e)^{V_{\mathcal{F}}}\left(\frac{1}{{\epsilon}}\right)^{2V_{\mathcal{F}}-1}\,$
(3.4)

for some constants $C$. We will use this in our proofs.

###  3.2 Estimating $f_{0}$ when $m_{0}$ is given

In this subsection, we assume to know $m_{0}$ and consequently (2.2)
simplifies to:

$\displaystyle\min_{f\in\mathcal{F}}$ $\displaystyle\ \ \ \
\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$ (3.5) $\displaystyle{\sf subject\ to}$
$\displaystyle\ \ \ \ f(x_{i})\geq(y_{i}-m_{0}(x_{i}))^{2}\ \ \ \forall\ 1\leq
i\leq n\,,$

with the corresponding population version:

$\displaystyle\min_{f\in\mathcal{F}}$ $\displaystyle\ \ \ \ \mathbb{E}[f(X)]$
(3.6) $\displaystyle{\sf subject\ to}$ $\displaystyle\ \ \ \
\mathbb{P}\left(f(X)\geq(Y-m_{0}(X))^{2}\right)=1.$

If $f_{0}\in\mathcal{F}$, then by Lemma 2.1, it is the optimal solution of
(3.6). If $f_{0}\notin\mathcal{F}$, then the model mis-specification error,
denoted by $\Delta(\mathcal{F})$, is defined as follows:

$\Delta(\mathcal{F})=\min_{\begin{subarray}{c}f\in\mathcal{F}\\\ f\geq
f_{0}\end{subarray}}\mathbb{E}[f(X)-f_{0}(X)]\,.$ (3.7)

The following theorem presents a bound on the difference between the expected
width of the estimated prediction interval and its corresponding population
minimizer in terms of $\Delta(\mathcal{F})$ and the complexity of the
underlying function class.

######  Theorem 3.5.

Let $\hat{f}$ be the solution of (3.5). Then for any $t\geq 0$, with
probability at least $1-2e^{-t}$ the following bound holds:

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq\Delta(\mathcal{F})+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}\,.$
(3.8)

Moreover, if $\mathcal{F}$ is a VC class, then with probability at least
$1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{\mathrm{VC}(\mathcal{F})}{n}}-c\sqrt{\frac{t}{n}}\,$ (3.9)

for some positive constant $c$. If $\mathcal{F}$ is a non-VC class, then with
probability at least $1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-\frac{2}{\delta}\cdot\mathcal{R}_{n}(\mathcal{F})-\sqrt{\frac{2t}{n}}\,$
(3.10)

for any $\delta>0$. Here $\widehat{\mathrm{PI}}_{100\%}(X)$ and
$\widehat{\mathrm{PI}}_{100\%}(X,\delta)$ are defined in (2.3) and (2.4),
respectively.

######  Remark 3.6.

In Theorem 3.5, for non-VC class, we have obtained a coverage guarantee in
terms of $\mathcal{R}_{n}(\mathcal{F})$ and $\delta$. Suppose the covering
number of $\mathcal{F}$ grows exponentially with some $\alpha>0$. For
$0<\alpha\leq 2$, we know that
$\mathcal{R}_{n}(\mathcal{F})\leq\mathcal{O}(n^{-1/2})$ and all we need is
$\delta\gg n^{-1/2}$ for its contribution to be negligible, whereas for
$\alpha>2$, we know that
$\mathcal{R}_{n}(\mathcal{F})\leq\mathcal{O}(n^{-1/\alpha})$ and $\delta\gg
n^{-1/\alpha}$ is needed. For the former case, if we choose $t=c_{1}\log{n}$
and $\delta=c_{2}(\log{n})^{-1/2}$ for some positive constants $c_{1}$ and
$c_{2}$, we have

$\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-c_{3}\sqrt{\frac{\log{n}}{n}}\ \ \text{w.p. greater than }1-n^{-c_{1}}\,,$

for a positive constant $c_{3}$. Similar choices of $(t,\delta)$ for
$\alpha>2$ yields

$\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-c_{3}\frac{\sqrt{\log{n}}}{n^{\frac{1}{\alpha}}}\ \ \text{w.p. greater than
}1-n^{-c_{1}}\,.$

###  3.3 Estimating $f_{0}$ when $m_{0}$ is unknown

In this section, we present our analysis when we estimate both $(m_{0},f_{0})$
using either i) one-step estimation (Algorithm 1) which estimates both the
functions simultaneously or ii) two-step estimation (Algorithm 2), where we
first estimate $m_{0}$ via regression and then plugin that estimator in (2.2)
to estimate $f_{0}$. We divide our analysis for these methods into two
subsections; 3.3.1 studies Algorithm 1 and subsection 3.3.2 analyzes Algorithm
2.

####  3.3.1 One-step estimation

In one-step procedure, we estimate $(m_{0},f_{0})$ simultaneously by solving
(2.2). Analogous to (3.7), we define $\Delta(\mathcal{M},\mathcal{F})$ as:

$\Delta(\mathcal{M},\mathcal{F})=\min_{\begin{subarray}{c}m\in\mathcal{M},f\in\mathcal{F}\\\
f\geq f_{m}\end{subarray}}\mathbb{E}[f(X)-f_{0}(X)]\,.$ (3.11)

Furthermore, as evident from the constraint of (2.2), the complexity of
$(\mathcal{M},\mathcal{F})$ appears through the function
$g(X,Y)=f(X)-(Y-m(X))^{2}$ as we need $g(X_{i},Y_{i})\geq 0$ for all $1\leq
i\leq n$. Therefore, we define a function class $\mathcal{G}$ as:

$\mathcal{G}=\left\\{g(x,y)=f(x)-(y-m(x))^{2}:f\in\mathcal{F},m\in\mathcal{M}\right\\}\,.$
(3.12)

Our next theorem provides a bound for the width of the estimated prediction
interval obtained via Algorithm 1 as well as a coverage guarantee in terms of
the model mis-specification error and the complexity.

######  Theorem 3.7.

Let $\hat{f}$ be the solution of (2.2). Then for any $t\geq 0$, with
probability at least $1-2e^{-t}$ the following bound holds:

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq\Delta(\mathcal{M},\mathcal{F})+4\mathcal{R}_{n}(\mathcal{F})+4{B_{\mathcal{F}}}\sqrt{\frac{t}{2n}}\,.$
(3.13)

Moreover, if $\mathcal{G}$ is a VC class, then with probability at least
$1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{\mathrm{VC}(\mathcal{G})}{n}}-c\sqrt{\frac{t}{n}}\,$ (3.14)

for some positive constant $c$. If $\mathcal{G}$ is a non-VC class, then with
probability at least $1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-\frac{2}{\delta}\cdot\mathcal{R}_{n}(\mathcal{G})-\sqrt{\frac{2t}{n}}\,$
(3.15)

for any $\delta>0$.

Note that, the upper bound (3.13) incorporates the complexity of $\mathcal{M}$
only through the bias term $\Delta(\mathcal{M},\mathcal{F})$. Therefore, if
$\mathcal{M}$ is highly complex (say set of all measurable functions), then
$\Delta(\mathcal{M},\mathcal{F})$ is small, and the bound on (3.13) improves.
However, we pay the price in the coverage, which explicitly depends on the
complexity of $\mathcal{G}$. Here the complexity of $\mathcal{G}$ is typically
the maximum of complexity of $\mathcal{M}$ and $\mathcal{F}$, i.e.
$\mathcal{R}_{n}(\mathcal{G})=\mathcal{O}(\max\\{\mathcal{R}_{n}(\mathcal{M}),\mathcal{R}_{n}(\mathcal{F})\\})$.
To understand this trade-off, consider a simple example where $\mathcal{M}$ is
a collection of all measurable functions. Then we can always choose $\hat{m}$
such that $Y_{i}=\hat{m}(X_{i})$ for all $1\leq i\leq n$ (note that a $n$
degree polynomial is sufficient to perfectly fit the data). Consequently, for
any $\mathcal{F}$, we choose $f\in\mathcal{F}$ which has the smallest
empirical mean (It is $0$ if $0\in\mathcal{F}$), without any regard to the
prediction interval leading to poor coverage. That is why we need to control
the complexity of $\mathcal{G}$ for adequate coverage as illustrated in (3.14)
and (3.15). The choice of $\delta$ aligns with the discussion in Remark 3.6.

####  3.3.2 Two-step procedure

In this subsection, we analyze the two-step procedure as outlined in Algorithm
2. As the one-step procedure optimizes over two function classes $\mathcal{M}$
and $\mathcal{F}$ simultaneously, it can be computationally challenging. For
example, $\mathcal{M}$ is a class of smooth $k$ times differentiable functions
(e.g. Hölder class, Sobolev class, or Besov class of functions) then it is not
immediate how to optimize $m$ and $f$ efficiently in one shot. To overcome
this issue, we resort to Algorithm 2, where we first learn $\hat{m}$ and then
solve (2.2) with $\hat{m}$ fixed, i.e.:

$\displaystyle\min_{f\in\mathcal{F}}$ $\displaystyle\ \ \ \
\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$ (3.16) $\displaystyle{\sf subject\ to}$
$\displaystyle\ \ \ \ f(x_{i})\geq(y_{i}-\hat{m}(x_{i}))^{2}\ \ \ \forall\
1\leq i\leq n\,.$

Since $\hat{m}$ is fixed and we only optimize over the function class
$\mathcal{F}$, this two-step procedure is more computationally efficient
compared to the one-step procedure.

Suppose that $\hat{m}$ is a consistent estimator, i.e.
$\|\hat{m}(x)-m_{0}(x)\|_{\infty}$ is small with high probability. One may
expect that the optimal value of (3.16) is close to the optimal value of
(3.5), as the constraints in both problems are close to each other. However,
this is not always the case as it heavily relies on $\mathcal{F}$. Consider a
toy example where $\mathcal{F}=\\{f_{0}(x),f_{0}(x)+1\\}$ and
$\hat{m}(x)=m_{0}(x)+\varepsilon$. Here $\varepsilon$ is an arbitrarily small
constant. The optimal solution for (3.16) is $f_{0}(x)+1$ since $f_{0}(x)$
does not satisfy the constraint, whereas the optimal solution for (3.5) is
$f_{0}(x)$. Hence the difference between the two optimal values is $1$, which
does not depend on $\varepsilon$ due to the lack of richness of the class
$\mathcal{F}$ for a continuity. In other words, if $\mathcal{F}$ does not
admit some sort of _continuity_ with respect to the perturbation of the
constraints, then it is not possible to quantify the closeness of the optimal
value of (3.16) and (3.5) with respect to the closeness between $\hat{m}$ and
$m_{0}$. We now state our continuity assumption; define $f^{\star}$ to be the
_best approximator_ of $f_{0}$ in $\mathcal{F}$, i.e.
$f^{\star}={\arg\min}_{f\in\mathcal{F}:f\geq f_{0}}\mathbb{E}[f(X)]$. Our
continuity assumption requires that $f^{\star}$ should not be an isolated
point, i.e. given $c>0$, $f^{\star}+c\in\mathcal{F}$ for all _small_ $c$. In
other words, our assumption brings certain continuity to our optimization
problem with respect to the perturbation of $m_{0}$. We now formally state the
assumption:

######  Assumption 3.8.

Let $f^{\star}={\arg\min}_{f\in\mathcal{F}:f\geq f_{0}}\mathbb{E}[f(X)]$.
There exists $D>0$ such that for any $c\in[0,D]$, we have
$f^{\star}+c\in\mathcal{F}$.

Based on the above assumption, we now present a theoretical guarantee for the
two-step estimation below:

######  Theorem 3.9.

Let $\hat{f}$ be the solution of (3.16). Suppose that Assumption 3.8 holds and
$\|\hat{m}-m_{0}\|_{\infty}\leq c_{1}D$. Then for any $t>0$, with probability
at least $1-2e^{-t}$, we have

$\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq\Delta(\mathcal{F})+c_{2}\|\hat{m}-m_{0}\|_{\infty}+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}},$

where $c_{1}$ and $c_{2}$ are constants depending on the boundness of $Y$ and
$\mathcal{M}$. Moreover, if $\mathcal{F}$ is a VC class, then with probability
at least $1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{\mathrm{VC}(\mathcal{F})}{n}}-c\sqrt{\frac{t}{n}}\,$ (3.17)

for some constant $c$. If $\mathcal{F}$ is a non-VC class, then with
probability at least $1-e^{-t}$ the following coverage guarantee holds:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-\frac{2}{\delta}\cdot\mathcal{R}_{n}(\mathcal{F})-\sqrt{\frac{2t}{n}}\,$
(3.18)

for any $\delta>0$.

We now compare our theoretical findings for the one-step and the two-step
algorithm. In terms of the bound on the excess width, the key difference
between Theorem 3.7 and Theorem 3.9 is in the bias term; the bias term
$\Delta(\mathcal{M},\mathcal{F})$ in Theorem 3.7 is replaced by
$\Delta(\mathcal{F})+c_{2}\|\hat{m}-m_{0}\|_{\infty}$ in Theorem 3.9. As
$\hat{m}\in\mathcal{M}$ by definition, it is immediate that the bias term of
Theorem 3.9 is larger than that of Theorem 3.7 and consequently Algorithm 1
provides narrower width. This difference is negligible if
$\|\hat{m}-m_{0}\|_{\infty}\to 0$ and $m_{0}\in\mathcal{M}$.

On the other hand, the two-step estimation is more computationally efficient
than one-step since estimating $m_{0}$ using regression techniques is easy to
implement. Furthermore, the bias term in Theorem 3.9 is more interpretable as
it precisely quantifies how an estimation error in $m_{0}$ propagates through
the estimation error in the average width of the prediction intervals.
Furthermore, establishing guarantees for $\hat{m}$ obtained through two-step
estimation is also easier compared to that of one-step, as $\hat{m}$ depends
on both function classes $\mathcal{M}$ and $\mathcal{F}$ in one-step
estimation.

###  3.4 Coverage guarantee for $\alpha>0$

To achieve a user-specific coverage level $\alpha>0$, we shrink our $100\%$
prediction interval by a factor of $\hat{\lambda}_{\alpha}$ as outlined in
Algorithms 1 and 2. For the purpose of theoretical analysis, we split $n$ data
into two parts: (1) $\mathcal{D}_{1}$: $n_{1}$ samples that is used to obtain
the $100\%$ prediction interval; (2) $\mathcal{D}_{2}$: $n_{2}=n-n_{1}$
samples that is used to learn the shrink level $\hat{\lambda}_{\alpha}$. More
specifically, we first obtain the estimators $\hat{m}$ and $\hat{f}$ based on
$n_{1}$ samples via either one-step or two-step estimation and construct the
$100\%$ prediction interval as in (2.3) or (2.4) depending on the complexity
of the underlying function classes. We then select a subset of
$\mathcal{D}_{2}$ which lies inside the estimated $100\%$ prediction interval,
i.e.

$\mathcal{I}:=\left\\{1\leq i\leq
n_{2}:(y_{i}-\hat{m}(x_{i}))^{2}\leq\hat{f}(x_{i})+\delta\right\\}$

and estimate $\hat{\lambda}_{\alpha}$ based on this selected sub-samples, i.e.

$\hat{\lambda}_{\alpha}:=\inf_{\lambda}\left\\{\lambda:\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}\mathbf{1}_{\\{(y_{i}-\hat{m}(x_{i}))^{2}>\lambda(\hat{f}(x_{i})+\delta)\\}}\leq\alpha\right\\}\,.$

The choice of $\mathcal{I}$ is to avoid some _extreme observations_ in
$\mathcal{D}_{2}$; for example if $(Y_{i}-\hat{m}(X_{i}))^{2}$ is much larger
than $\hat{f}(X_{i})+\delta$ for some of the observations in
$\mathcal{D}_{2}$, then choosing $\hat{\lambda}_{\alpha}$ based on the entire
samples $\mathcal{D}_{2}$ may become unstable yielding a very large
$\hat{\lambda}_{\alpha}$, whereas intuitively we expect
$0\leq\hat{\lambda}_{\alpha}\leq 1$. To avoid this pitfall, we restrict
ourselves to the observations in $\mathcal{I}$. In the previous sections, we
have shown that with high probability
$(Y-\hat{m}(X))^{2}\leq\hat{f}(X)+\delta$ and therefore $\mathcal{I}$ is
expected to be close to $\mathcal{D}_{2}$.

Define $p_{n_{1}}$ to be the following conditional probability:

$p_{n_{1}}=\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\mid\mathcal{D}_{1}\right)\,.$

In the previous coverage guarantees and Remark 3.6, we have discussed the
choice of $\delta$ and quantified the rate of convergence of $p_{n_{1}}$ to
$1$ in terms of the complexity of the corresponding function class. We prove
the following theorem to establish a guarantee of our prediction interval for
$100(1-\alpha)\%$ coverage:

######  Theorem 3.10.

Let $\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)$ be the output of
afforementioned technically modified Algorithm 1. Then under some weak
continuity condition (i.e. $(y-\hat{m}(x))/\hat{f}(x)$ has continuous
distribution), with probability at least $1-2{n_{2}}^{-10}$, we have

$\displaystyle\left|\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-(1-\alpha)\right|\leq
2(1-p_{n_{1}})+c\cdot\sqrt{\frac{\log n_{2}}{n_{2}}}\,$

for some absolute constants $c$ and large enough $n_{1},n_{2}$.

The coverage guarantee of the preceding theorem relies on $p_{n_{1}}$, which
can be further bounded via Theorem 3.7. In particular, if both $\mathcal{M}$
and $\mathcal{F}$ are VC class of functions (e.g. for model aggreagation type
problems), then by Theorem 3.7, with probability at least $1-n^{-c_{1}}_{1}$,
we have

$\displaystyle
p_{n_{1}}=\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{\log n_{1}}{n_{1}}}\,.$ (3.19)

For non-VC class, choosing $\delta=c_{2}(\log n_{1})^{-1/2}$ yields with
probability greater than $1-n_{1}^{-c_{1}}$:

$p_{n_{1}}\geq 1-c_{3}\sqrt{\frac{\log n_{1}}{n_{1}}},\ \ \text{ for
$0<\alpha<2$}$

and

$p_{n_{1}}\geq 1-c_{3}\frac{\sqrt{\log n_{1}}}{n^{\frac{1}{\alpha}}_{1}},\ \
\text{ for $\alpha>2$}.$

Thus, for Donsker class, with probability at least
$1-2{n_{2}}^{-10}-n^{-c_{1}}_{1}$, we have

$\displaystyle\left|\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-(1-\alpha)\right|\lesssim\sqrt{\frac{\log
n_{1}}{n_{1}}}+\sqrt{\frac{\log n_{2}}{n_{2}}}$

for large enough $n_{1},n_{2}$. For the case $\alpha>2$, with probability at
least $1-2n^{-10}_{2}-n^{-c_{1}}_{1}$, we have

$\displaystyle\left|\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-(1-\alpha)\right|\lesssim\frac{\sqrt{\log
n_{1}}}{n^{\frac{1}{\alpha}}_{1}}+\sqrt{\frac{\log n_{2}}{n_{2}}}$

for large enough $n_{1},n_{2}$.

##  4 Applications

In this section, we present several examples where our theory can be applied
to obtain a prediction interval with small width as well as desired coverage
guarantee.

###  4.1 Aggregation of Predictive Intervals

In this section we expand upon the example of the model aggregation presented
in Section 2.3. Recall that we start with some candidate estimators
$\mathcal{M}_{0}=\\{m_{1},\dots,m_{L}\\}$ of $m_{0}$ and
$\mathcal{F}_{0}=\\{f_{1},\dots,f_{K}\\}$ of $f_{0}$. And we consider the
following linear span of $\mathcal{F}_{0}$ and $\mathcal{M}_{0}$:

$\displaystyle\mathcal{F}$
$\displaystyle:=\left\\{f(x)=\sum^{K}_{j=1}\alpha_{j}f_{j}(x)\,\bigg{|}\,\alpha_{j}\geq
0\right\\}$ $\displaystyle\mathcal{M}$
$\displaystyle:=\left\\{m(x)=\sum^{L}_{j=1}\beta_{j}m_{j}(x)\,\bigg{|}\,\beta_{j}\in\mathbb{R}\right\\}.$
(4.1)

One may further truncate the functions in $\mathcal{M}$ and $\mathcal{F}$ at a
desired level $B_{\mathcal{M}}$ and $B_{\mathcal{F}}$. We then optimize over
the function classes $\mathcal{M}$ and $\mathcal{F}$ to find a suitable linear
combination that obtains small width and adequate coverage by solving (2.8).
For establishing theoretical guarantees, note that $\mathcal{M}$ and
$\mathcal{F}$ are truncated linear vector spaces generated by the functions in
$\mathcal{M}_{0}$ and $\mathcal{F}_{0}$ and consequently they are VC class of
functions as a linear subspace spanned by finitely many functions is VC with
VC dimension less than or equal to the 2 + dimension of the subspace (see
Lemma 2.6.15 of Van Der Vaart and Wellner, (1996)). Therefore, by Theorem 3.7,
we have the following Corollary of whose proof is deferred to Appendix:

######  Corollary 4.1.

We consider Algorithm 1 with $(\mathcal{M},\mathcal{F})$ defined in (4.1).
With probability at least $1-e^{-t}$, the average squared width and the
coverage of the estimated prediction interval (2.9) satisfies:

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]$
$\displaystyle\leq\Delta(\mathcal{M},\mathcal{F})+c\sqrt{\frac{K}{n}}+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}$
$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle\geq 1-c\sqrt{\frac{K+L^{2}}{n}}-c\sqrt{\frac{t}{n}}\,,$

where $c$ is some constant.

Although the above corollary is based on the one-step method, the analysis of
the two-step method follows similarly. First, we estimate $m_{0}$ via finding
the optimal linear combination of the functions of $\mathcal{M}$ and then we
use that $\hat{m}$ to estimate $f_{0}$ by solving a linear programming
problem. Assumption 3.8 is satisfied if we include a constant estimator in our
candidate estimators. We skip the details here for the sake of brevity.

###  4.2 Nonparametric Construction of Predictive Intervals

In the previous application, we assume to have some candidate functions at our
disposal and then our goal was to find a suitable linear combination with
minimal width and coverage guarantee. Here we consider the general non-
parametric setting, i.e. both $m_{0}$ and $f_{0}$ (the maximum quantile of
$\sqrt{v_{0}(x)}\xi$ given $x$) are _smooth_ functions, say they belong to a
Hölder class or a Sobolev class. Suppose
$\mathcal{X}\times\mathcal{Y}\subseteq[0,1]^{d+1}$. Recall that, a Hölder
class of functions from $[0,1]^{d}$ to ${\mathbb{R}}$, denoted by
$C_{L}^{k}([0,1]^{d})$, is defined as:

$\displaystyle C_{L}^{k}([0,1]^{d})$
$\displaystyle=\left\\{f:[0,1]^{d}\to{\mathbb{R}}:\max_{\begin{subarray}{c}\mathbf{l}=(l_{1},\dots,l_{d})\\\
\sum_{i}l_{i}\leq\lfloor
k\rfloor\end{subarray}}\sup_{x}\left|\frac{\partial^{\mathbf{l}}}{\partial
x_{1}^{l_{1}}\dots\partial x_{d}^{l_{d}}}f(x)\right|\right.$
$\displaystyle\hskip
50.00008pt\left.+\max_{\begin{subarray}{c}\mathbf{l}=(l_{1},\dots,l_{d})\\\
\sum_{i}l_{i}=\lfloor
k\rfloor\end{subarray}}\sup_{x}\frac{\left|\frac{\partial^{\mathbf{l}}}{\partial
x_{1}^{l_{1}}\dots\partial
x_{d}^{l_{d}}}f(x)-\frac{\partial^{\mathbf{l}}}{\partial
x_{1}^{l_{1}}\dots\partial x_{d}^{l_{d}}}f(y)\right|}{\|x-y\|^{k-\lfloor
k\rfloor}}\leq L\right\\}$ (4.2)

Assume $\mathcal{F}=C_{F}^{\alpha}([0,1]^{d})$ and
$\mathcal{M}=C_{M}^{\beta}([0,1]^{d})$ and the function classes are well-
specified to avoid potential biases, i.e. $f_{0}\in\mathcal{F}$ and
$m_{0}\in\mathcal{M}$. Then, it is immediate that
$\Delta(\mathcal{M},\mathcal{F})=0$. For one-step procedure, we solve (2.2) by
optimizing over $\mathcal{M}=C_{M}^{\beta}([0,1]^{d})$ and
$\mathcal{F}=C_{F}^{\alpha}([0,1]^{d})$ simultaneously. Let $\hat{m}$ and
$\hat{f}$ be the solution. Theorem 3.7 yields the following Corollary of whose
proof is deferred to Appendix.

######  Corollary 4.2.

We consider Algorithm 1 with $\mathcal{F}=C_{F}^{\alpha}([0,1]^{d})$ and
$\mathcal{M}=C_{M}^{\beta}([0,1]^{d})$. Suppose that $f_{0}\in\mathcal{F}$ and
$m_{0}\in\mathcal{M}$. With probability at least $1-e^{-t}$, the following
holds:

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\lesssim\begin{cases}\sqrt{\frac{t}{n}}\,,&\text{
if }\alpha>d/2,\\\ \frac{1}{n^{\alpha/d}}+\sqrt{\frac{t}{n}}\,,&\text{ if
}\alpha<d/2\,.\end{cases}$

and

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle\geq 1-\frac{c}{\delta\sqrt{n}}-\sqrt{\frac{2t}{n}}\hskip
36.135pt\text{if }0<\gamma<2\,,$ (4.3)
$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle\geq 1-\frac{c}{\delta
n^{\frac{1}{\gamma}}}-\sqrt{\frac{2t}{n}}\hskip 36.135pt\text{if }\gamma>2\,.$
(4.4)

where $\gamma:=\max\\{d/\alpha,(d+1)/\beta\\}$ and $c$ is a positive constant.

For the two-step UTOPIA, the results are similar, with an additional factor of
$\|\hat{m}-m_{0}\|_{\infty}$ in Theorem 3.9, where $\hat{m}$ is first
estimated by non-parametric regression. When $\mathcal{M}$ is well specified,
then by standard techniques of non-parametric regression, one may obtain
$\hat{m}$ such that

$\|\hat{m}-m_{0}\|_{\infty}=O_{p}\left((n/\log{n})^{-\frac{\beta}{2\beta+d}}\right)$

For more details, see Theorem 4.3 of Belloni et al., (2015). Last but not the
least, we also have a guarantee for $100(1-\alpha)\%$ coverage for any
$\alpha>0$ which follows immediately from Theorem 3.10 with a bound on
$p_{n_{1}}$ obtained in (4.3) or (4.4) (depending on $\gamma$).

####  4.2.1 Estimation via deep neural network

In modern statistical applications, often a predictor is estimated via deep
neural networks for the ease of implementation and its adequate approximation
ability. Estimation of non-parametric mean function under various conditions
via fully/sparsely connected neural networks and its optimality has recently
been explored in several works including but not limited to Kohler and
Krzyżak, (2016); Bauer and Kohler, (2019); Schmidt-Hieber, (2020); Kohler and
Langer, (2021); Fan et al., (2022); Bhattacharya et al., (2023). In our
problem, we may also estimate the prediction interval, i.e. $f_{0}$ via deep
neural networks. Consider the two-step UTOPIA, namely, Algorithm 2. First, we
can estimate $m_{0}$ least-squares. We then estimate $f_{0}$ from (3.16) via
deep neural networks. Setting up the notations, define $\mathcal{F}$ as the
following function class:

$\mathcal{F}_{\mathrm{DNN}}(d,N,L,1):=\\{\text{DNNs with width $N$, depth $L$,
input dim $d$ and output dim $1$}\\}$

Here $\mathcal{F}_{\mathrm{DNN}}$ is the collection of all neural networks
with at most $L$ number of hidden layers and each layer has at most $N$
neurons. We further truncate the output of the neural network at level
$B_{\mathcal{F}}$. There are several possible choices of activation functions,
i.e. ReLU, leaky-ReLU, sigmoid, tanh etc. The complexity of the underlying
neural networks also depends on the activation function. Let $\hat{f}$ be the
solution of (3.16). By Bartlett et al., (2019) (Theorem 8),
$\mathcal{F}_{\mathrm{DNN}}(d,N,L,1)$ with piecewise linear activation
functions (e.g. ReLU or leaky-ReLU) is a VC class whose VC dimension is of the
order $WL\log{W}$, where $W$ is the total number of parameters of the neural
network and $L$ is the total number of hidden layers. As an example, if we use
fully connected neural networks, then $W\sim dN+N^{2}L$ and consequently VC
dimension is $\sim N^{2}L^{2}\log{NL}$ (if $d$ is much smaller than $NL$).
Thus, we do not need to extend the estimated prediction interval by $\delta$,
i.e. we consider the prediction interval given by (2.3). Theoerem 3.9 yields
the following Corollary of whose proof is deferred to Appendix:

######  Corollary 4.3.

Suppose the true $f_{0}$ belongs to Hölder class as defined in (4.2) with
smoothness $\alpha$, i.e. $f_{0}\in\mathcal{F}=C_{F}^{\alpha}([0,1]^{d})$. We
consider the two-step method with $\mathcal{F}$ being a set of DNNs with width
$N$, depth $L$ and piecewise linear activation functions. With probability at
least $1-e^{-t}$, the following holds:

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]$
$\displaystyle\leq
c\left(\frac{NL}{\log{N}\log{L}}\right)^{-\frac{2\alpha}{d}}$
$\displaystyle\qquad+c\frac{NL\sqrt{\log{NL}}}{\sqrt{n}}+c\|\hat{m}-m_{0}\|_{\infty}+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}$

and

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)\geq
1-c\sqrt{\frac{N^{2}L^{2}\log{NL}}{n}}-c\sqrt{\frac{t}{n}}.$

Choosing $(N,L)$ such that $NL\sim n^{\frac{d}{2(2\alpha+d)}}$, we obtain the
rate $n^{-\frac{\alpha}{2\alpha+d}}$, the standard non-parametric rate (upto
$\log$ factor). One advantage of using neural network method is its ability of
to learn unknown low-dimension composition structure. If the underlying
function is a composition of low-dimensional functions, then the convergence
rate is improved to the low-dimensional convergence rate (Kohler and Langer,,
2021; Schmidt-Hieber,, 2020; Fan et al.,, 2022).

###  4.3 Prediction via kernel basis

In this subsection, we elaborate on (2.12) from Section 2.4 that uses a
reproducing kernel Hilbert space to estimate $f_{0}$. The method in Section
2.4 can further be extended to the unknown mean case. We here illustrate the
one-step method via modifying (2.12) as follows:

$\displaystyle\min_{B\in\mathbb{S}^{n\times n},m\in\mathcal{M}}$
$\displaystyle{\sf tr}(KBK)$ (4.5) $\displaystyle{\sf subject\ to}$
$\displaystyle\left\langle
K_{i},BK_{i}\right\rangle\geq\left(y_{i}-m(x_{i})\right)^{2}\ \forall\ 1\leq
i\leq n$ $\displaystyle{\sf tr}(KB)\leq r,\ \ B\succeq 0\,.$

Here the class $\mathcal{M}$ can be anything a practitioner would like to
choose; it can a suitable non-parametric class (Hölder, Sobolev, Besov etc.)
or a different RKHS with another kernel $\tilde{K}$ (which can be the same as
$K$) or even aggregate various estimate of $m$. As an example, if the
practitioner decides to model $m_{0}$ by an RKHS with kernel $\tilde{K}$, then
we can modify (4.5) as:

$\displaystyle\min_{B,{\bm{\gamma}}\in{\mathbb{R}}^{n}}$ $\displaystyle{\sf
tr}(KBK)$ (4.6) $\displaystyle{\sf subject\ to}$ $\displaystyle\left\langle
K_{i},BK_{i}\right\rangle\geq\left(y_{i}-\sum_{j=1}^{n}\tilde{K}(x_{i},x_{j})\gamma_{j}\right)^{2}\
\forall\ 1\leq i\leq n$ $\displaystyle{\sf tr}(KB)\leq r,\ \ B\succeq 0\,.$

Similarly, if $\mathcal{M}=\\{m_{1},\dots,m_{L}\\}$ is a collection of some
estimates of $m_{0}$, then one may aggregate them:

$\displaystyle\min_{B,{\bm{\gamma}}\in{\mathbb{R}}^{L}}$ $\displaystyle{\sf
tr}(KBK)$ (4.7) $\displaystyle{\sf subject\ to}$ $\displaystyle\left\langle
K_{i},BK_{i}\right\rangle\geq\left(Y_{i}-\sum_{j=1}^{L}m_{j}(X_{i})\gamma_{j}\right)^{2}\
\forall\ 1\leq i\leq n$ $\displaystyle{\sf tr}(KB)\leq r,\ \ B\succeq 0\,.$

Let $\hat{B}$ and $\hat{m}$ be the solution of (4.5). We then define
$\hat{f}(x):=\langle K_{x},\hat{B}K_{x}\rangle$ with $K_{x}\in\mathbb{R}^{n}$
and $K_{x,i}=K(x_{i},x)$ and construct the prediction interval as in (2.4). We
now present the theoretical guarantee for the obtained prediction interval,
which is a corollary of Theorem 3.7:

######  Corollary 4.4.

We consider Algorithm 1 with
$\mathcal{F}=\\{\langle\Phi(x),\mathcal{A}(\Phi)(x)\rangle_{\mathcal{H}}:\|\mathcal{A}\|_{\star}\leq
r,\mathcal{A}\succeq 0\\}$. Suppose that $\sup_{x}K(x,x)\leq b$. Let
$(\hat{m},\hat{f})$ be obtained by solving (4.5). Then we have with
probability $\geq 1-e^{-t}$:

$\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq\Delta(\mathcal{M},\mathcal{F})+\frac{16r\sqrt{b}}{\sqrt{n}}\sqrt{\mathbb{E}\left[K(X,X)\right]}+4rb\sqrt{\frac{t}{2n}}\,.$

Furthermore we have

$\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\geq
1-\frac{2}{\delta}\cdot\left(\mathcal{R}_{n}(\mathcal{M})+\frac{16r\sqrt{b}}{\sqrt{n}}\sqrt{\mathbb{E}\left[K(X,X)\right]}\right)-\sqrt{\frac{2t}{n}}\,$

for any $\delta>0$.

##  5 Simulation studies

In this section, we present some simulation studies to underscore the
applicability of UTOPIA for constructing prediction intervals with a minimum
width and adequate coverage guarantee. Here we mostly focus on the model
aggregation, where we start with a few predictors of $f_{0}$ (details are
described below) and then aggregate by solving (2.8). As mentioned previously,
the aim of UTOPIA is to combine several initial estimators of $f_{0}$ (namely
the collection $\mathcal{F}$) to find a better one. Towards that direction, we
consider the following $\mathcal{F}$ which consists of six estimators of
$f_{0}$:

  1. 1.

Estimator 1$(f_{1})$: A neural network based estimator with depth=1, width=10
that estimates the 0.6 quantile function of $(Y-m_{0}(X))^{2}\mid X=x$.

  2. 2.

Estimator 2$(f_{2})$: A fully connected feed forward neural network with
depth=2 and width=50 that estimates the 0.7 quantile function of
$(Y-m_{0}(X))^{2}\mid X=x$.

  3. 3.

Estimator 3$(f_{3})$: A quantile regression forest estimating the 0.8 quantile
function of $(Y-m_{0}(X))^{2}\mid X=x$.

  4. 4.

Estimator 4$(f_{4})$: A gradient boosting model estimating the 0.9 quantile
function of $(Y-m_{0}(X))^{2}\mid X=x$.

  5. 5.

Estimator 5$(f_{5})$: A neural network based estimator with depth=1, width=10
that estimates $\mathbb{E}[(Y-m_{0}(X))^{2}\mid X=x]$.

  6. 6.

Estimator 6$(f_{6})$: The constant function $1$.

Here the quantile estimators are obtained by minimizing the corresponding
check loss. If the mean is known, then we solve (3.5), whereas for unknown
mean we resort to the two-step estimation procedure; first, we estimate the
conditional mean function $m_{0}$ using a three-layer neural network with 50
neurons in each hidden layer and then solve (3.16). For each of the simulation
setup, we first generate $\\{x_{i},y_{i}\\}^{n}_{i=1}$ i.i.d. from the data
generating process (described below) with $n=2200$. We then use $n_{pre}=1000$
data points to estimate $\\{f_{i}:1\leq i\leq 5\\}$ and $m_{0}$ (if unknown),
$n_{opt}=100$ data points to solve the optimization problem (3.5) for known
mean case and (3.16) for unknown mean case, and finally $n_{adj}=100$ data
points to adjust the interval to achieve $95\%$ coverage. Also, we leave
$n_{t}=1000$ data points as the testing data. We compare UTOPIA with three
baseline methods: (1) standard prediction band using linear quantile
regression (LQR); (2) split conformal prediction (SplitCF) (Lei et al.,,
2018); (3) semi-definite programming (SDP) method with Gaussian kernel
proposed by Liang, (2021). In the sequel, we describe our three different
setups, including (1) a location-scale model with a known mean function; (2) a
location-scale model with an unknown mean function; and (3) non-location-scale
model.

###  5.1 Setup 1

In our first simulation setup, we consider a simple location-scale model with
the known mean function $m_{0}$ (which is taken to be $0$ here for
simplicity). More precisely we assume $Y_{i}=\sqrt{v_{0}(X_{i})}\xi_{i}$ where
$X_{i}\sim{\rm Unif}([-1,1])$, $\xi_{i}\perp\\!\\!\\!\perp X_{i}$,
$\xi_{i}\sim{\rm Unif}([-1,1])$ and $v_{0}(x)=1+25x^{4}$. The results are
depicted in Figure 2 and Table 1.

(a) UTOPIA

(b) LQR

(c) SplitCF

(d) SDP

Figure 2: The blue dots, orange line and light blue band represent the testing data, the conditional mean function and the prediction band, respectively, for Setup 1. Table 1: Coverage and Average Width for Setup 1 Y-variable | Coverage | Width
---|---|---
UTOPIA | $95.4\%$ | $5.0837$
LQR | $93.8\%$ | $7.8812$
SplitCF | $94.5\%$ | $8.6148$
SDP | $96.6\%$ | $12.9346$

###  5.2 Setup 2

In our second setup for simulation, we extend our previous setup by adding an
unknown mean to the location-scale model. More precisely here our model is
$Y_{i}=m_{0}(X_{i})+\sqrt{v_{0}(X_{i})}\xi_{i}$, where as before
$\xi_{i}\perp\\!\\!\\!\perp X_{i}$ and both are generated from ${\rm
Unif}([-1,1])$. The true mean function $m_{0}(x)=1+5x^{3}$ and
$v_{0}(x)=1+25x^{4}$. The results are presented in Figure 3 and Table 2.

(a) UTOPIA

(b) LQR

(c) SplitCF

(d) SDP

Figure 3: The blue dots, orange line and light blue band represent the testing data, the conditional mean function, and the prediction band, respectively, for Setup 2. Table 2: Coverage and Average Width for Setup 2 Y-variable | Coverage | Width
---|---|---
UTOPIA | $95.0\%$ | $5.5559$
LQR | $95.8\%$ | $10.9486$
SplitCF | $95.4\%$ | $10.9755$
SDP | $98.2\%$ | $20.2919$

###  5.3 Setup 3

Our final setup is the most general setup, where we move beyond the location-
scale model. We generate the inputs $X_{i}$ from ${\rm Unif}([-1,1])$ and the
outputs $Y_{i}$ from a conditional Laplace distribution centered around its
mean $X^{2}_{i}+5X^{4}_{i}$, where $Y$ is conditioned to lie within two times
standard deviation $\sqrt{1+25X^{4}_{i}}$. More precisely, if
$Z\sim\mathrm{Lap}(m(x),\sigma(x))$ with
$m(x)=x^{2}+5x^{4},\sigma(x)=\sqrt{1+25x^{4}}$, then $Y\sim
Z|Z\in[-2\sigma(x),2\sigma(x)]$. We summarize the results in Figure 4 and
Table 3.

(a) UTOPIA

(b) LQR

(c) SplitCF

(d) SDP

Figure 4: The blue dots, orange line and light blue band represent the testing data, the conditional mean function and the prediction band, respectively, for Setup 3. Table 3: Coverage and Average Width for Setup 3 Y-variable | Coverage | Width
---|---|---
UTOPIA | $95.8\%$ | $14.3610$
LQR | $94.5\%$ | $20.7668$
SplitCF | $95.0\%$ | $21.1687$
SDP | $95.5\%$ | $23.7630$

All the methods yield coverage close to $95\%$. In each configuration, UTOPIA
delivers the lowest average width among the four approaches. This is due to
the misspecification of the quantile in the cases of LQR, SplitCF, and SDP. In
comparison, UTOPIA aggregates various estimators, enabling it to better
capture the shape of the quantile.

Although our primary focus has been on univariate inputs in the aforementioned
setups, it’s worth noting that UTOPIA exhibits compelling results even when
applied to multivariate configurations. We invite the reader to refer to
Appendix B for some simulations under multivariate setup and the corresponding
outcomes.

##  6 Real data analysis

In this section, we apply UTOPIA to two datasets: the Fama-French financial
data (Fama and French,, 1993) and the FRED-MD macroeconomic datas (McCracken
and Ng,, 2016). We compare UTOPIA with two baseline methods: (1) standard
prediction band using linear quantile regression (LQR); (2) split conformal
prediction (SplitCF) (Lei et al.,, 2018).

###  6.1 Fama-French data

This dataset contains annual and monthly measurements of four market risks
from July 1926 until February 2023. These variables include (a) Risk-free
return rate (RF), which refers to the one-month Treasury bill rate or interest
rate, (b) Market factor (MKT), indicating the additional return on the market
(i.e., the difference between market return and interest rate), (c) Size
factor (SMB, Small Minus Big), the mean returns for small and large portfolios
based on their market capitalization, and (d) Value factor (HML, High Minus
Low), expressing the divergence in returns for value and growth portfolios.

We use MKT (as $x$) to predict the other three variables HML, SMB, and RF (as
$y$), though it is known that they are weakly dependent. We follow the two-
step estimation procedure, i.e. we first pre-train the condition mean function
and the candidate estimators and then find the best linear combination of the
candidate estimators (i.e. model aggregation) via solving the optimization
problem. Here the conditional mean function is estimated using a two-layer
neural network with 10 neurons in the hidden layer. We consider six candidate
estimators which are constructed based on:

  * (1)

a two-layer neural network with width $10$ that estimates 0.05 quantile of
$Y\mid X=x$;

  * (2)

a three-layer neural network with width $50$ that estimates 0.35 quantile of
$Y\mid X=x$;

  * (3)

a quantile regression forest that estimates 0.65 quantiles of $Y\mid X=x$;

  * (4)

a gradient boosting model that estimates 0.95 quantiles of $Y\mid X=x$;

  * (5)

a two-layer neural network with width $10$ that estimates
$\mathbb{E}[(Y-m_{0}(X))^{2}\mid X=x]$;

  * (6)

a constant $1$.

We test the performance of UTOPIA on the testing data.

We consider standardized monthly data ($n=1160$, normalized to zero mean and
unit standard deviation). For each task (i.e. use $x$ = MKT to predict $y$=
HML, SMB or RF), we remove the outliers with respect to $y$ using standard
interquartile ranges (IQR) method (after which $n=1095,1126,1143$ for $y$=
HML, SMB, RF, respectively). We then (1) assign the first $80\%$ of the data
as the training data, which is used for pretraining estimators, solving the
optimization problem, and adjusting the interval based on the user-specified
coverage level ($95\%$ in our experiment), and (2) assign the remaining $20\%$
of the data as the testing data, which is used for testing the performance of
the obtained prediction band. Our results are shown in the following pictures
and tables.

UTOPIA | LQR | SplitCF
---|---|---
|  |
|  |
|  |

Figure 5: The blue dots, orange line and light blue band represent the testing data, the estimated conditional mean function and the prediction band, respectively, for the Fama-French dataset. Table 4: Coverage and Average Width for Fama-French Data | HML | SMB | RF
---|---|---|---
| Coverage | Width | Coverage | Width | Coverage | Width
UTOPIA | $96.35\%$ | $1.7187$ | $96.02\%$ | $2.1033$ | $100.00\%$ | $2.3186$
LQR | $97.26\%$ | $1.9835$ | $98.23\%$ | $2.5460$ | $100.00\%$ | $3.9645$
SplitCF | $97.26\%$ | $1.9342$ | $97.79\%$ | $2.2629$ | $100.00\%$ | $2.8225$

Figure 5 and Table 4 summarizes the result of our experiments. All these three
methods achieve desired coverage $(95\%)$ for three different tasks. In
comparison to the linear quantile regression and split conformal method, our
approach consistently achieves minimal bandwidth.

###  6.2 FRED-MD data

FRED-MD (Federal Reserve Economic Data – Macro Data) is a large-scale dataset
that provides monthly economic and financial data for the United States from
January 1959 to January 2023. It contains over 130 monthly macroeconomic time
series, spanning a wide range of different economic sectors and sub-sectors.
These variables are grouped into 8 groups: (1) output and income, (2) labor
market, (3) housing, (4) consumption, orders and inventories, (5) money and
credit, (6) interest and exchange rates, (7) prices, and (8) stock market.

In this experiment, we choose two variables from the group money and credit,
namely: (1) REALLN (i.e. real estate loans at all commercial banks) and (2)
DTCOLNVHFNM (i.e. consumer motor vehicle loans outstanding). We aim to predict
the above two variables (as $y$) using the rest of the variables. For data
pre-processing we mostly follow the experimental setup of (Fan and Gu,, 2022,
Appendix C), which, in turn, is based on McCracken and Ng, (2016). Breifly
speaking, this first transforms each variable according to the transformation
code (namely TCODE) as elaborated in the Appendix of McCracken and Ng, (2016)
and then removes missing values. We used the transformed data from the github
of Fan and Gu, (2022), namely `transformed_modern.csv`. After the data pre-
processing, we have $n=453$ observations for $p=126$ variables.

To predict a particular variable from the rest of the variables, we use the
FARM-predict technique (e.g, see Zhou et al., (2022)). Briefly speaking FARM-
predict assumes that the covariates $X$ is generated from a factor model
$X=Bf+u$ where $f$ is the underlying factor, $B$ is the loading matrix and $u$
is the idiosyncratic component. The idea is to first estimate
$(\hat{f},\hat{B})$ and set $\hat{u}=X-\hat{B}\hat{f}$. The select a subset of
$u$ which has strong conditional correlation with the residual after fitting
on $\hat{f}$, namely, $Y-\hat{\mathbb{E}}[Y\mid\hat{f}]$, where we use linear
model to estimate $\mathbb{E}[Y\mid f]$.

In our experiment we use two factors of $X$ and top three components of
$\hat{u}$ in terms of the absolute correlation with
$Y-\hat{\mathbb{E}}[Y\mid\hat{f}]$. Upon selecting $(\hat{f},\hat{u})$ we use
those as our predictors and employ three methods as before. We also use the
same two-step estimation process and data splitting method outlined in Section
6.1. The results of our experiments are presented in Figure 6 and Table 5.

(a) UTOPIA for REALNN

(b) UTOPIA for DTCOLNVHFNM

Figure 6: The blue dots, orange line and light blue band represent the testing data, the conditional mean function and the prediction band, respectively, for the FRED-MD data. Table 5: Coverage and Average Width for FRED-MD Data | REALLN | DTCOLNVHFNM
---|---|---
| Coverage | Width | Coverage | Width
UTOPIA | $94.51\%$ | $0.00049$ | $100.00\%$ | $0.00185$
LQR | $76.92\%$ | $0.00056$ | $97.80\%$ | $0.00325$
SplitCF | $25.27\%$ | $0.01097$ | $34.07\%$ | $0.01130$

As evident from Table 5, our approach attains the desired coverage ($\approx
95\%$) in both the prediction tasks along with smallest mean-squared width
compared to linear quantile regression and split conformal methods.

##  7 Conclusion and future research

In this paper, we propose a novel technique to aggregate several candidates
for prediction interval of response variable $Y$ given some observable
covariates $X$ with minimal bandwidth and adequate coverage guarantee. Our
technique relies on simple linear programming which is easy to implement.
Besides applying our method on both synthetic and real data, we have also
established theoretical guarantees. The problem of minimizing average
bandwidth is an interesting problem and requires further in-depth research. As
a starter, there is no known characterization of the optimal prediction band
with $100(1-\alpha)\%$ coverage and minimal average bandwidth, even for the
location scale model. Although our method asymptotically achieve minimal
bandwidth for $100\%$ coverage, we don’t have theoretical guarantee on the
bandwidth for general $\alpha>0$. It is expected that if the _shape_ of the
prediction interval remains changed for small $\alpha$, then our method should
also achieve minimal bandwidth for $\alpha>0$. Another interesting extension
of our approach would be to obtain heterogeneous linear combination of the
predictors (i.e. the coefficients also depend on $x$), which is expected to be
more spatially adaptive than our current method. We leave these directions for
future research.

## References

  * Barber et al., (2022) Barber, R. F., Candes, E. J., Ramdas, A., and Tibshirani, R. J. (2022).  Conformal prediction beyond exchangeability.  arXiv preprint arXiv:2202.13415.
  * Bartlett et al., (2019) Bartlett, P. L., Harvey, N., Liaw, C., and Mehrabian, A. (2019).  Nearly-tight vc-dimension and pseudodimension bounds for piecewise linear neural networks.  The Journal of Machine Learning Research, 20(1):2285–2301.
  * Bauer and Kohler, (2019) Bauer, B. and Kohler, M. (2019).  On deep learning as a remedy for the curse of dimensionality in nonparametric regression.  The Annals of Statistics, 47(4):2261–2285.
  * Belloni et al., (2019) Belloni, A., Chernozhukov, V., Chetverikov, D., and Fernández-Val, I. (2019).  Conditional quantile processes based on series or many regressors.  Journal of Econometrics, 213(1):4–29.
  * Belloni et al., (2015) Belloni, A., Chernozhukov, V., Chetverikov, D., and Kato, K. (2015).  Some new asymptotic theory for least squares series: Pointwise and uniform results.  Journal of Econometrics, 186(2):345–366.
  * Bhattacharya et al., (2023) Bhattacharya, S., Fan, J., and Mukherjee, D. (2023).  Deep neural networks for nonparametric interaction models with diverging dimension.  arXiv preprint arXiv:2302.05851.
  * Bousquet, (2002) Bousquet, O. (2002).  A bennett concentration inequality and its application to suprema of empirical processes.  Comptes Rendus Mathematique, 334(6):495–500.
  * Chernozhukov et al., (2010) Chernozhukov, V., Fernández-Val, I., and Galichon, A. (2010).  Quantile and probability curves without crossing.  Econometrica, 78(3):1093–1125.
  * Duchi, (2009) Duchi, J. (2009).  Probability bounds.  URL: http://www. cs. berkeley. edu/jduchi/projects/probability bounds. pdf.
  * Fakoor et al., (2021) Fakoor, R., Kim, T., Mueller, J., Smola, A. J., and Tibshirani, R. J. (2021).  Flexible model aggregation for quantile regression.  arXiv preprint arXiv:2103.00083.
  * Fama and French, (1993) Fama, E. F. and French, K. R. (1993).  Common risk factors in the returns on stocks and bonds.  Journal of financial economics, 33(1):3–56.
  * Fan and Gijbels, (1996) Fan, J. and Gijbels, I. (1996).  Local polynomial modelling and its applications: monographs on statistics and applied probability 66, volume 66.  CRC Press.
  * Fan and Gu, (2022) Fan, J. and Gu, Y. (2022).  Factor augmented sparse throughput deep relu neural networks for high dimensional regression.  arXiv preprint arXiv:2210.02002.
  * Fan et al., (2022) Fan, J., Gu, Y., and Zhou, W.-X. (2022).  How do noise tails impact on deep relu networks?  arXiv preprint arXiv:2203.10418.
  * Hu and Lei, (2020) Hu, X. and Lei, J. (2020).  A distribution-free test of covariate shift using conformal prediction.  arXiv preprint arXiv:2010.07147.
  * Kearns and Schapire, (1994) Kearns, M. J. and Schapire, R. E. (1994).  Efficient distribution-free learning of probabilistic concepts.  Journal of Computer and System Sciences, 48(3):464–497.
  * Koenker and Bassett Jr, (1978) Koenker, R. and Bassett Jr, G. (1978).  Regression quantiles.  Econometrica: journal of the Econometric Society, pages 33–50.
  * Koenker and Hallock, (2001) Koenker, R. and Hallock, K. F. (2001).  Quantile regression.  Journal of economic perspectives, 15(4):143–156.
  * Kohler and Krzyżak, (2016) Kohler, M. and Krzyżak, A. (2016).  Nonparametric regression based on hierarchical interaction models.  IEEE Transactions on Information Theory, 63(3):1620–1630.
  * Kohler and Langer, (2021) Kohler, M. and Langer, S. (2021).  On the rate of convergence of fully connected deep neural network regression estimates.  The Annals of Statistics, 49(4):2231–2249.
  * Kolmogorov and Tikhomirov, (1959) Kolmogorov, A. N. and Tikhomirov, V. M. (1959).  $\varepsilon$-entropy and $\varepsilon$-capacity of sets in function spaces.  Uspekhi Matematicheskikh Nauk, 14(2):3–86.
  * Ledoux and Talagrand, (1991) Ledoux, M. and Talagrand, M. (1991).  Probability in Banach Spaces: isoperimetry and processes, volume 23.  Springer Science & Business Media.
  * Lei et al., (2018) Lei, J., G’Sell, M., Rinaldo, A., Tibshirani, R. J., and Wasserman, L. (2018).  Distribution-free predictive inference for regression.  Journal of the American Statistical Association, 113(523):1094–1111.
  * Lei et al., (2013) Lei, J., Robins, J., and Wasserman, L. (2013).  Distribution-free prediction sets.  Journal of the American Statistical Association, 108(501):278–287.
  * Lei and Wasserman, (2014) Lei, J. and Wasserman, L. (2014).  Distribution-free prediction bands for non-parametric regression.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, pages 71–96.
  * Liang, (2021) Liang, T. (2021).  Universal prediction band via semi-definite programming.  arXiv preprint arXiv:2103.17203.
  * Lu et al., (2021) Lu, J., Shen, Z., Yang, H., and Zhang, S. (2021).  Deep network approximation for smooth functions.  SIAM Journal on Mathematical Analysis, 53(5):5465–5506.
  * Marteau-Ferey et al., (2020) Marteau-Ferey, U., Bach, F., and Rudi, A. (2020).  Non-parametric models for non-negative functions.  Advances in neural information processing systems, 33:12816–12826.
  * McCracken and Ng, (2016) McCracken, M. W. and Ng, S. (2016).  Fred-md: A monthly database for macroeconomic research.  Journal of Business & Economic Statistics, 34(4):574–589.
  * Moon et al., (2021) Moon, S. J., Jeon, J.-J., Lee, J. S. H., and Kim, Y. (2021).  Learning multiple quantiles with neural networks.  Journal of Computational and Graphical Statistics, 30(4):1238–1248.
  * Pollard, (1990) Pollard, D. (1990).  Empirical processes: theory and applications.  Ims.
  * Schmidt-Hieber, (2020) Schmidt-Hieber, J. (2020).  Nonparametric regression using deep neural networks with relu activation function.  The Annals of Statistics, 48(4):1875–1897.
  * Sen, (2018) Sen, B. (2018).  A gentle introduction to empirical process theory and applications.  Lecture Notes, Columbia University, 11:28–29.
  * Stankeviciute et al., (2021) Stankeviciute, K., M Alaa, A., and van der Schaar, M. (2021).  Conformal time-series forecasting.  Advances in Neural Information Processing Systems, 34:6216–6228.
  * Tibshirani et al., (2019) Tibshirani, R. J., Foygel Barber, R., Candes, E., and Ramdas, A. (2019).  Conformal prediction under covariate shift.  Advances in neural information processing systems, 32.
  * Vaart and Wellner, (1997) Vaart, A. v. d. and Wellner, J. A. (1997).  Weak convergence and empirical processes with applications to statistics.  Journal of the Royal Statistical Society-Series A Statistics in Society, 160(3):596–608.
  * Van Der Vaart and Wellner, (1996) Van Der Vaart, A. W. and Wellner, J. A. (1996).  Weak convergence.  Springer.
  * Vapnik and Chervonenkis, (2015) Vapnik, V. N. and Chervonenkis, A. Y. (2015).  On the uniform convergence of relative frequencies of events to their probabilities.  Measures of complexity: festschrift for alexey chervonenkis, pages 11–30.
  * Vovk et al., (2005) Vovk, V., Gammerman, A., and Shafer, G. (2005).  Algorithmic learning in a random world, volume 29.  Springer.
  * Vovk et al., (2009) Vovk, V., Nouretdinov, I., and Gammerman, A. (2009).  On-line predictive linear regression.  The Annals of Statistics, pages 1566–1590.
  * Xu and Xie, (2020) Xu, C. and Xie, Y. (2020).  Conformal prediction for time series.  arXiv preprint arXiv:2010.09107.
  * Zhou et al., (2022) Zhou, Y., Xue, L., Shi, Z., Wu, L., and Fan, J. (2022).  Measuring housing vitality from multi-source big data and machine learning.  Journal of the American Statistical Association, 117(539):1045–1059.

##  Appendix A Proofs

###  A.1 Proof of Lemma 2.1

First we show that $(l^{\star},u^{\star})$ is a feasible solution of (2.1).

$\displaystyle\mathbb{P}\left(m_{0}(X)+\sqrt{v_{0}(X)}F_{\xi\mid
X}^{-1}(0)\leq Y\leq m_{0}(X)+\sqrt{v_{0}(X)}F_{\xi\mid X}^{-1}(1)\right)$
$\displaystyle=\mathbb{P}\left(F_{\xi\mid
X}^{-1}(0)\leq\frac{Y-m_{0}(X)}{\sqrt{v_{0}(X)}}\leq F_{\xi\mid
X}^{-1}(1)\right)$ $\displaystyle=\mathbb{P}\left(F_{\xi\mid
X}^{-1}(0)\leq\xi\leq F_{\xi\mid X}^{-1}(1)\right)=1\,.$

For optimality, we prove even a stronger statement: For any other feasible
pair of functions $(l,u)$, we have $l^{\star}(X)\geq l(X)$ and
$u^{\star}(X)\leq u(X)$ almost surely. If $(l,u)$ is feasible, then:

$\displaystyle\mathbb{P}(Y\leq u(X))=1\implies$
$\displaystyle\mathbb{E}_{X}(\mathbb{P}(Y\leq u(X)\mid X))=1$
$\displaystyle\implies$ $\displaystyle\mathbb{P}(Y\leq u(X)\mid X)=1\text{
almost surely }$ $\displaystyle\implies$ $\displaystyle u(X)\geq F_{Y\mid
X}^{-1}(1)=u^{\star}(X)\text{ almost surely }$

A similar calculation obtains that $l(X)\geq l^{\star}(X)$ almost surely. This
completes the proof.

###  A.2 Proofs of Lemma 2.4

For any feasible solution $f\in\mathcal{F}$, we have

$\displaystyle\mathbb{E}\left[P\left(f(X)\geq\left(Y-m(X)\right)^{2}\mid
X=x\right)\right]=P\left(f(X)\geq\left(Y-m(X)\right)^{2}\right)=1.$

Thus, it holds almost surely for all $x\in\mathcal{X}$ that

$\displaystyle P\left(f(X)\geq\left(Y-m(X)\right)^{2}\mid X=x\right)=1.$

By the definition of $f_{m}$, we have $f(x)\geq f_{m}(x)$ for all
$x\in\mathcal{X}$. Because of the minimization criterion, (2.6) has a unique
solution $f_{m}$. In the sequel, we prove (2.7). By the definition of $f_{0}$,
for all $x\in\mathcal{X}$, we have

$\displaystyle P\left(\left(Y-m_{0}(X)\right)^{2}\leq f_{0}(X)\mid
X=x\right)=1.$

Also, note that

$\displaystyle\left(Y-m(X)\right)^{2}\leq\left(|Y-m_{0}(X)|+|m_{0}(X)-m(X)|\right)^{2}.$

Thus, we have for all $x\in\mathcal{X}$

$\displaystyle
P\left(\left(Y-m(X)\right)^{2}\leq\left(\sqrt{f_{0}(X)}+|m_{0}(X)-m(X)|\right)^{2}\,\Big{|}\,X=x\right)=1$

By the definition of $f_{m}$, it holds that

$\displaystyle f_{m}(x)\leq\left(\sqrt{f_{0}(x)}+|m_{0}(x)-m(x)|\right)^{2},\
\forall x\in\mathcal{X}.$ (A.1)

Next, we show equality holds for (A.1). By the definition of $f_{0}$, for any
$x\in\mathcal{X}$ and ${\epsilon}>0$, we have

$\displaystyle
P\left(\left(Y-m_{0}(X)\right)^{2}>f_{0}(X)-{\epsilon}\,\Big{|}\,X=x\right)>0.$

By Assumption 2.2, we have for any $x\in\mathcal{X}$

$\displaystyle
P\left(Y-m_{0}(X)>\sqrt{f_{0}(X)-{\epsilon}}\,\Big{|}\,X=x\right)=P\left(Y-m_{0}(X)<-\sqrt{f_{0}(X)-{\epsilon}}\,\Big{|}\,X=x\right)>0.$

For any fixed $x\in\mathcal{X}$, we define the following event

$\Omega_{x}:=\begin{cases}\left\\{Y-m_{0}(x)>\sqrt{f_{0}(x)-{\epsilon}}\right\\}&\text{if
}m_{0}(x)\geq m(x)\\\
\left\\{Y-m_{0}(x)<-\sqrt{f_{0}(x)-{\epsilon}}\right\\}&\text{if
}m_{0}(x)<m(x).\end{cases}$ (A.2)

Thus, $P(\Omega_{x}\mid X=x)>0$. Under the event $\Omega_{x}$, we have for any
${\epsilon}>0$

$\displaystyle\left(Y-m(x)\right)^{2}$
$\displaystyle=\left(Y-m_{0}(x)\right)^{2}+\left(m_{0}(x)-m(x)\right)^{2}+2\left(Y-m_{0}(x)\right)\left(m_{0}(x)-m(x)\right)$
$\displaystyle\geq
f_{0}(x)-{\epsilon}+\left(m_{0}(x)-m(x)\right)^{2}+2\sqrt{f_{0}(x)-{\epsilon}}\cdot|m_{0}(x)-m(x)|$
$\displaystyle=\left(\sqrt{f_{0}(x)-{\epsilon}}+|m_{0}(x)-m(x)|\right)^{2}.$

Let ${\epsilon}\rightarrow 0$. Then, under the event $\Omega_{x}$, we have

$\displaystyle\left(\sqrt{f_{0}(x)}+|m_{0}(x)-m(x)|\right)^{2}\leq\left(Y-m(x)\right)^{2}.$

Consequently, it holds that

$\displaystyle
P\left(\left(\sqrt{f_{0}(X)}+|m_{0}(X)-m(X)|\right)^{2}\leq\left(Y-m(X)\right)^{2}\,\Big{|}\,X=x\right)$
$\displaystyle=P\left(\left(\sqrt{f_{0}(x)}+|m_{0}(x)-m(x)|\right)^{2}\leq\left(Y-m(x)\right)^{2}\,\Big{|}\,X=x\right)$
$\displaystyle\geq P(\Omega_{x}\mid X=x)>0$

By the definition of $f_{m}$, we have

$\displaystyle P\left(\left(Y-m(X)\right)^{2}\leq
f_{m}(X)\,\big{|}\,X=x\right)=1$

Thus, we have

$\displaystyle P\left(\left(\sqrt{f_{0}(X)}+|m_{0}(X)-m(X)|\right)^{2}\leq
f_{m}(X)\,\Big{|}\,X=x\right)>0,$

which implies

$\displaystyle\left(\sqrt{f_{0}(x)}+|m_{0}(x)-m(x)|\right)^{2}\leq f_{m}(x).$
(A.3)

Combining (A.1) and (A.3), we finish the proofs.

###  A.3 Proof of Proposition 3.3

The proof of the corollary relies on the bound of the Rademacher complexity
under various growth conditions.

Polynomial growth and Exponential growth with $0<\alpha<2$: As
$\|f\|_{\infty}\leq B$, we can take $F(x)=B$ to be the envelope function.
Applying (Vaart and Wellner,, 1997, Lemma 2.14.1) we obtain:

$\displaystyle\mathcal{R}_{n}(\mathcal{F})$
$\displaystyle=\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}f(X_{i})\right]$
$\displaystyle\leq\frac{B}{\sqrt{n}}\sup_{Q}\int_{0}^{1}\sqrt{1+\log{\mathcal{N}({\epsilon}B,L_{2}(Q),\mathcal{F})}}\
d{\epsilon}$
$\displaystyle\leq\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\log{\mathcal{N}({\epsilon},L_{\infty},\mathcal{F})}}\
d{\epsilon}$
$\displaystyle\leq\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\alpha\log{\frac{1}{{\epsilon}}}}\
d{\epsilon}\leq C\sqrt{\frac{\alpha}{n}}\,.$

Now for polynomial growth we have:

$\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\log{\mathcal{N}({\epsilon},L_{\infty},\mathcal{F})}}\
d{\epsilon}\leq\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\alpha\log{\frac{1}{{\epsilon}}}}\
d{\epsilon}\leq C\sqrt{\frac{\alpha}{n}}\,.$

And for exponential growth with $0\leq\alpha<2$:

$\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\log{\mathcal{N}({\epsilon},L_{\infty},\mathcal{F})}}\
d{\epsilon}\leq\frac{1}{\sqrt{n}}\int_{0}^{B}\sqrt{1+\left(\frac{1}{{\epsilon}}\right)^{\alpha}}\
d{\epsilon}\leq\frac{C}{\sqrt{n}}$

Here the constant $C$ depends on $(B,\alpha)$.

Exponential growth with $\alpha>2$:  For $\alpha>2$, i.e. for non-Donsker
class, the calculation is different. Denote by
$\|f\|_{n}=\sqrt{(1/n)\sum_{i}f^{2}(X_{i})}$. Fix $\tau>0$. Let
$N_{\tau}=\mathcal{N}(\tau,L_{\infty},\mathcal{F})$ and
$\mathcal{F}_{\tau}=\\{f_{1},\dots,f_{N_{\tau}}\\}$ be the center of covers.
Therefore for each $f$, there exists $\tau(f)\in\\{1,\dots,N_{\tau}\\}$ such
that $\|f-f_{\tau(f)}\|_{\infty}\leq\tau$. Therefore, we have:

$\displaystyle\mathcal{R}_{n}(\mathcal{F})$
$\displaystyle=\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}f(X_{i})\right]$
$\displaystyle=\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}\left(f(X_{i})-f_{\tau(f)}(X_{i})+f_{\tau(f)}(X_{i})\right)\right]$
$\displaystyle\leq\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}\left(f(X_{i})-f_{\tau(f)}(X_{i})\right)\right]+\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}_{\tau}}\frac{1}{n}\sum_{i}{\epsilon}_{i}f(X_{i})\right]$
$\displaystyle\leq\tau+\frac{1}{\sqrt{n}}\mathbb{E}_{{\epsilon},X}\left[\max_{f\in\mathcal{F}_{\tau}}\frac{1}{\sqrt{n}}\sum_{i}{\epsilon}_{i}f(X_{i})\right]\,.$
(A.4)

Now by Dudley’s theorem, we have

$\displaystyle\mathbb{E}_{{\epsilon},X}\left[\max_{f\in\mathcal{F}_{\tau}}\frac{1}{\sqrt{n}}\sum_{i}{\epsilon}_{i}f(X_{i})\right]$
$\displaystyle\leq
C\left(\int_{\frac{\tau}{4}}^{B}\sqrt{\log{\mathcal{N}({\epsilon},L_{\infty},\mathcal{F}_{\tau})}}\
d{\epsilon}+\int^{\frac{\tau}{4}}_{0}\sqrt{\log{\mathcal{N}({\epsilon},L_{\infty},\mathcal{F}_{\tau})}}\
d{\epsilon}\right)$ $\displaystyle\leq
C_{1}\left(\int_{\frac{\tau}{4}}^{B}{\epsilon}^{-\frac{\alpha}{2}}\
d{\epsilon}+\int^{\frac{\tau}{4}}_{0}\tau^{-\frac{\alpha}{2}}\
d{\epsilon}\right)$ $\displaystyle\leq C_{2}\tau^{1-\frac{\alpha}{2}}\,.$

Here $C_{2}$ depends on $(B,\alpha)$. Using this in (A.3) we obtain:

$\mathcal{R}_{n}(\mathcal{F})\leq\tau+\frac{C_{2}}{\sqrt{n}}\tau^{1-\frac{\alpha}{2}}$

As this is true for all $\tau>0$, we have:

$\displaystyle\mathcal{R}_{n}(\mathcal{F})$
$\displaystyle\leq\inf_{\tau>0}\left\\{\tau+\frac{C_{2}}{\sqrt{n}}\tau^{1-\frac{\alpha}{2}}\right\\}\,.$

Choosing $\tau=n^{-1/\alpha}$, we conclude the proof.

###  A.4 Proof of Theorem 3.5

Firstly, we prove the bandwidth guarantee (3.8), which mainly relies on the
following proposition:

######  Proposition A.1.

Let $\mathcal{F}$ be a class of $B_{\mathcal{F}}$-bounded functions and
$\\{x_{i}\\}^{n}_{i=1}$ be i.i.d samples drawn from some distribution. Then
for any $t>0$, with probability at least $1-2e^{-t}$, we have

$\displaystyle\left|\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\right|\leq
2\mathcal{R}_{n}(\mathcal{F})+2B_{\mathcal{F}}\sqrt{\frac{t}{2n}}$

###### Proof of Proposition A.1.

By McDiarmid’s inequality, we know with probability at least $1-e^{-t}$,

$\displaystyle\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\leq\mathbb{E}_{X}\left[\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\right]+2B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

By Giné-Zinn symmetrization, we further have

$\displaystyle\mathbb{E}_{X}\left[\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\right]\leq
2\mathbb{E}_{{\epsilon},X}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}f(X_{i})\right]=2\mathcal{R}_{n}(\mathcal{F}).$

Consequently, we have with probability at least $1-e^{-t}$,

$\displaystyle\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\leq
2\mathcal{R}_{n}(\mathcal{F})+2B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Similarly, we have with probability at least $1-e^{-t}$,

$\displaystyle\sup_{f\in\mathcal{F}}\left(\mathbb{E}[f(X)]-\frac{1}{n}\sum_{i}f(x_{i})\right)\leq
2\mathcal{R}_{n}(\mathcal{F})+2B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Thus, we prove Proposition A.1. ∎

With Propostion A.1 in hand, we now prove (3.8) in the following. Suppose
$f^{\star}$ is population optimizer over $\mathcal{F}$, i.e.

$f^{\star}={\arg\min}_{f\in\mathcal{F}:f\geq f_{0}}\mathbb{E}[f(X)]\,.$

Then we have
$\Delta(\mathcal{F})=\mathbb{E}[f^{\star}(X)]-\mathbb{E}[f_{0}(X)]$. Moreover,
$f^{\star}$ is a feasible solution of (3.5). By the optimality of $\hat{f}$,
it holds that

$\displaystyle\frac{1}{n}\sum_{i}\hat{f}(x_{i})\leq\frac{1}{n}\sum_{i}f^{\star}(x_{i}),$

which implies

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\quad=\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\frac{1}{n}\sum_{i}\hat{f}(x_{i})+\frac{1}{n}\sum_{i}\hat{f}(x_{i})-\frac{1}{n}\sum_{i}f^{\star}(x_{i})+\frac{1}{n}\sum_{i}f^{\star}(x_{i})-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\quad\leq\sup_{f\in\mathcal{F}}\left(\mathbb{E}[f(X)]-\frac{1}{n}\sum_{i}f(x_{i})\right)+\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)$
$\displaystyle\quad\leq
4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Here the last inequality follows from Propostion A.1. Consequently, we have

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]$
$\displaystyle=\Delta(\mathcal{F})+\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\leq\Delta(\mathcal{F})+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Thus, we prove (3.8). The coverage guarantee (3.10) can be seen as a special
case of Theorem 3.7 with $\mathcal{M}=\\{m_{0}\\}$. Thus we do not go into
details here.

###  A.5 Proofs of Theorem 3.7

###### Proof.

Let $\hat{m}$ and $\hat{f}$ be the solution of (2.2). And we denote

$\displaystyle
m^{\star},f^{\star}={\arg\min}_{m\in\mathcal{M},f\in\mathcal{F}:f\geq
f_{m}}\mathbb{E}[f(X)].$

Then we have
$\Delta(\mathcal{M},\mathcal{F})=\mathbb{E}[f^{\star}(X)]-\mathbb{E}[f_{0}(X)]$.
Moreover, $(m^{\star},f^{\star})$ is a feasible solution of (2.2). By the
optimality of $(\hat{m},\hat{f})$, we know that

$\displaystyle\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})\leq\frac{1}{n}\sum^{n}_{i=1}f^{\star}(x_{i}).$

It then holds that

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\quad=\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})+\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})-\frac{1}{n}\sum^{n}_{i=1}f^{\star}(x_{i})+\frac{1}{n}\sum^{n}_{i=1}f^{\star}(x_{i})-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\quad\leq\sup_{f\in\mathcal{F}}\left(\mathbb{E}[f(X)]-\frac{1}{n}\sum_{i}f(x_{i})\right)+\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)$
$\displaystyle\quad\leq
4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Here the last inequality follows from Propostion A.1. Consequently, we have

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]$
$\displaystyle=\Delta(\mathcal{M},\mathcal{F})+\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\leq\Delta(\mathcal{M},\mathcal{F})+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

We then prove the coverage guarantees (3.14) and (3.15) in the following.

Case 1: $\mathcal{G}$ is VC class. We denote
$\hat{g}(x,y)=\hat{f}(x)-(y-\hat{m}(x))^{2}$. It then holds that

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle=\mathbb{E}\left[\mathbf{1}_{\\{\hat{g}(X,Y)\geq
0\\}}\,\big{|}\,\mathcal{D}\right]$
$\displaystyle\geq\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{\hat{g}(x_{i},y_{i})\geq
0\\}}-\sup_{g\in\mathcal{G}}\left|\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{g(x_{i},y_{i})\geq
0\\}}-\mathbb{E}\left[\mathbf{1}_{\\{g(X,Y)\geq 0\\}}\right]\right|$
$\displaystyle=1-\sup_{g\in\mathcal{G}}\left|\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{g(x_{i},y_{i})\geq
0\\}}-\mathbb{E}\left[\mathbf{1}_{\\{g(X,Y)\geq 0\\}}\right]\right|.$

Since $\mathcal{G}$ is VC class, we know $\\{\mathcal{G}\geq 0\\}=\\{\\{g\geq
0\\}:g\in\mathcal{G}\\}$ is also VC class with VC dimension
$\leq\mathrm{VC}(\mathcal{G})$ (see Lemma 2.6.18 in (Van Der Vaart and
Wellner,, 1996)). By symmetrization, Haussler’s bound and Proposition 3.3 with
polynomial growth, we have

$\displaystyle\mathbb{E}\left[\sup_{g\in\mathcal{G}}\left|\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{g(x_{i},y_{i})\geq
0\\}}-\mathbb{E}\left[\mathbf{1}_{\\{g(X,Y)\geq 0\\}}\right]\right|\right]\leq
c\sqrt{\frac{\mathrm{VC}(\mathcal{G})}{n}}$

for some constant $c$. Then Theorem 2.3 in Bousquet, (2002) yields with
probability at least $1-e^{-t}$:

$\displaystyle\sup_{g\in\mathcal{G}}\left|\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{g(x_{i},y_{i})\geq
0\\}}-\mathbb{E}\left[\mathbf{1}_{\\{g(X,Y)\geq 0\\}}\right]\right|\leq
c\sqrt{\frac{\mathrm{VC}(\mathcal{G})}{n}}+c\sqrt{\frac{t}{n}}.$

Therefore we show that with probability at least $1-e^{-t}$,

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X)\,\big{|}\,\mathcal{D}_{1}\right)\geq
1-c\sqrt{\frac{\mathrm{VC}(\mathcal{G})}{n}}-c\sqrt{\frac{t}{n}}$

for some constants $c$.

Case 2: $\mathcal{G}$ is non-VC class. To prove the coverage guarantee for
non-VC class $\mathcal{G}$, we define the following function

$\displaystyle h_{\delta}(t):=\begin{cases}1&\text{if }t\leq-\delta\\\
-\frac{t}{\delta}&\text{if }-\delta<t<0\\\ 0&\text{if }t\geq 0.\end{cases}$

Note that

$\displaystyle\mathbf{1}_{\\{t<-\delta\\}}\leq
h_{\delta}(t)\leq\mathbf{1}_{\\{t<0\\}},\ \forall t\in\mathbb{R}.$

We denote $\hat{g}(x,y)=\hat{f}(x)-(y-\hat{m}(x))^{2}$. It then holds that

$\displaystyle\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle=\mathbb{E}\left[\mathbf{1}_{\\{\hat{g}(X,Y)<-\delta\\}}\,\big{|}\,\mathcal{D}\right]$
$\displaystyle\leq\mathbb{E}\left[h_{\delta}(\hat{g}(X,Y))\,\big{|}\,\mathcal{D}\right]$
$\displaystyle\leq\frac{1}{n}\sum^{n}_{i=1}h_{\delta}(\hat{g}(x_{i},y_{i}))+\sup_{g\in\mathcal{G}}\left(\mathbb{E}\left[h_{\delta}(g(X,Y))\right]-\frac{1}{n}\sum^{n}_{i=1}h_{\delta}\left(g(x_{i},y_{i})\right)\right).$

Since $(\hat{m},\hat{f})$ is feasible for (2.2), we have

$\displaystyle\frac{1}{n}\sum^{n}_{i=1}h_{\delta}(\hat{g}(x_{i},y_{i}))\leq\frac{1}{n}\sum^{n}_{i=1}\mathbf{1}_{\\{\hat{g}(x_{i},y_{i})<0\\}}=0.$

By the definition of $h_{\delta}$, we know $\|h_{\delta}\|_{\infty}\leq 1$ and
$h_{\delta}$ is $\delta^{-1}$-Lipschitz. Thus, we have with probability at
least $1-e^{-t}$,

$\displaystyle\sup_{g\in\mathcal{G}}\left(\mathbb{E}\left[h_{\delta}(g(X,Y))\right]-\frac{1}{n}\sum^{n}_{i=1}h_{\delta}\left(g(x_{i},y_{i})\right)\right)$
$\displaystyle\quad\leq
2\mathbb{E}_{\mathcal{D}}\mathbb{E}_{\mathbf{{\epsilon}}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum^{n}_{i=1}{\epsilon}_{i}h_{\delta}\circ
g(x_{i},y_{i})\right|]+\sqrt{\frac{t}{2n}}$
$\displaystyle\quad\leq\frac{2}{\delta}\cdot\mathbb{E}_{\mathcal{D}}\mathbb{E}_{\mathbf{{\epsilon}}}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum^{n}_{i=1}{\epsilon}_{i}g(x_{i},y_{i})\right]+\sqrt{\frac{t}{2n}}$
$\displaystyle\quad=\frac{2}{\delta}\cdot\mathcal{R}_{n}(\mathcal{G})+\sqrt{\frac{t}{2n}},$

where the first inequality follows from the McDiarmid’s inequality and the
symmetrization argument and the second inequality follows from Talagrand’s
contraction lemma (Ledoux and Talagrand,, 1991). Consequently, we have with
probability at least $1-e^{-t}$,

$\displaystyle\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)\leq\frac{2}{\delta}\cdot\mathcal{R}_{n}(\mathcal{G})+\sqrt{\frac{t}{2n}}.$

This completes the proof. ∎

###  A.6 Proofs of Theorem 3.9

###### Proof.

Suppose that $Y\in\mathcal{B}(0,R)$ and $\mathcal{M}$ is
$B_{\mathcal{M}}$-bounded. Note that

$\displaystyle\left|\left(y_{i}-\hat{m}(x_{i})\right)^{2}-\left(y_{i}-m_{0}(x_{i})\right)^{2}\right|$
$\displaystyle=\left|\left(2y_{i}-\hat{m}(x_{i})-m_{0}(x_{i})\right)\left(\hat{m}(x_{i})-m_{0}(x_{i})\right)\right|$
$\displaystyle\leq(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}.$

Thus, we have

$\displaystyle\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})$
$\displaystyle=\inf_{\begin{subarray}{c}f\in\mathcal{F}\\\
f(x_{i})\geq(y_{i}-\hat{m}(x_{i}))^{2}\end{subarray}}\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$
$\displaystyle\leq\inf_{\begin{subarray}{c}f\in\mathcal{F}\\\
f(x_{i})\geq(y_{i}-m_{0}(x_{i}))^{2}+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}\end{subarray}}\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$

Recall that

$\displaystyle f^{\star}={\arg\min}_{f\in\mathcal{F}:f\geq
f_{0}}\mathbb{E}[f(X)],$

which implies $f^{\star}(x_{i})\geq f_{0}(x_{i})\geq(y_{i}-m_{0}(x_{i}))^{2}$
for any $i=1\ldots n$. If
$(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}\leq D$, then by
Assumption 3.8, we have
$f^{\star}+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}\in\mathcal{F}$.
Furthermore, we have

$f^{\star}(x_{i})+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}\geq(y_{i}-m_{0}(x_{i}))^{2}+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}.$

Consequently, we have

$\displaystyle\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})$
$\displaystyle\leq\inf_{\begin{subarray}{c}f\in\mathcal{F}\\\
f(x_{i})\geq(y_{i}-m_{0}(x_{i}))^{2}+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}\end{subarray}}\frac{1}{n}\sum^{n}_{i=1}f(x_{i})$
$\displaystyle\leq\frac{1}{n}\sum^{n}_{i=1}f^{\star}(x_{i})+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}.$

Thus, it holds with probability at least $1-2e^{-t}$ that

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\leq\frac{1}{n}\sum^{n}_{i=1}\hat{f}(x_{i})-\frac{1}{n}\sum^{n}_{i=1}f^{\star}(x_{i})+2\left|\sup_{f\in\mathcal{F}}\left(\frac{1}{n}\sum_{i}f(x_{i})-\mathbb{E}[f(X)]\right)\right|$
$\displaystyle\leq(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

By the definition of $f^{\star}$, we have

$\displaystyle\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]$
$\displaystyle=\Delta(\mathcal{F})+\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f^{\star}(X)]$
$\displaystyle\leq\Delta(\mathcal{F})+(B_{\mathcal{M}}+3R)\cdot\|\hat{m}-m_{0}\|_{\infty}+4\mathcal{R}_{n}(\mathcal{F})+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

The coverage guarantee can be seen as a special case of Theorem 3.7 with
$\mathcal{M}=\\{\hat{m}\\}$. Thus we do not go into details here. ∎

###  A.7 Proof of Theorem 3.10

Recall that

$\displaystyle
p_{n_{1}}:=\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1}\right)=\mathbb{P}\left(\left(Y-\hat{m}(X)\right)^{2}\leq\hat{f}(X)+\delta\,\big{|}\,\mathcal{D}_{1}\right)$

and

$\displaystyle\mathcal{I}:=\left\\{i\in\\{1,\dots,n_{2}\\}:(y_{i}-\hat{m}(x_{i}))^{2}\leq\hat{f}(x_{i})+\delta\right\\}.$

The following lemma establishes a lower bound on $|\mathcal{I}|$:

######  Lemma A.2.

For any $\beta>0$, with probability at least
$1-e^{-\frac{p_{n_{1}}{n_{2}}}{2n^{2\beta}_{1}}}$, we have
$|\,\mathcal{I}|\geq(1-{n_{1}}^{-\beta})p_{n_{1}}{n_{2}}$.

###### Proof.

We define a sequence of ${n_{2}}$ random variables
$\\{Z_{1},\dots,Z_{n_{2}}\\}$ as follows

$Z_{i}=\mathds{1}((X_{i},Y_{i})\in\mathcal{I})=\mathds{1}(Y_{i}\in\widehat{\mathrm{PI}}_{100\%}(X_{i},\delta)).$

Conditional on $\mathcal{D}_{1}$,
$Z_{i}\overset{i.i.d.}{\sim}\mathrm{Ber}(p_{n_{1}})$ and consequently
$|\,\mathcal{I}|=\sum^{n_{2}}_{i=1}Z_{i}\mid\mathcal{D}_{1}\sim\mathrm{Bin}({n_{2}},p_{n_{1}})$.
By the Chernoff bound for binomial distribution, we have

$|\,\mathcal{I}|\geq(1-{n_{1}}^{-\beta})p_{n_{1}}{n_{2}},\ \ \text{w.p.
greater than }1-e^{-\frac{p_{n_{1}}{n_{2}}}{2n^{2\beta}_{1}}}\,.$

∎

With Lemma A.2 in hand, we are ready to prove Theorem 3.10.

###### Proof of Theorem 3.10.

For notation simplicity, we define the function $g_{\lambda}(x,y)$ as:

$\displaystyle
g_{\lambda}(x,y):=\mathbf{1}_{\\{(y-\hat{m}(x))^{2}>\lambda(\hat{f}(x)+\delta)\\}}.$

Then we have

$\displaystyle\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)$
$\displaystyle\quad=\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]$
$\displaystyle\quad=\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})+\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]-\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i}).$

Note that by the choice of $\hat{\lambda}_{\alpha}$, we have

$\displaystyle\alpha-\frac{1}{|\mathcal{I}|}\leq\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})\leq\alpha.$

Consequently, we have

$\displaystyle\left|\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-\alpha\right|$
$\displaystyle\quad\leq\frac{1}{|\mathcal{I}|}+\left|\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|.$

By Lemma A.2, we know the first term is roughly $1/{n_{2}}$ when ${n_{1}}$ is
large enough. In the sequel, we upper bound the second term. Note that

$\displaystyle\left|\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|$
$\displaystyle\leq\underbrace{\left|\frac{{n_{2}}}{|\mathcal{I}|}\cdot\frac{1}{n_{2}}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\frac{1}{n_{2}}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})\right|}_{\rm(i)}+\underbrace{\left|\frac{1}{n_{2}}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})\right|}_{\rm(ii)}$
$\displaystyle\quad+\underbrace{\left|\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|}_{\rm(iii)}.$

Upper bound of (i): For (i), we have

$\displaystyle{\rm(i)}=\left|\frac{{n_{2}}}{|\mathcal{I}|}-1\right|\left|\frac{1}{n_{2}}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})\right|\leq\left|\frac{{n_{2}}}{|\mathcal{I}|}-1\right|\cdot\frac{|\mathcal{I}|}{{n_{2}}}=1-\frac{|\mathcal{I}|}{{n_{2}}},$

where the inequality follows from the fact that
$\|g_{\hat{\lambda}_{\alpha}}\|_{\infty}\leq 1$.

Upper bound of (ii): For (ii), we have

$\displaystyle{\rm(ii)}=\frac{1}{n_{2}}\sum_{i\notin\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})\leq\frac{1}{n_{2}}\cdot({n_{2}}-|\mathcal{I}|)=1-\frac{|\mathcal{I}|}{{n_{2}}}.$

Upper bound of (iii): For (iii), we have w.p. $\geq 1-e^{-t}$

$\displaystyle{\rm(iii)}$ $\displaystyle\leq\sup_{\lambda\geq
0}\left|\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}g_{\lambda}(x_{i},y_{i})-\mathbb{E}[g_{\lambda}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|$
$\displaystyle=\sup_{\lambda\geq
0}\left|\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}g_{\lambda}(x_{i},y_{i})-\mathbb{E}[g_{\lambda}(X,Y)\mid\mathcal{D}_{1}]\right|$
$\displaystyle\leq
2\mathbb{E}_{{\epsilon},\mathcal{D}_{2}}\left[\sup_{\lambda\geq
0}\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}{\epsilon}_{i}\cdot
g_{\lambda}(x_{i},y_{i})\,\bigg{|}\,\mathcal{D}_{1}\right]+\sqrt{\frac{t}{2n_{2}}},$

where the last inequality follows from a simple application of McDiarmid’s
inequality and Giné-Zinn symmetrization. Applying (Vaart and Wellner,, 1997,
Lemma 2.14.1) we obtain:

$\displaystyle\mathbb{E}_{{\epsilon},\mathcal{D}_{2}}\left[\sup_{\lambda\geq
0}\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}{\epsilon}_{i}\cdot
g_{\lambda}(x_{i},y_{i})\,\bigg{|}\,\mathcal{D}_{1}\right]\lesssim\frac{1}{\sqrt{{n_{2}}}}\sup_{Q}\int_{0}^{1}\sqrt{\log{N({\epsilon},L_{2}(Q),\mathcal{G}_{\lambda})}}\
d{\epsilon},$

where $\mathcal{G}_{\lambda}:=\\{g_{\lambda}\mid\lambda\geq 0\\}$. Note that
we define the function class $\mathcal{G}_{\lambda}$ conditioning on
$\mathcal{D}_{1}$, i.e. fixing $(\hat{m},\hat{f})$. By (Vaart and Wellner,,
1997, Theorem 2.6.7), we have for any $0<{\epsilon}<1$ that

$\displaystyle\sup_{Q}N({\epsilon},L_{2}(Q),\mathcal{G}_{\lambda})\leq
cV(\mathcal{G}_{\lambda})(16e)^{V(\mathcal{G}_{\lambda})}\left(\frac{1}{{\epsilon}}\right)^{V(\mathcal{G}_{\lambda})-1},$

where $c$ is an universal constant and $V(\mathcal{G}_{\lambda})$ is the VC
dimension of the function class $\mathcal{G}_{\lambda}$. Notice that

$\mathcal{G}_{\lambda}=\\{\mathbf{1}_{\\{\lambda(\hat{f}(x)+\delta)-(y-\hat{m}(x))^{2}<0\\}}\,|\,\lambda\geq
0\\},$

which implies $V(\mathcal{G}_{\lambda})\leq 2$. Consequently, we have

$\displaystyle\mathbb{E}_{{\epsilon},\mathcal{D}_{2}}\left[\sup_{\lambda\geq
0}\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}{\epsilon}_{i}\cdot
g_{\lambda}(x_{i},y_{i})\,\bigg{|}\,\mathcal{D}_{1}\right]$
$\displaystyle\lesssim\frac{1}{\sqrt{{n_{2}}}}\sup_{Q}\int_{0}^{1}\sqrt{\log{N({\epsilon},L_{2}(Q),\mathcal{G}_{\lambda})}}\
d{\epsilon}$
$\displaystyle\leq\frac{1}{\sqrt{{n_{2}}}}\int_{0}^{1}\sqrt{\log\left\\{cV(\mathcal{G}_{\lambda})(16e)^{V(\mathcal{G}_{\lambda})}\left(\frac{1}{{\epsilon}}\right)^{V(\mathcal{G}_{\lambda})-1}\right\\}}\
d{\epsilon}$
$\displaystyle\leq\frac{1}{\sqrt{{n_{2}}}}\int_{0}^{1}\sqrt{\log\left\\{2c(16e)^{2}\left(\frac{1}{{\epsilon}}\right)\right\\}}\
d{\epsilon}\lesssim\frac{1}{\sqrt{{n_{2}}}}\,.$

As a result, we have w.p. $\geq 1-e^{-t}$

$\displaystyle{\rm(iii)}\leq
2\mathbb{E}_{{\epsilon},\mathcal{D}_{2}}\left[\sup_{\lambda\geq
0}\frac{1}{n_{2}}\sum^{n_{2}}_{i=1}{\epsilon}_{i}\cdot
g_{\lambda}(x_{i},y_{i})\,\bigg{|}\,\mathcal{D}_{1}\right]+\sqrt{\frac{t}{2n_{2}}}\leq
c\cdot\sqrt{\frac{t}{{n_{2}}}}$

for some constants $c$.

Combining the upper bounds for (i), (ii) and (iii), we have w.p. $\geq
1-e^{-t}$

$\displaystyle\left|\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|\leq
2\left(1-\frac{|\mathcal{I}|}{{n_{2}}}\right)+c\cdot\sqrt{\frac{t}{{n_{2}}}},$

which implies

$\displaystyle\left|\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-\alpha\right|$
$\displaystyle\quad\leq\frac{1}{|\mathcal{I}|}+\left|\frac{1}{|\mathcal{I}|}\sum_{i\in\mathcal{I}}g_{\hat{\lambda}_{\alpha}}(x_{i},y_{i})-\mathbb{E}[g_{\hat{\lambda}_{\alpha}}(X,Y)\mid\mathcal{D}_{1},\mathcal{D}_{2}]\right|$
$\displaystyle\quad\leq\frac{1}{|\mathcal{I}|}+2\left(1-\frac{|\mathcal{I}|}{{n_{2}}}\right)+c\cdot\sqrt{\frac{t}{{n_{2}}}}.$

By Lemma A.2, with probability at least
$1-e^{-\frac{p_{n_{1}}{n_{2}}}{2n^{2\beta}_{1}}}$, we have
$|\,\mathcal{I}|\geq(1-{n_{1}}^{-\beta})p_{n_{1}}{n_{2}}$ for any $\beta>0$.
Consequently, for any $t,\beta>0$, with probability at least
$1-e^{-\frac{p_{n_{1}}{n_{2}}}{2n^{2\beta}_{1}}}-e^{-t}$, we have

$\displaystyle\left|\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-\alpha\right|$
$\displaystyle\quad\leq\frac{1}{(1-{n_{1}}^{-\beta})p_{n_{1}}{n_{2}}}+2(1-p_{n_{1}}+{n_{1}}^{-\beta}p_{n_{1}})+c\cdot\sqrt{\frac{t}{{n_{2}}}}.$

We further choose $(\beta,t)$ such that $p_{n_{1}}{n_{2}}/2n^{2\beta}_{1}=t$
and $t=10\cdot\log{n_{2}}$. Then with probability at least $1-2{n_{2}}^{-10}$,
we have

$\displaystyle\left|\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-\alpha\right|$
$\displaystyle\quad\leq\frac{1}{\left(1-\sqrt{\frac{20\log{n_{2}}}{p_{n_{1}}{n_{2}}}}\right)p_{n_{1}}{n_{2}}}+2(1-p_{n_{1}})+c_{1}\cdot\sqrt{\frac{\log{n_{2}}}{{n_{2}}}}.$

When ${n_{1}},{n_{2}}$ are large enough such that

$\displaystyle
p_{n_{1}}\geq\frac{1}{2},\quad\sqrt{\frac{40\log{n_{2}}}{{n_{2}}}}\leq\frac{1}{2},$

we have

$\displaystyle\left|\mathbb{P}\left(Y\notin\widehat{\mathrm{PI}}_{100(1-\alpha)\%}(X,\delta)\,\big{|}\,\mathcal{D}_{1},\mathcal{D}_{2}\right)-\alpha\right|$
$\displaystyle\leq\frac{4}{{n_{2}}}+2(1-p_{n_{1}})+c_{1}\cdot\sqrt{\frac{\log{n_{2}}}{{n_{2}}}}$
$\displaystyle\leq
2(1-p_{n_{1}})+c_{2}\cdot\sqrt{\frac{\log{n_{2}}}{{n_{2}}}}.$

Thus, we finish the proof. ∎

###  A.8 Proof of Corollary 4.1

As $\mathcal{F}$ is the linear span of $K$ many functions, we have
$\mathrm{VC}(\mathcal{F})\leq K+2$. Then by Haussler’s bound and Proposition
3.3 with polynomial growth, we have

$\displaystyle\mathcal{R}_{n}(\mathcal{F})=\mathbb{E}\left[\sup_{f\in\mathcal{F}}\frac{1}{n}\sum_{i}{\epsilon}_{i}f(X_{i})\right]\lesssim\sqrt{\frac{\mathrm{VC}(\mathcal{F})}{n}}\lesssim\sqrt{\frac{K}{n}}\,.$

As

$\mathcal{G}=\left\\{\sum^{K}_{j=1}\alpha_{j}f_{j}(x)-\left(y-\sum^{L}_{j=1}\beta_{j}m_{j}(x)\right)^{2}\,\bigg{|}\,\alpha_{j}\geq
0,\beta_{j}\in\mathbb{R}\right\\}$

is the linear span of $\mathcal{O}(K+L^{2})$ many functions, we have
$\mathrm{VC}(\mathcal{G})\lesssim K+L^{2}$. Thus we have

$\displaystyle\mathcal{R}_{n}(\mathcal{G})$
$\displaystyle=\mathbb{E}\left[\sup_{g\in\mathcal{G}}\frac{1}{n}\sum_{i}{\epsilon}_{i}g(X_{i})\right]\lesssim\sqrt{\frac{\mathrm{VC}(\mathcal{G})}{n}}\lesssim\sqrt{\frac{K+L^{2}}{n}}\,.$

Then Coroallry 4.1 directly follows from Theorem 3.7.

###  A.9 Proof of Corollary 4.2

By Lemma 6.6 of (Sen,, 2018, 2022 version) we have:

$\log{N({\epsilon},C_{L}^{k}([0,1]^{d}),\|\cdot\|_{\infty})}\leq
c\cdot{\epsilon}^{-\frac{d}{k}}\,,$

where the constant $c$ depends on $(L,d,k)$. If $k>d/2$ then the function
class is Donsker (i.e. it satisfies (3.3) with $0<\alpha<2$), otherwise it is
non-Donsker. By Proposition 3.3, the Rademacher complexity of $\mathcal{F}$
can be bounded as:

$\mathcal{R}_{n}(\mathcal{F})\leq\begin{cases}Cn^{-\frac{1}{2}}\,,&\text{ if
}\alpha>d/2,\\\ Cn^{-\frac{\alpha}{d}}\,,&\text{ if }\alpha<d/2.\end{cases}$

Therefore, for the one-step procedure, Theorem 3.7 yields the following bound:

$\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\lesssim\begin{cases}\sqrt{\frac{t}{n}}\,,&\text{
if }\alpha>d/2,\\\ \frac{1}{n^{\alpha/d}}+\sqrt{\frac{t}{n}}\,,&\text{ if
}\alpha<d/2\,.\end{cases}$

To ensure the coverage guarantee, we need to find the complexity of
$\mathcal{G}$ as defined in (3.12). It is immediate that $\mathcal{G}\subseteq
C_{F}^{\alpha}([0,1]^{d})\oplus C_{G}^{\beta}([0,1]^{d+1})$ for some constant
$G$. Therefore, the covering number of $\mathcal{G}$ is bounded by the product
of the covering numbers of $C_{F}^{\alpha}([0,1]^{d})$ and
$C_{G}^{\beta}([0,1]^{d+1})$ and consequently $\mathcal{G}$ is a Donsker class
if $0<\gamma:=\max\\{d/\alpha,(d+1)/\beta\\}\leq 2$ and non-Donsker otherwise.
This observation, along with Proposition 3.3 and Theorem 3.7 yields with
probability at least $1-e^{-t}$:

$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle\geq 1-\frac{C}{\delta\sqrt{n}}-\sqrt{\frac{2t}{n}}\hskip
36.135pt\text{if }0<\gamma<2\,,$ (A.5)
$\displaystyle\mathbb{P}\left(Y\in\widehat{\mathrm{PI}}_{100\%}(X,\delta)\,\big{|}\,\mathcal{D}\right)$
$\displaystyle\geq 1-\frac{C}{\delta
n^{\frac{1}{\gamma}}}-\sqrt{\frac{2t}{n}}\hskip 36.135pt\text{if }\gamma>2\,.$
(A.6)

###  A.10 Proof of Corollary 4.3

Define $\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}})$ is the
truncated version of the functions of $\mathcal{F}_{\mathrm{DNN}}$ at level
$B_{\mathcal{F}}$. Here we assume fully connected neural networks unless
mentioned otherwise. A bound on the VC dimension of
$\mathcal{F}_{\mathrm{DNN}}$ can be found in (Bartlett et al.,, 2019, Lemma),
i.e.

$\mathrm{VC}(\mathcal{F}_{\mathrm{DNN}})\leq CWL\log{W}\,,$

where $W$ is the total number of parameters. Furthermore, as truncation does
not change the VC dimension, the same bound also holds for
$\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}})$. For a fully
connected neural network with width $N$ and depth $L$, the total number of
parameters is $N(d+1)+(L-2)(N^{2}+N)+1$ (including bias) which is of order
$N^{2}L$ for fixed $d$. Therefore, the VC dimension of
$\mathcal{F}_{\mathrm{DNN}}$ is of the order $N^{2}L^{2}\log{NL}$.
Furthermore, we may also truncate $f$ at level $B_{\mathcal{F}}$ to refrain it
from being unbounded. Hence, from Theorem 3.9 and using the fact the that
Rademacher complexity of a VC class is bounded by $\sqrt{V/n}$, we obtain:

$\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq\Delta(\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}}))+C\frac{NL\sqrt{\log{NL}}}{\sqrt{n}}+C\|\hat{m}-m_{0}\|_{\infty}+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}$
(A.7)

Now focus on
$\Delta(\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}}))$, we
have:

$\displaystyle\Delta(\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}}))$
$\displaystyle=\min_{\begin{subarray}{c}f\in\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}})\\\
f\geq f_{0}\text{ a.s.}\end{subarray}}\mathbb{E}[f(X)]-\mathbb{E}[f_{0}(X)]$
$\displaystyle\leq\min_{\begin{subarray}{c}f\in\mathrm{trunc}_{B_{\mathcal{F}}}(\mathcal{F}_{\mathrm{DNN}})\\\
f\geq f_{0}\end{subarray}}\|f-f_{0}\|_{\infty}$
$\displaystyle\leq\min_{\begin{subarray}{c}f\in\mathcal{F}_{\mathrm{DNN}}\\\
f\geq f_{0}\end{subarray}}\|f-f_{0}\|_{\infty}$ $\displaystyle\leq
2\min_{f\in\mathcal{F}_{\mathrm{DNN}}}\|f-f_{0}\|_{\infty}\,.$

To understand how the last inequality follows, define
$f^{*}={\arg\min}_{f\in\mathcal{F}_{\mathrm{DNN}}}\|f-f_{0}\|_{\infty}$ and
define $\delta=\|f^{*}-f_{0}\|_{\infty}$. Even if it may be possible
$f^{*}<f_{0}$, we have $f^{*}+\delta\geq f_{0}$ almost surely. As
$f^{\star}+\delta\in\mathcal{F}_{\mathrm{DNN}}$ (adding a constant does not
change the architecture), we have:

$\min_{\begin{subarray}{c}f\in\mathcal{F}_{\mathrm{DNN}}\\\ f\geq
f_{0}\end{subarray}}\|f-f_{0}\|_{\infty}\leq\|f^{*}+\delta-
f_{0}\|_{\infty}\leq
2\delta=2\min_{f\in\mathcal{F}_{\mathrm{DNN}}}\|f-f_{0}\|_{\infty}\,.$

Now, if $f_{0}$ is in Hölder/Sobolev class of functions, i.e. $\alpha$ times
differentiable with bounded derivatives, then (Lu et al.,, 2021, Theorem 1)
obtains the following approximation error:

$\min_{f\in\mathcal{F}_{\mathrm{DNN}}}\|f-f_{0}\|_{\infty}\leq
C_{1}\left(\frac{NL}{\log{N}\log{L}}\right)^{-\frac{2\alpha}{d}}$

Therefore, from (A.7) we have:

$\mathbb{E}[\hat{f}(X)\mid\mathcal{D}]-\mathbb{E}[f_{0}(X)]\leq
2C_{1}\left(\frac{NL}{\log{N}\log{L}}\right)^{-\frac{2\alpha}{d}}+C\frac{NL\sqrt{\log{NL}}}{\sqrt{n}}+C\|\hat{m}-m_{0}\|_{\infty}+4B_{\mathcal{F}}\sqrt{\frac{t}{2n}}.$

Since the VC dimension of $\mathcal{F}_{\mathrm{DNN}}$ is of the order
$N^{2}L^{2}\log{NL}$, by Theorem 3.9, we have