Conformal PID Control for Time Series Prediction
Authors: Anastasios Nikolas Angelopoulos
Emmanuel Candes
Ryan Tibshirani
Venue: NeurIPS 2023 poster (2023)
Keywords: conformal prediction
time series
uncertainty quantification
distribution shift
TL;DR: conformal prediction for time series that anticipates and forecasts future distribution shifts

Abstract: We study the problem of uncertainty quantification for time series prediction, with the goal of providing easy-to-use  algorithms with formal guarantees. The algorithms we present build upon ideas from conformal prediction and control theory, are able to prospectively model conformal scores in an online setting, and adapt to the presence of systematic errors due to seasonality, trends, and general distribution shifts. Our theory both simplifies and strengthens existing analyses in online conformal prediction. Experiments on 4-week-ahead forecasting of statewide COVID-19 death counts in the U.S. show an improvement in coverage over the ensemble forecaster used in
official CDC communications. We also run experiments on predicting electricity demand, market returns, and temperature using autoregressive, Theta, Prophet, and Transformer models. We provide an extendable codebase for testing our methods and for the integration of new algorithms, data sets, and forecasting rules at [this link](http://github.com/aangelopoulos/conformal-time-series).

Conformal PID Control for Time Series Prediction
Anastasios N. Angelopoulos University of California, Berkeley mggeloppulls@berkklleyy edu
Emmanuel J. Candès Stanford University candes@stanfore edu
Ryan J. Tibshirani University of California, Berkeley ryantibs@berkelee .dd edu
Abstract
We study the problem of uncertainty quantification for time series prediction, with the goal of providing easy-to-use algorithms with formal guarantees. The algo- rithms we present build upon ideas from conformal prediction and control theory, are able to prospectively model conformal scores an online setting, and adapt to the presence of systematic errors due seasonality, trends, and general distribution shifts. Our theory both simplifies and strengthens existing analyses in online confor- mal prediction. Experiments 4-week-ahead forecasting of statewide COVID-19 death counts in the U.S. show an improvement in coverage over the ensemble forecaster used in official CDC communications. We also experiments predicting electricity demand, market returns, and temperature using autoregres- sive. Theta, Prophet, and Transformer models. We provide an extendable codebase for testing our methods and for the integration of new algorithms, data sets, and forecasting rules
Introduction
Machine learning models run in production systems regularly encounter data distributions that change over time. This can due to factors such as seasonality and time-of-day, continual updating and re- training of upstream machine learning models, changing user behaviors, and These distribution shifts can degrade a model's predictive performance. They also invalidate standard techniques for uncertainty quantification, such as conformal prediction 35).
To address the problem of shifting distributions, we consider the task of prediction an adversarial online setting, as in [16]. In this problem setting, we observe a (potentially) adversarial time series of deterministic covariates 2, E X and responses y/ V. for & E are free to define any conformal score foreefunniin 8; X x ) - IR, which we can view as measuring the accuracy our forecast time L. We will assume with loss of generality that is negatively oriented (lower values mean greater forecast accuracy). For example, we may use the absolute error st (1 y) = ly f((z),, where fi is a forecaster trained on data up to but not including data at time t.
{1, 2, 3 ......,, As in standard
conformal prediction,
The challenge in the sequential setting is as follows. We seek to invert the score function to construct a conformal prediction set,
C. = {y E y: S1 y) - 9.)
(1)
where qt is an estimated 1 a-a quantile for the distribution of the score (x Vr.) at time f. In standard conformal prediction, we would take q, be level 1 sample quantile (up finite-sample correction) of s, s, B.). i < L; if the data sequence (I,, y1). i E N were i.i.d. or exchangeable, then this would yield 1 coverage 335 at each time However, in the sequential setting, which does not assume exchangeability (or any probabilistic model for the data for that matter), choosing in (T) yield coverage a formidable task. In fact, if are not willing to make any assumptions about the data sequence, then a coverage guarantee at time t would only be possible with trivial methods, which construct prediction intervals of infinite sizes.
httpe///ggtthh..ccm/mggggggggoos time-series
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Therefore, our goal is to achieve long-rum run coverage in time. That is, letting ert, = 1 for $ Ct). we would like to achieve, for large integers
1 T Σ
err, err, = =
under few assumptions, where o(1) denotes a quantity that tends to zero as T - 0c. We note that (2) is not probabilistic at all, and every theoretical statement we will make this paper holds deterministically. Furthermore, going beyond (2), also seek to design flexible strategies to produce the sharpest prediction sets possible, which not only adapt to, but also anticipate distribution shifts.
We call our proposed solution conformal PID control. It treats the system for producing prediction (PID) controller. In the language of control, the prediction sets take a control variable, Gt, and then produce process variable, ert, We seek to anchor err, to a set point, a. To do so, we apply corrections to g based on the error of the output, 9r = ert, - a. By reframing the problem in this language, we are able to build algorithms that have more stable coverage while also prospectively adapting to changes the score sequence, much in the same style as control system. See the diagram in Figure [
sets as a cyppotional-integaa---eeee
79t
quantiie tracker
a
gt
U++l
Σ
I
<>.9
g.)
yt-1 & C++1
integrattor
D
g
eIT,
eff
scorecaster
Figure Conformal PID Control, expressed as block diagram.
1.1
Peek results: methods
Three design principles underlie our methods:
Quantile tracking (P control). Running online gradient descent on the quantile loss (summed over all past scores) gives rise to method that we call quantile tracking, which achieves long-run coverage (2) under no assumptions except boundedness of the scores. This bound can be unknown. Unlike adaptive conformal inference (ACI) [16), quantile tracking does not return infinite sets after sequence of miscoverage events. This can be seen as equivalent to proportional (P) control.
2. Error integration (I control). By incorporating the running sum E (err; a:) of the coverage errors into the online quantile updates, we can further stabilize the coverage. This error integration scheme achieves long-run coverage (2) under no assumptions whatsoever on the scores (they can be unbounded). This can be seen as equivalent to integral (I) control.
3. Scorecasting (D control). To account for systematic trends in the scores this may be due to aspects of the data distribution, fixed or changing, which are not captured by the initial forecaster- train second model, scorecaster, to predict the quantile of the next score. While quantile tracking and error integration are merely reactive, scorecasting is forward-looking systematic trends in the errors and lead practical advantages in terms of coverage and efficiency (set sizes). This can be seen as equivalent to derivative (D) control. Traditional control theory would suggest using linear approximation of 9t It-1, but in our problem, will typically choose more advanced scorecasting algorithms that go well beyond simple difference schemes.
can residualize out
These three modules combine to make our final iteration, the conformal PID controller:
4++1
19tt + Tr
gt
+ gl
(3)
D

---

In traditional PID control, one would take rz(z) to be linear function of r. Here, we allow for nonlinearity and take r, to be saturation function obeying
(4) for constants c > O, and a sublinear, nonnegative, nondecreasing function h-we call function h satisfying these conditions admissible. An example the tangent integrator (xr) K1 tan( i log(t)/(tssm), where we set tan(z) = sign(x)) oo for x & I-T/22, x /2,, and Csat: K1 are constants. The choice of integrator r, design decision for the user, as is the choice of scorecaster g/.
I = h(t)
re(x) 2 b, and I x -c h(t) r((x)
We find it convenient to reparametrize (3).t to produce sequence of quantile estimates qt, & N used in the prediction sets ① as follows:
let q be any function of the past: Is. yi, qs for i S t,
(5)
then update 4++1 = Gc+1 + rr. Σ err, a)
Taking 9:+1 19. + g. recovers but find it generally useful to instead consider the formu- lation in (5) which will be our main focus in the exposition henceforth. Now we view gt+1 as the scorecaster, which directly predicts 4 using past data. A main result of this paper, whose proof is given in Appendix A is that the conformal PID controller (5) yields long-run coverage for any choice of integrator Te that satisfies the appropriate saturation condition, and any scorecaster âc+1-
Theorem 1. Let { heN be any sequence of numbers in b- b/2] and let {st be any sequence of score functions with outputs in b-b/2. 222 b/2). Here b > O, and may be infinite. Assume that rt satisfies for an functiio h. Then the iterations in (5) achieve long-run coverage (
To emphasize, this result holds deterministically. with no probabilistic model on the data (zt, ye)). t E N. (Thus in the case that the sequence is random, the result holds for all realizations of the random variables.) As we will soon see, this theorem can be seen as a generalization of existing results in the online conformal literature.
1.2 Peek results: experiments
COVID-19 death forecasting. To demonstrate conformal PID in practice, we consider 4-weck-ahead forecasting of COVID-19 deaths in California, from late 2020 through late 2022. The base forecaster ft that we use is the ensemble model from the COVID-19 Forecast Hub, which the model used for official CDC communications each time t we actually seek to predict the observed death count 31+4 at time t + 4.
COVID-19 forecasting [I0] 29]. In this forecasting problem, at
Figure 2 shows the central 80% prediction sets from the Forecast Hub ensemble model on the left panel, and those from our conformal PID method on the right. We use a quantile conformal score function, as in conformalized quantile regression [30], applied asymmetrically (i.c., separately) the lower and upper quantile levels). We use the tan integrator, with constants chosen heuristically (as described in Appendix C). and f1-regularized quantile regression as the scorecaster- particular, the scorecasting model at time predicts the quantile of the score at time + based on all previous forecasts, cases, and deaths, from all 50 US states. The main takeaway is that conformal PID control is able to correct for consistent underprediction of deaths in the winter wave of late 2020/early 2021. We can from the figure that the original ensemble fails to cover 8 times in stretch of 10 weeks, resulting in a coverage of 20%; meanwhile, conformal PID only fails to cover 3 times during this stretch, restoring the coverage to 70% (recall the nominal level is 80%).
How is this possible? The ensemble is mainly comprised of constituent forecasters that ignore geographic dependencies between states [II] for the sake of simplicity computational tractability. But COVID infections and deaths exhibit strong spatiotemporal dependencies, and most US states experienced the winter wave of late 2020/early 2021 at similar points in time. The scorecaster is thus able to learn from the mistakes made on other US states in order to prospectively adjust the ensemble's 's forecasts for the state of Califomia. Similar improvements can be seen for other states, and we include experiments for New York and Texas as examples in Appendix which also gives more details on the scorecaster and the results.
Electricity demand forecasting. Next we consider a data set on electricity demand forecasting in New South Wales [18), which includes half-hourly data from May 7, 1996 to December 5, 1998.
3

---

