Review

Conformal forecasting RNNs

Stankevičiūtė, Alaa and van der Schaar train an LSTM to issue \(H\)-step point forecasts and add to each horizon one calibration quantile of the absolute residual, Bonferroni-corrected so that the \(H\) intervals cover jointly with probability at least \(1-\alpha\). The half-width is one number per horizon for every input. Two of the objections below are already in the literature, and we set them out with their sources before adding anything. What we have not found stated anywhere is the flaw in the comparison itself. The baselines are calibrated for per-horizon coverage and then graded on a joint criterion that only CF-RNN was widened to meet, and the paper’s own uncorrected results show the gap is that widening rather than the method. Neither coverage nor width is a proper score, and rerunning the paper’s COVID-19 experiment with one reverses its conclusion: the quantile RNN it reports as failing has the best interval on every split, and the dropout RNN that covers half the time still beats the uncorrected CF-RNN on two splits of three. Demonstration 35 shows the width mechanism on a synthetic series.

What is proposed

The paper is Stankevičiūtė, Alaa and van der Schaar, “Conformal Time-Series Forecasting,” NeurIPS 2021. The data are many exchangeable sequences, each a history \(y_{1:T}\) with an \(H\)-step continuation \(y_{T+1:T+H}\). An LSTM \(M\) trained on some of the sequences issues the \(H\) point forecasts \(\hat y_{T+1:T+H} = M(y_{1:T})\) at once. On \(m\) held-out calibration sequences the per-horizon absolute residuals are recorded, and Algorithm 1 sets

$$\hat\varepsilon_h \;=\; \text{the } \big\lceil (m+1)(1-\alpha/H) \big\rceil\text{-th smallest of } \big\{\, |y^{(i)}_{T+h} - \hat y^{(i)}_{T+h}| : i = 1, \dots, m \,\big\}, \qquad h = 1, \dots, H,$$

and the interval at horizon \(h\) for a new sequence is \(\hat y_{T+h} \pm \hat\varepsilon_h\). Theorem 1 states that all \(H\) intervals contain the truth with probability at least \(1-\alpha\), by split conformal validity at each horizon and Boole’s inequality across horizons. The half-width \(\hat\varepsilon_h\) does not depend on the new sequence.

The experiments report two numbers per method: the fraction of test sequences whose truth lies inside all \(H\) intervals, and the mean interval width. The baselines are a multi-quantile RNN (MQ-RNN, called QRNN in the code) trained by pinball loss at the 5% and 95% levels, and a Monte Carlo dropout RNN (DP-RNN) whose interval is the sample mean plus or minus \(1.645\) sample standard deviations. The paper writes that the baselines “both fail to achieve target coverage, sometimes reporting coverage rates as low as zero,” that “lack of coverage guarantees makes them less useful in high stakes real-world applications,” and that “CF-RNN adapts its prediction interval widths to reliably match the required target coverage.”

What has already been said

Lin, Trivedi and Sun raise the vacuous Bonferroni quantile twice. The appendix of Conformal Prediction with Temporal Quantile Adjustments (NeurIPS 2022) states that under the correction “the PI will always by [sic] infinitely wide if \(\alpha/T \le 1/(N+1)\),” and that “the original paper implemented the split-conformal incorrectly by ignoring \(v_{N+1}\), which is why this issue did not appear.” Their companion paper, Conformal Prediction Intervals with Temporal Dependence (TMLR 2022), names the dataset: “The authors performed an incorrect split-conformal experiment, which is why the COVID19 dataset still has finite width in Stankevičiūtė et al. (2021).” Everything in our section on the COVID-19 guarantee below is a confirmation of theirs, with two details added: the cap that rescues it is visible in the released code, and under the paper’s own split sizes the vacuity hits the COVID-19 experiment alone, since MIMIC-III at \(H = 2\) and EEG at \(H = 10\) have calibration sets large enough for the corrected index to exist.

