Copula Conformal prediction for multi-step time series prediction

We thank the reviewer for the detailed review.

Our response to the weaknesses are:

1. Theoretical analyses: the i.i.d. assumption is imposed onHowever, this eliminates the dependency of data over time, which is typical and expected for time-series.

This is a major misunderstanding. The superscript i is the index of the data sample in the dataset; we are assuming each time series sample is independent, not the timesteps (denoted by underscript). Please refer to our notation in the second paragraph of section 4. Therefore, our method maintains the temporal dependency. Our figure 1 helps illustrate our setting of dependent time steps and independent data samples of whole trajectories. We understand that this could be confusing, so we added a parenthesis in the proof. Hope this helps to clear things out.

--

2. The simulation performance of CopulaCPTS seems very similar to Copula [Messoudi et al. 2021], which notably was not developed for this setting.

We respectfully disagree and argue that the difference in both the algorithm and empirical performance is significant. Note that in all of the Copula [Messoudi et al.] results in table 1, the coverage is invalid (below 90%), whereas all of the CopulaCPTS coverages are valid. Because the main appeal of conformal methods is the guaranteed validity, our algorithm brings significant theoretical and practical improvement.

Messoudi et al.’s work also assumes the multitarget data samples (equivalent to our $z^i$) are independent. You are right in that one can treat multi-step prediction as a multitarget problem and use [Messoudi et al.] directly, but I hope our strong experiment results in table 1 show the reason why we need to use CopulaCPTS instead of [Messoudi et al].

--

3. Performance of Copula on real-data examples is not reported.

Please see the last two rows of table 1 for quantitative results, figure 4 for qualitative examples, and figure 7,8, and 9 in the appendix for visualizations on real data.

--

4. [Messoudi et al. 2021] does not underperform much; Need more baselines.

In our response to weakness #2, we have shown that CopulaCPTS significantly improves upon [Messoudi et al. 2021] for our setting. We have added the Locally ellipsoidal CP as a baseline in table 7 in appendix C.7. The two methods achieve comparable results.

We thank the reviewer for your thoughtful and generous review.

To answer your questions:

Equation (4) different from prior work.

Apologies for this typo, the $\forall j$ should be inside the probability, as we are certifying the joint coverage of all timesteps. We have updated our draft for the expression both in equation (4) and theorem 4.1 in the new version.

--

In experiments, did you use the same training dataset for all methods or did you adapt it depending on whether the method need a calibration set?

Thank you for helping to clarify this! For baselines that do not require calibration, the calibration split is used for training the model. We have added this sentence in section 5.1.

We thank the reviewer for reading our paper carefully and writing a helpful review. Appreciate it!

In response to the weaknesses (W) and questions (Q):

W2: In the experiment, computational times are not given. The standard deviations of the experiments are made with only 3 runs of the algorithm. Maybe the algorithm is very time-consuming (this should be more discuss).

Thank you for pointing out this potential source of confusion! CopulaCPTS is not time-consuming to run. It does not grow exponentially with prediction horizon, because we search for $\mathbf{u}^*$ by SGD (appendix B.2). We ran the experiment for 3 runs because it is common practice in the machine learning literature, and that our runtime bottleneck was the baseline BJRNN. We have added this discussion paragraph in the appendix B.2 on SGD:

" … The optimization process to find $\mathbf{s}^*$ typically takes a few seconds on CPU. For each run of our experiments, the calibration and prediction steps of CopulaCPTS combined took less than 1 minute to run on an Apple M1 CPU. "

--

W3: typos

We have fixed the typos in our draft. Thank you for catching them!

--

Q1: Particular definition for the multivariate empirical quantile function (Eq. 13)

We constructed this vector empirical quantile function for the algorithm. Since comparisons between vectors are ambiguous, one cannot directly use a canonical empirical quantile function such as equation 10. Hence, we introduce the partial order and define a specific quantile function based on the partial order.

--

Q2: In Lemma A.2, how is it possible to have an equality? (Do scores not need to be continuous?).

The scores are continuous. The probability in lemma A.2 is a cumulative probability of the value of $\hat{F}_j(s_j)$, so the scores can be continuous and the equality can hold. Please let me know if this answers your question.

--

Q3: Isn't the definition of "exchangeability" in the paper rather a consequence of the "true" definition?

This definition of exchangeability is commonly used in papers on conformal prediction beginning from [Vovk et al. 2005]. We use this formulation as the “exchangeability assumption” in our paper for clarity and brevity.

--

Q4: is the "∀j" in the probability or outside ? In the proof, Eq. (20), this is inside but in Theorem 4.1, this is outside.

The $\forall j$ should be inside the probability as in Eq 20, as we are trying to certify the joint probability here. Thank you for catching this important typo, we have updated Theorem 4.1 in the new version.

--

Q5: In the experiments, the score is chosen to be an L2 norm (see, for example, step 9 of Algorithm 1). What are the implications of this choice on the results? For example, are the results very different if we use another norm?

