Distributional conformal prediction

Distributional conformal prediction

\linespread{1.15} \parskip 0in

\linespread{1.25} \parskip .05in

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) koenker1978, distribution regression (DR) foresi1995conditional,chernozhukov2013inference, as well as non-parametric and high-dimensional machine learning methods such as quantile neural networks taylor2000quantile and quantile trees and random forests 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 cattaneo2019prediction,chernozhukov2021conformal and to the problem of predicting individual treatment effects 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

for a given miscoverage level $\alpha\in (0,1)$. A prediction interval satisfying this property is said to be 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 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 romano2019malice. Second, as in Section (ref), 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), 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 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:

A prediction interval satisfying property (ref) as $T\rightarrow \infty$ is said to be approximately conditionally valid.\footnote{See, for example, lei2014distribution,sesia2019comparison,foygel2021limits for a further discussion of the difference between conditional and unconditional validity.}

While the requirement in (ref) 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 vovk2012conditional for a detailed discussion. By Proposition 2 of 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 vovk2012conditional. In 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 (ref) cannot be achieved in a distribution-free way for non-trivial predictions. By {vovk2012conditional,lei2014distribution,foygel2021limits}, any prediction set satisfying (ref) for every probability distribution of $(X_t,Y_t)$ has infinite Lebesgue measure with non-trivial probability. Therefore, we only aim to achieve (ref) 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 (ref) 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), 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), 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

where $Q(\tau, x)$ is the $\tau$-quantile of $Y_t$ given $X_t=x$. To implement the prediction interval (ref), a plug-in approach would replace $Q$ with a consistent estimator $\hat{Q}$

This approach exhibits two well-known drawbacks. First, it will often exhibit undercoverage in finite samples romano2019conformalized. Second, it is neither conditionally nor unconditionally valid under misspecification.

We build upon conformal prediction 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 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 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 (ref). 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 (ref). 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 optimal DCP. In Section (ref), 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.

Motivating Example

We illustrate the advantages of DCP relative to mean-based conformal prediction lei2018distributionfree based on the following simple analytical example.

Our motivating example draws on 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

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$ lei2018distributionfree,romano2019conformalized. This feature implies that $\mathcal{C}^{\rm reg}_{(1-\alpha)}$ is not adaptive to the heteroskedasticity in the location-scale model (ref) 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

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 (ref) and is conditionally valid.

Figure (ref) 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) 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.

Related Literature

We build on and contribute to the literature on conformal prediction vovk2005algorithmic,vovk2012conditional,vovk2009online,lei2013distribution,lei2014distribution,lei2018distributionfree,chern2018COLT,romano2019conformalized, the literature on model-free prediction politis2013modelfree,politis2015modelfree, as well as the literature on quantile prediction methods komunjer2013handbook.

Within the conformal prediction literature, our paper is most closely related to lei2014distribution, lei2018distributionfree, and romano2019conformalized. 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. lei2018distributionfree develop conformal prediction methods based on regression models for conditional expectations. However, as discussed in Section (ref), 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.

romano2019conformalized propose a split conformal approach based on QR models, which they call conformalized quantile regression (CQR). See also sesia2019comparison,kivaranovic2020adaptive for related approaches and 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 romano2019conformalized and we view our proposal as a (fully quantile-rank based) refinement of 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 lei2013distribution,lei2014distribution,eck2019efficient,izbicki2019flexible,izbicki2020cd and nearest-neighbor estimators gyorfi2020nearest. It is most closely related and can be viewed as an alternative to conformal histogram regression 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 politis2013modelfree and 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.

Distributional Conformal Prediction

Here we introduce DCP. We present a full and a split sample version of our method.

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)$ chen2016trimmed; see also the discussion in the conformalInference R-package (\url{https://github.com/ryantibs/conformal}).} Define the augmented data $Z^{(y)}=\{Z_t^{(y)}\}_{t=1}^{T+1}$, where

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 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

and obtain $p$-values as

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), 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).

Split Distributional Conformal Prediction

An important drawback of full DCP (Algorithm (ref)) 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 lei2018distributionfree,romano2019conformalized.

In Algorithm (ref), 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).

Theoretical Performance Guarantees

In this section, we establish the theoretical properties of our procedure. We focus on full-sample DCP (Algorithm (ref)). For the split-sample approach (Algorithm (ref)), we provide a modified version (Algorithm (ref)) in the SI Appendix and present its theoretical properties in Section (ref).

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 vovk2005algorithmic,hoeffding1952large.

The proof of Theorem (ref) is standard and omitted. Theorem (ref) 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) 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.

Assumption (ref) 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) below. Note that we can replace the consistency requirement in Assumption (ref) 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) is to guarantee the coverage probability for $Y_{T+1}$, the conditions in Assumption (ref) 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) 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.

By Assumption (ref), $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.

Theorem (ref) establishes the asymptotic unconditional validity of the procedure. Since Theorem (ref) 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) 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) for the time series case; see also 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) implies that the prediction intervals are conditionally valid when $F^*=F$.

Theorems (ref)--(ref) 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), we assume $F^*=F$. Since the first version of this paper was written, candes2021conformalized have provided more general results where $F^*\approx F$.}

Extension: Optimal DCP

In Section (ref), 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

where the functions $r_{1}(\cdot,\cdot),r_{2}(\cdot,\cdot)$ satisfy that for any $x\in\mathcal{X}$,

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), 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}$

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$.

Lemma (ref) 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 (ref). Compared to the choice of $\psi(y,x)=|F(y,x)-1/2|$ mentioned in Section (ref), 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), 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 (ref) and setting it to zero.} Therefore, the choice of $\psi(y,x)=|F(y,x)-1/2|$ mentioned in Section (ref) is optimal in these cases. However, Lemma (ref) 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

We illustrate this in Figure (ref) with $\alpha=0.1$. (ref) 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$.

The first result in Lemma (ref) 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) implies that the conformal procedure can recover the uniquely defined optimal interval, not just achieving the optimal length.

Lemma (ref) also confirms the insight of 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 (ref) 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) in the SI Appendix. We can provide the following guarantee.

The main requirements in Assumption (ref) 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 (ref) 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 (ref).

The following theorem provides a stronger result about $\widehat{\mathcal{C}_{(1-\alpha)}^{{\rm conf}}}(X_{T+1})$ based on stronger assumptions.

The key component of Assumption (ref) in the SI Appendix is consistent estimation of $b$. Theorem (ref) 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.

Empirical 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.

All computations were carried out in R R2021. Code and data for replicating the empirical results are deposited on Github (\url{https://github.com/kwuthrich/Replication_DCP}).

Predicting Stock Market Returns

Here we consider the problem of predicting stock market returns, which are known to exhibit substantial heteroskedasticity; see Chapter 13 in 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 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 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) 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.

Figure (ref) 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.}

Predicting Wages Using CPS Data

We consider the problem of predicting wages using individual characteristics. We use the 2012 CPS data provided in the R-package hdm 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 romano2019conformalized and 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) 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) 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) in the SI Appendix plots histograms of the predicted coverage probabilities.

Table (ref) 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.