Thanks a lot for reading our paper and for your insightful comments.

1. Lack of intuitive understanding: Thank you for your concern. To achieve a minimal-width prediction interval, it is crucial to accurately capture the shape of the interval. However, a single predictor may fail to do so due to its specific structural limitations. For example, the correct shape might be a polynomial function, but the predictor could output a linear function, resulting in model misspecification. By aggregating several predictors and considering a convex hull of these base predictors, we reduce model misspecification. This approach leads to an ensemble predictor that more accurately captures the interval's shape, resulting in a smaller width. The second step (shrinkage) is equivalent to conformal prediction, which maintains the coverage guarantee.

2. Comparison with Hosen et al. (2014): Thank you for pointing out this interesting paper. We will include a comparison with the previous methods from Hosen et al. (2014) in our revised version. Hosen et al. (2014) focus on ensemble prediction intervals obtained using neural networks through ranking and weighted averaging. While they consider both coverage and width, they combine these metrics into a single optimization objective (coverage width criterion). In contrast, our method focuses on finding the interval with the smallest width among those that provide sufficient coverage. Additionally, Hosen et al. (2014) do not provide theoretical guarantees for coverage and width, nor do they address domain shift scenarios. Our work, on the other hand, emphasizes providing a theoretically robust methodology specifically designed to handle domain shift challenges.

3. Limited experiments: Thank you for raising this point. We have now conducted three additional real-data experiments to evaluate our method. We have added one more baseline, the weighted quantile conformal prediction method. Please see the global response for details on these experiments. Last but not least, we have also conducted a simulation experiment quantifying the robustness of our method against conformal prediction. Please see our response to Reviewer 6ToN for details of that experiment.

4. What does $f \in \\mathcal{F}$ signify?: Sorry for the typo, it should be $\min_{l, u}$.

5. What does $m_0$ represent?: $m_0$ is the conditional mean function, as defined in line 137. We will clarify it in our revision.

6. Why not apply hinge to $\hat w(x)$?: The purpose of using a hinge function is to ensure the convexity of the problem. Specifically, we replace the indicator function for $(Y - m_0(x))^2 - f(x)$ with a convex surrogate (hinge function) to facilitate the use of standard convex optimization techniques.

7. Notation $\\mathcal{B}_{\mathcal{F}}$: $\\mathcal{B}\_{\mathcal{F}}$ is defined to be an upper bound of $\sup_{f\in\mathcal{F}}\|f\|_{\infty}$ (as shown in line 181 in Theorem 3.2). We will make it more clear in our revision.

8. Experiment on Algorithm 2: We have conducted an experiment on the Airfoil data using our second algorithm, which requires estimating the optimal transport map. Please see our global response for the details.

Thanks a lot for reading our paper and for your insightful comments.

1. The paper assumes familiarity with Fan et al. (2023) Thank you for raising this concern! We will carefully proofread the revised version of the manuscript according to your suggestions.

2. Mention optimal aggregation Thank you for pointing this out. Indeed, our current approach presents the general problem and then casts prediction interval aggregation as a special case. We will highlight the aggregation perspective more in the revised version of the manuscript.

3. Similarities and dissimilarities to Fan et al. (2023) Please see our global response. We selected the quantile levels of 0.85, 0.95, 0.9, and 0.9 because they are close to the desired coverage and we believe they effectively capture the correct shape of the interval. It is important to note that our method is not restricted to these specific quantiles; we can also apply the quantile levels used in Fan (2023) to our approach if needed.

4. Change of shape of prediction band Thank you for your feedback. We will clarify this in line 147 in the revision. Specifically, when we refer to the "the shape of the prediction band does not change much", we mean that the form of the minimal average width prediction band covering the data at a certain coverage level should not change much if we slightly change the coverage level, i.e., the form is stable with respect to the perturbation of the coverage level. For instance, if the optimal prediction band that covers $98\\%$ of the data is quadratic in $X$, then the minimal-width interval that covers $95\\%$ of the data is also a quadratic function of $X$, rather than abruptly changing to a linear function. Under this consistency assumption, a constant shrinkage factor $\lambda$ in the second step (i.e., (3.3) and (4.2)) is sufficient to adjust from $100\\%$ to $95\\%$ coverage. However, if this is not the case, we can modify our method to allow the shrinkage level $\lambda$ to depend on the input $X$, ensuring more accurate adaptation to the data distribution. For the multimodal distributions, we acknowledge that the optimal region is not necessarily a prediction band (i.e., $Y|X=x$ is not necessarily a single interval) and will include this in the limitation of the work.