The efficiency objection is also established. The same TMLR paper writes that CF-RNN “leads to unbalanced coverage (i.e. some TS receives poor coverage longitudinally while others high) and inefficient PIs.” Sun and Yu make it quantitative in Copula Conformal Prediction for Multi-step Time Series Forecasting (ICLR 2024): “We show that their method of applying Bonferroni correction produces inefficient confidence regions, especially for multidimensional data or long prediction horizons.” Replacing the correction with an estimated copula cuts the width by 33% on the same COVID-19 data and gives “a 40% sharper confidence region” on trajectory prediction, where they observe that “the confidence region produced CF-RNN is uninformatively large, as it covers all the lanes.” Cleaveland and co-authors (AAAI 2024) and at least two later papers repeat the conservativeness point in passing. Demonstration 35 illustrates the same conclusion by a different route, holding the point forecast and the calibration window fixed and varying only whether the width may depend on the state.

We have not found the next two objections anywhere in the citing literature.

The baselines are scored on a target they did not aim at

MQ-RNN is trained to produce, at every horizon, an interval whose per-horizon coverage is 90%. DP-RNN’s interval is a 90% normal band at each horizon. Neither is corrected for \(H\) horizons. Ten per-horizon 90% intervals cover jointly with probability \(0.9^{10} = 35\%\) if the misses are independent and somewhat more if they cluster, and the paper’s own uncorrected CF-RNN, Table 6, has per-horizon coverage between 85.5% and 91.6% on EEG and joint coverage of 59.4%. MQ-RNN’s 48% joint coverage on the same task, Table 5, is what a per-horizon 90% method looks like when judged on a joint criterion. The paper reads it as failure. The fair reading is that the baselines were never asked for joint coverage and CF-RNN was, and that widening any per-horizon interval by the same Bonferroni factor would produce the same comparison in reverse.

Confirming the COVID-19 guarantee is vacuous

The COVID-19 experiment uses \(m = 100\) calibration sequences and \(H = 50\) horizons at \(\alpha = 0.1\). The Bonferroni index is

$$\big\lceil (m+1)(1-\alpha/H) \big\rceil \;=\; \lceil 101 \times 0.998 \rceil \;=\; 101 \;>\; m ,$$

so the quantile in Algorithm 1 does not exist and, by the standard convention and by the proof of Theorem 1, \(\hat\varepsilon_h = +\infty\). The released code caps the quantile level at one and returns the largest calibration residual instead. We confirmed on a rerun that every Bonferroni half-width on this task equals the sample maximum. The intervals in Table 5, with a mean width of \(733.95 \pm 582.52\) and joint coverage of \(89.7 \pm 5.3\%\), are the largest error seen in calibration, added to and subtracted from the point forecast. That is a legitimate heuristic. It carries no finite-sample guarantee, and the paper presents it as the only method that “empirically achieve[s] the target joint coverage rates.” The same index calculation gives \(1733 \le 1823\) on MIMIC-III and \(3803 \le 3840\) on EEG, so the guarantee is real on those two.

Coverage and width are not a score

The paper states its selection rule: “as long as the coverage rate surpasses target coverage, the intervals should be as efficient (narrow) as possible.” Two numbers cannot rank methods on their own. A method can hit the coverage target with a width that is far too large on most inputs and too small on the rest, and a second method can hit the same target with a smaller average width and worse intervals. The proper score for an interval at level \(a\) is the interval score of Winkler,

$$\mathrm{IS}_a(l, u; y) \;=\; (u - l) + \frac{2}{a}(l - y)_+ + \frac{2}{a}(y - u)_+ ,$$

which is minimized in expectation by the true \(a/2\) and \(1-a/2\) quantiles and by nothing else. It charges width, and it charges each miss in proportion to its size at the rate \(2/a\). The paper’s rule is what the interval score becomes if the miss penalty is replaced by a threshold on the average miss rate, and the replacement is what lets a constant width pass.

The constant width is the point. The information a forecast interval should carry is how far the next value can wander from the point forecast given the state, and on real series that distance moves with the state. CF-RNN’s width is one quantile of the pooled residual distribution, so it is the pooled shape, and the loss from pooling is the quantity our information gap pages measure. The paper’s conclusion names the problem: future work “would involve making them more adaptive to individual observations.”

Rerun on their COVID-19 data