Base Forecaster
Conformal PID Control
Coverage
50%
6000
4000
2000
2020-00
2071-01
2011-34
2023-87
202-10
2022-01
2023.04
2032
2020-10
203361
20
2021-07
202-10
2072-07
2023
2022-77
Time
Figure 2: Results for 4-week ahead COVID-19 death forecasting in California. The left column shows the COVID-19 Forecast Hub ensemble model, and the right column shows conformal PID control using the tan integrator, and scorecaster given by l fi-penalized quantile regression past forecasts, cases, and deaths from all 50 states. The top row plots coverage, averaged over a trailing window of 10 weeks. The nominal coverage level is 1 line. The bottom row plots the prediction sets in gold along with the ground-trute truth times series (death counts). Miscoverage events are indicated by red dots. Summary statisties such the coverage and average set size are available Table
all
0.8 and marked by gray dotted
For the base forecaster we use Transformer model 34 as implemented in darts 19]. This is only re-trained daily, to predict the entire day's demand in one batch; this is a standard approach with Transformer models due to their high computational cost. For the conformal score, we use the asymmetric (signed) residual score. We use the tan integrator as before, and we use a lightweight Theta model re-trained at every time point (half-hour), as the scorecaster.
ACI
Conformal PID Control
Coverage
A of man Wiitt
75%
1.0
1996-94133
1996 296-0522
2006 05.23
1938 20-00110 09--03329
200s r0s6-eal deal
1938 o003
1998 -te--
1996 20-22
1096 1897-0177
1905 13
1996 45547
25.22
1906 25
1996 cl-7
1996 deus
1996 -e--ddoe
1996 26
1996
1908
1996
Time
Figure 3: Results for electricity demand forecasting. The left column shows adaptive conformal inference (ACT), and the right column shows conformal PID control. The base forecaster Transformer model, and we use tan integrator and a Theta scorecaster. The format of the figure follows that of Figure 2, except the nominal coverage now 1 a 0.9, and the coverage is averaged over trailing window of 50 points (we also omit the red dots which mark miscoverage events). Summary statistics are available in Table 2
The results are shown in the right panel of Figure 3 where adaptive conformal inference (ACI) [6 also compared in the left panel. In short, conformal PID control is able anticipate intraday variations in the scores, and produces sets that "hug" the ground truth sequence tightly; it achieves tight coverage without generating excessively large or infinite sets. The main reason why this improved is that the scorecaster has seasonality component built into its prediction model; in general, large improvements such as the one exhibited Figure should only expected when the base forecaster is imperfect, is the case here.

---

1.3 Related work
The adversarial online view of conformal prediction was pioneered by [66 in the same paper that first introduced ACI. Since then, there has been significant work towards improving ACI, primarily by setting the learning rate adaptively [T7 40 77 and incorporating ideas from multicalibration to improve conditional coverage [ It is worth noting that 7 also makes the observation that the ACI iteration can be generalized track the quantile of the score sequence, although their focus is on adaptive regret guarantees. Because the topic of adaptive learning rates for ACI and related algorithms has already been investigated heavily, we do not consider in the current paper. Any such method, such as those of [7 7
should work well in conjunction with proposed algorithms.
related but distinct line of work surrounds online calibration in the adversarial sequence model, which dates back to 14 [51, and connects in interesting ways to both game theory and online learning. We will not attempt to provide comprehensive review of this rich and sizeable literature, but simply as a few interesting examples of recent work.
highlight 25 24 23
Lastly, outside the online setting, note that several researchers have been interested in generalizing conformal prediction beyond the i.i.d. (or exchangeable) data setting: this includes 33 28 26
12 g
and for time series prediction, in particular, 9 31 39 3) The focus of all of these papers is quite
different, and they all rely on probabilistic assumptions on the data sequence to achieve validity; we include further discussion in Appendix B
2 Methods
We describe the main components of our proposal one at a time, beginning with the quantile tracker.
2.1 Quantile tracking
The starting point for quantile tracking is to consider the following optimization problem:
T
minimize 9EB
Σ
Pi-a - q).
(6)
for large T, where we abbreviate s, = (x vr) for the score of the test point, and P1-a denotes the
quantile loss at the level a, i.e., Pr (2) = T/2| for >> 0 and (1 r)]2| for 2 S 0. The latter the standard loss used in quantile regression [22 21). Problem (6) is thus simple convex (linear) program that tracks the 1 a quantile of the score sequence st, t E N. To see this, recall that for
continuously distributed random variable Z, the expected loss B[pi-a -24) q)] is uniquely minimized (over E IR) at the level 1 - a quantile of the distribution of Z.
In the sequential setting, where we receive one score at a time, a natural and simple approach is to apply online gradient descent (6). with a constant learning rate n O. This results in the update
g++1 = Qt + V/ ge) 9: + 1)(err, - a),
where the second line follows Vpi-a(St 91) 1 if 8t
9t eIr, 1, and Vp - qt)= = -a if s; I 9t < eIT, = O. Note that the update in 7 is highly intuitive: miscovered (committed
error) at the last iteration then increase the quantile, whereas if we
covered (did not commit an error) then we decrease the quantile.
Even though it is extremely simple, the quantile tracking iteration can achieve long-run coverage own its own, provided the scores are bounded.
Proposition 1. Les (s: hen be any sequence of numbers in b.. < b < Then the quantile tracking iteration 7 satisfies
T
b+ b+n nT
Σ err, a)
for any learning rate 1) > 0 and T 2 1. In particular, yields long-run coverage as in (2)
Teechnically, this is the online subgradient method; in slight abuse of notation, we write Vpr (0) to denote subgradient of pi-a O which can take any value aa,, -a).
5

---

A few remarks are in order. First, although Proposition 1 assumes boundedness of the scores, we do not need to know this bound in order to run 7 and obtain long-run coverage. As long as the scores lie in [-b, bj for any finite b, the guarantee goes through- agnostically and performs the same updates in any case. Notably, the adaptive conformal inference algorithm can be expressed as special case of the quantile tracker; see Appendix B.1 for details.
the quantile tracker proceeds
Second, for the learning rate, in practice we typically set v) heuristically, as some fraction of the highest score over trailing window B, = ma usually gives good results, and we use it in all experiments unless specified otherwise (we also set the window length A to be the same as the length of the burn-in period for training the initial base forecaster and scorecaster) R Extremely high learning rates result in volatile sets, while very low ones may fail to keep up with rapid changes in the score distribution.
5t On this scale, setting O.1
2.2 Error integration
Error integration is generalization of quantile tracking that follows the iteration:
4++1 Te Σ (err; - a1) )
(8)
where ri is a saturation function that satisfies 4 for an admissible function h; recall that we use admissible to mean nonnegative, nondecreasing, and sublinear. As we saw in [33), the quantile tracker uses a constant threshold function h. whereas h is now permitted to grow, as long as it grows sublinearly, i.c., h(t)/t 0 as + OC. A non-constant threshold function h can be desirable because it means that r. will "saturate" (will hit the conditions on the right-hand sides in less often, so corrections for coverage error will occur less often, and in this sense, a greater degree of coverage error can be tolerated along the sequence.
The next proposition, in particular its proof, makes the role of h precise. Importantly, Proposition suffices to prove Theorem
Proposition 2. Let (EE)IeN be any sequence of numbers in [-b, b). where b > O, and may be infinite. Assume that satisfies for an admissible function h. Then the error integration iteration (8) satisfies
1 T
h(T) -1
eIT, a)
(9)
T
for any T 1, where is the constant in In particular, this means
yields long-rut coverage
())
The choice of saturation function essentially corresponds to a choice of adaptive learning rate; see Appendix Dorendiifforrro details.
2.3 Scorecasting
The final piece to discuss is scorecasting. scorecaster attempts to forecast g++1 directly, taking advantage of any leftover signal that is not captured by the base forecaster. This is the role played by 4:+1 in (5). Scorecasting may be particularly useful when it difficult to modify or re-train the base forecaster. This can occur when the base forecaster is computationally costly to train (e.g., in Transformer model); or it can occur complex operational prediction pipelines where frequently updating forecasting implementation is infeasible. Another scenario where scorecasting may be useful is one in which the forecaster and scorecaster have access to different levels of data. For example, if public health agency collects epidemic forecasts from external groups, and forms an ensemble forecast, then the agency may have access to finer- teer-ggriiedd data that it can use to recalibrate the ensemble's 's prediction sets (compared the level of data granularity granted to the forecasters originally).
This motivates the need for scorecasting modular layer that "sits on top" of the base forecaster and residualizes out systematic errors in the score distribution. This intuition is made more precise
Teechniclll, this learning rate is not fixed, so Proposition does not directly apply. However, we can view it as special case of error integration and an application of Proposition] 2 thus provides the relevant coverage guarantee.

---

by recalling, as described above (following Proposition 2), that scorecasting combined with error integration as in (5) is just a reparameterization of error integration (8), where q q: q and s 2t are the new quantile and new score, respectively. A well-executed scorecaster could reduce the variability in the scores and make them more exchangeable, resulting more stable coverage and tighter prediction sets, as seen in Figure 3 On the other hand, an aggressive scorecaster with little or no signal can actually hurt by adding variance to the new score sequence s/ which could result in more volatile coverage and larger sets.
There is no limit to what can choose for the scorecasting model. We might like use model that can simultaneously incorporate seasonality, trends, and exogenous covariates. Common choices would be SARIMA (seasonal autoregressive integrated moving average) and ETS (error-trend- seasonality) models, but there are many other available methods, such as the Theta model [2]. Prophet model 32 and neural network forecasters; see (20 for a review.
2.4 Putting it all together
Briefly, we revisit the PID perspective, to recap how quantile tracking, error integration, and score- casting fit in and work in combination. It helps to return to 33. which copy again here:
a++1 = g nnlerr, a) + Σ CIT; a)
(10)
Quantile tracking is precisely given by taking g/. = Oe and r; = O. This can seen as equivalent to control: subtract de from both sides in (TO) and treat the increment u4+1 9t+11 - qi the process variable; then in this modified system, quantile tracking is exactly P control. For this reason, we use "conformal P control" to refer to the quantile tracker in the experiments that follow. Similarly, "conformal PI control" to refer to the choice g/ gt. and rt F 0 as a generic integrator (for us, tan is the default). Lastly, "conformal PID control" refers to letting of be generic scorecaster, and r, * be a generic integrator.
Experiments
In addition to the statewide COVID-19 death forecasting experiment described in the introduction, we run experiments on all combinations of the following data sets and forecasters.
Data sets:
Forccasters (all via darts [1]]
Electricity demand 18
Autoregressive (AR) model with 3 lags
Return (log price) of Amazon, Google, and Microsoft stock (27
Theta model with 0 = 2
Prophet model 32
Temperature in Delhi 37
Transformer model
34
In all cases except for the COVID-19 forecasting data set, we: re-train the base forecaster each time point; construct prediction sets using the asymmetric (signed) residual score; and use Theta model for the scorecaster. For the COVID-19 forecasting setting, we: use the given ensemble model as the base forecaster (no training at all); construct prediction sets using the asymmetric quantile score; and use an l di-penalized quantile regression as the scorecaster. fit on features derived from previous forecasts, cases, and deaths, as described in the introduction. And lastly, in all cases, we use a tan function for the integrator with constants chosen heuristically, as described in Appendx C
The results that we choose to show in the subsections below are meant to illustrate key conceptual points (differences between the methods). Additional results presented in Appendix Our GitHub repository, https:///gi //githu//.t/...... cctmm//anggelooolloooooooooooooo http://alooooomooooooooooooolllm time--ime------------- series provides the full suite of evaluations.
3.1 ACI versus quantile tracking
We forecast the daily Amazon (AMZN) opening stock price from 2006-2014. We this in log-space (hence predicting the return of the stock). Figure compares ACI and the quantile tracker, each
7

---

ACI
Conformal Control
Coverage
1000
Setts
500
2010
2014
2018
2008
2012
2016
2018
Time
AR
Prophet
Theta
Transformer
ACI
P Ctrl
ACI
P Ctll
ACI
P Ctrl
ACI
P Cull
Marginal coverage
0.894
0.894
0.904
0.888
0.895
0.894
0.906
0.887
Longest err sequence
6
3
13
6
5
21
9
Average set size Median size
00
17.6
00
51.7
06
17.8
70.4
13.4
50.7
13.1
61.9
75% quantile set size
28.4
22.3
117
72.1
20.9
22.6
179
98.4
90% quantile set size
44.1
37.7
236
114
38.6
38.4
302
153
95% quantile set size
49.9
46.2
140
48.2
46.9
196
Figure 4: Results for forecasting Amazon stock return, comparing ACI and quantile tracking (P control). The plots show AR the base forecaster; the table summarizes the results of all four base forecasters. We use the default learning rates for both ACI and quantile tracking: 7) = 0.005 and r 0.1 respectively.
ACI
Conformal Control
100%
Coverage
wo hwo jur an H
50%
1000
500
2302
2012
2014
2018
2014
2010
2017
2014
2010
Time
AR
Prophet
Theta
Transformer
ACI
P
Ctrl
ACI
P Ctrl
ACI
Ctrl
ACI
P Ctrl
Marginal coverage
0,9
0.898
0.901
0.897
0,9
0.896
0.902
0.897
Longest sequence
2
2
3
2
7
3
Average set size Median set size
28.4
42.6
27.8
54.2
00
18.6
60.6
31.9
00
18.7
35.5
75% quantile set size
00
37.4
0c
56.5
00
37.8
00
69.7
90% quantile set size 0c
66.3
93.4
0C
63.4
20
123
95% quantile set size
81.8
116
78.5
164
Figure 5; As in Figure 4 but with larger learning rates for ACI and quantile tracking: = 0.1 and 7 0.5B, respectively.