5. Discussion on $\lambda \le 1$ Thank you for the comments. We will clarify this discussion in the revised version of the paper. In the first step, we typically choose $\epsilon$ to be small; therefore, the average width of $\hat f_{\rm init}$ can be larger than necessary, and coverage is also typically larger than $95\\%$. Therefore, in the second step, we aim to shrink the interval (i.e., ensure $\lambda \leq 1$). However, it is not immediate that $\hat \lambda$, obtained by solving equation (3.3) is indeed $\le 1$. To ensure this, we prove in Lemma 3.3 and Lemma 3.4 that with high probability, we indeed have $\hat \lambda(\alpha) \le 1$.

6. Research on prediction interval estimation Thank you for pointing this out. In our revision, we will add a detailed related work section that covers significant research on prediction interval estimation, including methods applicable to both shifted and non-shifted data. We will discuss approaches such as conformal prediction, quantile regression, and RKHS-based methods, highlighting their contributions and relevance to our work.

7. Limitations of the paper are not discussed Thank you very much for highlighting this limitation! It is indeed true that if the conditional distribution of $Y$ given $X$ is multimodal, the prediction interval may not be the best thing to do. We will mention this in the limitation.

8. Comparison with conformal

(1) Methodology relation and difference: Step 2 of Algorithm 1, i.e., shrinkage (equation (3.3) and (4.2)) in our method is equivalent to (weighted) conformal prediction, using the score function $(y - \hat{m}(x))/\hat{f}(x)$, where $\hat{f}$ is obtained from the first step of our method. By efficiently aggregating various predictors, our $\hat{f}$ accurately captures the shape of the interval, resulting in a smaller bandwidth. In contrast, variance-adjusted conformal prediction selects $f$ as the variance function, which may not necessarily capture an optimal shape for a small bandwidth. Additionally, our method explicitly aims to minimize the bandwidth in the first step. In contrast, conformal prediction does not explicitly minimize bandwidth; it implicitly addresses this by setting the cutoff at the $1-\alpha$ quantile of the score function.

(2)Robustness: Our method is more robust compared to variance-adjusted conformal prediction. Please see our response to Reviewer 6ToN for the details of our experiment regarding robustness.

(3)Theoretical guarantees: We provide finite sample guarantees for both the width and coverage of our prediction intervals. In contrast, while typical conformal literature offers finite-sample marginal or asymptotic conditional guarantees for coverage, they often do not include finite-sample guarantees for width, as our paper does. For distribution shift [1] provides coverage guarantee assuming the density ratio $w(\cdot)$ is known, whereas our Theorem 3.4 explicitly highlights the effect of the estimation error of $w(\cdot)$.
We acknowledge that our method is particularly focused on addressing marginal coverage, and therefore, its pointwise guarantee is not immediately clear. We will include this as a limitation in our discussion and consider it a topic for future research.

[1] Tibshirani, R., et al., Conformal prediction under covariate shift.

9. Typo Thank you very much for carefully reading our paper and pointing out the typos! We will take care of it in the revised version of the manuscript.

Thanks a lot for reading our paper and for your insightful comments.

1. Extension of empirical evaluation: Thank you for your comment, please see our global response for the details of the additional experiments.

(i). We now have two more real-data experiments on our first method, Algorithm 1(which relies on the estimation of the density ratio of the source and the target covariates), one using the real estate valuation data and the other using the energy efficiency data. Furthermore, we have an additional experiment on our second method, Algorithm 2 (which relies on estimating the optimal transport map), using the Airfoil data.

(ii). We have now added one more baseline, namely the weighted quantile-adjusted conformal method, i.e. we now compare our against both variance-adjusted weighted conformal and weighted quantile-adjusted conformal method.

(iii). We have also conducted a simulation experiment to showcase the robustness of our algorithm in comparison to weighted conformal prediction. Kindly see our response to Reviewer 6ToN for the details of this experiment.

2. Figure 9: Thank you for pointing out this problem. We will take care of it in the revised version of the paper.

(I will quickly write a reply now, to try to give you some time to respond before the deadline, but I will also think more carefully about this later before making my final recommendation. I will still consider raising my score)

(notably, the AssetWealth dataset you referenced is 13GB)

• There are 8 different datasets in the TMLR paper i pointed to, I thought that maybe it could be easy to apply to at least one of those?

...suggesting its significant potential impact across diverse real-world scenarios

• The method has only been applied to a single dataset with a real-world distribution shift though, right? Are there some other tabular datasets with realistic distribution shifts?

These experiments consistently demonstrate the advantages of our method over WVAC and WQC

• Expect for "Airfoil data with Algorithm 2" where the coverage drops below 0.95? Also, is this the only evaluation of Algorithm 2? Why? (could perhaps be made more clear in general what the two different algorithms should be utilized for, how do I know which one to use in a real-world application?)

...performing dimensionality reduction. Once these preprocessing steps are completed and we obtain vector representations of the images, our method can be implemented efficiently and executed rapidly

• Why is dimensionality reduction a required step? Isn't this a quite significant limitation of your method then, which would limit the real-world utility also for non-image data? No dataset in the experiments has more than 8(?) features, what's the maximum number your method can handle?

Thanks a lot for reading our paper and for your insightful comments. We answer your questions and concerns in the following.

1. Novelty compared to Fan et al. (2023): Please see our global response.

2. Assumptions on the covariate shift: Covariate shift is a common assumption, even for high-dimensional data, e.g., see [1] or [2]. The impact of domain shift in our analysis is primarily through the estimation error of the density ratio and the optimal transport map. We acknowledge that high-dimensional settings can pose challenges, such as the curse of dimensionality. However, these issues can often be mitigated by incorporating more structure; without any structural assumption, we do not expect to gain any theoretical guarantee under covariate shift in high-dimension. An example of a structured distribution shift is the exponential tilt assumption (the one we used in our experiments), or that only a subset of variables undergoes the shift (sparsity type assumption). Similarly, we can assume some structure in the optimal transport map to avoid the curse of dimensionality. While our analysis is quite general, we believe that with specific structural assumptions, our method and analysis can be appropriately adjusted to incorporate these structures. Therefore, we agree that we need some structural assumptions in high-dimension, but without any assumption, the curse of dimensionality is unavoidable. Once we have that, our method is applicable to obtain a valid prediction band with a small width.

[1] Tripuraneni, N., et al., Covariate shift in high-dimensional random feature regression.

[2] Tripuraneni, N., et al., Overparameterization improves robustness to covariate shift in high dimensions.

3. Similarity and dissimilarity with conformal prediction: We agree with you that there is a similarity between our method and the conformal prediction in the second stage, i.e., the shrinkage (i.e., (3.3) and (4.2)), which is equivalent to (weighted) conformal prediction as both of these methods aim to derive weighted quantile of the sample score $s(x, y) = (y - \hat{m}(x))/\hat{f}(x)$. However, the key difference lies in the construction of $\hat f$; by efficiently aggregating various predictors, $\hat f$ accurately captures the shape of the interval, resulting in a smaller bandwidth.

Robustness check

Another advantage of our method is its robustness. Our method aggregates various predictor intervals, including an estimator for the conditional variance function (Estimator 5). In contrast, the weighted variance-adjusted conformal (WVAC) prediction interval [3] heavily relies on accurately estimating this conditional variance, while conditional quantile regression relies on accurate estimation of the conditional quantile function. We did a small robustness check using the following model:

$

& X \sim{\sf Unif}([-1, 1]), \ \ \xi \sim{\sf Unif}([-1, 1]),\ \  X,\xi \text{ independent}, \\
& Y = \sqrt{1 + 25 X^4} \xi .

$

We estimated the conditional variance using a random forest with varying depths {3, 5, 7, 15}. We generated n = 2500 samples, keeping $75\\%$ as source data and resampling the remaining 25% with weighted samples proportional to $w(x) \propto (1 + \exp{(-2x)})^{-1}$. As depth increases, overfitting leads to poor out-of-sample variance predictions. The following table is the result of 100 Monte Carlo experiments. The number inside the parenthesis is the median of coverage over these Monte Carlo iterations.

The above table implies our method is more robust to the misspecification of some model components and remains stable.

In a nutshell, our method has a similar coverage guarantee as the conformal prediction, which does not require many assumptions; our coverage guarantee (Theorem 3.4) requires a reasonable estimation of $w(\cdot)$, which is necessary for both our method and that of [3]; [3] even assumed $w(\cdot)$ is known, whereas our bounds quantify the estimation error of $w(\cdot)$. When we require a guarantee on the width, we require some assumptions on the predictor class.

[3] Tibshirani, R., et al., Conformal prediction under covariate shift.

4. Extensive empirical evaluation Please see our global response.

As the rebuttal period is coming to a close, we wanted to check if there are any remaining concerns that might be preventing you from adjusting your score. If there is any additional clarification we can provide, please let us know. We would be more than happy to address any questions you might have. Thank you again for your time and effort in reviewing our paper.