We reran the COVID-19 experiment with the authors’ code and hyperparameters, an LSTM with embedding size 20 trained for 1,000 epochs, on daily cases for the 307 English lower-tier local authorities from the UK Health Security Agency API over the same 150-day window, 100 observed days and 50 to forecast, split 162 training, 81 calibration and 64 test areas in the paper’s proportions, over 3 random splits. All three methods share the architecture, the hyperparameters and the training areas, so the only difference between them is how the interval is built. The Bonferroni index is \(\lceil 82 \times 0.998 \rceil = 82 > 81\), and on every split every Bonferroni half-width equalled the largest calibration residual. Every interval is scored at the level it claims. Means over splits, with standard deviations in parentheses.

The first row is the paper’s method and it is the calibration maximum at every horizon. It reaches 96% joint coverage against a 90% target, so by the paper’s criterion it wins, and it costs an interval score of 3,437 to get there.

The lower three rows all target 90% per horizon and can be compared with each other directly. MQ-RNN, which the paper reports as failing, has the best interval on 3 of the 3 splits: 238 on average against 394 for the uncorrected CF-RNN. Its coverage is 84% at a 90% target and the score still prefers it, because its misses are small. DP-RNN covers 51% of the time, which the paper reads as total failure, and it still scores better than the uncorrected CF-RNN on 2 of 3 splits. The uncorrected CF-RNN interval adds a constant to the same LSTM point forecast that MQ-RNN refines with two trained quantiles, and the constant is what costs it.

The same construction on the FRED daily universe

The skaters benchmark harness runs laplace, a zero-dependency online distributional forecaster, on the FRED daily universe of change series. For this study every series gets laplace at horizon \(H = 10\), and at each origin two sets of intervals are built from the same point forecasts, laplace’s predictive means. The CF intervals add the CF-RNN quantile of the last 500 resolved absolute residuals at each horizon. The laplace intervals are the quantiles of laplace’s own predictive distribution at the same levels. Both are scored at \(a = \alpha/H = 0.01\), the Bonferroni level with joint target 90%, and at \(a = \alpha = 0.1\) per horizon. 569 series, a median of 256 scored origins per series, scored on the most recent 2,000 changes after a 300-step warm-up. The universe is every cached daily series that clears the harness’s two standing screens, at least 500 changes and no single move dwarfing the series’ own scale, which is 569 of the 701 cached. Scores are on each series’ own scale, so the table reports medians over series and the share of series on which each interval scores better.

Per horizon at the Bonferroni level, the median score ratio is 1.151 at \(h = 1\) and 1.015 at \(h = 10\). The constant width over-covers at the joint level, 93.0% against a 90% target, because ten Bonferroni intervals that never share a state cannot use the correlation between horizons. The same run also scores laplace as a density, at a median held-out log score of 2.87 nats per observation. The constant-width interval has no density to score at all.

What survives

Theorem 1 is correct where its index exists. Split conformal quantiles per horizon with a Bonferroni correction give joint coverage over exchangeable sequences, and the paper is right that this holds for any point forecaster. The framing is also a real contribution: a set of independent sequences is exchangeable where the time steps within one series are not, and that is the right way to get a finite-sample guarantee out of a forecasting RNN. What does not follow is that the constant-width construction is a good interval, or that methods without the guarantee are worse. The guarantee is one number, the average coverage, and once it is met the comparison is a comparison of proper scores. On that comparison the constant width is beaten by the quantile RNN on every split of the paper’s own data.

The paper is Stankevičiūtė, Alaa and van der Schaar, Conformal Time-Series Forecasting, NeurIPS 2021, code at github.com/kamilest/conformal-rnn. The prior criticism is Lin, Trivedi and Sun, Conformal Prediction with Temporal Quantile Adjustments (NeurIPS 2022, arXiv:2205.09940) and Conformal Prediction Intervals with Temporal Dependence (TMLR 2022, arXiv:2205.12940), and Sun and Yu, Copula Conformal Prediction for Multi-step Time Series Forecasting (ICLR 2024, arXiv:2212.03281). The interval score is Winkler (1972) and Gneiting and Raftery (2007), equation (43). The FRED study is benchmarks/cfrnn_study.py in the skaters repository, and the COVID-19 rerun script is data/reviews/run_covid.py here. The UK Health Security Agency dashboard API supplies daily cases by specimen date for English lower-tier local authorities, where the paper’s code reads a 24 May 2021 snapshot of cases by publication date for United Kingdom local authorities, so the areas and the counts differ from the paper’s and the split sizes are scaled to the 307 areas available.