---

with its default learning rate: 1) = 0.005 for ACI, and 0.1B, for quantile tracking. We see that the coverage from each method decent, but oscillates nontrivially around the nominal level of 1 = 0.9 (with ACI generally having larger oscillations). Figure thuus increases the learning rate for each method: "j = 0.1 for ACI, and 1) = 0.5 B, for the quantile tracker. We now see that both deliver very tight coverage. However, ACT does so by frequently returning infinite sets, meanwhile, the corrections to the sets made by the quantile tracker are nowhere near as aggressive.
As final comparison, in Appendix from ever being infinite. This heuristic may be used by practitioners that want to guard against infinite sets, but it no longer has a validity guarantee for bounded or unbounded scores. The results in the appendix indicate that the quantile tracker has similar coverage this procedure, and usually with smaller sets.
we modify ACI to clip the sets in way that disallows them
3.2 The effect of integration
Next we forecast the daily Google (GOOGL) opening stock price from 2006-2014 (again done in log-space). Figure S compares the quantile tracker without and without additional integrator component (P control versus control). We purposely choose very small learning rate, 7) = 0.01 B₁, in order to show how the integrator series. The coverage of PI control begins to oscillate more towards the end the sequence, which we attribute at least in part to the fact that the integrator measures coverage errors accumulated over all time-and by the end of long sequence, the marginal coverage can still be close to 1 a even if the local coverage deviates more wildly. This be addressed by using local version of the integrator, an idea we return to in the discussion.
stabilize coverage, which it does nicely for most of the time
Conformal control
Conformal Control
Coverage
50%
Sets
2008
2010
2012
2014
2010
2018
2008
201e
2012
2012
2016
2018
Time
AR
Prophet
Theta
Transformer Ctrl
P Ctrl
PI Ctrl
P Ctrl
PI Ctrl
P Ctrl
PI Ctrl
PI Ctrl
Marginal coverage
0.895
0.896
0.882
0.888
0.891
0.896
0.817
0.897
Longest сr sequence
6
31
13
58
6
106
10
Average set size Median size
19.2
24.7
95.6
68.2
42.9
33.5
281
108
16.4
22.8
91.5
60.7
44.3
32.8
219
73.8
75% quantile set size
21.3
30.1
120
42.6
386
121
90% quantile set size
30.7
38.4
150
118
66.2
51.3
482
244
95% quantile set size
32.8
44.7
156
139
69.1
57.5
524
341
Figure 6: Results for forecasting Google stock return, comparing quantile tracking with and without the integrator (P control versus PI control). The plots show Prophet the base forecaster; the table summarizes the results of all four base forecasters. We purposely use a very small learning rate, 0.01 B,, in order to show how the integrator can stabilize coverage.
3.3
The effect of scorecasting
Figures aanddalreade Jarredy showcase examples which scorecasting offers significant improvement in coverage and set sizes. Recall that these were settings in which the base forecaster produces errors (scores) that have predictable trends. Further examples in the COVID-19 forecasting setting, which display similar benefits to scorecasting, are given in Appendix F

---

We emphasize that it is not always the case that scorecasting will help. In some settings, scorecasting may introduce enough variance into the new score sequence that the coverage or sets will degrade in stability. (For example, this will happen if we run a highly complex scorecaster on a sequence of i.i.d. scores, where there are no trends whatsoever.) In practice, scorecasters should be designed with care, just as one would design base forecaster; it is unlikely that using "out of the box" techniques for scorecasting will be robust enough, especially in high-stakes problems. A p examples in which scorecasting, run across all settings using generic Theta model, can hurt (for example, it adds noticeable variance to the coverage and sets in some instances within the Amazon data setting).
Discussion and extensions
Discussion. Our work presents a framework for constructing prediction sets in time series that is analogous (and indeed formally equivalent) to PID control, consisting of quantile tracking (P control), which is simply online gradient descent applied to the quantile loss; error integration (I control) to stabilize coverage; and scorecasting (D control) to remove systematic trends in the scores (errors made by the base forecaster).
We found that the combination of quantile tracking and integration consistently yields robust and favorable performance our experiments. Scorecasting provides additional benefits if there are trends left in scores that are predictable (and the scorecaster is well-designed), as is the case in some of our examples. Otherwise, scorecasting may add variability and make the coverage and prediction sets more volatile. Overall, designing the scorecaster (which includes the choice to even use one at all) is an important modeling step, just like the design of the base forecaster.
It is worth emphasizing that, with the exception of the COVID-19 forecasting example, our experi- ments are intended to be illustrative, and we did not look to use state-of of-the-art forecasters, or include any and all possibly relevant features for prediction. Further, while we found that using heuristics to set constants (such as the learning rate 17, and constants Csat, K₁ for the tan integrator) worked decently well, we believe that more rigorous techniques, along the lines of [1 7 tune these adaptively an online fashion.
can be used to
In this for all and (0, y) = for all The goal is to bound the deviation of the average risk t= L. LC Bt) from a. We state a result for the integrator below, and give its proof in Appendix A
Extensions. We now present an extension of our analysis to conformal risk control 1
6 13
setting, we are given sequence of loss functions L; 2n x x - [0, 1] satisfying L,
Proposition 3. Consider the iteration qi+1 = Ei that r, satisfies (4).fo an admissible function h. Also assume that C, ye) = 0 if S -b bndd y if q: 2 b, where b > O, and may be infinite. Then for all T Z 1,
vi) -a)), with L, as above. Assume
1 T
c. (T) +1
(L.
- a) <
(11)
T
for any T ' 1, where is the constant in
We briefly conclude by mentioning that we believe many other extensions are possible, especially with respect to the integrator. Broadly, we can choose to integrate in a kernel-weighted fashion,
Tt Σ [err; a) K ((i, (1, (1, (1, vil)
(12)
As special case, the kernel could simply assign weight I lift-i w, and weight otherwise, which would result in an integrator that aggregates coverage over a trailing window of length w. This can help consistently sustain better local coverage, for long sequences. As another special case, the kernel could assign weight based on whether x;; and r, lie in the same bin in some pre-defined binning of X' space, which may be useful for problems with group structure (where we want group-wise coverage). Various other choices and forms of kernels are possible, and it would be interesting to consider adding together a number of such choices (12) in combination, in a multi-resolution flavor, for the ultimate quantile update.
10

---

Acknowledgments
We would like to thank Tijana Zrnić, Amit Kohli, and Jordan Lekeufack for their valuable feedback. This work was supported by the National Science Foundation (NSF) Graduate Research Fellowship Program under grant no. 2146752, and the Office Naval Research (ONR) Multi- Univeritiv Research Initiative (MURI) Program under grant no.
References
[1] Anastasios N. Angelopoulos, Stephen Bates, Adam Fisch, Lihua Lei, and Tal Schuster. Confor- mal risk control. arXiv preprint arXiv:2208. 02804, 2022.
[2] Vassilis Assimakopoulos and Konstantinos Nikolopoulos. The theta model: A decomposition
approach to forecasting. International Journal of Forecasting, 16(4):521-330, 2000.
[3] Andreas Auer, Martin Gauch, Daniel Klotz, and Sepp Hochreiter. Conformal prediction for time series with modern Hopfield networks. arXiv preprint arXiv:2303. 12783, 2023.
[4] Rina Foygel Barber, Emmanuel J. Candès, Aaditya Ramdas, and Ryan J. Tibshirani. Conformal prediction beyond exchangeability, arXiv preprint arXiv:2202. 1415, 2022.
[5] Osbert Bastani, Varun Gupta, Christopher Jung. Georgy Noarov, Ramya Ramalingam, and Aaron Roth. Practical adversarial multivalid conformal prediction. In Advances in Neural Information Processing Systems, 2022.
[6] Stephen Bates, Anastasios Angelopoulos, Lihua Lei, Jitendra Malik, and Michael Jordan. Distribution-free, risk-controllingg prediction sets. Journal of the ACM, 68(6);1-34, 2021.
[7] Aadyot Bhatnagar, Huan Wang, Caiming Xiong, and Yu Bai. Improved online conformal prediction via strongly adaptive online learning. arXiv preprint arXiv:2302.07869, 2023.
[8] Emmanuel J. Candès, Lihua Lei, and Zhimei Ren. Conformalized survival analysis. Journal of the Royal Statistical Sociery: Series B, 85(1):24 2023.
[9] Victor Chernozhukov, Kaspar Wuthrich, and Zhu Yinchu. Exact and robust conformal inference methods for predictive machine learning with dependent data. In Proceedings of the Annual Conference Learning Theory, 2018
[10] Estee Cramer, Yuxin Huang, Yijin Wang, Evan L. Ray, Matthew Cornell, Johannes Bracher, Andrea Brennen, Alvaro J. Castro Rivadeneira, Aaron Gerding, Katie House, Dasuni Jayawar- dena, Abdul Hannan Kanji, Ayush Khandelwal, Khoa Le, Vidhi Mody, Vrushti Mody, Jarad Niemi, Ariane Stark, Apurv Shah, Nutcha Wattanchit, Martha W. Zorn, Nicholas G. Reich, and US COVID-19 Forecast Hub Consortium. The United States COVID-19 Forecast Hub dataset. Scientific Data, 9, 2022.
[11] Estee Y. Cramer, Evan L. Ray. Velma K. Lopez, Johannes Bracher, Andrea Brennen, Alvaro J. Castro Rivadeneira, Aaron Gerding, Tilmann Gneiting, Katie H. House, Yuxin Huang, Dasuni Jayawardena, Abdul H. Kanji, Ayush Khandelwal, Khoa Le, Anja Mühlemann, Jarad Niemi, Apury Shah, Ariane Stark, Yijin Wang, Nutcha Wattanachit, Martha W. Zorn, Youyang Gu, Sansiddh Jain, Nayana Bannur, Ayush Deva, Mihir Kulkarni. Srujana Merugu, Alpan Raval, Siddhant Shingi, Avtansh Tiwari, Jerome White, Spencer Woody, Maytal Dahan, Spencer Fox, Kelly Gaither, Michael Lachmann, Lauren Ancel Meyers, James G. Scott, Mauricio Tec, Ajitesh Srivastava, et al. Evaluation of individual and ensemble probabilistic forecasts of COVID-19 mortality in the United States. Proceedings of the National Academy of Sciences, 119(15)ee2 136221133111111 119, 2022.
[12] Clara Fannjiang, Stephen Bates, Anastasios N. Angelopoulos, Jennifer Listgarten, and Michael I. Jordan. Conformal prediction under feedback covariate shift for biomolecular design. Proceed- ings of the National Academy of Sciences, 119(43):c2204566111 19, 2022.
[13] Shai Feldman, Liran Ringel, Stephen Bates, and Yaniv Romano. Achieving risk control in online learning settings. arXiv preprint arXiv:2205. 09095, 2022.
[14] Dean Foster. proof of calibration via Blackwell's approachability theorem. Games and Economic Behavior, 29(1-2):73-78, 1999.
[15] Dean P. Foster and Rakesh V. Vohra. Asymptotic calibration. Biometrika, 85(2):379-33 -2200,9--99, 1998.

---

[16] Isaac Gibbs and Emmanuel J. Candès. Adaptive conformal inference under distribution shift. In Advances in Neural Information Processing Systems, 2021.
[17] Isaac Gibbs and Emmanuel J. Candès. Conformal inference for online prediction with arbitrary distribution shifts. arXiv preprint arXiv:2208.0844 2022.
[18] Michael Harries. Splice-2 comparative evaluation: Electricity pricing. Technical report. University of New South Wales, 1999.
[19]
Julien Herzen, Francesco Lässig. Samuele Giuliano Piazzetta, Thomas Neuer, Léo Tafti, Guil- laume Raille, Tomas Van Pottelbergh, Marek Pasieka, Andrzej Skrodzki, Nicolas Huguenin, Maxime Dumonal, Jan Košcisz, Dennis Bader, Frédérick Gusset, Mounir Benheddi, Camila Williamson, Michal Kosinski, Matej Petrik, and Gaël Grosch. Darts: User-friendly modern machine learning for time series. Journal of Machine Learning Research, 23(1):5442-5447, 2022.
[20] Rob J. Hyndman and George Athanasopoulos. Forecasting: Principles and Practice. OTexts, 2018.
[21] Roger Koenker. Quantile Regression. Cambridge University Press, 2005.
[22] Roger Koenker and Gilbert Bassett. Regression quantiles. Econometricu, 46(1):33-50, 1978.
[23]
Volodymyr Kuleshov and Shachi Deshpande. Online calibrated regression for adversarially robust forecasting. arXiv preprint arXiv:2302.121199 2023.
[24] Volodymyr Kuleshov and Stefano Ermon. Estimating uncertainty online against an adversary. In Association for the Advancement of Artificial Intelligence, 2017.
[25] Volodymyr Kuleshov and Percy Liang. Calibrated structure prediction. In Advances in Neural Information Processing Systems, 2015.
[26] Lihua Lei and Emmanuel J. Candès. Conformal inference of counterfactuals and individual treatment effects. Journal of the Royal Statistical Society: Series B, 83(5):911-938 2021.
[27] Cam Nguyen. S&P 500 stock data. http:://www....ggllllee.. htt.a//glle..cmmtttttttttsaaaaaa sandp500 2018.
[28]
Aleksandr Podkopaev and Aaditya Ramdas. Distribution-ffes uncertainty quantification for classification under label shift. In Uncertainty in Artificial Intelligence, 2021.
[29] Evan L. Ray, Logan C. Brooks, Jacob Bien, Matthew Biggerstaff, Nikos I. Bosse, Johannes Bracher, Estee Y. Cramer, Sehastian Funk, Aaron Gerding, Michael A. Johansson, Aaron Rumack, Yijin Wang, Martha Zorn, Ryan J. Tibshirani, and Nicholas G. Reich. Comparing trained and untrained probabilistic ensemble forecasts of COVID-19 cases and deaths in the United States. International Journal of Forecasting, 3963:: 2023.
[30] Yaniv Romano, Evan Patterson, and Emmanuel J. Candès, Conformalized quantile regression. In Advances in Neural Information Processing Systems, 2019.
[31] Kamile Stankeviciute, Ahmed M. Alaa, and Mihaela van der Schaar. Conformal time-series forecasting. Advances in Neural Information Processing Systems, 2021.
[32] Sean J. Taylor and Benjamin Letham. Forecasting at scale. The American Statistician, 72(1):37-
45, 2018.
[33] Ryan J. Tibshirani, Rina Foygel Barber, Emmanuel J. Candès, and Aaditya Ramdas. Conformal prediction under covariate shift. In Advances in Neural Information Processing Systems, 2019.
[34] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In Advances in Neural Informa- tion Processing Systems, 2017.
[35] Vladimir Vovk, Alex Gammerman, and Glenn Shafer. Algorithmic Learning in a Random World. Springer, 2005.
[36] Vladimir Vovk, Alexander Gammerman, and Craig Saunders. Machine-learning applications of algorithmic randomness. In International Conference on Machine Learning, 1999.
[37] Sumanth Vrao. Daily climate time series data. https: http:///wwwwkll.e.../// com ctom/dataees// sumanthvrao / ddidd-alllyylodd climate- time- series- data, 2017.
[38] Chen Xu and Yao Xie. Conformal prediction interval for dynamic time-series. International Conference Machine Learning, 2021.
12

---

[39] Chen Xu and Yao Xie. Sequential predictive conformal inference for time series. In International Conference on Machine Learning. 2023.
[40] Margaux Zaffran, Olivier Féron, Yannig Goude, Julie Josse, and Aymeric Dieuleveut. Adaptive conformal predictions for time series. In International Conference on Machine Learning, 2022.
13

---

A Proofs
We begin with the proof of Proposition The proof very simple, and we derive it as a corollary of Proposition given next, because the proof reveals something perhaps unforeseen about the quantile tracker: acts as an error integrator, despite only adjusting the quantile based on the most recent time step.
Proof of Proposition Without loss of generality, set q = 0. Unraveling the iteration
yields
4++1
err; - a).
(13)
For (x)))) 172. h(t) b, (4) holds with 1/n.. Proposition
2 now applies.
Proof of Proposition 2 Abbreviate Ey = Ei=i a). We will prove one side of the absolute
inequality in (9) namely, Ex S h(T) I, and the other side follows similarly. We use induction.
The base case, for T = holds trivially. Now assume the result true up to T - 1. We divide the
argument into two cases: either h(T 1) E S h(T 1) 1 or Fr S h(T 1).
In the first case, note that that (4) implies qr = (Er-1) 2 b and therefore 8r S Qr and err, = 0. This means that
ET F a (TT 1) +
((T) +1,
h is nondecreasing, which the desired result T. In the second case, just use errr S 1,
ET S ET + a S c h(T 1) + 1 a ase e- h(T) +1.
This again gives the desired result at T, and completes the proof.
Proof of Theorem 7 We can transform (5) by setting 1411 = a++1 û.+1, and this becomes an update of the form (8 with respect to 4'+1- Further, the score sequence this new parametertizion is s/ = ÿ.. which is in [-b, b] because both and qe are in [-b/2, b/2]. Applying Proposition gives the result.
that,proof we only prove one side of the absolute inequality (II), and use induction. Abbreviate E7 > (Li y1) a)).
Proof of Proposition 3 The proof is similar to that of Proposition 2 as in that proof,
The base case holds trivially. For the inductive step, either c h(T 1) < ET S c h(T - +1
or Er-1 S c- h(T 1). In the first case, we have saturated, so (C,, yt) = 0, and
Fr- a e N(T 1) 1) 1- S e M(T) +
as is nondecreasing, which is the desired result at T. In the second case, we just use boundedness of the loss T,, (C,, yt) 1, so
MT ET < Er + 1)
+1.
This again gives the desired result at T, and completes the proof.
B Comparisons to existing methods
This section contains descriptions of the relationship of our method to pre-existing methods for conformal prediction under distribution shift.
B.1 ACI as a special case
Though it may not be immediately obvious, adaptive conformal inference (ACT) actually a special case of the quantile tracker. ACI follows the iteration:
0++1 a₁ 7 a),
14