The beauty of conformal prediction is that the validity guarantee is independent of the choice of nonconformity scores. Our proof is also general with regard to which score is used. Because we model the copula directly on the empirical CDF and hence the rank of nonconformity scores, our efficiency argument is also generalizable to any score choice. There are works exploring domain-specific nonconformity scores to improve efficiency, but we regard that as an orthogonal line of work to ours.

We thank the reviewer for providing helpful feedback.

In response to the weaknesses (W) and questions (Q):

Q1/W1: What distinguishes the proposed method from Messoudi's work, apart from the introduction of multivariate entries?

Our copula-based method is significantly different from Messoudi et al. Our contribution is threefold:

• the two-step approach (section 4.1): Messoudi’s algorithm uses the same dataset to calibrate the individual scores and the copula, which results in inaccuracy in the estimation. We improve the process by (1) using separate calibration datasets for the individual time steps, and (2) constructing conformal predictive distributions for each timesteps. Our method leads to un-confounded calibrations and allows us to prove validity.
• validity guarantee (Theorem 4.1 and appendix A): In their paper, Messoudi et al. only showed empirical results and didn’t provide theoretical support for their performance. As we have shown in table 1 in the paper, directly applying Copula often results in invalidity, which is a fatal flaw because guaranteed coverage is the main appeal of using conformal prediction. Our algorithm, on the other hand, is provably valid.
• New algorithm for finding $s^*$ (Eq. 8, and appendix B.2): Copula calibration in Messsoudi et al is implemented as a grid search (grows exponentially with the dimension of the copula) or by assuming that all $\alpha_j$ are the same (sacrifices efficiency). We mitigate this by using SGD for optimization, improving search time from exponential time to constant time.

In conclusion, our algorithm resolved the shortcomings of Messoudi et al (invalidity and slowness), and presented a practical, sound, and fast algorithm for the big-data time series setting that does not yet exist in the conformal prediction literature.

--

W1.1 While the authors assert that their contribution is the multivariate nature of every entry, it's important to note that the nonconformity score remains scalar. Therefore, the process of conformal prediction remains largely the same as in the univariate case.

Our algorithm is not the same as the univariate case. Our contribution lies in combining multiple conformal predictions (can be multivariate or not, yes, with a scaler nonconformity score) together while maintaining coverage for all and having good efficiency. We argue that the second part departs from canonical CP and is novel.

--

W2: data type. While this paper primarily focuses on time series, upon closer examination, it becomes apparent that the temporal aspect of the data may not be as critical as initially assumed….

It is true that our algorithm is not limited to time series data. However, our work is motivated by many practical problems in time series that (1) has many exchangeable time series sequences, and (2) requires a confidence region coverage guarantee over a long prediction horizon, as exemplified in our choice of real-world datasets, anonymized medical data and autonomous driving. In addition, the structure of time series allows us to expand our algorithm to the autoregressive prediction setting, as explored in Appendix B.4.

--

W2.2: Furthermore, it's crucial to differentiate this work from other studies that delve into conformal prediction methods that extend beyond exchangeability.

We differentiated our setting with the CP-beyond-exchangeability method explicitly in our related works section (paragraph 2 under conformal prediction in section 2).

--

Q2: In comparison to another work [1] focusing on the multi-output version of conformal prediction that uses quantiles and can be extended to the multivariate case with the same score function, what distinguishes and contributes to this paper? I also recommend adding [1] as another baseline in the experimental part to prove the efficiency improvement.

Feldman et al. ’s work [1] on multivariate conformal quantile regression can be used in conjunction with our copula methods, hence not a direct baseline. Specifically, we can use them for the individual timesteps where the dimensionality is lower, and then use copula to jointly calibrate the time steps, to further improve efficiency.

We have attempted to implement this idea during the rebuttal. However, [1] requires calculating a grid over the output space, which grows exponentially with the output dimension. The original paper only evaluated performance for data with response dimension fewer than 4. In a time series setting, however, it is usual to have long prediction horizons. Running [1]’s algorithm on our particle simulation with output dimension = 50 (25 time steps x 2d output) requires 1048576.00 GiB of GPU memory.

Q3: In Equation 8, when aiming to enhance efficiency by minimizing the L1 norm of u, what is the rationale for choosing this particular norm?

Our algorithms holds as long as the constraint of $C_{empirical}(\mathbf{u}*) \geq 1-\alpha$ holds, and this constraint is independent of the minimization goal. Therefore, the choice of this loss function is not very critical. In our experiments, we found no significant difference between using L1 or L2 norms. We provide results on our first two synthetic datasets as example:

--

Q4. Typos.

We have fixed the loss function in B.2; thank you for pointing it out. Can you please elaborate on what typos you are referring to in equation (16) and (19)?

--

We thank the reviewer again for helping us improve our paper, and hope that our responses have cleared up the ambiguities. If so, we will appreciate it greatly if you can increase the score; If not, please let us know if you have any further questions.