The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
54,980 characters
Distributional conformal prediction
\title{Distributional conformal prediction\footnote{We are grateful to the Editor, two anonymous referees, Dimitris Politis, and Allan Timmermann for valuable comments. W\"uthrich is also affiliated with CESifo and the Ifo Institute. Chernozhukov gratefully acknowledges funding by the National Science Foundation. The usual disclaimer applies.}}
\author{Victor Chernozhukov\thanks{Massachusetts Institute of Technology; 50 Memorial Drive, E52-361B, Cambridge, MA 02142, USA; Email: \url{[email removed]}} \qquad Kaspar W\"uthrich\thanks{Department of Economics, University of California San Diego, 9500 Gilman Dr., La Jolla, CA 92093, USA; Email: \url{[email removed]}} \qquad Yinchu Zhu\thanks{Brandeis University; 415 South Street, Waltham, MA 02453, USA; Email: \url{[email removed]}}}
\date{First version on arXiv: September, 17 2019 \quad This version: \today}
\maketitle
\linespread{1.15} \parskip 0in
\begin{abstract}
We propose a robust method for constructing conditionally valid prediction intervals based on models for conditional distributions such as quantile and distribution regression. Our approach can be applied to important prediction problems including cross-sectional prediction, $k$-step-ahead forecasts, synthetic controls and counterfactual prediction, and individual treatment effects prediction. Our method exploits the probability integral transform and relies on permuting estimated ranks. Unlike regression residuals, ranks are independent of the predictors, allowing us to construct conditionally valid prediction intervals under heteroskedasticity. We establish approximate conditional validity under consistent estimation and provide approximate unconditional validity under model misspecification, overfitting, and with time series data. We also propose a simple ``shape'' adjustment of our baseline method that yields optimal prediction intervals.
\medskip
\noindent \textbf{Keywords:} prediction intervals, quantile regression, distribution regression, conditional validity, model-free validity
\end{abstract}
\maketitle
\newpage
\linespread{1.25} \parskip .05in
\section{Introduction}
We develop a robust approach for constructing prediction intervals based on models for conditional distributions. The proposed method is generic and can be implemented using a great variety of flexible and powerful methods, including conventional quantile regression (QR) \citep{koenker1978}, distribution regression (DR) \citep[e.g.,][]{foresi1995conditional,chernozhukov2013inference}, as well as non-parametric and high-dimensional machine learning methods such as quantile neural networks \citep[e.g.,][]{taylor2000quantile} and quantile trees and random forests \citep[e.g.,][]{chaudhuri2002,meinshausen2006quantile}.
We observe data $\left\{(Y_t,X_t)\right\}_{t=1}^T$, where $Y_t$ is a continuous outcome of interest and $X_t$ is a $p\times 1$ vector of predictors. Our task is to predict $Y_{T+1}$ given knowledge of $X_{T+1}$. This setting encompasses many classical cross-sectional and time series prediction problems. Moreover, our approach can be applied to synthetic control settings where the goal is to predict counterfactuals in the absence of a policy intervention \citep[e.g.,][]{cattaneo2019prediction,chernozhukov2021conformal} and to the problem of predicting individual treatment effects \citep[e.g.,][]{kivaranovic2020conformal,lei2020conformal}.
With iid (or exchangeable data), standard conformal prediction methods, which are based on modeling the conditional mean, yield prediction intervals $\widehat{\mathcal{C}_{(1-\alpha)}}$ that satisfy
\begin{eqnarray}
P\left(Y_{T+1}\in \widehat{\mathcal{C}_{(1-\alpha)}} \left(X_{T+1}\right)\right)\ge 1-\alpha \label{eq:uncond_coverage}
\end{eqnarray}
for a given miscoverage level $\alpha\in (0,1)$. A prediction interval satisfying this property is said to be \emph{unconditionally} valid. Unconditionally valid prediction intervals guarantee accurate coverage on average, treating $(Y_{T+1},X_{T+1})$ and $\left\{(Y_t,X_t)\right\}_{t=1}^T$ as random.
However, in many applications, unconditional validity may be unsatisfactory. Let us consider three examples; see \cite{romano2019malice,foygel2021limits} for further examples and discussions. First, from a fairness perspective, data-driven recommendation systems should guarantee equalized coverage across protected groups, in which case the goal is to construct prediction intervals that are valid conditional on a protected attribute such as race or gender \citep[][]{romano2019malice}. Second, as in Section \ref{sec:predicting_stock_returns}, consider the problem of predicting stock returns given the realized volatility. Since the distribution of returns is more dispersed when the variance is higher, a natural prediction algorithm should yield wider prediction intervals for higher values of volatility. That is, the prediction interval should be valid conditional on the known value of realized volatility rather than on average. Third, as in Section \ref{sec:predicting_wages}, suppose our goal is to predict wages based on an individual's education and experience. An unconditionally valid prediction interval exhibits coverage 90\% on average across all individuals but may contain the true wage of high-school dropouts with no work experience with probability zero. A more useful prediction interval should exhibit correct coverage conditional on an individual's observed education and experience and contain the true wage with 90\% probability for every single individual.
Motivated by this discussion, we develop a \emph{distributional conformal prediction} (DCP) method for constructing prediction intervals that are approximately valid conditional on the full vector of predictors $X_{T+1}$, while treating $Y_{T+1}$ and $\left\{(Y_t,X_t)\right\}_{t=1}^T$ as random:
\begin{eqnarray}
P\left(Y_{T+1}\in \widehat{\mathcal{C}_{(1-\alpha)}} \left(X_{T+1}\right)\mid X_{T+1}\right) \ge 1-\alpha + o_P(1). \label{eq:cond_coverage}
\end{eqnarray}
A prediction interval satisfying property \eqref{eq:cond_coverage} as $T\rightarrow \infty$ is said to be approximately \emph{conditionally} valid.\footnote{See, for example, \cite{lei2014distribution,sesia2019comparison,foygel2021limits} for a further discussion of the difference between conditional and unconditional validity.}
While the requirement in \eqref{eq:cond_coverage} is natural in many applications, there are also other notions of conditional validity. Instead of conditioning on $X_{T+1}$ (object conditional), one can also study the conditional coverage probability given the training sample $\left\{(Y_t,X_t)\right\}_{t=1}^T$ (training conditional) or given $Y_{T+1}$ (label conditional) or combinations of them; see \cite{vovk2012conditional} for a detailed discussion. By Proposition 2 of \cite{vovk2012conditional}, inductive conformal predictions (also known as split-sample conformal predictions) automatically achieve training conditional validity as long as the training sample is large enough. In classification problems (the support of $Y_{T+1}$ is a finite set), label conditional validity is often of great interest as it is important to know the error rates for different categories and provides useful information on false positive and false negative rates \citep{vovk2012conditional}. In \cite{vovk2012conditional}, label conditional validity is achieved by forming the conformity score within each category. Both training and label conditional validity can be achieved in a distribution-free way, i.e., for a given procedure, the conditional validity holds for any distribution of the data.
However, object conditional validity in the sense of \eqref{eq:cond_coverage} cannot be achieved in a distribution-free way for non-trivial predictions. By \cite{{vovk2012conditional,lei2014distribution,foygel2021limits}}, any prediction set satisfying \eqref{eq:cond_coverage} for every probability distribution of $(X_t,Y_t)$ has infinite Lebesgue measure with non-trivial probability. Therefore, we only aim to achieve \eqref{eq:cond_coverage} for a limited class of probability distributions. The construction of the proposed prediction set $ \widehat{\mathcal{C}_{(1-\alpha)}}$ relies on learning the conditional distribution $Y_t\mid X_t$ and we only hope for conditional validity in \eqref{eq:cond_coverage} in the class of distributions that can be learned well. In particular, this class of distributions are those satisfying our regularity conditions.
Our empirical results demonstrate the importance of using DCP instead of standard conformal prediction methods based on modeling the conditional mean. When predicting daily stock returns in Section \ref{sec:predicting_stock_returns}, the coverage probability of the 90\% mean-based conformal prediction interval can drop to around 50\% when the realized volatility is high. By contrast, DCP provides a coverage probability close to 90\% for all values of realized volatility. This finding is important since volatility tends to be high during periods of crisis when accurate risk assessments are most needed. When predicting wages in Section \ref{sec:predicting_wages}, we find that the DCP prediction intervals contain the true wage with probability close to 90\% for most individuals, whereas standard mean-based conformal prediction intervals either substantially under- or overcover.
To motivate DCP, note that a conditionally valid prediction interval is given by
\begin{equation}
\left[Q\left(\frac{\alpha}{2} , x\right),Q\left(1-\frac{\alpha}{2}, x \right)\right], \label{eq:ideal_prediction_interval}
\end{equation}
where $Q(\tau, x)$ is the $\tau$-quantile of $Y_t$ given $X_t=x$. To implement the prediction interval \eqref{eq:ideal_prediction_interval}, a plug-in approach would replace $Q$ with a consistent estimator $\hat{Q}$
\begin{equation}
\left[\hat{Q}\left(\frac{\alpha}{2},x \right),\hat{Q}\left(1-\frac{\alpha}{2},x \right)\right]. \label{eq:naive_prediction_interval}
\end{equation}
This approach exhibits two well-known drawbacks. First, it will often exhibit undercoverage in finite samples \citep[e.g.,][]{romano2019conformalized}. Second, it is neither conditionally nor unconditionally valid under misspecification.
We build upon conformal prediction \citep{vovk2005algorithmic,vovk2009online} and use the conditional ranking as a conformity score. This choice is particularly useful when working with regression models for conditional distributions such as QR and DR.\footnote{This transformation is also very useful in other prediction problems \citep[e.g.,][]{politis2015modelfree}.} Our method is conditionally valid under correct specification, while the construction of the procedure as a conformal prediction method guarantees the unconditional validity under misspecification. Let $F(y,x)=P(Y_t\leq y \mid X_t=x)$ denote the conditional cumulative distribution function (CDF) of $Y_t$ given $X_t=x$. Throughout the paper, we assume that $F(\cdot,X_t)$ is a continuous function almost surely. Our method is based on the probability integral transform, which states that the \emph{conditional rank}, $U_t:=F\left( Y_t,X_t\right)$, has the uniform distribution on $(0,1)$ and is independent of $X_t$.
To construct the prediction interval, we test the plausibility of each $y \in \mathbb{R}$. By the probability integral transform, conditional on $X_{T+1}$, $F(Y_{T+1},X_{T+1}) $ belongs to $[\alpha/2,1-\alpha/2] $ with probability $1-\alpha$. Thus, collecting all values $y\in \mathbb{R}$ satisfying $F(y,X_{T+1})\in[\alpha/2,1-\alpha/2] $ yields a conditionally valid prediction interval in the sense of \eqref{eq:cond_coverage}. We operationalize this idea by proposing a conformal prediction procedure based on the estimated ranks, $\hat{U}_{t}^{(y)}:=\hat{F}^{(y)}(Y_t,X_t)$. For each $y \in \mathbb{R}$, $\hat{F}^{(y)}$ is an estimator of $F$ obtained based on the augmented data, $\{(Y_t,X_t)\}_{t=1}^{T+1}$, where $Y_{T+1}=y$. Data augmentation is a key feature of conformal prediction. It implies the model-free unconditional exact finite-sample validity with iid (or exchangeable) data and, thus, guards against model misspecification and overfitting. Without data augmentation, the resulting prediction intervals are not exactly valid, not even with correct specification and iid data.
Our baseline method asymptotically coincides with the oracle interval in \eqref{eq:ideal_prediction_interval}.
This oracle interval may not be the shortest possible prediction interval in general. Therefore, we also develop a simple and easy-to-implement adjustment of our baseline method for improving efficiency, which we refer to as \emph{optimal DCP}. In Section \ref{sec:predicting_wages}, we show empirically that optimal DCP yields shorter prediction intervals than baseline DCP when the conditional distribution is skewed.
We establish the following theoretical performance guarantees for the baseline and optimal DCP.
\begin{enumerate}\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item[(i)] Asymptotic conditional validity under consistent estimation of the conditional CDF
\item[(ii)] Unconditional validity under model misspecification:
\vspace{-5pt}
\begin{itemize}\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item[(a)] Finite-sample validity with iid\ (or exchangeable) data
\item[(b)] Asymptotic validity with time series data
\end{itemize}
\item[(iii)] For optimal DCP:
\begin{itemize}\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item[(a)] Under weak conditions: asymptotic conditional validity and optimality (shortest length)
\item[(b)] Under strong conditions: asymptotic convergence to the optimal prediction interval
\end{itemize}
\end{enumerate}
\subsection{Motivating Example}
\label{sec:intro_conformal_prediction}
We illustrate the advantages of DCP relative to mean-based conformal prediction \citep[e.g.,][]{lei2018distributionfree} based on the following simple analytical example.
\begin{equation}
Y_t=X_t+X_t\varepsilon_t, \quad X_t\overset{iid}\sim \text{Uniform}(0,1), \quad \varepsilon_t\overset{iid}\sim N(0,1).\label{eq:dgp_motivating_example}
\end{equation}
Our motivating example draws on \cite{koenker1982robust,koenker2005book,lei2018distributionfree,romano2019conformalized}. We focus on the population conformal prediction (or oracle) problem under correct specification and abstract from finite sample issues.
Mean-based conformal prediction is based on the residuals $R_t=Y_t-E(Y_t\mid X_t)=Y_t-X_t=X_t\varepsilon_t$. The mean-based prediction interval is
\begin{equation}
\mathcal{C}^{\rm reg}_{(1-\alpha)}(x)=\left[x-Q_{|R|}(1-\alpha),x+Q_{|R|}(1-\alpha)\right], \label{eq:interval_regression}
\end{equation}
where $Q_{|R|}(1-\alpha)$ is the $(1-\alpha)$-quantile of the distribution of $|R_{t}|$. An important property and drawback of $\mathcal{C}^{\rm reg}_{(1-\alpha)}$ is that its length, $2 \cdot Q_{|R|}(1-\alpha)$, is fixed and does not depend on $X_{T+1}=x$ \citep{lei2018distributionfree,romano2019conformalized}. This feature implies that $\mathcal{C}^{\rm reg}_{(1-\alpha)}$ is not adaptive to the heteroskedasticity in the location-scale model \eqref{eq:dgp_motivating_example} and not conditionally valid.
DCP is based on the ranks $U_t=\Phi\left(\varepsilon_t \right)$, where $\Phi(\cdot)$ is the CDF of $N(0,1)$. The DCP prediction interval is
\begin{eqnarray}
\mathcal{C}^{\rm dcp}_{(1-\alpha)}(x)
=\left[x-x\cdot Q_{|\varepsilon|}(1-\alpha),x+x\cdot Q_{|\varepsilon|}(1-\alpha)\right], \label{eq:interval_regression}
\end{eqnarray}
where $Q_{|\varepsilon|}(1-\alpha)=\Phi^{-1}(1-\alpha/2)$ is the $(1-\alpha)$-quantile of $|\varepsilon_t|$. Unlike $\mathcal{C}^{\rm reg}_{(1-\alpha)}$, the length of $\mathcal{C}^{\rm dcp}_{(1-\alpha)}$, $2x\cdot Q_{|\varepsilon|}(1-\alpha)$, depends on $X_{T+1}=x$. Our construction automatically adapts to the heteroskedasticity in model \eqref{eq:dgp_motivating_example} and is conditionally valid.
Figure \ref{fig:intro} provides an illustration. Panel (a) shows that the conditional length of $\mathcal{C}^{\rm reg}_{(0.9)}$ is constant, whereas the length of $\mathcal{C}^{\rm dcp}_{(0.9)}$ varies as a function of $x$. $\mathcal{C}^{\rm dcp}_{(0.9)}$ is shorter than $\mathcal{C}^{\rm reg}_{(0.9)}$ for low values and wider for high values of $x$. Panel (b) shows that $\mathcal{C}^{\rm dcp}_{(0.9)}$ is valid for all $x$, whereas $\mathcal{C}^{\rm reg}_{(0.9)}$ overcovers for low values and undercovers for high values of $x$. Figure \ref{fig:intro} illustrates the advantage of our method. For predictor values where the conditional variance is low, it yields shorter prediction intervals, while ensuring conditional coverage for values where the conditional dispersion is large by suitably enlarging the prediction interval.
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.4\textwidth,trim=0 1.25cm 0 1cm]{graphics/picture_intro_length.pdf}
\includegraphics[width=0.4\textwidth,trim=0 1.25cm 0 1cm]{graphics/picture_intro_coverage.pdf}
\caption{Motivating example}
\label{fig:intro}
\end{center}
\end{figure}
\subsection{Related Literature}
\label{sec:literature}
We build on and contribute to the literature on conformal prediction \citep[e.g.,][]{vovk2005algorithmic,vovk2012conditional,vovk2009online,lei2013distribution,lei2014distribution,lei2018distributionfree,chern2018COLT,romano2019conformalized}, the literature on model-free prediction \citep{politis2013modelfree,politis2015modelfree}, as well as the literature on quantile prediction methods \citep[see, e.g.,][for a review]{komunjer2013handbook}.
Within the conformal prediction literature, our paper is most closely related to \cite{lei2014distribution}, \cite{lei2018distributionfree}, and \cite{romano2019conformalized}. \cite{lei2014distribution} propose conditionally valid and asymptotically efficient conformal prediction intervals based on estimators of the conditional density. We take a different and complementary approach, allowing researchers to leverage powerful regression methods for modeling conditional distributions, including QR and DR approaches. \cite{lei2018distributionfree} develop conformal prediction methods based on regression models for conditional expectations. However, as discussed in Section \ref{sec:intro_conformal_prediction}, this approach is not conditionally valid under heteroskedasticity. They also propose a locally weighted conformal prediction approach, where the regression residuals are weighted by the inverse of a measure of their variability. This approach can alleviate some of the limitations of mean-based conformal prediction but is motivated by and based on restrictive locations-scale models. By contrast, our approach is generic and exploits flexible and substantially more general models for the whole conditional distribution.
\cite{romano2019conformalized} propose a split conformal approach based on QR models, which they call conformalized quantile regression (CQR). See also \cite{sesia2019comparison,kivaranovic2020adaptive} for related approaches and \cite{vovk2020conformal} for a general approach to adaptive conformal prediction. Their approach is based on splitting the data into two subsets, $\mathcal{T}_1$ and $\mathcal{T}_2$. Based on $\mathcal{T}_1$, they estimate two separate quantile functions $\hat{Q}(\alpha/2,x)$ and $\hat{Q}(1-\alpha/2,x)$ and construct the prediction intervals as
\[
\left[ \hat{Q}(\alpha/2,x)-Q_{E}(1-\alpha),\hat{Q}(1-\alpha/2,x)+Q_{E}(1-\alpha)\right],
\]
where $Q_{E}(1-\alpha)$ is the $(1-\alpha)(1+1/|\mathcal{T}_2|)$-th empirical quantile of $$E_t=\max\left\{\hat{Q}(\alpha/2,X_t)-Y_t,Y_t- \hat{Q}(1-\alpha/2,x)\right\}$$ in $\mathcal{T}_2$. Constructing prediction intervals based on deviations from quantile estimates is similar to working with deviations from mean estimates, as the deviations are measured in absolute levels. By contrast, exploiting the probability integral transform, our approach is generic and relies on permuting ranks, which naturally have the same scaling on $(0,1)$. Note, however, that our paper was inspired by \cite{romano2019conformalized} and we view our proposal as a (fully quantile-rank based) refinement of \cite{romano2019conformalized}.
Our adjustment for constructing efficient prediction intervals is related to and inspired by conformal prediction literature on minimum-volume prediction sets based on density estimators \citep[e.g.,][]{lei2013distribution,lei2014distribution,eck2019efficient,izbicki2019flexible,izbicki2020cd} and nearest-neighbor estimators \cite{gyorfi2020nearest}. It is most closely related and can be viewed as an alternative to conformal histogram regression \citep{sesia2021conformal}. The main differences between our approach and conformal histogram regression are the following. First, our method is based on an optimization problem formulated in terms of estimated quantile functions and does not require estimating a conditional density or histogram. Second, we do not work with nested sets but instead use a simple adjustment of our baseline conformity score. Finally, our approach works for general outcome distributions and does not rely on assuming unimodal distributions.
Conceptually, our paper is further related to the transformation-based model-free prediction approach developed in \cite{politis2013modelfree} and \cite{politis2015modelfree} in that we rely on transformations of the original setup into one that is easier to work with (i.e., ranks which are uniformly distributed) and study the properties of our approach in a model-free setting. An important difference is the implementation of the resulting procedure. The transformation-based approach is based on the bootstrap, whereas our approach is based on permuting ranks. Permuting ranks estimated based on the augmented data guarantees the model-free finite sample validity of our method with exchangeable data. To our knowledge, no exact finite-sample validity results have been developed for the bootstrap-based approach.
\section{Distributional Conformal Prediction}
Here we introduce DCP. We present a full and a split sample version of our method.
\subsection{Full Distributional Conformal Prediction}
Let $y$ denote a test value for $Y_{T+1}$. We test plausibility of each value $y \in \mathbb{R}$, collect all plausible values, and report them as the prediction set. In practice, we consider a grid of test values $\mathcal{Y}_{\rm trial}$.\footnote{For example, we can choose $\mathcal{Y}_{\rm trial}$ to be a fine grid between $-\max_{1\le t \le T} |Y_t|$ and $\max_{1\le t \le T} |Y_t|$. This choice has a theoretical justification since, under exchangeability, $P\left(|Y_{T+1}|\ge \max_{1\le t \le T} |Y_t|\right)\le 1/(1+T)$ \citep{chen2016trimmed}; see also the discussion in the \texttt{conformalInference} \texttt{R}-package (\url{https://github.com/ryantibs/conformal}).} Define the augmented data $Z^{(y)}=\{Z_t^{(y)}\}_{t=1}^{T+1}$, where
\begin{equation}\label{eq: def Z}
Z_t^{(y)}=\begin{cases}
(Y_t,X_t) & \textrm{if}\ 1\leq t\leq T \\
(y,X_t) & \textrm{if}\ t=T+1
\end{cases}
\end{equation}
Based on the augmented dataset $Z^{(y)}$, we estimate the conditional CDF using a suitable method such as QR and DR, which are discussed in more detail in the SI Appendix. Let $\hat{F}^{(y)}$ denote the estimator for $F$ based on the augmented sample. If the original estimate is not monotonic, we rearrange it \citep[e.g.,][]{chernozhukov2009improving,chernozhukov2010quantile} so that $\hat{F}^{(y)}(\cdot,x)$ is always monotonic. To simplify the exposition, we keep these rearrangements implicit.
We compute the ranks $\{\hat{U}^{(y)}_t\}_{t=1}^{T+1}$, where
\begin{equation}\label{eq:augmented_data}
\hat{U}^{(y)}_t=\begin{cases}
\hat{F}^{(y)}(Y_t,X_t)& \textrm{if}\ 1\leq t\leq T \\
\hat{F}^{(y)}(y, X_t)& \textrm{if}\ t=T+1
\end{cases}
\end{equation}
and obtain $p$-values as
\begin{equation}
\hat{p}(y)=\frac{1}{T+1}\sum_{t=1}^{T+1}\mathbf{1}\left\{\hat{V}^{(y)}_{t}\ge \hat{V}^{(y)}_{T+1} \right\}, \label{eq:p_value}
\end{equation}
where $\hat{V}^{(y)}_{t}:=\psi(\hat{U}^{(y)}_{t})$, and $\psi(\cdot)$ is a deterministic function. For our baseline method, we use $\psi(x)=|x-1/2|$.
In Section \ref{sec:extension}, we show how to choose $\psi$ optimally to ensure efficiency. Prediction intervals are computed as $\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})=\left\{y\in \mathcal{Y}_{\rm trial}:\hat{p}(y)>\alpha \right\}$.\footnote{Instead of $\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})$ we typically report the closed interval
$\widetilde{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})=\left[\min\left(\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})\right),\max \left( \widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})\right)\right]$.}
We summarize our approach in Algorithm \ref{algo:conformal}.
\begin{algo}[Full DCP]
\label{algo:conformal}
\text{ }
\vspace{-2mm}
\begin{itemize} \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item[] \textbf{Input:} Data $\left\{(Y_t,X_t)\right\}_{t=1}^T$, miscoverage level $\alpha \in (0,1)$, a point $X_{T+1}$, test values $\mathcal{Y}_{\rm trial}$
\item[] \textbf{Process:}
For $y\in \mathcal{Y}_{\rm trial}$,
\begin{enumerate} \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item define the augmented data $Z^{(y)}$ as in \eqref{eq:augmented_data}
\item compute $\hat{p}(y)$ as in \eqref{eq:p_value}
\end{enumerate}
\item[] \textbf{Output:} Return $(1-\alpha)$ prediction set $\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}(X_{T+1})=\left\{y\in \mathcal{Y}_{\rm trial}:\hat{p}(y)>\alpha \right\}$
\end{itemize}
\end{algo}
\subsection{Split Distributional Conformal Prediction}
An important drawback of full DCP (Algorithm \ref{algo:conformal}) is its computational burden due to the grid search. Since $\hat{F}^{(y)}$ is obtained based on the augmented data, one has to choose $\mathcal{Y}_{\rm trial}$ and re-estimate the entire conditional distribution for all $y\in \mathcal{Y}_{\rm trial}$. Therefore, we propose a split conformal procedure that exploits sample splitting, avoids grid search, and only requires estimating $F$ once. Sample splitting is a popular approach for improving the computational performance of conformal prediction methods \citep[e.g.,][]{lei2018distributionfree,romano2019conformalized}.
\begin{algo}[Split DCP]
\label{algo: split dcp}
\text{ }
\begin{itemize} \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item[] \textbf{Input:} Data $\left\{(Y_t,X_t)\right\}_{t=1}^T$, miscoverage level $\alpha \in (0,1)$, point $X_{T+1}$
\item[] \textbf{Process:}
\begin{enumerate} \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item Split $\{1,\dots,T\}$ into $\mathcal{T}_1:=\{1,\dots,T_0\}$ and $\mathcal{T}_2:=\{T_0+1,\dots,T\}$
\item Obtain $\hat{F}$ based on $\{Z_t\}_{t\in \mathcal{T}_1}$
\item Compute $\{\hat{V}_t\}_{t\in \mathcal{T}_2}=\{\psi(\hat{U}_t)\}_{t\in \mathcal{T}_2}$, where $\hat{U}_t=\hat{F}(Y_t,X_t)$.
\item Compute $\hat{Q}_{\mathcal{T}_2}$, the $(1-\alpha)(1+1/|\mathcal{T}_2|)$ empirical quantile of $\{\hat{V}_t\}_{t\in \mathcal{T}_2}$.
\end{enumerate}
\item[] \textbf{Output:} Return $(1-\alpha)$ prediction set $\widehat{\mathcal{C}^{\rm split}_{(1-\alpha)}}(X_{T+1})=\left\{y:\psi\left( \hat{F}(y,X_{T+1})\right)\le \hat{Q}_{\mathcal{T}_2}\right\}$.
(Since $\hat{F}(\cdot, X_{T+1})$ is monotonic, $\widehat{\mathcal{C}^{\rm split}_{(1-\alpha)}}(X_{T+1})$ is an interval.)
\end{itemize}
\end{algo}
In Algorithm \ref{algo: split dcp}, we split $\{1,\dots,T\}$ into $\{1,\dots,T_0\}$ and $\{T_0+1,\dots,T\}$. With iid data, one can also consider random splits.
Split DCP lends itself naturally to simple in-sample validity checks with both cross-sectional and time series data as illustrated in Section \ref{sec:applications}.
\section{Theoretical Performance Guarantees}
\label{sec:theory}
In this section, we establish the theoretical properties of our procedure. We focus on full-sample DCP (Algorithm \ref{algo:conformal}). For the split-sample approach (Algorithm \ref{algo: split dcp}), we provide a modified version (Algorithm \ref{algo: optimal}) in the SI Appendix and present its theoretical properties in Section \ref{sec:extension}.
When the data are iid\ (or exchangeable), our method achieves finite-sample unconditional validity in a model-free manner, as a consequence of general results on conformal inference and permutation inference more generally \citep[e.g.,][]{vovk2005algorithmic,hoeffding1952large}.
\begin{thm}[Finite sample unconditional validity]\label{thm:finite_sample} Suppose that the data are iid or exchangeable and that the estimator of the conditional distribution is invariant to permutations of the data. Then
\[
P\left(Y_{T+1}\in \widehat{\mathcal{C}_{(1-\alpha)}^{\rm full}}\left(X_{T+1} \right)\right)\geq 1-\alpha.
\]
\end{thm}
The proof of Theorem \ref{thm:finite_sample} is standard and omitted. Theorem \ref{thm:finite_sample} highlights the strengths and drawbacks of conformal prediction methods. Most commonly-used estimators of the conditional CDF such as QR and DR are invariant to permutations of the data. As a result, Theorem \ref{thm:finite_sample} provides a model-free unconditional performance guarantee in finite samples, allowing for arbitrary misspecification of the model of the conditional CDF. On the other hand, it has a major theoretical drawback. Even with iid data, it provides no guarantee at all on conditional validity.
Our next theoretical results provide a remedy. We impose the following weak regularity conditions.
\begin{assumption}
\label{assu: basic} Suppose that there exists a non-random function $F^*(\cdot,\cdot)$ such that the following conditions hold as $T\rightarrow \infty$. Define $V_t:=\psi(F^*(Y_t,X_t))$ for $1\leq t\leq T+1$.
\begin{enumerate}\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}
\item There exists a strictly increasing continuous function $\phi:[0,\infty)\rightarrow[0,\infty)$
such that $\phi(0)=0$ and $(T+1)^{-1}\sum_{t=1}^{T+1}\phi(|\hat{V}_{t}-V_{t}|)=o_{P}(1)$
and $\hat{V}_{T+1}=V_{T+1}+o_{P}(1)$, where $\hat{V}_{t}:=\hat{V}_{t}^{(Y_{T+1})}=\psi(\hat{F}^{(Y_{T+1})}(Y_t,X_t))$ for $1\leq t\leq T+1$.
\item $\sup_{v\in\mathbb{R}}|\tilde{G}(v)-G(v)|=o_{P}(1)$, where $\tilde{G}(v)=(T+1)^{-1}\sum_{t=1}^{T+1}\mathbf{1}\{V_{t}< v\}$ and $G(\cdot)$ is the distribution function of $V_{T+1}$.
\item $\sup_{x_{1}\neq x_{2}}|G(x_{1})-G(x_{2})|/|x_{1}-x_{2}|$ is bounded.
\end{enumerate}
\end{assumption}
Assumption \ref{assu: basic} allows for some flexibility with respect to the model estimator. Here, we only require $F^*$ to be a non-random function, which may or may not be $F$. The interpretation is straight-forward when $F^*=F$ since this simply means that the estimator $\hat{F}$ is consistent for $F$. We discuss the case of $F^*\neq F$ after Theorem \ref{thm: unconditional validity asy} below. Note that we can replace the consistency requirement in Assumption \ref{assu: basic} with a stronger uniform consistency requirement, $\sup_{x,y}|\hat{F}(y,x)-F^\ast(y,x)|=o_P(1)$.
We also notice that the quantities $\hat{V}_{t}$ and $V_{t}$ are defined under the true $Y_{T+1}$. This means that $\hat{F}^{(y)}$ uses $y=Y_{T+1}$. In other words, the estimator $\hat{F}$ based on the sample $\{(X_t,Y_t) \}_{t=1}^{T+1}$ would be consistent for some $F^*$ if $Y_{T+1}$ were observed.\footnote{This is not really much different from assuming that $\hat{F}$ based on the sample $\{(X_t,Y_t) \}_{t=1}^{T}$ is consistent for some $F^*$.} Since the goal of Assumption \ref{assu: basic} is to guarantee the coverage probability for $Y_{T+1}$, the conditions in Assumption \ref{assu: basic} only need to hold for $y=Y_{T+1}$.
Notice that $\hat{F}$ is consistent for $F^*$ under a very weak norm, and no rate condition is required.
When $\psi(x)=|x-1/2|$, a simple example of $\phi(\cdot) $ in Assumption \ref{assu: basic} is $\phi(x)=x^q $ for some $q>0$; in other words, a sufficient condition is $(T+1)^{-1}\sum_{t=1}^{T+1}|\hat{F}(Y_t,X_t)-F^\ast(Y_t,X_t)|^{q}=o_{P}(1)$, which can be verified for many existing estimators with $q=2$.
The following lemma gives the basic consistency result.
\begin{lem}
\label{lem: basic}Let Assumption \ref{assu: basic} hold. Then $\hat{G}(\hat{V}_{T+1})=G(V_{T+1})+o_{P}(1)$, where $\hat{G}(v)=(T+1)^{-1}\sum_{t=1}^{T+1}\mathbf{1}\{\hat{V}_{t}< v\}$.
\end{lem}
By Assumption \ref{assu: basic}, $G(\cdot)$ is uniformly continuous and thus continuous. Since $G(\cdot)$ is the distribution function of $V_{T+1}$, we have that $G(V_{T+1})$ has the uniform distribution on $(0,1)$, i.e., $P(G(V_{T+1})\leq\alpha)=\alpha$. This implies the unconditional asymptotic validity.
\begin{thm}[Asymptotic unconditional validity]
\label{thm: unconditional validity asy}Let Assumption \ref{assu: basic}
hold. Then
\[
P\left(Y_{T+1}\in\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}\left(X_{T+1} \right)\right)=1-\alpha+o(1).
\]
\end{thm}
Theorem \ref{thm: unconditional validity asy} establishes the asymptotic unconditional validity of the procedure. Since Theorem \ref{thm:finite_sample} already establishes the unconditional validity in finite-samples for iid or exchangeable data without assuming any consistency of $\hat{F}$, the main purpose of Theorem \ref{thm: unconditional validity asy} is to address the case of non-exchangeable data (e.g., time series data with ergodicity), especially when the model is misspecified (i.e., if $F^*\neq F$).
To illustrate model misspecification, consider the popular linear QR model, which assumes $Q(\tau,x)=x^{\top}\beta(\tau)$ and thus $F(y,x)=F(y,x;\beta)=\int_{0}^{1} \mathbf{1}\{x^{\top}\beta(\tau) \leq y \} d\tau$. This model is typically estimated by $\hat{\beta}(\tau)=\arg\min_{\beta}\sum_{t=1}^{T+1}\rho_{\tau}(Y_t-X_t^{\top} \beta)$ with $\rho_{\tau}(a)=a(\tau-\mathbf{1}\{a< 0\})$. Under misspecification ($Q(\tau,x)\neq x^{\top}\beta(\tau)$), $\hat{\beta}(\tau)$ is still estimating $\beta^*(\tau)=\arg\min_{\beta}\sum_{t=1}^{T+1}E\rho_{\tau}(Y_t-X_t^{\top} \beta)$ and $F^*$ is defined using $\beta^*(\cdot)$, e.g., $F^*(y,x)=\int_{0}^{1} \mathbf{1}\{x^{\top}\beta^*(\tau) \leq y \} d\tau$. For parametric models, $F^*$ is usually the probability limit of $\hat{F}$. In general, we can consider a model $\mathcal{F}$ and minimize the empirical risk $\hat{F}=\arg\min_{g\in\mathcal{F}} \sum_{t=1}^{T+1} L(Y_t,X_t,g) $ for some loss function $L$. Even if the model is misspecified ($F\notin \mathcal{F}$), it is still possible to show that $\hat{F}$ is close (in some norm) to $F^*=\arg\min_{g\in\mathcal{F}} \sum_{t=1}^{T+1}E[ L(Y_t,X_t,g) ]$. In the SI Appendix, we provide a more detailed discussion of this and some theoretical results verifying the consistency requirement in Assumption \ref{assu: basic} for the time series case; see also \cite{chern2018COLT} for a general discussion of conformal prediction in time series settings.
The cost of allowing for misspecification is that one cannot guarantee conditional validity when $F^*\neq F$. On the other hand, Lemma \ref{lem: basic} implies that the prediction intervals are conditionally valid when $F^*=F$.
\begin{thm}[Asymptotic conditional validity]
\label{thm: conditional validity asy}Let Assumption \ref{assu: basic}
hold with $F^*=F$.
Then
\[
P\left(Y_{T+1}\in\widehat{\mathcal{C}^{\rm full}_{(1-\alpha)}}\left( X_{T+1}\right)\mid X_{T+1}\right)=1-\alpha+o_{P}(1).
\]
\end{thm}
Theorems \ref{thm: unconditional validity asy}--\ref{thm: conditional validity asy} establish the asymptotic validity of our procedure under weak and easy-to-verify conditions. They formalize the key intuition that conditional validity hinges on the quality of the estimator $\hat{F}$ of the conditional CDF.\footnote{In Theorem \ref{thm: conditional validity asy}, we assume $F^*=F$. Since the first version of this paper was written, \cite{candes2021conformalized} have provided more general results where $F^*\approx F$.}
\section{Extension: Optimal DCP}
\label{sec:extension}
In Section \ref{sec:theory}, we have seen that a generic conformity score $\psi(y,x)=|F(y,x)-1/2|$ leads to conditional validity if the conditional distribution $F$ can be estimated consistently. We now characterize an optimal choice of conformity score that results in the shortest prediction interval. Detailed implementation algorithms, technical assumptions, and proofs are provided in the SI Appendix.
Let $\mathcal{Z}$ and $\mathcal{X}$ denote the support of $Z_t=(Y_t,X_t)$ and $X_t$, respectively. The optimal prediction interval is
\begin{equation}
\mathcal{C}^{{\rm opt}}_{(1-\alpha)}(x)=[r_{1}(x,\alpha),\ r_{2}(x,\alpha)],\label{eq: opt con int Y}
\end{equation}
where the functions $r_{1}(\cdot,\cdot),r_{2}(\cdot,\cdot)$ satisfy that for any $x\in\mathcal{X}$,
\begin{equation}\label{eq: key opt}
r_{2}(x,\alpha)-r_{1}(x,\alpha)=\underset{ F(z_2,x)-F(z_1,x)\geq1-\alpha}{\min}\ z_2-z_1.
\end{equation}
The question is whether it is possible to design a conformity score that achieves the above optimal prediction interval. To answer this question formally, we consider a generic conformity score $\psi(y,x)$, which might contain components that need to be estimated.
Permuting a large number of values of $\psi(Y_{t},X_{t})$ in conformal predictions amounts to taking the sample $(1-\alpha)$-quantile of $\psi(Y_{t},X_{t})$; for example, following Algorithm \ref{algo: split dcp}, one would $(1-\alpha)(1+1/|\mathcal{T}_2|)$ empirical quantile of $\psi(Y_{t},X_{t})$. Assuming a law of large numbers, this empirical quantile would be close to the population $(1-\alpha)$-quantile of $\psi(Y_{t},X_{t})$, leading to the asymptotic conformal prediction interval for $Y_{T+1}$
\begin{equation}
\mathcal{C}^{{\rm conf}}_{(1-\alpha)}(X_{T+1})=\{y:\ \psi(y,X_{T+1})\leq Q_\psi(1-\alpha)\},\label{eq: conf int}
\end{equation}
where $Q_\psi(1-\alpha)$ is the $(1-\alpha)$-quantile of $\psi(Y_{t},X_{t})$. The following result shows how to construct the optimal conformity score $\psi$.
\begin{lem}
\label{lem: opt score}
Let $\psi_*(y,x)=|F(y,x)-b(x,\alpha)-(1-\alpha)/2|$, where $b(\cdot,\cdot)$ is a function satisfying that for any $x\in\mathcal{X}$,
\begin{equation}\label{eq: b fun condition}
b(x,\alpha)\in\underset{z\in[0,\alpha]}{\arg\min}\ Q(z+1-\alpha,x)-Q(z,x).
\end{equation}
Let $\mathcal{C}^{{\rm conf}}_{(1-\alpha)}(X_{T+1})$ be defined as in \eqref{eq: conf int} with the above conformity score $\psi_*$.
Assume that $F(\cdot,x)$ is a continuous function for any $x\in\mathcal{X}$. Then $Q_\psi(1-\alpha)=(1-\alpha)/2$ and
$$
\mu\left(\mathcal{C}^{{\rm opt}}_{(1-\alpha)}(X_{T+1})\right)=\mu\left(\mathcal{C}^{{\rm conf}}_{(1-\alpha)}(X_{T+1})\right)\ {\rm almost\ surely},
$$
where $\mu(\cdot)$ denotes the Lebesgue measure.
If the optimization problem in \eqref{eq: opt con int Y} has a unique solution for any $x\in\mathcal{X}$, then
$$
\mathcal{C}^{{\rm opt}}_{(1-\alpha)}(X_{T+1})=\mathcal{C}^{{\rm conf}}_{(1-\alpha)}(X_{T+1})\ {\rm almost\ surely}.
$$
\end{lem}
Lemma \ref{lem: opt score} motivates conformity scores of the form $\psi_*(y,x)=\left|F(y,x)-[b(x,\alpha)+(1-\alpha)/2]\right|$, where $b(\cdot,\cdot)$ solves \eqref{eq: b fun condition}. Compared to the choice of $\psi(y,x)=|F(y,x)-1/2|$ mentioned in Section \ref{sec:theory}, we can view $\psi_*$ as having a ``shape'' adjustment $b(x,\alpha)-\alpha/2$. Since $F(Y_t,X_t)$ is independent of $X_t$, the optimal conformity score measures the distance between two independent components: $F(Y_t,X_t)$ and $1/2+(b(X_t,\alpha)-\alpha/2)$. Hence, by Lemma \ref{lem: opt score}, in order to take into account the shape of the conditional distribution $F(\cdot,x)$, it suffices to consider the scalar quantity $1/2+(b(x,\alpha)-\alpha/2)$.
In some special cases, the ``shape'' adjustment can be shown to be zero, i.e., $b(x,\alpha)=\alpha/2$. One typical example is when $F(\cdot,x)$ is a symmetric uni-modal distribution with a well-defined conditional density.\footnote{In this case, $Q(1/2+\delta,x)-Q(1/2,x)=Q(1/2,x)-Q(1/2-\delta,x)$ and the conditional density is increasing on $(-\infty, Q(1/2,x))$ and decreasing on $(Q(1/2,x),\infty)$. One can show $b(x,\alpha)=\alpha/2$ by taking the first-order derivative for the optimization problem in \eqref{eq: b fun condition} and setting it to zero.} Therefore, the choice of $\psi(y,x)=|F(y,x)-1/2|$ mentioned in Section \ref{sec:theory} is optimal in these cases. However, Lemma \ref{lem: opt score} provides a construction that achieves optimality more generally. By the definition of $\psi_*$ and $Q_\psi(1-\alpha)=(1-\alpha)/2$, the prediction interval is
\begin{equation}\label{eq: interval form conformal}
\mathcal{C}^{{\rm conf}}_{(1-\alpha)}(x)=[Q(b(x,\alpha),x), Q(b(x,\alpha)+1-\alpha,x)].
\end{equation}
We illustrate this in Figure \ref{fig:intro} with $\alpha=0.1$. \eqref{eq: interval form conformal} implies that $b(x,\alpha)$ is the quantile-index of the lower bound of the interval. For the symmetric distribution in the left panel, we see $b(x,\alpha)=0.05$, which is $\alpha/2$. For the asymmetric distribution in the right panel, we see that $b(x,\alpha)=0.007$, which is far away from $\alpha/2=0.05$.
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.4\textwidth,trim=0 1.5cm 0 1cm]{graphics/illu_optimal_normal.pdf}
\includegraphics[width=0.4\textwidth,trim=0 1.5cm 0 1cm]{graphics/illu_optimal_chi2.pdf}
\caption{Optimal prediction intervals}
\label{fig:intro}
\end{center}
\end{figure}
The first result in Lemma \ref{lem: opt score} is general and allows for the lack of uniqueness of the optimal prediction interval. For example, if $F$ is the uniform distribution on a certain interval, then all conditionally valid prediction intervals have the same length. Clearly, in this case, achieving the optimal length is the only goal one can hope for.
When we can uniquely define the optimal prediction interval, Lemma \ref{lem: opt score} implies that the conformal procedure can recover the uniquely defined optimal interval, not just achieving the optimal length.
Lemma \ref{lem: opt score} also confirms the insight of \cite{lei2014distribution}: the optimal confidence set for $X_{T+1}=x$ should take the form $ \{y:\ f(y,x)\geq c(x)\}$ for some $c(x)>0$, where $f(y,x)=\partial F(y,x)/ \partial y$. Assume that $F(\cdot,x)$ is a uni-modal distribution and $f(\cdot,x)$ is a continuous function for any $x\in\mathcal{X}$. Then this confidence set is an interval. This means that $ \{y:\ f(y,x)\geq c(x)\}=[c_1(x),c_2(x)]$ and $f(c_1(x),x)=f(c_2(x),x)=c(x)$. We notice that $c_1(x),c_2(x)$ are related to our results in that $c_1(x)=Q(b(x,\alpha),x)$ and $c_2(x)=Q(b(x,\alpha)+1-\alpha,x)$. To see this, simply observe that the first-order condition of the optimization problem in \eqref{eq: b fun condition} is
$1/f(Q(z+1-\alpha,x),x)-1/f(Q(z,x))=0 $, which implies that
$$
f(Q(b(x,\alpha)+1-\alpha,x))=f(Q(b(x,\alpha),x)).
$$
To make the procedure operational, we provide the conformal prediction interval $\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})$ defined
in Algorithm \ref{algo: optimal} in the SI Appendix. We can provide the following guarantee.
\begin{thm}\label{thm: effigency gen 2}
Let Assumption \ref{assu: efficiency gen} in the SI Appendix hold. Then
\[
P\left(Y_{T+1}\in\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})\mid X_{T+1}\right)=1-\alpha+o_{P}(1)
\]
and
\[
\mu\left(\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})\right)\leq\mu\left(\mathcal{C}_{(1-\alpha)}^{{\rm opt}}(X_{T+1})\right)+o_{P}(1).
\]
\end{thm}
The main requirements in Assumption \ref{assu: efficiency gen} in the SI Appendix are consistency of $\hat F$ and that the density $f$ bounded below on its support. This is quite mild in the sense that it does not imply that the optimal prediction interval in \eqref{eq: opt con int Y} is uniquely defined. For example, it allows $f$ to be a uniform distribution. Therefore, as discussed above, the conformal prediction interval would have approximately the shortest length but might not converge to $\mathcal{C}^{{\rm opt}}_{(1-\alpha)}(X_{T+1})$ in \eqref{eq: opt con int Y}.
The following theorem provides a stronger result about $\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})$ based on stronger assumptions.
\begin{thm}
\label{thm: efficiency}Let Assumption \ref{assu: efficiency } in the SI Appendix hold.
Consider the conformal prediction interval $\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})$ defined
in Algorithm \ref{algo: optimal} in the SI Appendix. Then
\[
\mu\left(\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})\triangle\mathcal{C}_{(1-\alpha)}^{{\rm opt}}(X_{T+1})\right)=o_{P}(1),
\]
where $\triangle$ denotes the symmetric difference of sets (i.e.,
$A\triangle B=(A\backslash B)\bigcup(B\backslash A)$), $\mathcal{C}_{(1-\alpha)}^{{\rm opt}}(X_{T+1})$
is defined in \eqref{eq: opt con int Y}.
\end{thm}
The key component of Assumption \ref{assu: efficiency } in the SI Appendix is consistent estimation of $b$. Theorem \ref{thm: efficiency} shows that $\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})$ is close to $\mathcal{C}_{(1-\alpha)}^{{\rm opt}}(X_{T+1})$ in the sense that the symmetric difference between these two sets has vanishing Lebesgue measure.
\section{Empirical Applications}
\label{sec:applications}
We illustrate the performance of DCP in two empirical applications and provide a comparison to alternative approaches. We consider eight different conformal prediction methods.
\begin{enumerate}\setlength\itemsep{0pt}
\item \textbf{DCP-QR:} DCP with QR (Algorithm \ref{algo: split dcp})
\item \textbf{DCP-QR$^\ast$:} Optimal DCP with QR (Algorithm \ref{algo: optimal} in SI Appendix)
\item \textbf{DCP-DR:} DCP with DR (Algorithm \ref{algo: split dcp})
\item \textbf{CQR:} CQR with QR \citep{romano2019conformalized}
\item \textbf{CQR-m:} CQR variant \citep{sesia2019comparison,kivaranovic2020adaptive} with QR
\item \textbf{CQR-r:} CQR variant \citep{sesia2019comparison} with QR.
\item \textbf{CP-OLS:} Mean-based split conformal prediction with OLS
\item \textbf{CP-loc:} Locally-weighted conformal prediction \citep{lei2018distributionfree} with OLS
\end{enumerate}
All computations were carried out in \texttt{R} \citep{R2021}. Code and data for replicating the empirical results are deposited on Github (\url{https://github.com/kwuthrich/Replication_DCP}).
\subsection{Predicting Stock Market Returns}
\label{sec:predicting_stock_returns}
Here we consider the problem of predicting stock market returns, which are known to exhibit substantial heteroskedasticity; see Chapter 13 in \cite{elliott2016economic} for a recent review and the references therein. We use data on daily returns of the market portfolio (CRSP value-weighted portfolio) from July 1, 1926, to June 30, 2021.\footnote{The CRSP data are constructed from the Fama/French 3 Factors data \citep{data_returns} available from \href{https://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html}{Kenneth R. French's data library} (accessed August 17, 2021).} We use lagged realized volatility $X_t$ to predict the present return $Y_t$.\footnote{We compute realized volatility as the square root of the sum of squared returns over the last 22 days.} Daily returns are not iid and exhibit time series dependence. In the SI Appendix, we show that the key conditions underlying our theoretical results hold when the data are $\beta$-mixing. Several stochastic volatility models for asset returns, including the popular GARCH models, can be shown to be $\beta$-mixing \citep[e.g.,][]{boussama1998ergodicite,carrasco2002mixing,francq2006mixing}.
We evaluate the performance of the different methods by splitting the data into a holdout and a test sample. To account for the dependence in the data, we present results averaged over five consecutive prediction exercises. In the first exercise, we apply split conformal prediction with an equal split ($|\mathcal{T}_1|=|\mathcal{T}_2|$) to the first 50\% of observations and use the next 10\% for testing. In the second exercise, we drop the first 10\% of the observations, apply split conformal prediction to the next 50\% of observations, and use the next 10\% for testing and so on.
Figure \ref{fig:conditional_coverage_returns} plots the empirical coverage probabilities for 20 bins obtained by dividing up the support of $X_t$ based on equally spaced quantiles. DCP-QR and DCP-QR$^\ast$ yield prediction intervals with coverage levels that are almost constant across all bins and close to the nominal level. They outperform DCP-DR, which undercovers in high-volatility regimes. The conditional coverage properties of DCP-QR and DCP-QR$^\ast$ are very similar to CQR, CQR-m, CQR-r, and CP-loc. This suggest that location-scale models, which are nested by QR, provide a good approximation of the conditional distribution. CP-OLS exhibits overcoverage under low-volatility regimes and substantial undercoverage under high-volatility regimes. This finding has important practical implications since the volatility tends to be high during periods of crisis, which is precisely when accurate risk assessments are most needed.
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.8\textwidth,trim=0 1.5cm 0 0.5cm]{graphics/conditional_coverage_returns.pdf}
\caption{Conditional coverage 90\% prediction intervals by realized volatility}
\label{fig:conditional_coverage_returns}
\end{center}
\end{figure}
Figure \ref{fig:conditional_length_returns} shows the conditional length of the prediction intervals. DCP-QR, DCP-QR$^\ast$, CQR, CQR-m, CQR-r, and CP-loc yield prediction intervals of similar length.
The DCP-DR prediction intervals are somewhat shorter than those of the QR-based methods at the upper tail. Finally, CP-OLS yields prediction intervals that are almost constant across all values of realized volatility; they are longer at the lower tail and shorter at the upper tail.\footnote{The CP-OLS prediction intervals are not exactly constant because we are reporting results averaged over five experiments.}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.8\textwidth,trim=0 1.5cm 0 0.5cm]{graphics/conditional_length_returns.pdf}
\caption{Conditional coverage 90\% prediction intervals by realized volatility}
\label{fig:conditional_length_returns}
\end{center}
\end{figure}
\subsection{Predicting Wages Using CPS Data}
\label{sec:predicting_wages}
We consider the problem of predicting wages using individual characteristics. We use the 2012 CPS data provided in the \texttt{R}-package \texttt{hdm} \citep{hdm2016}, which contains information on $N=29217$ observations. Here we use the index $i$ instead of $t$. To illustrate the impact of skewness on the performance of the different prediction methods, we use the hourly wage as our dependent variable $Y_i$.\footnote{We obtain the hourly wage by exponentiating the log hourly wage provided in the dataset.} Predictors $X_i$ include indicators for gender, marital status, educational attainment, region, experience, experience squared, and all two-way interactions such that $\dim(X_i)=100$ after removing constant variables.
Following \cite{romano2019conformalized} and \cite{sesia2019comparison}, we evaluate the performance of the different methods by randomly holding out 20\% of the data for testing, $\mathcal{I}_{\rm test}$, and applying split conformal prediction with an equal split to the remaining 80\% of the data. We repeat the whole experiment 20 times.
Panel (a) of Table \ref{tab:coverage} shows that all conformal prediction methods exhibit excellent unconditional coverage properties, confirming the theoretical finite sample guarantees. To assess and compare the conditional coverage properties, for each method, we compute conditional coverage probabilities as the predictions from logistic regressions of $\{Y_i\in \widehat{\mathcal{C}^{\rm split}_{(1-\alpha)}}(X_i)\}_{i\in \mathcal{I}_{\rm test}}$ on $\{X_i\}_{i\in \mathcal{I}_{\rm test}}$, where
$\widehat{\mathcal{C}^{\rm split}_{(1-\alpha)}}$ is the split conformal prediction interval obtained by the corresponding method. The less dispersed the predicted coverage probabilities are around the nominal level $1-\alpha=0.9$, the better the overall conditional coverage properties of a method. Panel (b) of Table \ref{tab:coverage} plots the standard deviation of the predicted coverage probabilities.\footnote{Using $\sqrt{1/| \mathcal{I}_{\rm test}|\sum_{i \in \mathcal{I}_{\rm test}}(\widehat{\texttt{Coverage}}_i-0.9)^2}$, where $\widehat{\texttt{Coverage}}_i$ is the predicted coverage probability, instead of the standard deviation yields very similar results.} DCP-QR$^\ast$ yields the lowest dispersion of all methods. The predicted coverage probabilities based on DCP-QR are less dispersed than those obtained from CQR, CQR-m, CQR-r. CP-loc yields a higher dispersion than the methods based on QR and DR, which demonstrates the value-added of using flexible models of the conditional distribution. Overall, DCP performs much better than CP-OLS for which the predicted coverage probabilities exhibit a very high dispersion. Figure \ref{fig:conditional_coverage_cps} in the SI Appendix plots histograms of the predicted coverage probabilities.
Table \ref{tab:length} shows the average length of the prediction intervals. DCP-QR$^\ast$ produces the shortest prediction intervals among of all methods. This demonstrates the practical advantage of the shape adjustment when the conditional distribution is skewed. The results also suggest a trade-off between conditional coverage accuracy and average length. For example, CP-OLS and CP-loc, which both exhibit poor conditional coverage properties, yield shorter prediction intervals than DCP-QR.
\begin{table}[ht]
\centering
\caption{Coverage 90\% prediction intervals}
\setlength{\tabcolsep}{2pt}
\begin{tabular}{cccccccc}
\toprule
\midrule
DCP-QR&DCP-QR$^\ast$&DCP-DR&CQR&CQR-m&CQR-r&CP-OLS&CP-loc\\
\midrule
\multicolumn{8}{c}{(a) Unconditional coverage} \\
\cmidrule(l{5pt}r{5pt}){1-8} \
0.90 & 0.90 & 0.90 & 0.90 & 0.90 & 0.90 & 0.90 & 0.90 \\
\midrule
\multicolumn{8}{c}{(b) Dispersion of predicted conditional coverage ($\times 100$)} \\
\cmidrule(l{5pt}r{5pt}){1-8} \
1.80 & 1.71 & 3.08 & 2.21 & 2.36 & 2.30 & 11.13 & 4.11 \\
\midrule
\bottomrule
\end{tabular}
\label{tab:coverage}
\end{table}
\begin{table}[ht]
\centering
\caption{Average length 90\% prediction intervals}
\setlength{\tabcolsep}{2pt}
\begin{tabular}{cccccccc}
\toprule
\midrule
DCP-QR&DCP-QR$^\ast$&DCP-DR&CQR&CQR-m&CQR-r&CP-OLS&CP-loc\\
\midrule
34.22 & 29.61 & 33.69 & 34.52 & 34.84 & 34.63 & 33.84 & 32.66 \\
\midrule
\bottomrule
\end{tabular}
\label{tab:length}
\end{table}
\bibliographystyle{apalike}
\bibliography{pnas-biblio}
\newpage