---

which is equivalent to
- a1+1 = 1 at + (er,, a)
= ar 1/ (1 ar)).
for 3 =nfff 8t S Quantileg (fs1...., St-1)) This shows that ACI a particular instance of the quantile tracker that uses a secondary score s 3t and quantile g. quantile tracking
Thus, because
the same as a linear coverage integrator (T3) is ACI.
We can see here that ACI transforms unbounded score sequences into bounded ones, which then implies long-run coverage for any score sequence. This may, however, come at cost: ACI can sometimes output infinite or null prediction sets (when 0g drifts below O or above 1, respectively). Direct quantile tracking on the scale of the original score sequence does not have this behavior.
B.2 Weighted conformal prediction
A line of work by 33 and 4 deals with variants of weighted conformal prediction. Specifically, the latter variant allows the incorporation of sample weights for each data point. When the data come from probabilistic model, the sample weights allow for probabilistic bounds on the coverage of new test point.
These methods, and others (e.g.. (9) in the probabilistic setting, differ greatly from the approach herein. When there no probabilistic model, for example, all such approaches give vacuous coverage guarantees. The reason is that coverage on a new test point is impossible without probabilistic model. However, our setting we target long-run coverage, a completely different (and achievable) goal. Thus, neither the methods nor the proof techniques and theorem statements share any technical similarity to speak of.
For good measure, we include below in Appendix G evaluations against a version of 4 with trailing window the size of the burn-in period. This approach generally performs worse than both our approach and ACI, and does not have guarantees in our setting, so we do not include it in main text.
B.3 Sequential predictive conformal inference (SPCI)
recently proposed line of work exemplified by Ensemble Batch Prediction Intervals (EnbPI) [3 and its improved version, SPCI (39], seeks to provide coverage guarantees for black-box time-series predictors. However, the methods therein are not distribution-free, and SPCI rely on the following assumptions for validity:
1. The use of the quantile random forest (QRF) algorithm.
2. Model assumptions, such as consistency.
3. Distributional assumptions, such as smoothness.
Thus, the validity SPCI is significantly weaker than the methods herein, perhaps along the lines of a standard nonparametric analysis. In contrast, our method has finite-sample validity averaged over time for any input sequence. However, it is possible to use the QRF/SPCI algorithm as the scorecaster, which would endow it with distribution-fiee guarantee. We do not do this because the QRF is not a particularly effective time series forecasting method, and was outperformed by quantile LASSO (which use for the COVID experiments) and the Theta model. The performance of the QRF tends to be "blocky" and high-variance; for an example, see plots on the SPCI GitHub page
However, we note that the EnbPI and SPCI algorithms can be seen specific choices of scorecasters, and one consequence of our paper that these algorithms can be endowed distribution-fer guarantees when used as scorecasters in
C Heuristics for setting constants
Consider the tan integrator rz (ar) K₁ Lan(x (og(()((/CC))) where we set tan(x) sign(x) oo whenever & [-m/2, 2/2). and Cent: K₁ are constants. The constant Can is primarily in charge of guaranteeing that by time T, we want to have absolute guarantee of at least 1 - coverage.
15

---

Then we can set
Cant (1og(TT) 1/log(T))
ensure the function has an asymptote the correct point. The purpose of the constant K1 is to place the integrator on the same scale as the scores. So if B' is hypothesized bound on the magnitude of the scores, then one can set K₁ B'. In practice, these heuristics can be taken as starting place, and then the numbers can be fine-tuned during a burn-in period by hand or algorithmically. As alluded to previously, we believe there room for work in the style of 17 to rigorously tune these parameters online, but it is not the focus of our paper.
D Coverage integration as adaptive learning rate
As already mentioned the introduction, in our experiments we use a nonlinear saturation function (x(x)) K2 tan(x log(t)/ (e((t)))(sa))) where we set tan(x:) sign(x) oo for [-T/2,, xt//2 x/2]. and Csat: K1 > are constants that we choose heuristically (described in Appendix C In a sense, this tan integrator is akin to quantile tracker whose learning rate adapts to the current coverage gap. To this, we can use a first-order Taylor approximation, which shows (ignoring constants):
log(t)
(err:-a)
log(t) sec?
log(t 1) (-1 E(err, (err, - a) Vp1-a(5i-9)).
Σ
0++1 tan
~ git
1-1
effective learning rate
Above, sec(2) 1/ cos(2)) the secant function, which has U-shape; thus we can see from the above that the effective learning rate is higher for larger errors. Similar analyses for different integrators will give different adaptive learning rates; see below for another example with decaying leamingg rate.
Consider rz (x:) n L.. (This will give long-run coverage only for bounded scores, because " V Eizi(err, a).
condition (4) is only met for finite and not infinite b.) Then (8) becomes 4++1 which can be rewritten
Vl
Σ (er, - a) + "
(er: - a) -0)
" A (err, - ax).
@++1
V
V
V
qe +
This is approximately the quantile tracker with decaying learning rate, on the order of 1/ 1/vL.
E
Comparison to clipped ACI
Figures and & compare the quantile tracker to a clipped version of ACI which disallows infinite-sized sets by clipping the sets to the largest score far.
F More details on COVID-19 forecasting
In this experiment, the scorecaster receives as input the three most recent scores (i.e., quan- tile errors) of the ensemble forecaster, as well the three most recent case and death counts, from all 50 states. The scorecaster is an fi-penalized quantile regression as implemented by sklearn. ddinear_modhe at 10; in our experience, the performance of the scorecaster was fairly robust to this choice. Auto- matic selection (e.g., using coos-validation) could be the topic of future study. Figures and 10 shows the analogous experiments but for forecasting in New York and Texas.
Quanttleeeegressor. We fixed tuning parameter for the l penalty
G Further experiments
give a more comprehensive view of our results, examing all data sets, and a range of tuning parameters for each method. We restrict our attention to AR as the base forecaster; for the rest of the base forecasters, we refer the GitHub repository: hhtpss: : //tgttttttu/btt/ com/ http//aaaglooooo///////oooo
conformal seme--mmemeeeeee
16

---

ACI (clipped)
Conformal Control
Coverage
75%
1000
Sets
500
2043
2010
2014
2008
2017
2014
2018
Time
AR
Prophet
Theta
Transformer
clip ACT
P Ctrl
clip ACT
P Ctrl
clip ACI
P Ctrl
clip ACT
P Ctrl
Marginal coverage
0.898
0.898
0.884
0.897
0.898
0.896
0.884
0.897
Longest sequence
8
3
7
3
Average size
44,9
28.4
52.6
42.6
43.3
27.8
60.5
54.2
Median set size
41.8
18.6
38.8
31.9
27.3
18.7
36.6
35.5
75% quantile set size
58.9
37.4
66.9
56.5
59.5
37.8
85.5
69.7
90% quantile set size
93.9
66.3
137
93.4
94.7
63.4
148
123
95% quantile set size
136
81.8
166
116
136
78.5
182
164
Figure As Figure but with clipped ACI.
ACI (clipped)
Conformal Control
100%
Coverage
75%
50%
Sets
2002
2010
2012
2014
2016
2004
2010
2012
2017
2014
2013
AR
Prophet
Theta
Transformer clip ACI
clip ACI
Ctrl
clip ACI
Ctrl
clip ACI
Ctrl
P Ctrl
Marginal coverage
0.894
0.894
0.896
0.888
0.895
0.894
0.89
0.887
Longest em sequence
13
6
Average set size
19.5
17.6
69.5
51.7
17.9
17.8
115
70.4
Median set size
13.4
13.3
48.7
37.3
12.8
13.1
61.7
44,3
75% quantile set size
27.9
22.3
91.1
72.1
20.7
22.6
165
98.4
90% quantile set size
44
37.7
168
114
38.5
38.4
248
153
95% quantile set size
48.7
46.2
195
47.2
46.9
196
Figure 8: As in Figure but with clipped ACI.
17

---

Base Forecaster
Conformal PID Control
100%.
Coverage
Sets
1000
2020-10
2011-44
192107
201140
2922.1
2002 47
7000-00
291-100
2972-11
2022
2012-97
Time
Base Forecaster
Conformal PTD Ctrl
Marginal coverage
0.78
0.85
Longest sequence
8
3
Average size
310
472
Median set size
249
440
75% quantile set size
388
608
90% quantile set size
667
821
95% quantile size
754
1.03e+03
Figure 9: Results for 4-week ahead COVID-19 death forecasting in New York.
Base Forecaster
Conformal PID Control
100%
Coverage
15%
4000
Sets
2020-0
2021-1
2011-44
2021.7
202130
2022-01
2022
2002
2010-10
2027e1
202147
2021-10
2072-11
2012-07
Base Forecaster
Conformal PID Ctrl
Marginal coverage
0.75
0.85
Longest sequence
7
2
Average size Median set size 75% quantile set size
556
827
482
706
759
1.03e+03
90% quantile set size
993
1.36e+03
95% quantile size
1.1
1.57c+03
Figure 10: Results for 4-week ahead COVID-19 death forecasting in Texas.
18

---

Base Forecaster
Conformal PID Ctrl
Marginal coverage
0.82
0.86
Longest sequence
6
2
Average set size
625
Median set size
$12
688
75% quantile set size
754
1.01e+03
90% quantile size
1.12e+03
1.47e403
95% quantile size
1.45e+03
1.8c+03
Table 1: Summary statisties for COVID-19 death forecasting in California, as Figure
AR
Transformer
ACI
Conformal PID Control
ACI
Conformal PID Ctrl
Marginal coverage
0.899
0,9
0.899
0.901
Longest sequence
3
2
Average set size Median set size
0.177
0.174
0.406
0.178
0,426
0.175
75% quantile set size
0.484
0.21
0.574
0,206
90% quantile set size
0.672
0.236
0.233
95% quantile size
0.252
0.249
Results for the Prophet and Theta models are not available because darts does not support intermittent retraining for these algorithms.
Table 2: Summary statistics for electricity forecasting, as in Figure 3
For each experiment, we describe the data set in new subsection, and two plots are included: one for the coverage, and one for the prediction sets. Fach column in the plots represents different method, and each row is different learning rate. For the quantile tracker, the learning rate is to be interpreted as the multiplier in front of Br. Each method given a different color, which stays consistent throughout the plots. We use tan integrator and Theta scorecaster throughout, just in the main text experiments. There is an additional method, 'Trail', corresponding to the proposal of Barber et al. 4 with trailing window.
G.1 Amazon/Google
These data sets are part of a multivariate time series consisting of thirty blue-chip stock prices, including those of Amazon (AMZN) and Google (GOOGL), from January 1, 2006 to December 31, 2014. We attempt to forecast the daily opening price of each of Amazon and Google stock, on a log scale. Available to the scorecaster are the previous open prices of all 30 stocks.
G.2 Microsoft
This data set is a univariate time series consisting of single stock open price, that of Microsoft (MSFT), from April 2015 to May 31, 2021.
G.3 Daily temperature in Delhi
This data set contains the daily temperature (averaged over 8 measurements in 3 hour periods), humidity, wind speed, and atmospheric temperature the city of Delhi from January 2003 April 24, 2017, scraped using the Weather Underground API.
G.4 Electricity demand forecasting
This data set measures electricity demand in New South Wales collected at half-hour increments from May 7th, 1996 to December 5th, 1998 (we zoom in on the first 2000 time points). There are also
19

---

several other variables collected, such as the demand and price in Victoria, the amount of energy transfer between New South Wales and Victoria, and so on. These are given as covariates to the scorecaster. The demand value is normalized by default to lie in [0, 1/.
G.5 Synthetic data sets
We perform some experiments two synthetic score sequences which include change points and other behaviors difficult to produce using real data. In this setting, there is no ground truth sequence, so we do not plot the sets. Instead, we plot the scores themselves in one column, and the quantiles qe produced by each algorithm in a different column (when 9t 2 we cover). The general goal is for far off. that corresponds to the "set being too large
q: to track the or too small" in a situation where we would be constructing sets out of these scores.
quantile of 9t. and it
We consider i.i.d. sequence of scores, noisy increasing sequence of scores, and mix of change points and trends. Our codebase describes the score generation procedure in more detail.
20

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
1.0
0.5
0xge08 1%
I- cugedb. 0%
I- nigett 8%%
0.0
1.0
0.5
cxgr30.0%
I=0.5, oxg-a9.0%
5, ovg-89.75
In0055, cege59.4%
0.0
1.0
Coverage
0.0
1.0
0.5
ovge09.5%
ImD
Iwo 10000 oxgrod.0
0.0
1.0
0.5
1-0.0001, oxg-80.15 %
I-00, 0xg-0.0%
[=0, exg-87.96
[-40, cxg-7..%%
0.0
3000
2000
3030
1030
2000
Time
Conformel
Contermal
300
I-001
1000
17
500
300
1003
2803
2800
1000
380
1000
2008
300
2000
Time
Figure 11: Results for the Amazon data set.
21

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
1.0
0.5
-xge07 7%
t.ggr.xg- 0%
I-0. oxpedtt.ss 0%
1-G.5, oig
[ 5, oxget88 8%%
0.0
1.0
0.5
cxg-30.0%
I=0.1, 0xg-89.7%
I=0.1, oog.a9.00
In:0.,, cege-90.79
0.0
1.0
Coverage
cugelll
ougedi 6%
0.0
1.0
0.5
1 evge0R.0
Iw0.01 oxge88.5
ImD 01.01
InOO 01. oxgeOR.0K
0.0
1.0
0.5
1-0.0001, oxg-87.15
I=00, 0vg-0.0%
[-0, oxg-89.6%
[-40, oxg-59.8%
0.0
3000
2000
3030
1030
2000
1000
3000
Time
Conformel
Canterral
500
500
bm0.00
1000
17
500
500
1003
2803
2800
1000
2008
Time
Figure 12: Results for the Google data set.
22

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
1.0
0.5
0xge01 55%
t.ggr.xg-0 0%
I- nige80 0%
(ww) 1og-18888 58.8%
1.0
cxg-39.0%
In0.5, 0xg-89.0%
oxg-89.0%
Oege-50 8%
0.0
1.0
AA deea Ae
Coverage
ougett
0.0
1.0
0.5
Imo
IniO 7%
0.0
1.0
0.5
1-0.0001, oxg-88. 2%
I=0, 0xg-0.0%
[-0, oxg-89.8%
[+40, oxg-5988%
0.0
3000
2000
3030
2000
3000
Time
Cantermal
100
23
1-00
100
20
1003
2803
2800
1000
200
300
Time
Figure 13: Results for the Microsoft data set.
23

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
1.0
0.5
0xge60 35%
cxge000 1%
[- nige00 0%
1.0
cxg:30.1%
0.g-90.1%
oxg-30.1%
0ege:90:11
0.0
1.0
Coverage
oge00.00
oegen00:
0.0
1.0
0.5
035 exge90.3%
Im0.05 05,
1m0.05 8%.8%
Iwo oxge00.00
0.0
1.0
0.5
1-0.0001, 0xg-92.5%
I-0, 0xg-0.0%
[-0, ogg-70.8%
[-40, ovg-78..%
0.0
1003
1430
1000
1200
1000
1800
Time
grond ENI
MW
WY
WW
M
MVV
WW
M
20
mnne
ImA DS
Time
Figure 14: Results for the Delhi temperature data set.
24

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
1.0
0.5
0xge00.350 55%
1.gg000 1%
I- nige00 4%
(we (0g-30 1%
0.0
1.0
0.5
0xgr30.1%
In0.5, 90g-90.15
I=0.5, oxg-30.1%
In0055, ovge30.2%%
0.0
1.0
Coverage
ouge00 1%
oogen00 D%%
0.0
1.0
0.5
025, exgr90.19 1N
Iwo 0S,
I-o oxge08.5%
Ini 0%
0.0
1.0
0.5
1-0.0001, oxg-99.8%
I-0, 0xg-0.0%
[=0, oxg-88.65
[+40, oxg-59.8%
0.0
2000
1500
2000
1030
1000
2000
1500
2000
Time
Cantermal
0.99
3.2
1.0
0.7
0.5
0.55
3.75
0.50
2.2
180
I-000001
075
02
1006 1228
1000
520
1228
1E28
Time
Figure 15: Results for the electricity demand forecasting data set.
25

---

Trail
ACI
Conformal
Conformal P1 Begee
Conformal PID Prightore "orme
1.0
0.5
0xge00..00 0%
1.0g...g0000.0 5%
[-|, nigett 8%%
[- mag-00.0s 0%
Irw), moget8 2%
0.0
1.0
0.5
r-0.05 oxgr30.4%
1=0.5, 00g-90.15
I=0.5, 0xg291.2%
In0055, degein .0%
0.0
1.0
Coverage
ouge90.00
oegett
0.0
1.0
0.5
r-0 exgr09 9%
Im0.0s 05, oxge80.55
1m0.05 oxge80.8%
I-O oxge81 ..2%
0.0
1.0
0.5
1-0.0001, 0xg-89.3% 3%
[-0, 0xg-0.0%
[-0, ogg-86.00
[+40, ovg-9888.. 28.2%
0.0
505
1000
753
1000
250
750
500
1030
900
1030
Time
Trai
ACI
Conforme
Cantermal
0.0
0.4
6
0.2
0.4
0.4
KHD
120
315
726
320
1786
1700
Time
Figure 16: Results for an i.i.d. score sequence.
26

---

Trail
ACI
Conformal
Conformal P1
Conformal PID
0.5
orgen
cxgettt 5%
I- niget8 2%
[- magett
0.0
0.5
cxg-30.8%
In0.5, 0.g-88.2%
I=0.5, ovg
In0.5s 0xge30.44%
0.0
Coverage
0.5
ougett 9%
rigelt
oigeli
Oegen00
0.0
0.5
1-0 cvg-87.3%
I-0
I-
I-00055, oxg-90.15
0.0
0.5
I-- 0001, ovg-1 1.0%
[-0, 0xg-0.0%
[-0, oxg-42.25
[-40, ovg-000111
1%
0.0
500
1000
250
750
1000
900
1030
Time
Score
Trai
Conforme
Contoral
WMW
Cantermal
Inda
MMM
03
1.1
IHEUT
1.2
0.4
1.1
1O00001
033
120
315
726
321
1786
1706
Time
Figure 17: Results for an increasing score sequence.
27

---

Trail
ACI
Conformal
Conformal P1
Conformal PID Ve Svilge
1.0
/ Not
Cxge01 00%
1 2%
I- rigett IS
b- magedt
(we) 1oge91755
0.0
1.0
A
cxg-92.0%
1:0.5. 09.-92.4%
I=0.5, o.99921%
In005, oxge02.01%
0.0
1.0
Coverage
ougedd
0.0
1.0
0.5
1- 005 exge92.7%
05,
Iwi 0.ge82.2%
0.0
1.0
mm, o1gr20.85
0.5
1-0.0001, 01g-87.8%
[--0, cxg-88.5% 5%
[-40 ogg-50.5%
0.0
2003
3030
4000
4000
2000
4030
4000
Time
Conferme
Canterral
54000011
2012
3023
3000
306
Time
Figure 18: Results for a score sequence that is a mix of change points and trends.
28


Reviews:
[{"review_id": "m8rkUVyWHd", "reviewer": "Reviewer_Cw13", "summary": "This paper aims to provide asymptotically valid conformal prediction (CP) regions for the time series prediction problem. Similar to the adaptive conformal inference (ACI) framework [11], the main idea is to choose the appropriate quantile of the non-conformity scores (or equivalently tuning the miscoverage rate $\\alpha$) at each time step to fulfill the required error rate. The framework performs an online quantile update while taking into account the sum of coverage errors overtime to stabilize the coverage better. They provided theoretical analysis that guarantees the asymptotic coverage of the predictions made by their framework. Numerical experiments on real datasets show that their method brings more stable and efficient prediction sets than ACI. ", "questions": "Isn't it a good idea to provide a concrete explanation of why the proposed framework never returns infinite sets while this can happen in ACI? ", "limitations": "The authors fairly mention the limitations of the proposed framework in the discussion section. ", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "In general, I liked the idea of designing a PID-like framework. Thanks to this structure, the prediction regions turn out to be more stable and efficient. Moreover, it ensures that their framework never outputs infinite size prediction sets, which is a big advantage over ACI. Also, they provided many numerical experiments in the paper and the Appendix to illustrate the power of their proposed framework highlighting the effect of each block and different choices of hyperparameters. ", "weaknesses": "1. The related work section is not well-written in general. Also, it does not represent the works in the intersection of conformal prediction and time series prediction in a fair manner. Indeed many good papers are missing, such as \n\t- Foygel Barber, R., Candes, E. J., Ramdas, A., and Tibshirani, R. J. Conformal prediction beyond exchangeability. arXiv preprint arXiv:2202.13415, 2022,\n\t- Sun, S. and Yu, R. Copula conformal prediction for multi-step time series forecasting. ArXiv, abs/2212.03281, 2022,\nto mention a few. \n\n2. The comparisons are limited to ACI. Though it is a natural baseline to compare with, many other CP frameworks are designed for the time series prediction problem. It is interesting to see how they behave compared to the methods proposed in this paper. \n\n3. The scorecasting block seems unnecessary and, to some extent, redundant. Accounting for a wrong choice of learning model (e.g., a model that does not account for the seasonality effect) can not be a part of the design of a CP framework. This point is also mentioned in the paper in lines 286-288. However, it seems that the authors tried to put this block to force the similarity of their framework and the (well-known) PID controller.  \n\n4. There are numerous design choices (such as the saturation function) and hyperparameters in the framework, and it's worth noting that several of them may not be immediately intuitive for people in the Machine Learning or Conformal Prediction community when it comes to making appropriate choices.\n\n5. Generally speaking, I am still unsure to what extent it is worth exploring frameworks with an asymptotic coverage property, while the popularity of conformal prediction lies in its ability to provide a finite-sample guarantee.\n\n6. Minor Comments:\n\t- There are some typos in Equation 5. \n\n\t- Line115: There is a typo in the definition of the pinball loss (\\tau should be replaced with \\alpha). \n\n\t- Line 178: The equation reference should be (10) instead of (11).  \n\n\t- In the proof of Proposition 2., the capital C has to be replaced with the small c. \n\n\t- Bringing the extensions of the proposed method (e.g., conformal risk control and Proposition 3. or conditional coverage) into the discussion section is unappealing.\n\t- In order to facilitate the comparisons for the readers,  it is a good idea to report the Avg. size of the prediction sets and Avg. coverage in the figures like what you have in the Appendix figures. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the accurate summary and the thoughtful review. We enjoyed engaging with your comments.\n\nThe main response (in a separate comment) contains a summary of the extensive improvements we have made to the paper, including new experiments, finite-sample theoretical statements, automatic parameter selections, and a much-expanded section on related work. Herein we respond to the reviewer's specific concerns, and sometimes reference the main response.\n\n>The related work section is not well-written in general. Also, it does not represent the works in the intersection of conformal prediction and time series prediction in a fair manner. Indeed many good papers are missing, such as _Conformal prediction beyond exchangeability_ and _Copula conformal prediction for multi-step time series forecasting_, to mention a few. The comparisons are limited to ACI. Though it is a natural baseline to compare with, many other CP frameworks are designed for the time series prediction problem. It is interesting to see how they behave compared to the methods proposed in this paper.\n\nThank you. We have improved the related work section, which we now realize was perhaps written too hastily. We added the papers you suggested as well as:\n\n“Exact and Robust Conformal Inference Methods for Predictive Machine Learning with Dependent Data” — Chernozhukov et al.\n“Conformal prediction with temporal quantile adjustments” — Lin et al.\n“Inductive conformal martingales for change-point detection” — Volkhonskiy et al.\n\nWe are happy to include further references at the reviewer’s suggestion.\n\nFurther, if you would find it valuable, we can add a new subsection titled “Comparisons to Existing Methods”, that contextualizes our work among the existing conformal time series papers. We also point out a comparison to NexCP in our original submission. Please see the main response (5) for more discussion about this. \n\n> The scorecasting block seems unnecessary and, to some extent, redundant. Accounting for a wrong choice of learning model (e.g., a model that does not account for the seasonality effect) can not be a part of the design of a CP framework. This point is also mentioned in the paper in lines 286-288. However, it seems that the authors tried to put this block to force the similarity of their framework and the (well-known) PID controller.\n\nIt was not our intention to force a similarity to PID; rather, we believe that scorecasting can be useful in various practical situations where the base forecaster cannot easily be retrained, or cannot be easily updated, whether for computational or other practical reasons that have to do with the complexity of the original prediction pipeline. The COVID forecasting example is one such real-data situation. We have continued to work on this example and now adding a custom scorecaster makes a considerable improvement in the final coverage result — even though a massive, multi-institution effort was put into the initial model. This should clarify the real-world need.\n\nWe have added text to explain better when scorecasting might be useful (and when it is not needed and may even hurt). Our initial version did not explain this sufficiently well, and we thank you for the critical feedback.\n\nThere are numerous design choices (such as the saturation function) and hyperparameters in the framework, and it's worth noting that several of them may not be immediately intuitive for people in the Machine Learning or Conformal Prediction community when it comes to making appropriate choices.\n\nAgreed. In the new version of the paper, we provide default choices for all these parameters; see the main response (3).\n\n> Generally speaking, I am still unsure to what extent it is worth exploring frameworks with an asymptotic coverage property, while the popularity of conformal prediction lies in its ability to provide a finite-sample guarantee.\n\nThank you. In the new version of the paper, we state finite-sample guarantees for all these algorithms; see (2) in the main response.\n\n> In order to facilitate the comparisons for the readers, it is a good idea to report the Avg. size of the prediction sets and Avg. coverage in the figures like what you have in the Appendix figures.\n\nFor every plot, we have included a table that includes seven numerical comparisons: \n\n* Coverage.\n* Length of the longest sequence of errors.\n* Average set size.\n* 50%, 75%, 90%, and 95% quantiles of set size.\n\nWe included these results for four different forecasting algorithms (ARIMA, FB Prophet, Theta model (theta = 2), and a Transformer).\n\nThis addition is a small piece of a much larger effort we made to improve our experiments: see (1) of our main response. \n\n> Bringing the extensions of the proposed method (e.g., conformal risk control and Proposition 3. or conditional coverage) into the discussion section is unappealing.\n\nWe have shortened the discussion section which was previously too long and distracting. We can add a separate subsection on extensions to the camera-ready version if that is deemed best.\n\n> Minor Comments\n\nThank you for the fixes."}, {"review_id": "KlDz0L4fbW", "reviewer": "Reviewer_aXqf", "summary": "This paper examines the issue of how to parameterize conformal prediction in the time series setting. It argues that standard conformal inference methods would not provide valid inferences in the sequential setting, which lacks exchangeability. To resolve this, the paper suggests using online quantile tracking, integrating errors for stabilizing coverage, and utilizing scorecasting in the presence of systematic trends. The authors establish the corresponding theoretical properties, such as asymptotic coverage and conduct several experiments to justify their approach.", "questions": "1. How the specific selection of the function $r_t$ and $\\hat{q}_t$ influence the\noverall performance of the proposed method in conformal prediction for time series?\n2. In the context of time series conformal prediction, how does this new method compare\nwith other popular methods in terms of robustness, and computational efficiency?\n3. How does a non-constant threshold function $h$ help with tolerating a greater degree of\ncoverage error throughout the sequence?", "limitations": "1. Despite demonstrating strong results across various experiments, the authors acknowledge\npotential limitations, including the challenge of proving asymptotic coverage, potential\ninapplicability to different time series, and limitations in scorecasting.\n2. The methods they develop are specific to cases where response $y_t$ is revealed at each time\npoint. However, there are many settings in which we receive the response in a delayed\nfashion or in large batches.", "rating": 6, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1. A new method for conformal prediction in the in the time series setting with a clear presentation of the details and background.\n2. The authors show the proposed method has a major practical benefit, that unlike ACI (Gibbs and Candes, 2021), no infinite sets are produced.\n3. The authors demonstrate the practical utility of their methods with an extensive series of experiments on a wide range of real-time series data sets.", "weaknesses": "1. Although this paper present a general theorem that their method achieves asymptotic coverage, there is insufficient discussion on the coverage guarantees such as the lower bounds and upper bounds on coverage.\n2. The authors also explicitly avoid addressing the choice of automatic algorithms for the\ntwo parameter $C_{sat}$ and $K_I$ in the nonlinear saturation function.\n3. In the sequential setting (7), the authors didn’t show how to select the learning rate\n$\\eta$, what factors might affect the selection of the learning rate $\\eta$, and what impact this\nselection have on the performance of the method.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> This paper examines the issue of how to parameterize conformal prediction in the time series setting. It argues that standard conformal inference methods would not provide valid inferences in the sequential setting, which lacks exchangeability. To resolve this, the paper suggests using online quantile tracking, integrating errors for stabilizing coverage, and utilizing scorecasting in the presence of systematic trends. The authors establish the corresponding theoretical properties, such as asymptotic coverage and conduct several experiments to justify their approach.\n\nThank you for the summary, and we appreciate your engagement. \n\nPlease see the main response (in a separate comment) for a summary of the extensive improvements we have made to the paper. Below we reference the specific points related to your critique. We hope our responses answer the questions sufficiently to earn the reviewer’s support. Please let us know how we can improve if you still feel we are missing something.\n\n> Although this paper presents a general theorem that their method achieves asymptotic coverage, there is insufficient discussion on the coverage guarantees such as the lower bounds and upper bounds on coverage.\n\nIn the new version of the paper, we have included finite-sample bounds that exactly quantify the lower and upper bounds on coverage; see the main response (2).\n\n> The authors also explicitly avoid addressing the choice of automatic algorithms for the two parameter and in the nonlinear saturation function.In the sequential setting (7), the authors didn’t show how to select the learning rate, what factors might affect the selection of the learning rate, and what impact this selection have on the performance of the method.\n\nWe fully agree. In the new version, we have added default parameter selections for all parameters; see the main response (3).\n\n> How the specific selection of the function r and h influence the overall performance of the proposed method in conformal prediction for time series? How does a non-constant threshold function help with tolerating a greater degree of coverage error throughout the sequence?\n\nWe can add a new paragraph in our methods section addressing this point. The function $h(t)$ (which is the main determinant of $r(t)$) parametrizes a tradeoff between coverage and smoothness. For any function $h(t)$ the coverage gap will be no more than $(h(t)+1)/T$. Here are some example choices:\nTaking $h(t) = c$ enforces that, for all $T$, the long-run coverage $1-\\frac{1}{T}\\sum_{t=1}^T \\mathrm{err}_{t}$ is within $1-\\alpha \\pm (c+1)/T$.\nIf the function $h(t)$ grows (super-)linearly, asymptotic coverage is not guaranteed, but it can still be achieved approximately. For example, taking $h(t) = cT$ enforces that the same average coverage as above is within $1-\\alpha \\pm (c + 1/T)$.\nHowever, if the user sets, say, $h(t) = \\sqrt{T}$, then the long-run coverage is within $1-\\alpha \\pm \\frac{\\sqrt{T} + 1}{T} \\approx 1-\\alpha \\pm \\frac{1}{\\sqrt{T}}$. \nAll above guarantees hold deterministically, and simultaneously for all $T$, by Proposition 2.\n\nOne may see the coverage gap $(h(t)+1)/T$ and think about setting $h(t)$ to $0$ or to set it as small as possible. The problem is that this leads to wild oscillations in the sets, like having too high a learning rate.\n\nIn other words, the function $r(t)$ essentially acts as an adaptive learning rate, where the learning rate gets larger if the historical coverage gap is bigger. We show this formally using a Taylor expansion. Larger learning rates mean sharper discontinuities between adjacent sets — an extreme case is ACI with a high learning rate, which oscillates between infinite-sized and zero-sized sets. The benefit of setting, say, $h(t) = \\sqrt{T}$, is that the effective learning rate is lower, which means smaller oscillations in the sets while still maintaining the desired coverage guarantees.\n\n> In the context of time series conformal prediction, how does this new method compare with other popular methods in terms of robustness, and computational efficiency?\n\nWe can add a new subsection, titled “Comparison to Existing Methods”. We also evaluated against NexCP in the appendix. Please see the main response comment (5).\n\n> Despite demonstrating strong results across various experiments, the authors acknowledge potential limitations, including the challenge of proving asymptotic coverage, potential inapplicability to different time series, and limitations in scorecasting. The methods they develop are specific to cases where response is revealed at each time point. However, there are many settings in which we receive the response in a delayed fashion or in large batches.\n\nThe COVID-19 experiment already showcases a setting with delayed feedback — the observations have a massive 4-week delay. This case is already covered by our framework. With a little bit more work, we have been able to significantly improve the COVID-19 results since the initial submission. \n\nOur results would be easily extended to the batched-sequential setting as well. This is a straightforward extension of our existing analysis. We can add a note about this fact in a new extensions section; see the main response comment.\n"}, {"review_id": "s7V6yzZnvh", "reviewer": "Reviewer_bvK1", "summary": "The paper tackles the problem of conformal UQ under distribution shifts. This is a very realistic setting where the upstream prediction model cannot be frequently updated, and we still want to certify some level of safety. The paper cleverly fames the problem as a PID control problem, and provides a practical algorithm that they show is effective on synthetic and real datasets.", "questions": "- Can you elaborate on the reason why Enbpi and SPCI from Xu & Xie are not good baselines for comparing your model? Many have used previous timesteps as covariate X, similar to your scorecasting approach.   ", "limitations": "The authors provided a good address of limitations in their discussion section. ", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- The setting addressed is impactful and practical. A major hindrance of using CP in ML systems deployment is that the exchangeability assumption is often violated. This work provides a practical adaptation while still maintaining some guarantees of the method. \n\n- The PID framing of the problem & solution is clever and insightful. (very Big Brain, if i may use gen-z lingo.) The two Gibbs & Candes paper from 2021 and 2022 for CP under distribution shifts was able to circumvent the exchangeability assumption and achieve some guarantees by transforming the distribution shift into a online optimization problem for a single parameter $\\alpha^*$, which they prove is similar to a P controller.  This work takes the idea a step further, and introduces the integrater (I) and forecasting (D) algorithms, completing the picture (fig 2 <- such a nice figure). \n\n- The method can utilize time-series modeling approaches to adapt to seasonality. \n\n- The experiment section is well-presented. It walks you through part-by-part of the PID controller and shows the impact of each. The examples chosen are clear and convincing. \n\n ", "weaknesses": "1. The paper did not define some concepts it used in its writing. For example: \n- line 29, exchangeability \n- line 25, sharpest\n- line 83, burn-in period\n- line 137 (not really a definition just... what are you trying to say here?) \"the quantile tracker proceeds agnostically and performs the same updates in any case.\"\n\nI think explaining them, even just in the appendix, will help with clarity. \n\n2. The paper only compared to ACI & ablations of its own method. A common choice of algorithm for the scenario introduced by the authors would be Enbpi and SPCI from Xu & Xie (references 29 and 30 in the paper).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> The paper tackles the problem of conformal UQ under distribution shifts. This is a very realistic setting where the upstream prediction model cannot be frequently updated, and we still want to certify some level of safety. The paper cleverly frames the problem as a PID control problem, and provides a practical algorithm that they show is effective on synthetic and real datasets.\n\nWe thank the reviewer for their kind comments and helpful suggestions. We are interested in the reviewer’s thoughts on our expansion/revision of the paper as well. \n\nPlease see the main response (in a separate comment) for a summary of the extensive improvements we have made to the paper. Below we reference the specific points related to your critique. \n\n> The setting addressed is impactful and practical. A major hindrance of using CP in ML systems deployment is that the exchangeability assumption is often violated. This work provides a practical adaptation while still maintaining some guarantees of the method.\nThe PID framing of the problem & solution is clever and insightful. (very Big Brain, if i may use gen-z lingo.) The two Gibbs & Candes paper from 2021 and 2022 for CP under distribution shifts was able to circumvent the exchangeability assumption and achieve some guarantees by transforming the distribution shift into a online optimization problem for a single parameter, which they prove is similar to a P controller. This work takes the idea a step further, and introduces the integrator (I) and forecasting (D) algorithms, completing the picture (fig 2 <- such a nice figure).\n\nWe thank the reviewer! We also enjoyed making Figure 2 :)\n\n> The paper did not define some concepts it used in its writing. For example... (lists examples)\n\nWe appreciate the fixes. We can and will make edits to these places for clarity.\n\n> The paper only compared to ACI & ablations of its own method. A common choice of algorithm for the scenario introduced by the authors would be Enbpi and SPCI from Xu & Xie (references 29 and 30 in the paper). Can you elaborate on the reason why Enbpi and SPCI from Xu & Xie are not good baselines for comparing your model? Many have used previous timesteps as covariate X, similar to your scorecasting approach.\n\nAbsolutely. We agree that the previous methods by Xu and Xie share a relationship to scorecasting. The main difference is that these methods do not have a distribution-free guarantee, while ours does. Indeed, the method of Xu and Xie can be used as a scorecaster — but the scorecasting framework works for a larger class of methods, including more advanced ones that might not have the guarantees from that paper. This is explained in more detail in _main response point 5._\n\nThat said, we can include a new subsection, titled “Comparison with Existing Methods.” This goes through the related literature in greater detail, and explains the exact mathematical connections — including which methods we chose to compare to, and why.\n\nWe hope this makes the relationship between our new results and related works clearer. "}, {"review_id": "lmedcUrJi0", "reviewer": "Reviewer_H9oS", "summary": "The paper proposes an adaptation to existing adaptive conformal prediction in two novel ways: (i) by tracking the quantile via online regression over the running sum of the errors; and (ii) by incorporating a second model to anticipate the quantile in the next instant (\"scorecasting\"). The authors include theoretical guarantees for the resulting method, as well as comparisons against a benchmark in the field (ACI, which is shown to be a a special case of the proposed algorithm) over many different datasets. ", "questions": "- Why is asymptotic coverage the best one can hope for in the adversarial setting? Is this a formal result? It would be helpful to develop this idea a bit more in the text.\n- What are advantages of using the scorecaster in the conformal prediction setup vs adapting the model instead?\n- Any guidance on how should one pick $r_t$ and the learning rate? How was this decided for the figures in the main paper? The paper seems to imply trial and error. The impact on figures in the appendix is extremely large!\n- When do you expect ACI to be better than Quantile+Integrator+Scorecaster?\n- What are advantages and disadvantages of your proposal vs SPCI, given both achieve coverage of time series without creating infinitely-wide intervals? How does the empirical performance compares? In what ways is your proposal better?\n- Can you explain what you mean by \"Our algorithms may not be universally useful on time series different from those we used for evaluation\" (line 283)? When should it be useful?\n- There is a method called \"Trail\" in the appendix plots, but no explanation. What is it?", "limitations": "Yes.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "- Clear: the paper is well-written and very clear regarding its contributions.\n- Limitations: the discussions on the method's limitations are very helpful.\n- Theoretical results: the mathematical proofs are easy to follow and the notation is clear.\n- Empirical results: the empirical studies are extensive and cover several different datasets. The experiments in the supplementary material are also very thorough.\n- Code: the code is well-organized and readable (but please include the `environments.yml` file, all the necessary datasets, and instructions on how to reproduce figures in the paper! See below).", "weaknesses": "- Theoretical guarantees: the guarantees of the proposed method are only asymptotic. Unfortunately, no finite-sample results are included. None of the theoretical developments in the paper showcase advantages that the proposed method has over leading alternatives. \n- Baselines: there are many recent proposals in the field (even in your related works section), but the only baseline is ACI. SPCI (Xu and Xie, 2022) also obtained asymptotic longitudinal coverage of time series without creating infinitely-wide intervals; can you underline in what ways your proposal is better? \n- Experimental results: Can you include more details on how the hyperparameters were chosen for the examples in the main paper (for example, $r_t$)? Having too many hyper parameters can be overwhelming to the end-user (I appreciate the discussion in appendices C and D!). It would be great to include at least one example of adaptive learning in the proposed algorithm. \n- Scalability: all the time series considered, both in the main text and in the supplementary material (SM), are small, with the largest ones comprising around 3000 data points. In scenarios where new data arrives every second, for example, it is not clear if the method would work as-is or if any modifications would be needed. A discussion about this and an experiment with a very long time series would be great.\n- References: the authors point out the field of online calibration, and mention connections to game theory and online learning, but never explicitly discuss the connection to their work. It would be nice to better understand what tools from the field you are using or building upon.\n- Relationship to control theory: seems slightly out-of-place. Either explain the connection better (perhaps in the SM) or develop it in another paper? As it is, it's hard for the CP audience understand what the PID is, or why this connection is relevant.\n- Reproducibility: adding the proposed method's pseudocode would help a lot, either to the main text or to the SM if space is lacking. Trying to reproduce the results following the README.md, I was instructed to install the dependencies listed on `environment.yml`, but no such file was provided and I had to track them down one by one. There also seems to be a dataset missing (`datasets/preds_cases.csv`) and not all plots from the paper were generated.\n- Minor writing suggestion: on one hand, the relation to control theory seems currently very incipient and more of a curiosity, but it is presented prominently in a section of its own. On the other hand, the interesting Proposition 3 is cramped in the final section. You might consider adding an independent section for extensions and incorporate the control theory analogy into the final discussion section.\n- Other minor suggestions: I don't see the point of including the zoom in the figures; it oftentimes hides parts of the figure and the underlying patterns are clear from the whole graph. Besides the many plots, adding a table with a summary of the results among all datasets could help a reader to draw conclusions faster. The meaning of sublinear function is given in line 164, but the term is used before in lines 55 and 67, so we suggest bringing the definition up. Display (5) has a floating \\cdot in the last inequality and inconsistent usage (c \\cdot h(t) vs -ch(t)). Figure 1 caption says prediction sets are in purple, but it is actually in orange.\n- Typos: line 115: $(1-\\alpha)$ should be $\\tau$; line 139: the em-dash is confusing.; line 151: \"an particular\" should be \"a particular; proof of proposition 2: $C$ should be $c$; display under line 192: the second $\\text{err}_i$ should be $\\text{err}_t$; line 201: \"can made\" should be \"can be made\"; line 295: \"saturation function as in 2\" should be \"saturation function that satisfies (5)\" or \"saturation function as in Proposition 2\"line 301: missing \\mathrm for err_i.; the bibliography is not printed in the SM and the references do not match the main text, so it ca be confusing.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the accurate summary, and the engagement (esp. on reproducibility!)\n\nSee the main response (in a separate comment) for a summary of the extensive improvements we have made to the paper. Below we reference the specific points related to your critique. (We have shortened/paraphrased some questions for space.)\n\n> Code and reproducibility concerns: please include the environments.yml file, datasets, and instructions.\n\nSo sorry for the omission! We have made the code much more usable (added requirements.txt file, instructions, all datasets), towards public release.\n\n> Guarantees: guarantees are only asymptotic. No finite-sample results are included. \n\nIn the new version of the paper, we have added finite-sample bounds for all algorithms; see the main response comment (2). In addition, existing work such as Xu and Xie does not have a distribution-free guarantee. We elaborate in the main response comment (5), and will add a new section to the camera-ready.\n\nThe Xu and Xie result requires the use of a specific model (QRF), and distributional smoothness assumptions leading to, e.g., model consistency. It is also only asymptotic.\n\nOur results have finite-sample validity for adversarial sequences.\n\nPlease see the main responses (2) and (5) for further details.\n\n> Experimental results: adaptive learning rate?\n\nYes, definitely! In the new paper, we have provided a default choice for all parameters; see the main response (3).\n\n> Scalability\n\nGreat question. The reviewer makes a great point that sometimes samples are collected every second or even faster. \n\nFor extremely long time series that require extremely fast predictions (i.e. with a sub-millisecond latency), such as those in high-speed-trading, we may need some modifications of our methods. Quantile tracking would work as-is, with no changes. The windowed version of integral control might be more appropriate for the extremely long sequence. Finally, efficient scorecasting algorithms would be needed to ensure fast prediction.\n\nWe will add a paragraph in the discussion section mentioning this as an opportunity for future work in which non-obvious limitations of our approach might appear. The challenge of building a system of this scale and speed is likely to be an interesting and independent project.\n\nThank you for the wonderful point. It is a really nice suggestion.\n\n> References: It would be nice to better understand what tools from online calibration you are using or building upon.\n\nThe online learning community works on problems in the adversarial sequence model, but they do not provide guarantees of the sort studied in our work. We do not build upon the technical tools from that area — and we have tried to survey that literature for close connections, to little avail.\n\nWe will mention this in the new “Comparison to Existing Methods” section. The quantile tracker is equivalent to online gradient descent, an algorithm well-studied in online optimization. However, regret guarantees cannot be used to prove guarantees of the same form as in our work. The idea of what we call “saturation” does not seem to appear in that literature (perhaps because it isn’t natural for systems that are not prediction sets), and it is the core idea of our proofs.\n\n> Relationship to control theory: seems slightly out-of-place.\n\nTotally agreed. In the new rewrite, we have made this connection far more precise and formal, starting in the introduction.\n\n> Add pseudocode \n\nThe idea of pseudocode is great! We will add it to the supplement, since we’re already running out of space, but we will reference it very clearly in the main text.\n\n> Add exensions section?\n\nWe’ll create a dedicated extensions section; see the main response (4).\n\nOn the note of control, we agree that having its own section was overkill. At the same time, we have made the connection more precise to the point that we feel it can be more of a motivational point, so we have added it to the introduction and removed the dedicated section.\n\n> Why is asymptotic coverage the best one can hope for in the adversarial setting?\n\nGreat question. The original statement is actually not fully true (c.f. finite sample results). For a proof sketch of the impossibility of coverage at time $t$ (as opposed to the average version), consider the following.\n\nConformal picks any quantile $q_{t+1}$.\n\nThe adversary sees $q_{t+1}$, and picks $s_{t+1} = q_{t+1} + \\epsilon$. \n\nThen conformal miscovers, unless $q_{t+1} = \\infty$ with some probability $\\alpha$.\n\nTherefore, conformal fails at any time $t$ without randomization. So we instead need to target an averaged form of coverage. This is not far from a formal argument. We would be happy to write it down at your request.\n\n> Any guidance on how should one pick and the learning rate? \n\nSee the main response comment (3).\n\n> What are advantages of using the scorecaster in the conformal prediction setup vs adapting the model instead?\n\nIf the base model can be adapted easily, it might be better to do that. But sometimes it is hard to adapt the base model (for engineering or political reasons), or systematic trends in the scores persist even after extensive engineering. E.g., the covid model was trained by the CDC, but still can be improved by our method!\n\n> When do you expect ACI to be better than Quantile+Integrator+Scorecaster?\n\nIf very bad choices of parameters are chosen for the integrator, or the scorecaster is really badly trained and actually adds error, then ACI can definitely be better. (We have added a sentence to this effect.) \n\n> Comparison to SPCI?\n\nSee the main response comment (5).\n\n> There is a method called \"Trail\" in the appendix plots, but no explanation. What is it?\n\nThis is Barber et al. with a trailing window; see main respons  (5).\n\n> Besides the many plots, adding a table with a summary of the results among all datasets could help a reader to draw conclusions faster.\n\nSee (1) of the response.\n\n> Minor points\n\nThank you! Very helpful.\n"}]

Author rebuttal:
We are grateful to the four expert and engaged reviewers, who took a clear interest in the paper and suggesting ways to improve it. Thank you!

All four said the practical benefits are substantial, as evidenced by comprehensive experiments. Three lean towards acceptance.

We hope to respond to the critical comments through the following extensive revisions:

**1. Major experimental revamp (see PDF attached)** We tripled the number of experiments and significantly extended our evaluations.

The COVID-19 experiment was completely reworked to look at state-level forecasts, and includes a custom scorecaster that learns geographic correlations between cases and deaths. This provides a large improvement over the CDC model used in official communications on COVID forecasting. See PDF Fig. 1.

We also included experiments with all combinations of the following.

_Forecasters:_ Transformer, Facebook Prophet, Theta(2), and AR(3).

_Datasets:_ Electricity demand in New South Wales, Returns of AMZN, MSFT, and GOOGL, and Temperature fluctuations in Delhi.

Each experiment now includes a table with these metrics:
* Coverage over the sequence.
* Longest continuous sequence of miscoverages.
* Mean set size.
* 50%, 75%, 90%, and 95% quantiles of set size. 
(See PDF Figs 5 and 6; we couldn't include all tables due to space, and had to abbreviate/shrink the text. Apologies for this.)
All experiments now use score functions that allow for asymmetric prediction sets. The results are much improved. (See PDF Figs 2-4.)

**2. All results have been rewritten in finite samples.**
Reviewers noted that asymptotic coverage is less interesting than finite-sample coverage. 

We rewrote the paper so _all results are finite-sample_, and give both upper- and lower- bounds on the coverage gap (the distance of the empirical coverage from alpha).

The bound in Proposition 1 is

$\frac{1}{T}\sum\limits_{t=1}^T\mathrm{err}_t=\alpha \pm \frac{b+\eta}{\eta T}$

And Proposition 2 is

$\frac{1}{T} \sum\limits_{t=1}^T\mathrm{err}_t=\alpha \pm \frac{ch(t)+1}{\eta T}$

Note: The bounds suggest that large $\eta$ is desirable from the perspective of coverage. However, when $\eta$ is massive, the gradient step $q_t + \eta g$ can over/under-shoot $s_{t+1}$, leading to sets that do not properly track the score quantile and instead oscillate.  See point (3) in this comment for an automated choice of $\eta$ that works well in all our experiments.

These propositions are not new content; the proofs are the same as the old ones. Previously, we stated the asymptotic version for ease of exposition, but the finite sample results were already proven.

Towards conditional coverage, Prop. 3 can be extended to provide simultaneous coverage over any covariate-dependent kernel-weighted time-windows. We can add this if reviewers would find it valuable.


**3. Automatic selection of the learning rate.**
We provided an automatic choice of the learning rate and default parameters for the integrator. This learning rate has good performance on all our datasets, which are relatively diverse. We recommend it for practical use. This avoids all hand-tuning. The setting is:

$\eta = 0.1 \hat{B}_t$, 

Where $\hat{B}_t$ is an estimate of the max score computed within a trailing window. (Because this learning rate is technically changing online, it gains its validity via the integration argument, not the standard quantile tracking one.)

_All figures in the PDF use this learning rate._

The procedure is less sensitive to the integration parameters. For these, we provide defaults in the Appendix.

**4. New subsection on extensions**
Reviewers asked for a separate “Extensions” subsection in the camera-ready version, which we can add.

**5. More precise related work, and comparisons to SPCI.**

We propose a new subsection titled “Comparisons to Existing Methods”, if the reviewers would find it valuable. Its purpose would be to draw connections between our paper and related methods. Note that at the request of a reviewer, we already added a more comprehensive lit review of conformal prediction for time series. 

We included evaluations against NexCP in our original appendix, but forgot to mention it in text. The “Trail” column is NexCP with a trailing window of weights. ACI typically performs better than NexCP, so we didn’t feel it was worth putting into the main text.

Regarding SPCI, our method is _valid in greater generality_ than SPCI. Importantly, SPCI relies on the following assumptions for validity:

* The use of the quantile random forest (QRF) algorithm.
* Both model assumptions, such as consistency, and distributional assumptions, such as smoothness.
* The result is intrinsically asymptotic due to the consistency requirement.

Thus, we believe that a direct comparison of SPCI and the method herein is not appropriate. The validity of SPCI is significantly weaker, perhaps along the lines of a standard nonparametric analysis. In contrast, our method has finite-sample validity averaged over time for any input sequence. However, it is possible to use the QRF/SPCI algorithm as the scorecaster, which would endow it with a distribution-free guarantee. 

In our original submission, we did not pursue the QRF/SPCI angle because, in our initial experiments, the QRF was not a particularly effective time series forecasting method. It was outperformed by quantile LASSO (which we use for the COVID experiments) and the Theta model. The performance of the QRF tends to be “blocky” and high-variance; for an example, see plots on the authors’ GitHub: https://github.com/hamrel-cxu/SPCI-code/blob/main/tutorial_electric_EnbPI_SPCI.ipynb

If, after reading our response, the reviewers would still like to see an explicit comparison to SPCI, they can inform us, and we will include one (this essentially means adding QRF as a scorecaster, which is easy to do).

Finally, we once more thank the reviewers for their comments. The paper has improved a lot because of their feedback. 

Meta-review: This paper has been well-received by the reviewers, who all agree that the authors make a substantial contribution to uncertainty quantification for time series forecasting. A couple of issues have nevertheless been raised in the reviews, but these could be clarified in the rebuttal. Eventually, all reviewers are in favour of accepting the paper.
Decision: Accept (poster)