Copula Conformal Prediction for Multi-step Time Series Forecasting
Authors: Sophia Huiwen Sun
Rose Yu
Venue: Submitted to ICLR 2023 (2023)
Keywords: Conformal Prediction
time series
uncertainty quantification
calibration
RNN

Abstract: Accurate uncertainty measurement is a key step to building robust and reliable machine learning systems. Conformal prediction is a distribution-free uncertainty quantification algorithm popular for its ease of implementation, statistical coverage guarantees, and versatility for underlying forecasters. However, existing conformal prediction algorithms for time series are limited to single-step prediction without considering the temporal dependency. In this paper we propose a Copula-based Conformal Prediction algorithm for multivariate, multi-step Time Series forecasting, CopulaCPTS. On several synthetic and real-world multivariate time series datasets, we show that CopulaCPTS produces more calibrated and sharp confidence intervals for multi-step prediction tasks than existing techniques.

Under review as a conference paper at ICLR 2023
COPULA CONFORMAL PREDICTION FOR MULTI-STEP TIME SERIES FORECASTING
Anonymous authors Paper under double-blind review
ABSTRACT
Accurate uncertainty measurement is a key step to building robust and reliable machine learning systems. Conformal prediction is distribution-fri -freee uncertainty quantification algorithm popular for its ease of implementation, statistical cov- erage guarantees, and versatility for underlying forecasters. However, existing conformal prediction algorithms for time series are limited to single-step predic- tion without considering the temporal dependency. In this paper we propose a Copula Conformal Prediction algorithm for multivariate, multi-step Time Series forecasting, CopulaCPTS. On several synthetic and real-world multivariate time series datasets, we show that CopulaCPTS produces more calibrated and sharp confidence intervals for multi-step prediction tasks than existing techniques.
INTRODUCTION
Deep learning models are becoming widely used in high-risk settings such as healthcare, transporta- tion, and finance. In these settings, it is important that model produces calibrated uncertainty to reflect its own confidence and to assist decision making. Confidence regions are a common approach to quantify forecast uncertainty (Khosravi et al., 2011). A (1 -cconidieece region p1 for random variable y is valid if it contains y/ s true value with high probability: PPyy E ( 2 - d. Note one can make r' infinitely wide satisfy validity. For the confidence region to be useful, we want to minimize its area while remaining valid; this is known as the efficiency of the region.
Conformal prediction is a powerful method that produces confidence regions with kinite-saple guarantees of validity (Vovk et al. 2005 Lei et al., 2018). Furthermore, it makes no assumptions about the forecast model or the underlying data distribution. Its generality, simplicity, and statistical guarantees have made conformal prediction popular for many real world applications including drug discovery Eklund et al. 2015) and robotics (Luo et al. 2021).
In this paper we present more calibrated and efficient conformal prediction algorithm for multi-step times series forecasting. This type of problem is ubiquitous- - examples include predicting hurricanc paths, vehicle trajectories, and financial epidemic forecasts. To quantify uncertainty, we want a "cone uncertainty" that covers the entirety of the forecast horizon. Sttankeviciite et al. (202T) presented algorithm, CF-RNN, for multi-step time series with the assumption that all time steps modeled independently. In practice, however, we found that CF-RNN produces confidence regions often too large be useful, especially when the problem multivariate has a long forecast horizon.
We introduce CopulaCPTS, a Copula-based based Conformal Prediction algorithm for multi-step Time Series forecasting. We improve efficiency by utilizing copulas to model the dependency between forecasted time steps. A copula is a multivariate cumulative distribution function that models the dependence between multiple random variables. is widely used in economic forecasting (Nessen 2007| Patton 2012) and has been introduced to conformal algorithms by Messoudi et al. (2021) to model correlation between multiple targets in non-temporal settings. Copula processes have been explored in generative models for time series (Salinas et al.) 20199; Drouin et al. (022). We found that copulas are effective in capturing uncertainty relations between time steps. Our contributions are:
CopulaCPTS is a general uncertainty quantification algorithm that can be applied to any multivariate multi-step forecaster, with statistical guarantees of validity.
CopulaCPTS produces significantly sharper and more calibrated uncertainty estimates than state-of-te the-an baselines on benchmark datasets, both synthetic and real.

---

Under review as a conference paper at ICLR 2023
We present a re-calibration technique such that CopulaCPTS can produce valid confidence intervals for time series forecasts of varying lengths.
2
RELATED WORK
Uncertainty Quantification for Deep Time-Series Forecasting The two major paradigms of Uncertainty Quantification (UQ) methods for deep neural networks are Bayesian and Frequentist. Bayesian approaches estimate a distribution over the model parameters given data, and then marginal- ize these parameters to form output distributions via Markov Chain Monte Carlo (MCMC) sampling
(Weling & Teh 2011 Neal 2012, Chen et al. 2014) or variational inference (VI) GGrvees 2011
Kingma et al. 2015 Blundell et al. 2015 Louizos & Welling, 2017). Wang et al. (2019) Wu et al.
(2021) propose Bayesian Neural Network (BNN)-based models for UQ of spatiotemporal forecasts. In practice, Bayesian UQ has two major drawbacks: (1) it requires significant modification to the training process of deep models, and (2) BNNs can be computationally expensive and difficult to optimize, especially for larger networks [Lakshminarayanal et al. 2017 Zadrozny & Elkan 2001)) UQ for deep neural network time series forecasts often adopt approximate Bayesian inference such as MC-dropout (Gal & Ghahramani 2016b Gal al. 2017).
Frequentist UQ methods emphasize robustness against variations in the data. These approaches either rely on resampling the data or learning interval bound to encompass the dataset. For time series forecasting UQ. frequentist approaches include ensemble methods such as bootstrap (Efron &
Hastie 2016 Alaa & van der Schaar 2020) and jackknife methods (Kim et al. 2020 Alaa & Van
Der [ehhaar|20200 interval prediction methods include interval regression through proper scoring
rules (Kivaranovic et al. 2020 Wu et al. 202T, and quantile regression (akeecci et 2006), with
many recent advances for time series UQ (Taaasovska & Lopez- Paz, 2 Gasthaus et al. 2019
Park et al. 2022 Kan et al. 2022). Many of the frequentist methods produces asymptotically valid
confidence regions and can be categorized as distribution-fies UQ techniques as they are (1) agnostic
to the underlying model and (2) agnostic to data distribution.
Conformal Prediction. Conformal prediction (CP) an important member of distribution-free UQ methods; we refer readers to Angelopoulos & Bates (222) for comprehensive introduction and survey. CP has become popular because of its simplicity, theoretical soundness, and low computa- tional cost. A key difference between CP and other UQ methods is that under the exchangeability assumption, conformal methods guarantee validity in finite samples (Vovk et al. 2005).
Most relevant to our work is recent progress expanding conformal prediction time-series forecasting. According to (Stankevičūūté et al. (202T) there are two settings: when the data is generated from one single time series from multiple independent time series. For the first setting,
ACI (Gibbs & Candes 202I) and EnbPI (Xu & Xie 2021] developed CP algorithms that relaxes the exchangeability assumption while maintaining asymptotic validity via online learning (former) and ensembling (later); Zaffran et al. (2022) further improves on the online adaptation, and Sousa al. (2022) combines EnbPI with conformal quantile regression (Romano et al.. (019)) to model heteroscedastic time series. These algorithms for single time series are not designed to model multiple independent time series where extracting common patterns can improve forecast and UQ results. The validity guaranteed by these methods are on average over time steps and are often asymptotic, rather than covering the full horizon as in our setting. Stankevičiuté et al. (202T) shares with us the multi-step forecasting for multiple time series setting, though they focuses univariate medical time series. We show that their method of applying Bonferroni correction produces inefficient confidence regions, especially when data is multidimensional. Neeven & Smirnoy (2018) Messoudi et al.) (2020) are CP algorithms for multi-target regression in the non-temporal setting, creating box-like regions to account for the correlations between the labels. In our work develop CoupulaCPTS for multivariate multi-step time series forecasting.
Neeven
BACKGROUND
3.1 INDUCTIVE CONFORMAL PREDICTION (ICP)
Let D = (Z 2n) be a dataset of size n. We denote 2; (æ, 9:) as a sample of an input and output pair that follows the distribution P. where is the data index. The input space X and
2

---

Under review as a conference paper at ICLR 2023
target space Y can be two arbitrary measurable spaces, their Cartesian product X Y the sample space. We begin the algorithm by splitting the dataset into a proper training set Dtrain with [Ptrain] m and a calibration Deal with Peal n m. The objective of conformal prediction is produce valid confidence region (Definition
Definition 1 1 (Validity). Given a new data pair (X, (X,Y) ~ P and a desired coverage rate 1-a E (0, 1).
region (x) is valid if P(X,N)~(Y E (l'"(())) 21-a
An important component of conformal prediction is the nonconformity score A Z T R, function to capture how well sample 2 conforms to the proper training set. For example, we may choose nonconformity score A(z, Dirain) to be the L2-loss on the new sample:
A(2, Dtrain) = l f())]
where f - is a forecasting model trained with the proper training set Dtrain. In the following sections, we refer the nonconformity score of sample = given forecasting model as Az, f) for simplicity. Let 8; = A(z f) denote the nonconformity score of sample 2, and let Scol = [A(2,, f) denote the set of nonconformity scores of all samples in Deat. We define the empirical p-quantile of the set of nonconformity scores S as:
Q(p, S) inf(f'
Σ 5.ES
1. 2 p}
S
Definition (Exchangeability). In dataset (zi, , 2n) of length n, any of its n. permulations are equally probable.
Under the exchangeablility assumption (definition) 2 conformal prediction is theoretically guaranteed to be valid (definition] TT. Let 2n+1 = (In+1 1n+11 Z be new sample, and let [0, 1] be chosen sigmificance level (and 1 the confidence level). Since Deat U {2++1} probability of A( A(2n+1)'s ''s rank among {A(=)),, is uniform. Therefore.
Uf J is exchangeable, the
P(A(zn+1, S Q(1 a, Seal (foo))) = [[(Deal] + 1)(1 - a)l 21-a
(3)
[Deail +1
The conformal confidence region is constructed in equation We say sample is covered if the
true value is in the confidence region Vn+1 E rl (In-1). Equation 3 means the probability of yn+1 being covered is greater than 1 a, hence the confidence region valid.
(In+1) {y A(zn+1, Q(( a, Scal U (co)))
(4)
The procedure introduced above is known as inductive or split conformal prediction, as it splits the dataset into training and calibration sets to reduce the amount of computation (Vovk et al. 2005 Lei & Wasserman, 2012). We introduce it as it is commonly used for uncertainty quantification for machine learning models (Papadopoulos 2008). Although our method in this paper, CopulaCPTS, is implemented based on inductive CP, it can easily adjusted for other CP variants.
3.2 COPULA AND ITS PROPERTIES
In this paper we will model the dependency between time step's uncertainty with copulas. Copula is a concept from statistics that describes the dependency structure a multivariate distribution. We use Copulas to capture the joint distribution for multiple future time steps.
This section will briefly introduce its notations and concepts.
Definition 3 (Copula). Given random vector (X,,... -x),, define the marginal cumulative den- sity function (CDF) as F((x) = P(X, S x), the copula of (X1, .x)xx)) is the joint CDF of (UI.... Uk) (F (X1), Fx(X))), meaning
C(ui. ... u+) P(U S u.....,,
... ,US S us)
Alternatively
C(u .u)) P(X, < (ur),,, X2 S (wx))
3

---

Under review as a conference paper at ICLR 2023
In other words, the Copula function captures the dependency structure between the variable Xs; we can view an & dimensional copula C' : 1] + [0, 1J as a CDF with uniform marginals. A fundamental result the theory of copula is Sklar's theorem.
Theorem 1 (Sklar's s theorem). Given a joint CDF F(X1,'.. X&) and the marginals Fz(x), there exists a copula such that
F(11,"" ,,) C(F(())) ,F((x))
for all x2 E [-00, oo] and 1,..,
Sklar's theorem states that for all multivariate distribution functions, there exists copula function such that the distribution can be expressed using the copula and multiple univariate marginal distributions. To give an example, when all the X are independent, the copula function is known the product
copula: C(u1. ... ux) = II 144.
Sample distribution CDF
Univariate Distributions
Copule Fundton Clti, My)
CDF
CDF of M2,022
-35
L11
148
4.5
Rg
Semple distribttin CDF
Figure I: An example copula, where we express a multivarite Gaussian with correlation p 0.8 with two univariate distributions and Copula function C'u1, u2).
4 COPULA CONFORMAL PREDICTION FOR TIME SERIES (COPULACPTS)
Many real world decision making tasks make use of time series forecasts that are multivariate and can predict multiple steps into the future. In these applications, we want a "cone of uncertainty" that covers the entire course of our forecasts. Existing time series conformal prediction methods either only provide coverage guarantce for individual time step forecasts (Gibbs & Candes) 2021 Xu & Xie 2021) produce confidence regions often too inefficient to be useful, especially in multivariate settings (Stankevičhüté et al. 2021). Hence, we propose a multi-step conformal prediction algorithm for time series, CopulaCPTS, that guarantees validity over multi-step forecasts. We improve efficiency of the confidence regions by modeling dependency of the time steps using a copula function.
We denote the time series dataset as D {(x
1=4+++)) 1 where X1at E gxxd is t time steps of imput which has dimension dy. Vt+1.++ E Rt is time steps of output of dimension dy. In the traditional time series forecasting setting we have d₂ = dy, but they do not have to be necessarily
(i)
equal. Given the dataset D, new test sample x (n+1) and confidence level 1 - a. our algorithm returns l confidence intervals, [[-- ..... To,rtt.r rl-a], one for each time step, such that:
IP[ Vh (I,..., Yeth
(5)
for any underlying predictive model. We define equation as the validity condition for the multi-step time series setting. We will use superscript x" to index data, and subscript y, to index the time steps within the multi-horizon y.
There are two characteristics desired in uncertainty quantification methods: calibration and efficiency. model calibrated when the predicted confidence level corresponds the probability of events falling into the predicted range. It is reflected when equality holds in the validity condition Pyy e TI-a)- 1 a. Efficiency, on the other hand, refers to the size of the confidence region. There is
4

---

Under review as a conference paper at ICLR 2023
trade-off between validity and efficiency, as one can always 1 be infinitely large to satisfy the validity condition. In practice, we want to achieve that the measure of the confidence region (e.g. its area length) to be as small possible, given that the validity condition holds. In the following sections we will introduce CopulaCPTS, conformal prediction algorithm that is both calibrated and efficient for multivariate time series forecasts.
4.1 INDUCTIVE CONFORMAL PREDICTION (ICP) FOR MULTIVARIATE FORECASTS
We view the multidimensional target of dimension d as point in R space, we simply have Y, IRd for each time step t. this paper we chose nonconformity score to be the L-2 distance A((æ, f) := yy y f(x)].. where f is the forecast model trained on the proper training set Drain The confidence region [ -nl-a( (D,x)) therefore is d-dimensional ball. We chose this metric because we are forecasting trajectories in space. Since the conformal prediction algorithm produces valid confidence regions regardless of the choice of nonconformity score, can choose other metrics such as Mahalanobis (Johnstone & Cox 2021) or L-1 (Messoudi et al. 2021 distance based on domain needs, and algorithm will still hold.
4.2 COPULA MULTI-STEP CONFORMAL REGRESSION
We showed in section 3.1 that given confidence level 1 - a and a forecast model f. the ICP algorithm finds nonconformity score 81-a such that for new sample drawn from the data distribution ~ P:
P(A(z, f) s s1-a)
We can use the algorithm to estimate an empirical cumulative distribution function (CDF) for the random variable A(2 f).
F(s): :- P(A(2, f 8)
Σ 1 A(s
(6)
Deal
2
For the multi-step conformal prediction algorithm, we will estimate this empirical CDF for each step of the time series, denoted
(sn) P(A(2,, f) f) 5 st) for h 1,...,k
Let on denote the probability produced by the each time step's CDF function Fa(ss), h E 1...., k. By Sklar's theorem, we have:
F F(s (k) 8k) C(Fi(s)() Fa(ss)) C(1 a1. oik)
(8)
We want to find the set of confidence levels 1 on, such that the entire predicted time series trajectory covered by the intervals with confidence 1 a, i.e. C(1 -a,,..,, anx) 2 1 a. The purpose of using copula is to model the dependency between the multiple predicted time steps, so we can better capture the confidence region for the joint probability. This allows our algorithm to produce more efficient confidence regions. We estimate the copula C' with the calibration set Deat. and then we can search for values a% for each time step to obtain the desired multi-step coverage.
We adopt the empirical copula (Ruschendori) 1976) as our default copula. One may pick other parametric copula functions, such the Gaussian copula, to introduce inductive bias and improve sample efficiency when calibration data is scarce. The empirical copula is non-parametric method of estimating marginals directly from observation, and hence does not introduce any bias. For the joint distribution of time series with & time steps, the copula of vector of probabilities [0, 1 defined as
Cempirscal (u)
Inicu
Σ
II
Luj <u.
i=++1
i=m+1
hml
Where m is the size of the calibration set. Here, the u,s the cumulative probabilities for each data in the calibration set D Deat with size n m.
u' (u ........ w (f)(s)) .... (s)))) i {m 1.....n)
5

---

end for end for Under review as a conference paper at ICLR 2023
Algorithm 1: Copula Conformal Time Series Prediction
Input: Dataset D (Fi, (e) ..... test input X1:t target significant level a. k pl-a
Output: Confidence regions [1 pl-a
// Training
Randomly split dataset D into training and calibration datasets D Deroin U Deat. where
Dirain| = and |Deat| = n
Train k-step forecasting model f with training set Dtrain-
// Calibration
Initialize of nonconformity scores for each time step Sn 0 for h ,,...,k..
for (xi, ((+x444) Deal do
&++L++R f f(x1.t)
for h = 1 to k do
Sn Sn Ulii4t vital
13 Construct CDFs M.... Fx as equation 6
Construct copula C(Fi())..... (-)) equation 9
18 Search for .i....,.$ such that (sfi ... Fs (s))
16 // Prediction
17 gn41 V++144* f(x]:tt
18 for h = 1 to lt do
19 IA (y lly
y <
20 end for
return 1
1 1-a
Note that to fulfill the validity condition of Eqqution 5 we only need to find
u (F)(()).... (()))) such that Cempinical (u") 2
(10)
We can find s7 & through any search algorithm. As exhaustive search is exponential to the prediction horizon k we implement the search with stochastic gradient descent using PyTorch. We study the effectiveness of this method in Appendix C.6. The confidence region for each time step is constructed as the set of all y/ E y, such that the nonconformity score less than s7. Algorithm [ summarizes the CoupulaCPTS procedure. We prove that CopulaCPTS produces valid confidence regions for multi-step time series forecasting (Proposition D in Appendix A
Proposition 1 (Validity of CopulaCPTS). The confidence regions provided by CopulaCPTS (algo- rithm are valid. i.e. IP[ vh E (1,..., k), Yt-h [ 1pll-1 -
4.3 COPULACPTS AUTO-REGRESSIVE FORECASTING
Auto-regressive forecasting is a common framework in time series forecasting. So far, we've been looking at forecasts for a predetermined number of time steps One can use a fixed length model to forecast for longer horizons k autoregressively by taking model output as part of the input. In the conformal prediction setting, we want to not only to autoregressively use the point value forecasts, but also propagate the uncertainty measurement.
we assume the time series to be stationary, then the copula remains the same for any sliding window of k steps, i.e. C(ui...., uk) = C'uu, (1++)) Given a model f to predict time steps. Hence, after we've ve found the wi, -i..... ux such that C(ui,-.., mx)
a,, simply have search for u such that C(u).... "t, 1)+1) ) > 1 a. The guarantee proven in Proposition [ still holds for the new estimation. On the other hand, if the time series non-stationary, we can fit copulas C1 (uI,..., (t),.. Cy & (ug-kr-.. uy ), one for each autoregressive prediction, given that we have k' steps of data our calibration set. This way, we transform the autoregressive problem

---

Under review as a conference paper at ICLR 2023
into k' multi-step problems that can be solved by CopulaCPTS. follows that each of the autoregressive predictions are valid. Appendix[B3. provides an example scenario where re-estimating the copula necessary and improves validity.
EXPERIMENTS
In this section, we show that CopulaCPTS produces more calibrated and efficient confidence regions compared to existing methods on two synthetic datasets and two real-world datasets. We also demonstrate that CopulaCPTS's advantage is more evident over longer prediction horizons. We also show its effectiveness in the autoregressive prediction setting.
Baselines. We compare our model with three representative works in different paradigms of uncertainty quantification for neural network time series prediction: the Bayesian- motivated Monte Carlo dropout RNN (MC-dropout) by Gal & Ghahramani (201 the frequentist blockwise jackknife RNN (BJRNN) by Alaa & Van Der Schaar (2020), and a conformal forecasting RNN (CF-RNN) by Stankevičiute et al. (2021). We use the same underlying prediction model for post-hoc uncertainty quantification methods BJRNN, CF-RNN, and CopulaCPTS. The MC-dropout RNN of the same architecture but is trained separately, as it requires extra dropout step during training and inference.
Metrics. We evaluate calibration and efficiency for each method. For calibration, we report the empirical coverage on the held-out test set. Coverage should be as close to the desired confidence level 1 - a possible. Coverage is calculated as coveragei-a y P(y E (x)) R Ei (y,, the confidence region. The metric should be as small as possible while being valid (coverage maintains above pre-specified confidence level). The area or volume is calculated as area -a (e>x||T>=<x))) x ΣT [[)-=((,))]
(xx)). For efficiency, we report the average area (2D) volume (3D) of
5.1 SYNTHETIC DATASETS
We first test the effectiveness of our models with two synthetic spatiotemporal datasets interacting particle systems (Kipf et al.) 2018), and drone trajectory following (simulated with Robtins
(Sakai et al. 2018)). For the particle simulation we predict Yr+1:++h where = 335 h = 25 and 60, h 10, and V/ R To add randomness to the tasks,
Yt R2; for drone simulation added Gaussian noise of a = .01 and .05 to the dynamics of the particle simulation, and a = .02 to drone dynamics. We generate 5000 samples for each dataset, and split the data by 45/45/10 for train, calibration, and test respectively. Please see Appendix C.1 for forecaster model details.
We visualize the calibration and efficiency of the methods in Figuure 2 for confidence levels 1-a 0.5 to 0.99. We can see that Copula- RNN, the red lines, are more calibrated and efficient compared other baseline methods, especially in tasks with large noise (particle dataset, 0.05). We can see that for harder tasks (particle a = 0.05, and drone trajectory prediction), MC-Dropout is overconfident, whereas BJ-RNN and CF-RNN produce very large (hence inefficient) confidence regions. This behavior of CF-RNN is expected because they apply Bonferroni correction to account for joint prediction for multiple time steps, which is an upper bound of copula functions. The numbers for confidence level 90% are presented in Table [ A quantitative comparison of confidence regions for drone simulation can be found in Figure in the appendix.
5.2 REAL WORLD DATASETS
COVID-19. We replicate the experiment setting of Stankevičiūtt et al. (2021) and predict new daily cases of COVID-19 in regions of the UK. The models take 100 days of data as input and forecast 50 days into the future. We used 200 time series for training, 100 for calibration, and 80 for testing.
Vehicle trajectory prediction. The Argoverse autonomous vehicle motion forecasting dataset (Chang al. 2019) is a widely used vehicle trajectory prediction benchmark. The task is predict 3 second trajectories based on all vehicle motion in the past 2 seconds sampled at 10Hz. Because trajectory prediction a challenging task, we utilize state-of-the-e ct al., 2020) as the underlying model for CF-RNN and Copula-RNN (details in Appendix C.I). Flexibility of underlying forecasting model is an advantage of host-hoc conformal prediction methods.
prediction algorithm LaneGCN (Liang
7

---

Under review as a conference paper at ICLR 2023
Coverage (particle .00
Coverage (particle 05)
Coverage (drone .02)
1.0
0.8
0.8
0.4
0.2
1- 8
0.0
1-c a
Area (particle .01)
Area (particle .05)
Volume (drone 00)
5
10
40
3
5
20
1
0.5
0.6
0.7
0.8
0.9
0.5
0.6
0.7
0.8
0.9
0.5
0.6
0.7
0.8
0.9
1
MC-Dropout
BJRNN
CF-RNN
Copula-RMi
Calibration
Figure 2: Calibration (upper row) and efficiency (lower row) comparison on different - a levels for simulated data sets. For calibration, the goal is coincide with the green dotted calibration line as closely as possible while staying valid (above the green linc). Note that copula methods are more calibrated across different significance levels. For efficiency, want the metric be low as possible. Copula-RNi -NNN outperforms the baselines consistently. (MC-dropout for the right two experiments produces invalid regions, so don't consider its efficiency.)
MC-dropout
BJRNN
CF-RNN
Copula-RNN
Particle Sim (a .01)
Coverage Area
63.9 + 16.1
98.9 0.2
97.0 2
91.5
221
0.81 2.00
2.24 0.59
3.13
13224
1.06
0.36
Particle Sim (a .05)
Coverage
16.1
1 4.3
100.0
0.00
94.5
1.5
90.3
0.7
Area
0.79 0.02
12.13
5.79
= 0.51
4.50
0.07
Drone Sim
Coverage
84.5
10.8
90.8
91.6
90.0
(a = .02)
Volume
9.64 2.133
49.57
3.77
32.18 + 13.66
16.52
+ 7.08
COVID-19 Daily Cases
Coverage Area
19.1 $.1
79.2
+ 30.8
95,4 1.9
92.1
+ 1.0
34.14 = 0.84
823.3
+ 529.7
610.2
96.0
429.0
15.1
Argoverse
Coverage
27.9
3.1
92.6
= 9.2
98.8 1.9
90.4 0.3
Trajectory
Area
127.6 + 20.9
880.8
156.2
396.9 + 18.67
126.8 + 12:22
Table 1: Performance on synthetic and real world datasets with target confidence - a = 0.9. Methods that are invalid (coverage below 90%) greyed out. CopulaCPTS achieves high level of calibration (coverage is close to 90%) while producing more efficient confidence regions (small area).
For model-dependent baselines MC-dropout and BJRNN, have to train an RNN forecasting model from scratch for each method, which is additional computational cost.
Results in Table I show that, like for synthetic datasets, CopulaCPTS is both more calibrated and efficient compared to baseline models for real world datasets. For the trajectory prediction task, learning the copula results in 77% sharper confidence region while still remaining valid for the 90% confidence interval. We visualize two samples from each dataset in Figure 3 The importance of efficiency in these scenarios is clear the confidence regions need to be narrow enough for them to be useful for decision making. Given the same underlying prediction model, we can see that CopulaCPTS produces much sharper region while still remaining valid.

---

Under review as a conference paper at ICLR 2023
MAA. hurccatt
1966
Cop 15-
groune truth
LB
30
chea
thre
Figure 3: Illustraiions of 90% confidence regions given by CF-RNN (blue) and CopulaCPTS (orange) on two real-worle datasets. For COVID-19 forecast (left 2) we see that CopulaCPTS produces sharper yet covering confidence regions. For Argoverse (right 2) we illustrate regions at time steps 1, 10, 20, and 30. Note that the confidence region produced CF-RNN uninformatively large, as covers all the lanes. Overall, CopulaCPTS is able to produce much more efficient confidence regions while maintaining valid coverage. These examples also illustrate the importance of efficiency.
Comparison of models at different horizon lengths. CopulaCPTS is algorithm designed to produce calibrated and efficient confidence regions for multi-step time series. When the predic- tion horizon long, CopulaCPTS's advantage more pronounced. Figure shows performance comparison across increasing time horizons on the particle dataset. Additional experiment results can found in 3 Appendix CopulaCPTSS achieves a 30% decrease area 20 time steps compared to CF-RNN, the best performing baseline; the decrease is above 50% at 25, This experiment shows the significant improvement using copula on modeling the joint distribution of future time steps.
CopulaCPTS for the Autoregressive prediction setting. We test the autoregressive setting (section 4.3) on the COVID-19 dataset. We train RNN model with and use it to autoregressively forecast the next 14 steps. Soble[compare the performance of re-estimating the copula for each 7-step forecasts versus using fixed copula calibrated using the first 7 steps. We compare the model to 14-step joint forecaster using CopulaCPTS as well.
Coverage (Partilee
Area (Particle a= _=.01] -0)
1.00
2.5
Method
Coverage
Area
0.95
AR re-estimate
90.7
377.9
AR fixed
88.2
340.5
0.85
1.0
Joint
90.8
429.4
0.80
0.3
0.75
10
15
20
0.0
5
10
15
20
25
Table 2: Performance of autoregres- sive (AR) CopulaCPTS. Re-estimating copula gives valid confidence re- gion over time and is more efficient than joint CopulaCPTS forecast.
prediction horizon
prediction horizon
MC-Oropoot
UJNN
CI-NNN
RMM
Calibretion
Figure 4: CopulaCPTS remains more calibrated and efficient than baselines over increasing forecast horizons.
CONCLUSION AND DISCUSSION
In this paper we present CopulaCPTS, a conformal prediction algorithm for multidimensions and multi-step time-serics prediction. CopulaCPTS significantly improves calibration and efficiency of multi-step conformal confidence intervals by incorporating copulas model the joint distribution of multiple timesteps. We prove that CopulaCPTS has finite sample validity guarantee. experiments show that the algorithm outperforms state-of-the-ane art models all 4 benchmark datasets and on varying prediction horizons. On the flip side, CopulaCPTS assumes that the copula estimation is accurate and have calibration data for the entirety of the prediction horizon, even in the autoregressive case. We leave it to future work to relax these assumptions (for example by adjusting u' online given prediction errors in Gibbs & Candes (20T)) make CopulaCPTS applicable to the online setting where distribution shift is present.
our

---

Under review as a conference paper at ICLR 2023
REPRODUCIBILITY STATEMENT
All experiments included in this paper are repeated over 3 runs to account for randomness in neural network training. We report the mean and the standard deviation of the experiment results in our tables. We have included code for constructing the synthetic datasets, and source code for our model implementation, experiments, and visualizations in the supplementary material. For real world datasets, we refer readers Appendix C for detailed descriptions of how to obtain and preprocess the data.
REFERENCES
Ahmed Alaa and Mihacla Van Der Schaar. Frequentist uncertainty in recurrent neural networks via blockwise influence functions. In International Conference Machine Learning, pp. 175-190. PMLR, 2020.
AM Alaa and M van der Schaar. Frequentist uncertainty in recurrent neural networks via blockwise influence functions. ICML, 2020.
Anastasios Angelopoulos and Stephen Bates. A gentle introduction conformal prediction and distribution-free uncertainty quantification. arXiv preprint arXiv:2107.075// 2021.
Charles Blundell, Julien Cornebise, Koray Kavukcuoglu, and Daan Wierstra. Weight uncertainty in neural network. International conference on machine learning, pp. 1613-1622. PMLR, 2015.
Ming-Fang Chang, John W Lambert, Patsorn Sangkloy, Jagjeet Singh, Slawomir Bak, Andrew Hartnett, De Wang, Peter Carr, Simon Lucey, Deva Ramanan, and James Hays. Argoverse: 3d tracking and forecasting with rich maps. In Conference on Computer Vision and Pattern Recognition (CVPR). 2019.
Tianqi Chen, Emily Fox, and Carlos Guestrin. Stochastic gradient hamiltonian monte carlo. In International conference on machine learning, pp. 1683-1691, PMLR, 2014.
Alexandre Drouin, Étienne Marcotte, and Nicolas Chapados. Tactis: Transformer-attentioonal copulas for time series. arXiv preprint 2rX2::2202.03528, 2022.
Bradley Efron and Trevor Hastie. Computer Age Statistical Inference, volume 5. Cambridge University Press, 2016.
Martin Eklund, Ulf Norinder, Scott Boyer, and Lars Carlsson. The application of conformal prediction to the drug discovery process. Annals of Mathematics and Artificial Inielligence, 74(1):117-132, 2015.
Yarin Gal and Zoubin Ghahramani. Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In international conference on machine learning, pp. 1050-1059.
PMLR, 2016a.
Yarin Gal and Zoubin Ghahramani. A theoretically grounded application of dropout in recurrent neural networks. Advances neural information processing systems, 29, 2016b.
Yarin Gal, Jiri Hron, and Alex Kendall. Concrete dropout. In NIPS, pp. 3581-35901 2017.
Jan Gasthaus, Konstantinos Benidis, Yuyang Wang, Syama Sundar Rangapuram, David Salinas, Valentin Flunkert, and Tim Januschowski. Probabilistic forecasting with spline quantile function RNNs. In AISTATS 22, pp. 1901-1910, 2019.
Isaac Gibbs and Emmanuel Candes. Adaptive conformal inference under distribution shift. Advances in Neural Information Processing Systems, 34, 2021.
Alex Graves. Practical variational inference for neural networks. Advances in neural information processing systems, 24, 2011.
Chancellor Johnstone and Bruce Cox. Conformal uncertainty sets for robust optimization. In Conformal and Probabilistic Prediction and Applications, pp. 72-90. PMLR, 2021.
10

---

Under review as a conference paper at ICLR 2023
Kelvin Kan, François-Xavice Aubet, Tim Januschowski, Youngsuk Park, Konstantinos Benidis, Lars Ruthotto, and Jan Gasthaus. Multivariate quantile function forecaster. In International Conference on Artificial Intelligence and Statistics, pp. 10603-10621. PMLR, 2022.
Abbas Khosravi, Saeid Nahavandi, Doug Creighton, and Amir F Atiya. Comprchensive review of neural network-based prediction intervals and new advances. IEEE Transactions on neural networks, 22(9): 2011.
Byol Kim, Chen Xu, and Rina Barber. Predictive inference free with the jackknife+-after-boood bootssssr. In H. Larochelle. M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing Systems, volume 33, pp. 4138-411 4149. Curran Associates, Inc., 2020.
Durk Kingma, Tim Salimans, and Max Welling. Variational dropout and the local reparameterization trick. Advances in neural information processing systems, 28, 2015.
Thomas Kipf, Ethan Fetaya, Kuan-Chich Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. arXiv preprint 180.vv:802.0440 2018.
Danijel Kivaranovic, Kory D Johnson, and Hannes Leeb. Adaptive, distribution-free prediction intervals for deep networks. In AISTATS, pp. 4346-4356, 2020.
Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. Simple and scalable predictive
uncertainty estimation using deep ensembles. Advances in neural information processing systems, 30, 2017.
Jing Lei and Larry Wasserman. Distribution free prediction bands. arXiv preprint arXiv: arXiv:120..5 1203.542. 2012.
Jing Lei, Max G'Sell, Alessandro Rinaldo, Ryan Tibshirani, and Larry Wasserman. Distribution-fre predictive inference for regression. Journal of the American Statistical Association, 113(523): 1094-11 1194-111,, 2018
Ming Liang, Bin Yang, Rui Hu, Yun Chen, Renjie Liao, Song Feng. and Raquel Urtasun. Learning lane graph representations for motion forecasting. In European Conference on Computer Vision, pp. 541-556. Springer, 2020.
Christos Louizos and Max Welling. Multiplicative normalizing flows for variational bayesian neural networks. In International Conference on Machine Learning, pp. 2218-22111211 2227. PMLR, 2017.
Rachel Luo, Shengjia Zhao, Jonathan Kuck, Boris Ivanovic, Silvio Savarese, Edward Schmerling, and Marco Pavone. Sample-efficient safety assurances using conformal prediction. arXiv preprint arXiv:2109. 14082, 2021.
Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Conformal multi-targel regression using neural networks. In Conformal and Probabilistic Prediction and Applications, pp. 63. PMLR, 2020.
Soundouss Messoudi, Sébastien Destercke, and Sylvain Rousseau. Copula- conformal predic- tion for multi-target regression. Pattern Recognition, 120: 108101 2021.
Radford M Neal. Bayesian learning for neural networks, volume 118. Springer Science & Business Media, 2012.
Jelmer Neeven and Evgueni Smirnov. Conformal stacked weather forecasting, In Conformal and Probabilistic Prediction and Applications, pp. 233-233. PMLR, 2018.
Roger B Nelsen. An introduction to copulas. Springer Science & Business Media, 2007.
Harris Papadopoulos. Inductive conformal prediction: Theory and application neural networks. In
Tools in artificial intelligence. Citeseer, 2008
Youngsuk Park, Danielle Maddix, François-Xave Xavier Aubet, Kelvin Kan, Jan Gasthaus, and Yuyang Wang. Learning quantile functions without quantile crossing for distribution-free time series forecasting. In International CaciinaaCoooofee Conference on Artificial Intelligence and Statistics, pp. 8127-8150.
PMLR, 2022.
11

---

Under review as a conference paper at ICLR 2023
Andrew Patton. A review of copula models for economic time series. Journal of Multivariate Analysis, 110:4-18, 2012.
Yaniv Romano, Evan Patterson, and Emmanuel Candes. Conformalized quantile regression. Advances in neural information processing systems, 32, 2019.
Ludger Ruschendorf. Asymptotic distributions of multivariate rank order statistics. The Annals of Statistics, pp. 912-923, 1976.
Atsushi Sakai, Daniel Ingram, Joseph Dinius, Karan Chawla, Antonin Raffin, and Alexis Paques. Pythonrobotics: a python code collection of robotics algorithms. CoRR, abs/1808. 18808.10000 2018.
David Salinas, Michael Bohlke-Schneideer Laurent Callot, Roberto Medico, and Jan Gasthaus. High- dimensional multivariate forecasting with low-rank gaussian copula processes. Advances in neural information processing systems, 32, 2019.
Martim Sousa, Ana Maria Tomé, and José Moreira. A general framework for multi-step ahead adaptive conformal heteroscedastic time series forecasting. arXiv preprint arXiv:2207 14219, 2022.
Kamilé Stankevičiuté, Ahmed Alaa, and Mihacla van der Schaar. Conformal time-series forecasting. In Advances in Neural Information Processing Systems, 2021.
Natasa Tagasovska and David Lopez-Paz. Single-model uncertainties for deep learning. In NeurlPS, pp. 6417-6428, 2019.
Ichiro Takeuchi, Quoc Le, Timothy Sears, Alexander Smola, et al. Nonparametric quantile estimation. MIT Press, 2006.
Vladimir Vovk, Alexander Gammerman, and Glenn Shafer. Algorithmic learning in random world. Springer Science & Business Media, 2005.
Bin Wang. Jie Lu, Zheng Yan, Huaishao Luo, Tianrui Li, Yu Zheng, and Guangquan Zhang. Deep uncertainty quantification: machine learning approach for weather forecasting. In Proceedings of the 25th ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pp. 2087-2095, 2019.
Max Welling and Yee W Teh. Bayesian learning via stochastic gradient langevin dynamics. In ICML, pp. 681-688, 2011.
Dongxia Wu. Liyao Gao, Matteo Chinazzi, Xinyue Xiong, Alessandro Vespignani, Yi-An Ma, and Rose Yu. Quantifying uncertainty in deep spatiotemporal forecasting, In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery Data Mining, pp. 1841-1851, 2021.
Chen Xu and Yao Xie. Conformal prediction interval for dynamic time-series. International Conference Machine Learning. PMLR, 2021.
Bianca Zadrozny and Charles Elkan. Obtaining calibrated probability estimates from decision trees and naive bayesian classifiers. In lcml, volume 1, pp. 609-616. Citeseer, 2001.
Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Dieuleveut. Adaptive conformal predictions for time series. In International Conference on Machine Learning, pp. 25834 25866. PMLR, 2022.
12

---

Under review as a conference paper at ICLR 2023
APPENDIX
A PROOF OF PROPOSITION 1
Proposition (Validity of CopulaCPTS). The confidence region provided by CopulaCPTS (algorithm 7 valid. ie. P( (l,...,k), k), Nrth E T,, al-ai
Proof. We estimated u² = (Fi F((st) such that
Compirical (u") Au 2
(1I)
Let (xf+t .8+1144) ~ X Y be a new data point. Denote ý = fx the forecast given by the trained model. Because the Fy functions and Cempirical monotonously increasing, have:
P[ E (l,...,k) k), Ye+h E ry T
C(B((stt))..,
2 (..((s)).....(((5))
2 a
B ADDITIONAL ALGORITHM DETAILS
B. UPPER AND LOWER BOUNDS COPULAS
To provide better understanding of the properties of Copulas, consider the Frechet-Hoeffding Bounds (Theorem) 2. In fact, the Trechet-Hoeffdiig upper- and lower- bounds are both copulas. The lower bound is percisely the Bonferroni correction used in (Stankevičiute al. (202T) therefore by estimating the copula more precisely instead of using lower bound, we have a guaranteed efficiency improvement for the confidence region.
Theorem (The Frechet- Hoeffding Bounds). Consider a copula C(ui,..., Then
max (1 - + Ui 0} s C(ui.... .4) & min (u.,..., uR
B.2 NUMERICAL. OPTIMIZATION WITH SGD FOR SEARCH
The empirical inverse CDF modeled the same way as Equation 6 is the same as the empirical quantile estimation (Equation 2)
F Fー) = inf(s:
Σ
A(a.f)so) 2 p}
(12)
Deat
We find the optimal sf in Equation 10 and Algorithm Ib minimizing the following loss:
L(s
(t)
Σ
I h=1
un < F (sw) -a)
We use the Adam optimizer and optimize for 500 steps. An alternative and faster way of search is to
assume $1 82 Sg which simplifies the search from exponential time constant time.
13

---

Under review as a conference paper at ICLR 2023
B.3 AUTOREGRESSIVE PREDICTION
In the context of this paper to forecast autoregressively is given input X1:t and l step forecasting model f, perform prediction
844144 f(x1a)
81+20+R+1 (x2:1;;
until all k time steps are predicted.
We now provide toy scenario illustrate when re-estimating the copula is necessary and improves validity. Consider a time series of three time steps (o, 41, l2. The two scenarios are illustrated in Figure 5 In both scenarios the mean and variance of all timeseps
0 and respectively. In scenario (a). to = & and hence their covariance is 1. The copula estimated on to and t is Co.1 (fo)), ((t,))
F(t) P(4)). This copula will significantly underestimate the confidence region of t2 where its covariance with & is -1. In fact the coverage of Co:1 (Fi (t1), F2(t2) 0.74. On the other hand, (b) illustrates a scenario where the copula for any consecutive time series remain the same Co Cj. In this case, applying Co directly to forecast C₁ achieves precisely 90% coverage.
0.5
325
3.55
209
522
0.50
315
308
(a) Time steps are positively correlated while and are negatively correlated
(b) Stationary series where each time steps are uncorrelated
Figure Two scenarios illustrate the autoregressive case
C EXPERIMENT DETAILS AND ADDITIONAL RESULTS
C.I UNDERLYING FORECASTING MODELS
Particle Dataset The underlying forecasting model for the particle experiments is an 1-layer LSTM network with embedding size 24. The hidden state is then passed through a linear network to forecast the timesteps concurrently (output has dimension k dyy). We train the model for 150 epochs with batch size 128. Hyperparameters of the network are selected through a model search by performance on 5-fold cross validation split of the dataset. The architecture and hyperparameters are shared for all baselines and CopulaCPTS in Table 1.
Drone For the drone trajectory forecasting task, particle dataset, but with hidden size increased to 128. We train the model for 500 epochs with batch size 128. The same architecture and hyperparameters are shared for all baselines and CopulaCPTS reported in Table 1.
he same LSTM forecasting network the
Covid-19 The base forecasting model for Covid-19 19 dataset is the same as the synthetic datasets, with hidden size 128 and were trained for 150 epochs with batch size 128. The same architecture and hyperparameters are shared for all baselines and CopulaCPTS reported in Table
14

---

Under review as a conference paper at ICLR 2023
Argoverse As highlighted in the main text, we utilize state-of-the-art prediction algorithm LaneGCN (Liang et al. 2020) as the underlying forcaster model for CF-RNN and Copula- RNN. We reter the readers their paper and code base for model details. The architec- ture of the RNN network used for MC-Dropout and BJRNN is Encoder-Decoder net- work. Both the encode and decoder contains LSTM layer with encoding size 8 and hid- den size 16. We chose this architecture because the is part of the official Argoverse base- lines (https: : / / thub. com/ seeeetsiiinngggggvrrrereeeee forecast igg) and demon- strates competitive performance.
C.2 COVID-19 DATASET
The COVID-19 dataset is downloaded directly from the official UK government website ht http:::/coronaviuuue lype and new BueeuusDaByyllliihDaett for metric. There are in total 380 regions and over 500 days of data, depending when it is downloaded. We selected 150-day time series from the collection construct our dataset.
data. wev.uk/detalll /dwwnload by selecting region for area
C.3 CALIBRATION AND EFFICIENCY CHART FOR COVID-19
Figure comparison of calibration and efficiency for the daily new COVID 19 cases forecasting.
Coverage (50 days)
Width (daily averaqe!
1.0
1500
0.8
1250
0.6
750
0.4
500
0.2
250
0.0
0.7
0.8
1-2
1-
MC dropout
CF-RNN
BJRNN
Copula RNN
Copula EncCec
Figure 6: Calibration and efficiency comparison on different level for COVID-19 19 Daily Forecasts. The copula methods (orange and red lines) are more calibrated (coinciding with the green doted line) and sharp (low width) compared baselines.
To see if the daily fluctuation due to testing behaviour disrupts other method, also the same experiment on weekly aggregated new cases forecast. We take 14 weeks data as input and output forecasts for the next weeks. The results are illustrated in Figure 7 The weekly forecasting scenario gives us similar insights as the daily forecasts.
Coverage (over weeks)
Width (weekly average)
1.0
3000
0.8
2500
0.6
1500
1000
500
0.5
0.6
0.7
0.8
0.9
0.5
0.6
0.7
0.8
0.9
1
1
mcdp
dmnn
hjmn
capula-mn
ropulla-EncDee
Figure 7: Covid Weekly Forecasts
15

---

Under review as a conference paper at ICLR 2023
C.4 ARGOVERSE
The Argoverse autonomous vehicle dataset contains 205,942 samples, consisting of diverse driving scenarios from Miami and Pittsburgh. The data can be downloaded from the official Argoverse dataset website. We split 90/10 into training set and validation set of size 185,348 and 20,594 respectively, The official validation set of size 39,472 is used for testing and reporting performance. We preprocess the scenes filter out incomplete trajectories and cap the number of vehicles modeled 60. If there are less than 60 cars in the scenario, we insert dummy cars into them to achieve consistent car numbers. For map information, we only include center lanes with lane directions as features. Similar to vehicles, we introduce dummy lane nodes into each scene to make lane numbers consistently equal 650.
C.5 ADDITIONAL EXPERIMENT RESULTS
We present in figure 8 and 9 some qualitative results for uncertainty estimation.
To test how the effects of copulaCPTS compare with baseline on other base forecasters, we also include encoder-decoder architecture with the same embedding size as the RNN models introduced in appendix C.I for each dataset. The results are presented Table 3 We omit these results in the main text because we found that they do did not bring significant improvement to time series forecasting UQ.
Table comparee model performance compared across different prediction horizons. We show that the advantage of our method is more pronounced for longer horizon forecasts.
45
(a)
(b) MC Dropout
(c) CF-RNN
Figure 8: 99% Confidence region produced by three methods for the drone dataset. Copula methods (a) produces more consistent, expanding cone of uncertainty compared to MC-Dropout (b) sharper one compared to CE-RNN (c).
Figure 9: Illustrations for confidence regions given by CF-RNN (blue) and CopulaCPTSS (orange) in at time steps O, 10, 20, and 30. Note that in order to achieve 90% coverage, the regions are larger than needed, especially straight-lane cases like the middle two. Using copulas couple together time steps results in much smaller region, while achieving similarly good coverage.
16

---

Under review as a conference paper at ICLR 2023
Particle Simulation (a = .01)
Coverage (90%)
Area (90%)
Coverage (99%)
Area (99%)
MC-dropout
63.9
0.81
i
73.1
1.14
BJRNN
98.9
+ 0.2
2.24
= 0.59
100 = 0.0
3.75
+ 0,71
CF-RNN
97.0
+ 2.3
3.13
3.24
100 0.21
4.34
Copula-RNN
91.5
2.1
1.06
= 0.36
97.4 2.4
4.47
+ 1.48
Copula-EncDec
90.2 + 2.3
1.19 + 0.17
99.6 0.7
1.83 0.09
Particle Simulation (a .05)
Coverage (90%)
Area (90%)
Coverage (99%)
Area (99%)
MC-dropout
16.1 4.3
0.79 = 0.02
33,9 5.1
2.12
0.03
BJRNN
100.0 +00
12.13 = 0.39
100.0 0.0
15.43
0.85
CF-RNN
94.5 1.5
5.79
= 0.51
99.8 2.2
19.21
8.19
Copula-RNi
90.3 1 0.7
4.50 + 0.07
99.1 1 0.8
12.82
+ 3.98
Copula-Ence EncDec
91.4 + 1.1
4.40 + 0.15
98.7 0.1
9.31 1.97
Drone Simulation
.02)
Coverage (90%)
Area (90%)
Coverage (99%)
Area (99%)
MC-dropout
84.5 + 10.8
9.64
2.13
90,0 7.8
16.02
3.62
BJRNN
90.8 22
49.57
1 3.77
100.0 1 4.0
65.77
4.56
CF-RNN
91.6 1.2
32.18
+ 13,66
100.0 0.0
36.79
1 14.03
CF-EneDec
100.0
1:
21.83
1 26.29
100.0
0.0
25.03
1 12.53
Copula-RNN
90.0 1.5
16.52
+ 15.08
98.5 + 0.5
21.48
8.91
COVID-19 Daily Cases Dataset
Coverage (90%)
Area (90%)
Coverage (99%)
Area (99%)
MC-dropout
19.1 5.1
34.14 14 = 0.84
100.0 = 0.0
1106.57
25.41
BJRNN
79.2 308
823.3
+:20 529.7
85.7 27.5
149187.
1:
CF-RNN
95.4 19
610.2
1 96.0
100.0
= 0.0
121435.
1: 26495.
CF-EneDec
91.7
= 1.4
570.3
22.1
100.0
108130
+ 10889.
Copula-RNN
92.1 1.0
429.0
15.1
100.0 0.0
88962.
+ 9643.
Dopula-EneDec
90.8 0.3
429.4
# 27.9
100.0 = 0.0
60852.
12263.
Argoverse Trajectory Prediction Dataset
Coverage (90%)
Area (90%)
Coverage (99%)
Area (99%)
MC-dropout BJRNN
27.9 13.1
127.6
20.9
31.5 3.9
242.1 54,0
92.6 9.2
880.8
+ 156.2
100.0
= 0.0
3402.8
= 268.
CF-LaneGCN
98.8 1.9
396.9 + 18.67
100.
0.2
607.2 8.67
Copula-LaneGENN
90.4 1: 0.3
126.8 1: 12.22
99.1 1: 0.4
324.1 1. 42.22
Table 3: Additional results. Copula methods achieves high level of calibration while producing sharper predictions regions. The sharpness gain is even more pronounced in at higher confidence levels (99%), where we want the prediction region be useful while remaining valid.
17

---

Under review as a conference paper at ICLR 2023
Step
Steps
15 Steps
Method
Coverage
Area
Coverage
Area
Coverage
Area
MC-Dropout
97.8 2.0
0.4 0.M4
88.0
0.69
0.25
52.3 14
0.94
BJRNN
45.3 = 39.4
0.27
= 0.18
97.7 2.1
2.69
+ 1.79
95.5 2.8
19.99
- 4.83
CE-RNN
100.0 +
0.01
+: 0.01
77.8 =
0.8
+:
66.7 0.0
18.82
3.73
De-EncDec
89.9 = 19.2
0.01
+ 0.01
100.0 +.00
0.75
+ 0.99
88.9
19.2
13.07
= 16.1
Copula-RNN
90.1 + 0.2
0.01
+ 0.01
89.8
0.54
+ 0:45
90.1 12
8.25
3.44
Copula- Dec
90.0 1 0.3 0.01 1 0.0
90.3
0.67 1. 1.01
90.5 1 0.5
7.13 1.55
Table Performance comparison across different horizons 90% confidence level on the drone simulation dataset. The improvement on efficiency is more pronounced when the horizon is longer.
C.6 STUDY ON 0% SEARCH
Figure 10 shows the on values for each an = (st) used in Copula CPTS outlined in line We present ot values searched using two methods of searching, with dichotomy search for constant value for the horizon as Messoudi et al. (202T), and by stochastic gradient descent as outlined in section 4.2.
15 of Algorithm[] [
The Or values are an indicator of how interrelated the uncertainty between each time step are: Bonferroni Correction used in SStankevicüüte et al. (202T) (grey dotted line in Figure To) assumes that the time steps are independent, with CopulaCPTS we have lower - at levels while having valid coverage (blue and orange lines in Figure I0). This shows that the uncertainty of the time steps not independent, and we are able to utilize this dependency to shrink the confidence region and still maintaining the coverage guarantee.
[able[]shows that there are no significant difference between coverage and area performance for the two search methods with in the scope of datasets we study in this paper. However, we want to highlight that SGD search is O(n) complexity to optimization steps, regardless to the prediction horizon. SGD also allows for varying O% which might be useful in some settings, for example capturing uncertainty spikes for some time steps as seen in the COVID-19 dataset of Figure I0 Dichotomy search, on the other hand, is (n)og(n)) complexity the search space depends on granularity, and will be O(kllog(kn) if we want to search for varying an.
Dataset
Coverage (90%)
Area
Fixed an
Varying 0+
Fixed on
Varying On
Particle (a .01)
91.7
1.9
91.5 4 2.1
1.13
+ 0.45
1.06
1 0.36
Particle (a = .05)
92.1
+ 1.3
90.3
1 0.7
4.89
1: 0.05
4.50
1: 0.07
Drone
90.3
1 0.5
90.0
1 1.5
15.92
1.98
16.52
1 7.08
Covid-19 19
92.9
# 0.1
92.1
4 1.0
498.44
6.36
429.0
+
Argoverse
90.2
0.1
90.4
+ 0.3
117.1
126.8
+ 12.2
Table 5: Coverage and area comparison between stochastic search for fixed and SGD for Varying ahe We not see significant difference between the performance two.
18

---

Under review as a conference paper at ICLR 2023
Particles
1.000
0.95
0585
0.905
0.80
0.980
0.75
0.955
0970
10
Prediction Horizon
Prediction Horizon
Drone Simulation
COVID-19
0.90
09975
0985
0989
N w
0.75
0970
0965
0:850
09825
0.800
10
Prediction Horizon
Prediction Horizon
Argoverse
0.995
0.900
0.900
Bonferroni Correction
Fixed
0.70
Varying
0960
Prediction Horizon
Figure 10: Comparison between dichotomy search for fixed ag values (blue) and stochastic gradient search for varying an (blue) through timesteps. Shaded regions are the standard deviation of the values over runs.
19

---

Under review as a conference paper at ICLR 2023
C.7 COMPARISON TO ADDITIONAL BASELINES
We include comparison to two additional simple UQ baselines on the particle simulation dataset.
L2-Conformal. L2-Cooformal uses the same underlying RNN forecaster as CF-RNN and Copula RNN. We use that nonconformity score of the vector norm of all timesteps concatenated together 18++124 - -141:++kl to perform ICP. As there are no analytic way to represent uncertainty region on 2-D space, we calculate the area and plot the region for L. Conformal baseline with the maximum deviation at each timestep such that the vector norm still stays within range.
x dyy-dimensional
Direct Gaussian. Direct Gaussian use the same model architecture and training hyperparameters. with the addition of a linear layer that outputs the variance for each timestep, and is optimized using negative log loss, proper scoring rule for probabilistic forecasting. We obtain the area by analytically calculating the 90% confidence interval for each variable.
Results in Table 6 shows that L2-conformal produces inefficient confidence area, and directly outputting variance under-covers test data. These results align with previous findings and motivates method, which is both more calibrated and sharper compared these baselines. We show visualization in Figure to to illustrate the different properties of the methods qualitatively.
Particle (a .01)
Particle (a = .05)
Method
Coverage (90%)
Area
Coverage (90%)
Area 4
L2- 1.2-Confrmma
88.5 = 0.4
7.21 + 0.35
89.7
7.21 0.35
Direct Gaussian
11.9 0.09
0.07 = 0.31
0.0 0.0
0.08 = 0.02
CF-RNN
97.0 = 2.3
3.13 + 3.24
97.0 + 2.3
5.79 = 0.51
Copula-RNN
91.5 + 2.1
1.06 + 0.36
90.3 + 0.7
4.50 + 0.07
Table 6: Comparison with two additional baselines on the particle dataset.
0.2
01
01
90
0.
-01
-011
022
-02-
-03
-03
-12
-1.0
-088
0.
-0.4
(a) L.2 Conformal
(b) Direct RNN Gaussian
01
01
00
00
-01
-0.1
-02
-02
-03
-033
-04
-0.8
-0.6
-0.4
-12
-10
-08
02
CF-RNN
Copula-RNN
Figure 11: Visualization of on sample from the Particle dataset's test set.
20


Reviews:
[{"review_id": "xNcozyBB1X", "reviewer": "Reviewer_aA8r", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "Summary of the paper\nIn this paper, the authors present a calibrated and efficient conformal prediction algorithm for multi-step time series forecasting and utilize copulas to model the dependency between forecasted time steps. The results show the proposed method performed well on synthetic and real-world datasets both in calibration and efficiency.\n\nContributions\n\n(1) The proposed method is a general uncertainty quantification algorithm that has wide applications.\n\n(2) CopulaCPTS shows competitive performance over other baselines and the experiment is sufficient.\nStrength and weaknesses", "summary_of_the_review": "In this paper, the authors present a method to quantify the uncertainty of the prediction model, which is novel and meaningful. The designed experiments are efficient and achieve the expected results. Hence, this is a good paper that could be accepted.\n", "strength_and_weaknesses": "Pros:\n\n(1) The proposed method is novel and the proof is solid.\n\n(2) The conducted experiments are comprehensive and show outstanding performance.\n\n(3) The concerned problem has practical applications in the prediction model.\n\nCons:\n\n(1) Some case studies can be added to interpret its detailed application better.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "fRmp_C4KT3", "reviewer": "Reviewer_FqUe", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes a novel conformal prediction method for time-series data, called Copula Conformal Prediction for multi-step time series forecasting (CopulaCPTS), for uncertainty quantification. In particular, the paper uses the concept of “copula” from statistics; the copula is a joint CDF on data, which fully captures the dependency structure of the data, like in time-series data. The proposed CopulaCPTS theoretically justified its validity (in Proposition 1), and empirically evaluated seven synthetic datasets and two real datasets for its efficacy. The claimed contribution are (1) CopulaCPTS is a general UQ algorithm for any multivariate multi-step forecaster with statistical validity, (2) CopulaCPTS produces sharper uncertainty estimates compared with other approaches, and (3) CopulaCPTS produces valid confidence intervals for time-varying time-series prediction.", "summary_of_the_review": "The core contribution of the paper is the coverage guarantee of the proposed approach, but to my understanding, the proof is not valid; so I would vote for rejection. However, I’m willing to adjust my understanding after discussion.\n\n========== after rebuttal\n\nThanks for the clarification and additional experiments. I believe that exploiting the concept of copula is interesting, but (1) the paper relies on a strong assumption (i.e., an accurate Copula is given, which is even not explicitly mentioned in the paper or written in a confusing way), (2) my main concerns were not addressed during the rebuttal, and (3) I think the paper is not ready to publish based on the current manuscript status; so, I maintain my score. \n\n\nThe following includes the remaining main concerns. \n\n* The probability in coverage is taken over a calibration set. Note that the probability in the full CP coverage is taken over a training set (aka a proper training set in CP), a calibration set (aka a training set in CP), and a new sample. The probability in inductive CP (ICP) coverage is taken over a calibration set and a new sample (where the training set is used for learning a score function). The probability in training-conditional ICP coverage is taken over a new sample. I’d recommend checking the difference between ICP and training-conditional ICP in [1]. \n* The above definition is crucial in proving Proposition 1 if the accurate Copula is unknown. If Proposition 1 is proven without this strong assumption, I’m wondering if the calibration set can be reused so Bonferroni correction needs to be used in this step. In the worst case, the rigorous Proposition 1 produces a more conservative algorithm such that a more rigorous proposed approach may not be clearly better than CF-RNN. In other words, by making an assumption that the accurate Copula is known, the proposed algorithm becomes heuristic, and it looks unfair to compare it with a more rigorous approach like CF-RNN — usually rigorous algorithms are more conservative (in achieving a desired coverage) than heuristic ones. In this sense, I think proving Proposition 1 by estimating a Copular from data is required as for the first step. \n", "strength_and_weaknesses": "**Strengths**:\n\n* Highlight the concept on Copula for time-series data\n* Extensive experiments on both simulated and interesting real datasets.\n\n**Weaknesses**:\n* Proposition 1 is invalid (or clarification is required)\n* A simple but important baseline is missing\n\n(Weakness 1) First of all, the ICP validity definition on coverage (Definition 1)  is not correct (unless there are reasons, which are not described). The probability needs to be taken over X, Y, and also a calibration set (e.g., [R1]). Also, this coverage guarantee is not used consistently. In particular, the validity statement in Proposition 1 contains the expectation over y, but to my understanding, the probability here is also taken over y. I guess these misunderstandings might contribute to the incorrect proof of Proposition 1. In particular, (11) assumes u* is given, which needs to be found in Line 15 of Algorithm 1 based on a calibration set. As u* depends on a calibration set, it is a random vector. However, this randomness is not taken into account in the proof as the probability in the proof is also taken only over a new data point. Finally, it is unclear how the exchangeability assumption is used in the proof. \n\nTo my understanding, the proof needs to be corrected based on the right coverage definition, and the proposed approach is reevaluated based on the correction.\n\n[R1] https://arxiv.org/pdf/1904.06019.pdf\n\n\n(Weakness 2) I think the beauty of conformal prediction is that it provides a coverage guarantee for any score function, as also highlighted in the paper. Thus, one straightforward way to handle time-series data is using a score function that aggregates scores of future time steps and simply runs the standard ICP on the aggregated score. For example, based on the paper’s notation, a RNN output and a label are $\\hat{\\mathbf{y}}, \\mathbf{y} \\in \\mathbb{R}^{k \\times d_y}$, respectively, then, a score function can be simply $|| \\hat{\\mathbf{y}} - \\mathbf{y} ||$, which is equivalent to the negative log likelihood of the Gaussian with the mean of $\\hat{\\mathbf{y}}$ and the identity covariance. This score function can be improved by estimating the diagonal of the covariance matrix (assuming zero correlation) by modifying RNN to estimate the standard deviations of $\\hat{\\mathbf{y}}$ (I believe this is a quite standard techniques to convert a point estimator to a uncertainty estimator, e.g.,[R2,R3]). Given any of these two score functions, we can run the standard ICP. \n\nI think the limitations of these two simple baselines need to be discussed and empirically evaluated to justify the efficacy of the proposed approach. Moreover, we can use more sophisticated forecasting model (e.g. Bayesian filtering with state evolution [R4,R5]) for a score function along with the standard ICP for the coverage guarantee. What’s the limitation of having a sophisticated score function for time-series data with the standard ICP for the coverage guarantee compared with the proposing new conformal prediction dedicated to time-series data along with a naive score function? \n\n[R2] https://openaccess.thecvf.com/content_CVPR_2019/papers/He_Bounding_Box_Regression_With_Uncertainty_for_Accurate_Object_Detection_CVPR_2019_paper.pdf \n[R3] https://arxiv.org/pdf/2001.00106.pdf\n[R4] https://arxiv.org/pdf/1805.11122.pdf\n[R5] https://arxiv.org/pdf/1605.07148.pdf\n\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "1: The main claims of the paper are incorrect or not at all supported by theory or empirical results.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ei7X0IYDqv", "reviewer": "Reviewer_gQHv", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper proposes conformal prediction method for multivariate multistep TS forecasting.", "summary_of_the_review": "This work proposed simple and yet effective method for multi-step forecasitng,  even with several unclear descriptions like general forecasting settings and experiment settings.", "strength_and_weaknesses": "## Strength\n- develop more efficient conformal prediction method with the simple yet effective idea of Coupla.\n- clear motivation and idea\n- extensive experiments are done\n\n\n\n## Weakness\n- multivariate parts that author claim to be new seems to not be really in the novel parts, but multi-step ones are main focus. Otherwise, kindly describe why existing methods by Stankevicˇiu ̄te ̇ et al. (2021) cannot handle this setting (with small modifications).\n- enhance UQ literatures and discussion on other classes. For example there are several advances in forecasting based on quantile regresison methods for univariate and multivariate https://arxiv.org/abs/2202.11316,https://arxiv.org/abs/2111.06581, and even conformal prediction enhance quantile regression method itself https://arxiv.org/abs/1905.03222.\n- more deep dive on behvarious ${\\alpha_t\\}$. For example, is this related to the author claim the capability of capturing dependency\n- the experiment setting are poorly described and difficult to examine univatie/multivate cases and which various base forecasters are used and why.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "e-YaZkFGksZ", "reviewer": "Reviewer_oH8D", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper describes a new conformal method to perform guaranteed prediction in multi-step multi-variate time series. In order to achieve better calibration and efficiency across dimensions, it proposes the use of copulas. ", "summary_of_the_review": "The paper introduces a new method to perform multi-dimensional, multi-step ahead time series prediction. Experiments performed on various data sets are convincing, but remain limited in the number of dimensions predicted (at most 3 dimensions in the case of drone trajectory prediction). This is however fine, as the focus is on multi-step time series prediction. \n\nAuthors may be interested in knowing that the approach of Johnstone and Cox has been recently adapted to cope with normalized scores, hence providing adaptive ellipsoids (see \"Messoudi, S., Destercke, S., & Rousseau, S. Ellipsoidal conformal inference for Multi-Target Regression.\").\n\n", "strength_and_weaknesses": "+: a new method to solve a difficult problem, easily applicable, which appears sound from what I understood. \n\n+: experimental results seem to be good, and are accompanied with code (that I could not run on my laptop in due time, but in reason of my install and the lateness of my review... I will update if I manage to run it in the next days)\n\n-: many typos in the paper, suggesting a rushed in submission. Among others:\n* exchangability should be exchangeability (numerous times)\n* serires --> series\n* equation 1 is an $L_1$ norm, not a $L_2$. Also, if I am correct non-conformity scores in split conformal should be computed on $\\mathcal{D}_{cal}$\n* equation 4, $A(x_{n+1},\\hat{f})$ should be $A((x_{n+1},y),\\hat{f})$ to be consistent with earlier notations\n* time-step is written time-step, timestep, time step at different places\n* joined CDF --> joint CDF\n* P5, there is an $\\epsilon_t$, which should probably be $\\alpha_t$ or $1-\\alpha_t$\n* Citations are sometimes inconsistent. Format Name (year) should be when the name is supposed to be read, while (name, year) should be used when the name should not be read. P5, Ruschendorf should not be read, hence (name, year) (\\citep) should be used. This happens in other places. \n* P7: forgotten parenthesis in the definition of coverage$_{1-\\alpha}$\n\n-: in the case of auto-regressive, it is unclear how is treated the newly produced point in order to produce further predictions. In particular, since the newly predicted point is a conformal sets, do authors run conformal prediction using weak-labels/intervals (as in \"Cauchois, M., Gupta, S., Ali, A., & Duchi, J. (2022). Predictive Inference with Weak Supervision.\"), or do they simply use the new point-wise prediction? What guarantees that the previous validity guarantees with respect to the ground-truth will hold (as the prediction is not the true value)? Also, how does one obtain the cumulative distribution for the non-conformity scores of the k+1 prediction? Does one use the same cumulative distribution, and if this is the case, why do we need to re-run an optimisation? \n\n-: While the optimisation program introduced to find suitable time-step wise confidence degrees seems fine, it is unclear what we can/should expect from it? In particular, shall we in practice have very similar confidence degrees or very imbalanced ones? What justify to have different confidence degrees at each run? Why not optimizing to minimize the obtained volume? How does it compare to the constant time dichotomic search assuming all confidence degrees being equal? ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject