Conformal prediction interval for dynamic time-series

Chen Xu Yao Xie
[c]School of Industrial
and Systems Engineering
Georgia Institute of Technology
Email<EMAIL_ADDRESS>[c]School of Industrial
and Systems Engineering
Georgia Institute of Technology
Email<EMAIL_ADDRESS>]

We develop a method to build distribution-free prediction intervals for time-series based on conformal inference, called |EnPI| that wraps around any ensemble estimator to construct sequential prediction intervals. |EnPI| is closely related to the conformal prediction (CP) framework but does not require data exchangeability. Theoretically, these intervals attain finite-sample, approximately valid average coverage for broad classes of regression functions and time-series with strongly mixing stochastic errors. Computationally, |EnPI| requires no training of multiple ensemble estimators; it efficiently operates around an already trained ensemble estimator. In general, |EnPI| is easy to implement, scalable to producing arbitrarily many prediction intervals sequentially, and well-suited to a wide range of regression functions. We perform extensive simulations and real-data analyses to demonstrate its effectiveness.

§ INTRODUCTION

Constructing prediction intervals (PIs) is a fundamental problem in statistics and machine learning –  given a sequence of scalar observations, along with exogenous variables as feature vectors, how to predict the future while quantifying the uncertainty? It is common to use regression functions for prediction; regression functions can be arbitrarily complex, such as random forest <cit.> and various deep neural network structures <cit.>. Multiple regression functions are often combined into an ensemble learner, to increase accuracy and decrease variance for point estimates <cit.>. It is informative to build PIs around point estimations to quantify where actual values lie. However, this task can be challenging without strong assumptions on the underlying data distribution. It is even more challenging for time-series fitted with arbitrary regression estimators, especially when test data come in sequence.

This paper directly addresses the challenges above by building distribution-free PIs for dynamic time-series data with a coverage guarantee. In particular, we efficiently build PIs for point estimates from ensemble estimators without refitting multiple ensemble models. We also construct PIs sequentially, so that our method can produce arbitrarily many PIs during test time.

Recently, a broad family of methods called conformal prediction (CP) is becoming popular. Formally introduced in <cit.>, this method assigns “conformity scores”[In this paper, conformity scores are calculated as residuals from fitting a regression algorithm $\mathcal{A}$ on the training data.] to training data and test data. Via permutation tests, one obtains $p$-values for hypothesis testing; inverting the test generates prediction intervals for test data. Under exchangeability in data, this procedure generates valid average coverage of the test point. Many works such as <cit.> operate under this logic. For comprehensive surveys and tutorials, we refer readers to <cit.>. In particular, our method is very close to the method by <cit.>, which constructs an efficient wrapper around ensemble learners to generate prediction intervals. However, all these methods above assume exchangeability in data, and for time-series, this assumption is hardly reasonable.

Recently, adapting CP methods under distribution shift has also been an important area. The work by <cit.> uses weighted conformal prediction intervals when the shifted distribution on test data is proportional to the pre-shift training distribution. Another recent work by <cit.> provides a coverage guarantee when the shifted distribution lies in an $f$-divergence ball around the training distribution. However, both works still assume i.i.d. or exchangeable training data, making them not directly suitable for time-series.

Nevertheless, the works by <cit.> shine light on how to design the permutation schemes carefully and impose certain consistency conditions on the estimator, so that the resulting coverage is approximately valid for time-series data with strongly-mixing errors. As a result, these works provide theoretical bases for the method in this paper to work. Note that these works also guarantee that the inductive conformal prediction method (ICP) by <cit.> has valid coverage for time-series data. However, ICP requires data splitting, making the intervals too wide; this paper's method does not require any data splitting.

Below is a tabular summary of the coverage guarantee of CP methods under different assumptions[Remark: different works proceed under these assumptions differently (due to varying purposes). Hence, the table presents the best attainable coverage guarantees. Yet, doing so may not be ideal in practice since intervals may be too wide under these guarantees.

The guarantee is proved in <cit.> under exchangeability, in <cit.> under covariate shift, and in <cit.> under strongly mixing time series.]:

Assumption     [l]Exchangeable     [l]Covariate
Shift     [l]Strongly

In Theory     $1-\alpha$     $1-\alpha$     $\approx 1-\alpha$

Empirically     $1-\alpha$    $1-\alpha$     $1-\alpha$

Furthermore, we mention several other methods outside the CP literature, which also allow predictive inference. Traditional time-series methods, such as ARIMA(p,d,q), exponential smoothing, state-space models (e.g., Kalman Filter), have been widely successful <cit.>. On the other hand, <cit.> use a discriminative learning framework to optimize the expected error rate under a budget constraint on interval sizes. Although it provides PAC-style, finite-sample guarantees, the training data are still assumed i.i.d.

Lastly, this work is not closely related to online learning. Two common online learning objectives are bounding expected mistakes until the first success <cit.> or maximizing cumulative reward over the test period under a reward function <cit.>. Although we build prediction intervals sequentially and may consider a coverage failure as a mistake, our primary focus is only upper bounding the percentage of a mistake by $\alpha$ over an arbitrarily long prediction period. We do not care when mistakes happen or penalize such mistakes.

We summarize the two main contributions as follows:

* We present a robust and computationally efficient predictive inference method around any ensemble estimator for constructing sequential prediction intervals, without any data splitting. It also works well for small-sample problems.

* We show that these prediction intervals enjoy approximately valid coverage under mild and verifiable assumptions on the underlying data and regression predictor. In particular, the method is suitable for non-stationary time-series.

§ PROBLEM SETUP

The problem is formally stated as follows. We assume $Y_t=f(X_t)+\epsilon_t$ for an unknown function $f:\mathbf{R}^d\rightarrow \mathbf{R}$. Note, $X_t$ here is expected to either be exogenous time-series sequences or be the history of $Y_t$. We implicitly assume that $f$ is independent of time, so $Y_t$ are all generated by the same function. Note, we do not assume that errors $\epsilon_t$ are independent over time; the precise assumption will be stated in Section <ref>. Notation-wise, upper case $X_i,Y_i$ are random variables and lower case $x_i,y_i$ are data.

Assume we have $\{y_i,x_i\}_{i=1}^T$ as realization of this random process up to time $T$. Our goal is to construct prediction intervals for $Y_t$ using $x_t$ when $t>T$. We do so sequentially: the algorithm produces a prediction interval $C_{T,T+1}$ for $Y_{T+1}$ and uses $y_{T+1}$, as soon as it becomes available, to subsequently produce prediction interval for $Y_{T+2}$.

Given a significance level $\alpha$, $C_{T,t}$ often depends on $\alpha$; we henceforth denote it as $C^{\alpha}_{T,t}$ (first $T$ denotes we use $T$ data points from the past to produce the prediction interval at $t$). In particular, the prediction interval $C^{\alpha}_{T,t}, t>T$  needs to satisfy the following average coverage guarantee:

\begin{equation}\label{eq:ave_cov}
P(Y_t \in C^{\alpha}_{T,t}) \geq 1-\alpha.
\end{equation}

In this paper, we require
\[C^{\alpha}_{T,t}:=\hat{f}_t(x_t) \pm (1-\alpha)\text{ quantile of }
% \footnote{For a distribution $F$, if $Z\sim F$,
% \vspace{-0.32cm}
% \begin{equation*}
%     \text{Quantile}(1-\alpha,F)=\text{inf}\{z:\mathbb{P}(Z\leq z)\geq 1-\alpha\}.
% \end{equation*}
% For a discrete distribution of values $a_1,\ldots,a_n$, where $\delta_{a_i}$ denotes a point mass distribution at $a_i$,
% \vspace{-0.32cm}
% \begin{equation*}
% \text{Quantile}(1-\alpha,a_{1:n})=\text{Quantile}(1-\alpha,\sum_{i=1}^n \delta_{a_i}/n)
% \end{equation*}
% \vspace{-0.6cm}}
\{\hat{\epsilon_i}\}_{i=t-T}^{t-1},
\]
with $\hat{\epsilon}_i:=y_i-\hat{f}_i(x_i)$. We remark that although the underlying data generating processing uses a fixed unknown $f$, our estimator $\hat{f}$ may be different for each data point. This happens because we use ensemble methods to estimate $f$, potentially creating different ensemble estimators for each data point. Nevertheless, we will see that these $\hat{f}_i$ are close to the true $f$ under appropriate model assumptions, hence guaranteeing the performance of our method.

We can equivalently conduct a hypothesis testing at time $t$ with the prediction interval, where the null hypothesis is that all values in $C^{\alpha}_{T,t}$ are plausible values for $Y_t$ under $\alpha$ significance level.

As noted earlier, the inequality <ref> holds exactly under exchangeable $(X_i,Y_i)$; such an assumption is clearly violated for time-series in general. However, by adopting results from earlier learning theory works <cit.>, we can ensure (<ref>) holds approximately, meaning that the probability of under-coverage can be bounded at any finite sample size $T$ and approaches zero as sample size reaches infinity.

From now on, a prediction interval is called valid if its average coverage is $1-\alpha$ under $\alpha$ significance level and is called efficient if it is shorter in width than other valid prediction intervals.

§ ENPI ALGORITHM

Given the setup above, we now present Algorithm <ref>, named |EnPI|; it uses ensemble estimators to approximate $f$ and calculates residuals $\hat{\epsilon}_i$ in an leave-$i$-out fashion.

We briefly discuss the inputs to this algorithm: we let $\mathcal{A}$ denote a particular regression algorithm that maps data into a regression estimator.[In general, $\mathcal{A}$ can consist of different types of regression algorithms (e.g., non-parametric regression, neural network, linear models, etc.), if we believe the true model belongs not to any individual family but their union.] The decision threshold $\alpha$ is user-defined, under which smaller thresholds yield wider intervals. We do not place any restriction on the aggregation function $\phi$, which can typically be chosen as mean, median, or truncated mean. The fixed positive integer $\widetilde{B}$ determines how many models trained under $\mathcal{A}$ are combined into an ensemble estimator. Theoretically, any $\widetilde{B}$ so that $\widetilde{B}e^{-1}$ is greater than 1 guarantees validity, and we expect a larger $\widetilde{B}$ to yield shorter and more stable intervals. Empirically, letting $\widetilde{B}e^{-1}\approx$ 10 to 20 is sufficient, especially for computationally intensive methods such as neural networks. Lastly, note that each response variable $y_t, t>T$ is revealed after we compute the prediction intervals at time $t$, so we never use residual $\epsilon_t$ for computing the prediction interval at time $t$.
Sequential Distribution-free Ensemble prediction interval  (|EnPI|)

[1]
Training data $\{(x_i, y_i)\}_{i=1}^T$, regression function $\mathcal{A}$, decision threshold $\alpha$, aggregation algorithm $\phi$, a fixed positive integer $\widetilde{B}$, and test data $\{(x_t,y_t)\}_{t=T+1}^{T+T_1} $, with each $y_t$ observable in sequence only after the prediction interval is constructed for time $t$.
Prediction intervals $\{C^{\phi,\alpha}_{T,t}(x_t)\}_{t=T+1}^{T+T_1}$
Draw $B\sim \mbox{Binomial}(\widetilde{B},e^{-1})$.
$b = 1, \dots, B$
Sample with replacement an index set $S_b=(i_1,\ldots,i_T)$ from training indices $(1,\ldots,T)$
Compute $\hat{f}^b=\mathcal{A}(\{(x_i,y_i) \mid i \in S_b \})$.
Initialize $\boldsymbol\epsilon=\{\}$
$\hat{f}_{-i}^{\phi}(x_i)=\phi(\{\hat{f}^b(x_i) \mid i \notin S_b \})$
Compute $\epsilon_i^{\phi}=|y_i-\hat{f}_{-i}^{\phi}(x_i)|$
$\boldsymbol\epsilon=\boldsymbol\epsilon \cup \{\epsilon_i^{\phi}\}$
Compute $\hat{y}_t^{\phi}=(1-\alpha) \text{ quantile of } \{\hat{f}^{\phi}_{-i}(x_t)\}_{i=1}^T$
Let $w^{\phi}= (1-\alpha) \text{ quantile of } \boldsymbol\epsilon$.
Return $C^{\phi,\alpha}_{T,t}(x_t)=\left[\hat{y}_t^{\phi} - w^{\phi},\hat{y}_t^{\phi} + w^{\phi}\right]$
Compute $\epsilon_t^{\phi}=|y_t-\hat{y}_t^{\phi}|$
$\boldsymbol\epsilon=(\boldsymbol\epsilon- \{\epsilon_1^{\phi}\}) \cup \{\epsilon_t^{\phi}\}$ and reset index of $\boldsymbol\epsilon$.

We make two important comments for |EnPI| in the following. First, we remark on the details of the construction. The choice of a random $B$ in Line 1 is required for exchangeable data, as discussed by <cit.>. Due to the concentration of binomial random variables, $B$ is typically close to $\widetilde{B}e^{-1}$. Despite being very similar to the jackknife+-after-bootstrap (J+aB) method by <cit.>, |EnPI| exploits time-series structure in Line 12-18, by producing prediction intervals sequentially using a sliding window of size $T$. This is often reasonable and necessary for time-series data since data are closer to the test point when $t$ is more predictive of future values than further ones. Meanwhile, appropriate aggregation $\phi$ can serve different purposes, such as reducing mean squared error (MSE) under mean, avoiding sensitivity to outliers under median, or achieving both under trimmed mean. Lastly, when each $\hat{f}^{\phi}_{-i}$ is close to $f$ in terms of MSE, refitting the ensemble estimators is not necessary, so residuals at time $t>T$ are all calculated using ensemble estimators trained on the first $T$ points. This claim will be more rigorously analyzed in Section <ref>.

Next, we remark on the computational cost. The most costly operations happen in Line 2-5 when training $B$ many $\hat{f}^b$, especially for a deep neural network regressor. Nevertheless, we do not need to retrain any estimator for $\hat{f}^{\phi}_{-i}$ in later lines. In comparison, if we naively build leave-$i$-out ensemble predictors by training $\mathcal{A}$ for $B$ many times on $\{x_j,y_j\}_{j\neq i}$, the computational cost in terms of fitting $\mathcal{A}$ drastically increases to $Bn$, which is infeasible in reality. The performance of these two implementations is similar under exchangeable data, as shown by <cit.>. At this stage, one may naturally ask why shouldn't we build $\mathcal{A}$ once on the whole training data and using the in-sample training residuals to compute prediction intervals (this method is titled Jackknife). However, due to overfitting, Jackknife tends to undercover even for exchangeable data; <cit.> presents a degenerate example in which Jackknife attains zero coverage. In the next section, we will show that using in-sample errors is prohibited theoretically for time-series data. Lastly, to avoid this overfitting issue and ensure computational efficiency, <cit.> proposes the ICP method to train on a subsample of the training data and use the rest to calibrate residuals. Although ICP prediction intervals are also approximately valid, we will demonstrate in the experiments that its intervals widths are often broader and less stable than those from |EnPI|. This situation happens because ICP faces data-splitting trade-off, yet |EnPI| cleverly uses all $T$ observations as both training data for ensemble estimators and calibration data for leave-$i$-out residuals.

§ THEORETICAL ANALYSIS

Without loss of generality and for notation simplicity, we only consider our algorithm's validity on the first test point with index $T+1$. We comment on why validity holds for all points from $T+2$ onward at the end of this section.

We first provide an elementary equivalence between the prediction interval $C^{\alpha}_{T,T+1}$ and the empirical $p$-value at $T+1$, where $\hat{p}_{T+1}:=T^{-1}\sum_{i=1}^T\mathbbm{1}\{\hat{\epsilon}_i>\hat{\epsilon}_{T+1}\}$:

\[
Y_{T+1} \in C^{\alpha}_{T,T+1} \mbox{ iff } \hat{p}_{T+1} \geq \alpha.
\]

Indeed, this holds due to basic algebraic computation. Intuitively, $Y_{T+1}$ lies in the prediction interval constructed using the $1-\alpha$ quantile of past $T$ training errors if and only if this empirical $p$-value is not significant under $\alpha$.

Therefore, our method covers $Y_{T+1}$ with probability  $1-\alpha$ approximately if the distribution of $\hat{p}_{T+1}$ is approximately uniform. More precisely, we aim to ensure that $|P(\hat{p}_{T+1}\leq \alpha)-\alpha|$ is small. Notation-wise, define the empirical distribution of $\hat{p}$ as $\hat{F}(x)=P(\hat{p}_{T+1}<x)$.

Earlier literature in learning theory provides such a bound. In particular, we remark that our Algorithm <ref>, which uses symmetric baseline predictors, can be fitted into Algorithm 1 by <cit.> with $p=1, T_1=1$, and non-overlapping block permutations. As a result, we borrow the same assumptions and theorem from <cit.>, whose Theorem 1 is a more concise version than Theorem 2 by <cit.>. Both theorems therein validate the approximately uniform distribution of $\hat{p}_{T+1}$.

The precise assumptions and theorem are as follows:

$\{\epsilon_i\}_{i=1}^T$ are stationary, strongly mixing, with sum of mixing coefficients bounded by $M$.

Let $\{\delta_{T},\gamma_{T}\}$ be sequences of numbers converging to zero. Assume that with probability $1-\gamma_{T}$, estimation errors are small:

* Mean squared error: $\sum_{i=1}^{T+1} (\hat{\epsilon}_i-\epsilon_i)^2/(T+1) \leq \delta_{T}^2$;

* Pointwise prediction error: $|\hat{\epsilon}_{T+1}-\epsilon_{T+1}| \leq \delta_{T}$;

* The pdf of $\epsilon_{T+1}$ is bounded above by a constant $D$.

The following theorem, with minor revisions, follows from these two assumptions:

For any $\alpha \in (0,1)$, the empirical $p$-value $\hat{p}_{T+1}$ obeys:

\begin{align*}
&|P(\hat{p}_{T+1}\leq \alpha)-\alpha |\\
&~~\leq C_1T^{-1/4}\log(T)+C_2(\delta_T+\sqrt{\delta_T})+\gamma_T,
\end{align*}

where constants $C_1=\max(6,\sqrt{M}),C_2=\max(2D+4,4D)$.

For proof of Theorem <ref>, we direct the readers to proofs of Lemma 11, 12, and 14 in <cit.>, after which Theorem <ref> follows.

We remark that although $T^{-1/4}\log(T)$ converges to zero really slowly, empirical results show that even at small $T$, $|P(\hat{p}_{T+1}\leq \alpha)-\alpha| \approx 0$, so the bound may be improved. On the other hand, $\delta_T=o(1)$ by assumption, so $\delta_T=O(\sqrt{\delta_T})$.

In the remainder of this section, we validate why the ensemble regressors in |EnPI| satisfy the assumptions above reasonably well. We also give examples and discuss the implications of these assumptions.

Assumption <ref> is a very mild assumption on the original process $\{(Y_t,X_t)\}$. The series can exhibit arbitrary dependence and be highly non-stationary, but still have strongly-mixing (or even i.i.d.) errors. In particular, any stationary Markov chains that are Harris recurrent and aperiodic are strongly mixing <cit.>. ARMA with i.i.d. innovations and GARCH processes are also strongly mixing. A more comprehensive list of mixing processes is given in <cit.>, including Gaussian random fields, Gibbs fields, continuous-time processes, etc.

Assumption <ref> requires scrutiny. The famous No Free Lunch Theorem <cit.> implies that all regression algorithms have equal average performance when the target function $f$ is completely random. In other words, we cannot expect any $\mathcal{A}$ to satisfy this assumption when $f$ is completely unrestricted. Nevertheless, this assumption has been verified for certain classes of $f$, as two examples below show:
[Neural Networks]
Assume $f$ is suitably smooth, in the sense that its Fourier transform has a bounded first moment of the magnitude distribution. <cit.> provide a root-mean-square convergence rate of $o_P(T^{-1/4})$ for general neural networks sieve estimators, by letting the number of hidden units $r_T$ increases with $T$ on the order of $O(T)$. Together with results in <cit.>, this result implies the desired root-$n$ asymptotic normality for plug-in estimators of smooth functionals using general neural networks sieve estimators.

[Sparse High-dimensional Linear Model]
Given $\mathcal{F}_M$, which is a dictionary of basis functions,
<cit.> assume that $f$ is well approximated by a linear combination of a small number of members of $\mathcal{F}_M$. In this case, they first establish approximate equivalence between the Lasso and the Dantzig selector in estimating $f$ and then show the distance between their prediction losses is of the same order as the prediction generalization error. In particular, with probability at least $1-M^{1-A^2/8}, A > 2\sqrt{2}$, the bound (Theorem 5.1) is a product of three factors: the ratio of the sparsity of $f$ with sample size $T$, how ill-posed the Gram matrix for sparse estimation is, and the log size of the dictionary $\mathcal{F}_M$.

Beyond specific classes of $f$, we note that this assumption closely relates to the risk (or generalization error) of $\mathcal{A}$. In particular, when $\mathcal{A}$ is uniformly stable, its generalization error can also be bounded using the empirical error and leave-one-out error, as in Theorem 2 of <cit.>. Examples of algorithms that are stable can be found in <cit.>, <cit.>,  <cit.>.

On the other hand, we remark that our leave-$i$-out predictors $\hat{f}^{\phi}_{-i}$ to $f$ are ensemble predictors, so we need to analyze how ensemble predictors can satisfy Assumption <ref>. Our analysis below works directly for $\phi$ as a weighted mean of individual regressors (such as mean or truncated mean); we provide some ideas for using general $\phi$ in the end.

We first note that the strongly-mixing assumption on error processes $\{\epsilon_t\}$ and $\{\hat{\epsilon}_t\}$ implies their Ergodicity. Ergodicity thus implies that

\begin{align*}
\sum_{i=1}^{T+1} (\hat{\epsilon}_i-\epsilon_i)^2/(T+1) & \rightarrow \mathbbm{E}[|\hat{\epsilon}-\epsilon|]\\
& =\mathbbm{E}[|\hat{f}^{\phi}(x)-f(x)|],
\end{align*}

which is the expected mean square error for each ensemble predictor. We drop the index because in expectation over the random binomial draw of $B$, each $\hat{f}^{\phi}$ is combined using the same number of individual predictors, all of which are computed using the same baseline regression algorithm. As a result, $\hat{f}^{\phi}_{-i}$ and $\hat{f}^{\phi}_{-j}$ are equal in expectation. Thus, we can directly analyze the generalization error of ensemble estimators. <cit.> provides a decomposition of the generalization error into bias, variance, and covariance. When each individual estimator $\hat{f}^b$ is trained under the same baseline model $\mathcal{A}$, every estimator has the same bias, variance, and the covariance between  $\hat{f}^{b}$ and  $\hat{f}^{b'}$ is bounded by their variance. As a result, we can see that the MSE of $\hat{f}^{\phi}$ is smaller than any of individual estimators. As the consistency assumption is imposed on $\mathcal{A}$, we can find $\delta_{T,b},\gamma_{T,b}$ under different assumptions of $f$ (see examples <ref> and <ref>). Thus, $\hat{f}^{\phi}$ satisfies the small MSE assumption. By a similar analysis, we can also show the prediction error for $\hat{f}^{\phi}$ is smaller than that of any $\hat{f}^{b'}$, so $\hat{f}^{\phi}$ also satisfies the small pointwise prediction error assumption. In the appendix, we provide more details on how to relate the generalization error of $\hat{f}^{\phi}$ to that of $\hat{f}^b$.

In general, combining a diverse set of regressors with small or even negative covariance is ideal for having small MSE <cit.>. This, in turn, reduces uncertainty in prediction and makes the bound in Theorem 1 converge faster. Furthermore, assume $\mathcal{A}$ being used satisfies Assumption <ref>. From the analyses above, so long as $\phi$ is not so ill-chosen that the generalization error and expected prediction error of $\hat{f}^{\phi}$ exceeds the maximum error of individual regressors, $\hat{f}^{\phi}$ will satisfy Assumption <ref>. In particular, boosting techniques such as AdaBoost <cit.>, which improves on weak learners, can be used. Nevertheless, how to efficiently produce leave-$i$-out boosting predictors with the same theoretical guarantees as <cit.> remains unclear.

Lastly, based on the previous arguments, we can see why the same approximate coverage guarantee holds at every point after $T+1$. The whole sequence is assumed to have strongly mixing errors, and when each ensemble predictor approximates the unknown $f$ well, the same consistency condition carries over to all the training residuals.

Remark: The above analyses are valid because $|\hat{\epsilon}_i|, i \geq 1$ are not in-sample error when $\hat{f}^{\phi}_{-i}$ does not involve $i$ in training, so that $|\hat{\epsilon}_i-\epsilon_i| \not\approx |\epsilon_i|$. In fact, the theory in <cit.> disallows the use of in-sample error, so Jackknife discussed earlier is prohibited.

§ SIMULATED AND REAL-DATA EXAMPLES

For experiments, we aim to show that (1) |EnPI| is approximately valid for different types of data, with or without feature time-series $X_t$. Moreover, the coverage is robust under different input parameters; (2) we compare |EnPI| with other methods to demonstrate that our method generally produces more efficient intervals.

Simulation. We simulate three different types of time-series and consider four real datasets as follows:

(1) First, we generate a multivariate stationary time-series (Multi):

\begin{equation} \label{eq:simul}
Y_{t}=X_{t}^{\prime} \beta+\varepsilon_{t}, 1 \leq t \leq T+T_1,
\end{equation}

where $X_t,\beta \sim N(0,I_p)$ and independent over time. We require $\|\beta \|_2=1$. To induce serial correlation, we let $\varepsilon_{0}=0$ and $\varepsilon_{t}=\rho \varepsilon_{t-1}+\xi_{t}, \xi_{t} \stackrel{i i d}{\sim} N\left(0,1-\rho^{2}\right)$.

(2) Second, we generate an 1D non-stationary simple random walk (Rand) with deterministic drift:

\begin{equation} \label{randwalk}
Y_t=\alpha+Y_{t-1}+\epsilon_t, \epsilon_t \sim N(0,1), 1 \leq t \leq T+T_1,
\end{equation}

where we choose $X_t=(Y_{t-1},\ldots,Y_{t-d})$ to remember the past $d$ values.

(3) Lastly, we construct a network model (Network) with random yet sparse connections among $k$ nodes. Let $Y_i=(Y_{0,i},Y_{1,i},\ldots,Y_{T+T1,i})$ be the data collected at node $i$. Initially, $Y_{0,i} \sim N(0,1)$ for all $i\in [k]$. To specify neighbors of $i$, we first randomly select an integer $N_i$ from $0.2k$ to $0.4k$ and then randomly pick $N_i$ indices from $[k]$ to be the $i$'s neighbors, with indices in the set $S_{N_i}$. For each $1\leq t\leq T+T_1$, let

\begin{equation} \label{network}
Y_{ti}=\sum_{j \in S_{N_i}} \alpha_{ij} Y_{t-1,j} +\epsilon_t, \epsilon_t \sim N(0,1),
\end{equation}

where $X_{ti}=(Y_{t-d,1},\ldots,Y_{t-d,k},\ldots,Y_{t-1,1},\ldots,Y_{t-1,k}) \in \mathbb{R}^{dk}$ contains past $d$ values from all locations.

For all three simulated series, we let $T=T_1=200$. For Multi, we let $p=300, \rho=0.75$, so features are high-dimensional and error correlations are high. For Rand, we let $\alpha=2, d=5$, since $Y_t$ drifts away really fast. We remark that ARIMA can model Rand well. For Network, we let $k=d=10,\alpha_{ij}=1$ for all $i,j$. We remark $\alpha_{ij}=1$ implies a simple linear relationship between neighbors and nodes, which can be modeled well by linear regression; non-linear regression can model general node-neighbor relationship.

Real data. We also consider four different real datasets. Three of them are from the UCI machine learning repository, and the fourth is retrieved from the NSRDB website[dataset available at <https://nsrdb.nrel.gov/>]. For computation efficiency, we choose at most the first 10000 points of each series, for greenhouse gas observation <cit.>, appliances energy usage <cit.>, Beijing air quality <cit.>, and solar radiation level measured in DHI at Downtown Atlanta. Detailed descriptions can be found in the appendix. We use both multivariate versions ($X_t$ is exogenous time-series sequences) and univariate versions ($X_t$ is past history of the response series) of these datasets.

Comparison methods. We compare |EnPI| with two other CP methods and with ARIMA(10,1,10). The first CP method is the inductive conformal prediction (ICP) by <cit.>, since its prediction intervals are also valid by the argument in <cit.>. The other CP method is the weighted ICP (WeightedICP) proposed by <cit.>, which is proven to work when the test distribution shifts in proportion to the training distribution; it is generally more applicable than ICP.

We use logistic regression to estimate the weights for WeightedICP. We do a 50:50 split into proper training set and calibration set for ICP and WeightedICP and modify them to construct sequential prediction intervals as line 12-18 in |EnPI|; the length of sliding window is equal to that of the calibration data. Lastly, we use ARIMA implemented in Python's |statsmodel| package  under default parameter specification.

Input Parameters. We choose four baseline regression functions for $\mathcal{A}$: ridge, lasso, random forest (RF), and neural networks (NN). The first three are implemented in the Python |sklearn| library, and the last is built using the |keras| library. Below are the parameter specifications unless otherwise specified:

* For ridge and lasso, the penalty parameter $\alpha$ is chosen with generalized cross-validation over ten uniformly spaced grid points between 0.0001 to 10 (the package default $\alpha$ is 1). Higher $\alpha$ means more robust regularization.

* For RF, we build ten trees under the mean-squared-error (MSE) criterion. We restrict the maximum search depth of each tree as 2 for faster training. We only allow each tree to split into features rather than samples, so that combining RFs trained on subsets of the training data is reasonable for |EnPI|.

* For NN, we add three hidden layers, each having 100 hidden nodes, and apply 20$\%$ dropout after the second hidden layer to avoid overfitting. We use the Relu activation between hidden layers. The optimizer is Adam with a fixed learning rate of $5\times 10^{-4}$ under the MSE loss. Batch size is fixed at 100 with a maximum of 1000 epochs. We also use early stopping if there is no improvement in training error after ten epochs.

Since the three CP methods train on random subsets of training data, we repeat all experiments below for 10 trials, where each trial splits training data independently. On the other hand, ARIMA is deterministic given training data, so we only train it once. We let $\alpha=0.1$ and use 30$\%$ total data for training unless otherwise specified.

Lastly, we use |EnPI| under $\widetilde{B}=100$ and $\phi$ as taking the sample mean. Thus, each of the ensemble predictors in  |EnPI| is a leave-$i$-out bagging predictor. Note, $\widetilde{B}$ implies |EnPI| on average combines 37 predictors from the same regression model.

Results. Figure <ref> and <ref> show average coverage and width vs. $1-\alpha$ for simulated and real-data series, respectively. Results on all real datasets except the solar radiation levels appear in the appendix. Each average value is taken over $T_1$ test points and 10 independent trials. Lines labeled “Ensemble" denotes |EnPI| and green dash-dotted diagonal lines show the ideal average coverage at each $1-\alpha$, which are equally spaced between [0.05,0.95]. From Figure <ref>, we can see that none of the CP methods loses coverage at any $1-\alpha$ level for Multi and Network, and they have almost the same widths. On the other hand, the middle figure on Rand shows that ARIMA and WeightedICP lose coverage at certain values of $1-\alpha$ despite having similar average interval widths as |EnPI| and ICP. Meanwhile, intervals are slightly shorter on average when using |EnPI|. Additional plots under other regression algorithms are in the appendix.

On the other hand, Figure <ref> also shows that CP methods are approximately valid across different $1-\alpha$, with nearly identical average widths. These methods overcover at smaller $1-\alpha$ but never undercover. Meanwhile, although ARIMA produces shorter intervals than CP methods on the multivariate solar data, it may produce wider intervals on the univariate data and lose coverage when $1-\alpha$ is near 1. The performance of ARIMA is often less stable than CP methods.

Simulation: average coverage and width vs. $1-\alpha$ using ridge regression.

Solar data: average coverage and width vs. $1-\alpha$ using random forest regression.

In addition, Figure <ref> shows grouped coverage and width boxplots using ridge, RF, and NN. The x-axis in each group is the percentage of total data used for training, and each box shows results over ten trials. The dashed-dotted black line is the target 0.9 coverage. All boxes tightly center around the target coverage in terms of coverage, with |EnPI| having less variance. For achieving validity, these methods are all suitable for small-sample problems. In terms of width, |EnPI| shows much less variance than the other two methods regardless of the regression algorithm used. On the multivariate version, |EnPI| under different amounts of training data may or may not produce shorter intervals than the other two if ridge or RF is used; however, they are shorter using NN. Interval widths from NN are much shorter than the widths by ridge or RF, so using a complex regression algorithm on this highly non-stationary data shows the efficiency of |EnPI|. On the univariate version, |EnPI| almost always produces shorter intervals than the other two methods, with widths similar to those on the multivariate version. Similar analyses apply to these plots. In general, we notice that |EnPI| produces more efficient intervals than the other two methods when trained on a small percentage of total data.

Solar data: coverage (top two) and width (bottom two) boxplots over 10 trials.

Lastly, we summarize average coverage, width, and Winkler score[Let the upper and lower end of the prediction interval at time $t$ under level $\alpha$ be $L_t(\alpha)$, $U_t(\alpha)$, so width is $W_t(\alpha)=U_t(\alpha)-L_t(\alpha)$. Then, Winkler score (WS) is:

\[
(WS)_t= \begin{cases}
W_t(\alpha), & \text{if } L_t(\alpha) \leq y_t\leq U_t(\alpha)\\
W_t(\alpha)+2\cdot \frac{L_t(\alpha)-y_t}{\alpha}, & \text{if } y_t < L_t(\alpha)\\
W_t(\alpha)+2\cdot \frac{y_t-U_t(\alpha)}{\alpha}, & \text{if } y_t > U_t(\alpha)
\end{cases}
\]
] from using the same regression function on all four datasets in multiple tables, which appear in the appendix. Winkler score was used in <cit.> as a quantitative score to capture both validity and efficiency, and is particularly useful when a good coverage rate does not indicate an optimal width. In general, we observe that all three CP methods are valid in almost all cases. Nevertheless, |EnPI| clearly outperforms the other two CP methods in higher average Winkler scores and shorter interval widths. Comparing to ARIMA, |EnPI| performs significantly worse on the multivariate Appliances Energy and Beijing air datasets, but is better on the other two datasets. In particular, the performance gaps become smaller or even disappear when |EnPI| is used on univariate versions of these datasets. We think such behaviors occur because the feature series collected for some datasets are less predictive of the response series than the response series's history. ARIMA's noticeable drawback is that it significantly undercovers on the Greenhouse Gas dataset and overcovers on the Beijing air dataset. Lastly, we remark that |EnPI| performs worse on the univariate than the multivariate version of the Greenhouse Gas dataset, likely because the response series is assumed to come from a weighted combination of its multivariate feature series.

From these results, we see that although theory guarantees valid coverage only when $\mathcal{A}$ is consistent, empirical results are valid even under potentially misspecified models. Nevertheless, regression algorithms closer to the true model produce more efficient intervals. Meanwhile, another crucial part for attaining efficient intervals is properly defining $X_t$ as a predictor for $Y_t$. Based on real-data results over multivariate and univariate versions of the same dataset, we sometimes see that a proper choice of $x_t$ affects the empirical performances more than the choice of $\mathcal{A}$. We think this is particularly appealing for the method to be applicable on a broader scale: even without $\mathcal{A}$ that closely represent the true data-generating model, researchers with suitable choices of feature $X_t$, which often require deep domain knowledge, can use |EnPI| to produce approximately valid and reasonably efficient intervals for their problems of interest.

§ CONCLUSION

In this paper, we present a predictive inference method titled |EnPI| for dynamic time-series. The method is approximately valid and asymptotically exact under mild assumptions on the estimation quality of the baseline regression algorithm and the error process's behavior. Computationally, |EnPI| is an efficient wrapper around any regression algorithm, especially for ensemble methods. Through extensive experiments, we demonstrate the versatility of |EnPI| on a wide range of time-series data, including non-stationary time-series and network data. These experiments also show that |EnPI| provides shorter intervals than other common predictive inference methods, making it preferable in reality.

In the future, we aim to extend the problem along these directions. Firstly, move from average validity to conditional validity, as in <cit.>. Secondly, bound the interval widths and quantify closeness to oracle prediction intervals, as in <cit.> for the iid case. Thirdly, modify the method for problems with missing data and quantify how far ahead the method can construct prediction intervals with validity guarantee. Fourthly, generalize this problem to high-dimensional response series. Lastly, connect conformal prediction with change-point detection, as in <cit.>, leading to robust and computationally efficient change-point detection techniques.

Conformal prediction interval for dynamic time-series:

Supplementary Materials

§ ADDITIONAL DERIVATIONS

Assume each $\hat{f}^{b}$ is trained under the same regression model $\mathcal{A}$. We show in detail two things. First, we show why the generalization error (measured in terms of mean square error) for the ensemble estimator $\hat{f}^{\phi}$ is always smaller than that of each $\hat{f}^{b}$. Second, we show why the expected prediction error of $\hat{f}^{\phi}$ is also smaller. For simplicity, we prove for $\phi$ as the sample mean with equal weights, but the proof idea easily generalizes to taking sample mean with unequal weights.

Consider a generic case, where we have $T$ training data $\{z_i\}_{i=1}^T, z_i=(x_i,y_i)$ and wish to evaluate the generalization error of a single estimator $\hat{f}^b, b=1,\ldots,B$ on a new point $Z_0$. We let $z^T=\{z_i\}_{i=1}^T$. In our case, each $\hat{f}^b$ is a pre-ensemled estimator trained on a subsample with replacement from $z^T$. In expectation over random draws of subsamples, each $\hat{f}^b$ is hence trained on the same number of the data from $z^T$. Therefore, we define $\text{GErr}(\hat{f}^b)=\mathbb{E}_{Z_0}\mathbb{E}_{Z^T}[(Y_0-\hat{f}^b(X_0)]$. Now, under our assumption that $Y_i=f(X_i)+\sigma^2$, <cit.> shows that

\[
\text{GErr}(\hat{f}^b)=\mathbb{E}_{X_0}\{\text{Var}(\hat{f}^b|X_0)+\text{Bias}(\hat{f}^b|X_0)^2\}+\sigma^2,\]

where $\text{Var}(\hat{f}^b|X_0)$ and $\text{Bias}(\hat{f}^b|X_0)$ are conditional variances and bias given $X_0=x_0$:

\begin{align*}
& \text{Var}(\hat{f}^b|X_0)= \mathbb{E}_{Z^T}\left[\left(\hat{f}^b(X_0|Z^T)-\mathbb{E}_{Z^T}[\hat{f}^b(X_0|Z^T)])^2\right)\right]\\
& \text{Bias}(\hat{f}^b|X_0)=\mathbb{E}_{Z^T}[\hat{f}^b(X_0|Z^T)]-f(X_0).
\end{align*}

It is easy to observe that when each $\hat{f}^b$ is trained under the same model on subsets of $z^T$, the bias and variance of different models are the same. Namely, for any $b\neq b'$ we have

\begin{align*}
&\text{Var}(\hat{f}^b|X_0)= \text{Var}(\hat{f}^{b'}|X_0) \\
& \text{Bias}(\hat{f}^b|X_0)= \text{Bias}(\hat{f}^{b'}|X_0)
\end{align*}

Thus, we call the conditional variance and bias for pre-ensemble estimators as $\text{Var}(\hat{f}|X_0)),\text{Bias}(\hat{f}|X_0)$, respectively. By Corollary 1 in <cit.>, we see that

\begin{align*}
\text{GErr}(\hat{f}^{\phi})=\mathbb{E}_{X_0}\{ \frac{1}{B}\text{Var}(\hat{f}|X_0)&+(1-\frac{1}{B})\text{Cov}(X_0)\\&+\text{Bias}(\hat{f}|X_0)^2 \}+\sigma^2,
\end{align*}

\[\text{Cov}(X_0)=\frac{1}{B(B-1)}\sum_{b}\sum_{b'\neq b} \text{Cov}(\hat{f}^b,\hat{f}^{b'}|X_0).
\]

Now, from the inequality that $\text{Cov}(\hat{f}^b,\hat{f}^{b'}|X_0) \leq \sqrt{\text{Var}(\hat{f}^b|X_0)\text{Var}(\hat{f}^{b'}|X_0)}=\text{Var}(\hat{f}|X_0),$  we see that

\begin{align*}
\text{GErr}(\hat{f}^{\phi}) & \leq \mathbb{E}_{X_0}\{ \frac{1}{B}\text{Var}(\hat{f}|X_0)+(1-\frac{1}{B})\text{Var}(\hat{f}|X_0)\\&+\text{Bias}(\hat{f}|X_0)^2 \}+\sigma^2 \\
& = \text{GErr}(\hat{f}).
\end{align*}

Thus, $\text{GErr}(\hat{f}^{\phi})\leq \text{GErr}(\hat{f}^{b})$ for any individual pre-ensemble estimator.

Now, we prove the expected error of $\hat{f}^{\phi}$ is smaller as well. Note that the expected error $\mathbb{E}_{X_0}[\mathbb{E}_{Z^T}[|Y_0-\hat{f}^b(X_0|Z^T)|]]$ for different $b$ is the same as well, so we call it $\text{PErr}(\hat{f})$.

\begin{align*}
& \text{PErr}(\hat{f}^{\phi}) \\
= & \mathbb{E}_{X_0}\mathbb{E}_{Z^T}\left[|Y_0-\hat{f}^{\phi}(X_0|Z^T)|\right] \\
= & \mathbb{E}_{X_0}\mathbb{E}_{Z^T}\left[|Y_0-\frac{1}{B}\sum_b\hat{f}^b(X_0|Z^T)|\right]\\
\leq & \mathbb{E}_{X_0}\mathbb{E}_{Z^T}\left[\frac{1}{B}\sum_b|Y_0-\hat{f}^b(X_0|Z^T)|\right]\\
= & \frac{1}{B}\sum_b \mathbb{E}_{X_0}\mathbb{E}_{Z^T}\left[|Y_0-\hat{f}^b(X_0|Z^T)|\right]\\
= & \mathbb{E}_{X_0}\mathbb{E}_{Z^T}\left[|Y_0-\hat{f}(X_0|Z^T)|\right],
\end{align*}

where the inequality follows from multi-dimensional triangle inequality.

§ ADDITIONAL EXPERIMENTS

§.§ Real Datasets

Descriptions:  The first series contains Greenhouse Gas observation (Greenhouse) <cit.> from 5.10 till 7.31, 2010, with four samples every day and 6 hours apart between data points. The goal is to find the optimal weights for the 15 observation series to match the synthetic control series. The second dataset contains appliances energy usage (Appliances) <cit.>. Consecutive data points are 10 minutes apart for about 4.5 months. We can use 27 different humidity and temperature indicators to predict the appliances' energy use in Wh. The third dataset on Beijing air quality (Beijing air) <cit.> contains air pollutants data from 12 nationally-controlled air-quality monitoring sites. The data is from 3.1, 2013 to 2.28, 2017. The goal is to predict PM2.5 air pollutant levels using 10 different air pollutants and meteorological variables. We use the data from the Tiantan district. The last dataset contains solar radiation level (Solar) measured in DHI at Downtown Atlanta, where consecutive data points are 30 minutes apart from the beginning of 2017 to the end of 2018. We have 7 other variables such as temperature, solar zenith angle, wind speed, etc. for predicting DHI levels. In general, these series can be very non-stationary, especially solar radiation and Beijing air quality.

Raw Data Plot:  Figure <ref> shows the raw data plot for Beijing Air and Solar radiation datasets, for the first 2000 data points of their response time series. We can see periodic fluctuations of different magnitude, which actually persist throughout the whole dataset. These fluctuations indicate changing mean and variance, with certain seasonal patterns that leads to clear non-stationary in the data.

Plot of first 2000 response data points for Beijing air and Solar radiation datasets.

§.§ Simulation Results

Figure <ref> shows additional average coverage and width vs $1-\alpha$ using other three regression algorithms: lasso, RF, and NN.

* For Multi (top row), we see that coverage is always valid at different $1-\alpha$ and |EnPI| returns significantly shorter width under random forest than the other two CP methods. It returns slightly shorter widths using NN as well. It is noticeable that average widths returned by RF and NN are one magnitude wider than those from lasso (or ridge). This happens because the underlying data generating process assumes $Y_t$ is a linear function of $X_t$, so linear regression functions model the data the best.

* For Rand (middle row), we see that coverage behavior using lasso is the same as that from ridge (Figure 1 in the main text). However, the average widths are all wider than those from ridge and ARIMA returns the shortest average widths. This is reasonable because the dataset can be well modeled by ARIMA and there is not a sparse relationship of $Y_t$ depending on its past. On the other hand, CP methods by RF and NN completely lose coverage, although |EnPI| always return shorter intervals than the other two CP methods. Such phenomenons further demonstrate that using complicated models to model data with simple data generating process can be very unsatisfactory.

* For Network, lasso and RF maintain valid coverage at all $1-\alpha$ and RF returns shorter widths for |EnPI|, confirming the efficiency of its prediction intervals. On the other hand, there is under-coverage using NN for all three CP methods, but those from |EnPI| still return the shortest intervals among all three CP methods.

Overall, we remark that CP methods under inappropriate choices of the regression algorithm lead to poor coverage, specially when the data are generated using simple models yet the regression models are complex. Nevertheless, at the same level of coverage, |EnPI| almost always return shorter intervals on average than the other two CP methods in comparion.

Simulation: average coverage and width vs. $1-\alpha$ using other regression algorithms.

§.§ Real-data Results

We first present additional grouped boxplots for the other three datasets (both multivariate and univariate versions) as in Figure 3 of the main text. Then, we show tables for using different regression algorithms on all real datasets (both univariate and multivariate versions), where the interpretation of these results is described near the end of the experiment section of the main text. Lastly, it is natural to use recurrent neural networks (RNN) for modeling time-series data, so we compare the performance between NN and RNN.

(1) We observe the following trends from the grouped boxplots:

* Figure <ref> (Greenhouse) shows slightly under-coverage at a low percentage of training data, but the coverage is approximately valid at higher percentages for all three regression algorithms on both multivariate and univariate versions. Coverage values for multivariate data are closer to the target 0.9 dash-dotted line because the goal of using the Greenhouse gas dataset is to find the optimal weights for predicting the response time-series using exogenous multivariate time-series. On the other hand, intervals are almost always shorter for |EnPI| using RF and NN. They may not be shorter using ridge but the boxes are much tighter, reflecting greater stability of |EnPI| than the other two

* Figure <ref> (Appliances) reveals similar but tighter coverage patterns. In particular, coverage is more valid on the univariate version, likely because the past history of energy use is more predictive of future energy use than the exogenous variables such as humidity and temperature of the ambient space (e.g. kitchen, bathroom, living room, etc.) On the other hand, intervals by NN and |EnPI| are the shortest among all combinations, especially at lower percentages of training data. This reflects that energy use is better modeled by complicated models than simpler models, especially when the sample size available for training is small relative to the test period. In general, |EnPI| returns more stable and shorter intervals than the other two methods, and more complex regression algorithms perform better on the univariate version of the dataset.

* Figure <ref> (Beijing air) shows similar but even tighter coverage patterns. There are no significant differences between using the multivariate or the univariate version, especially for complex regressors such as RF and NN. In terms of width, NN and RF for |EnPI| are much shorter than those for the other two CP methods and there are no major differences between using multivariate or univariate versions. All CP methods using ridge are more stable at higher percentages of training data on the univariate version, and |EnPI| is still the tightest among all three.

In general, we think |EnPI| is stable across different combinations of regression algorithms and dataset. For interval efficiency, datasets with possibly complicated relationships between exogenous feature series and the response series (e.g. Appliances and Beijing air datasets) are better modeled by complex regression algorithms such as RF and NN. However, using complex regression algorithms on datasets with simple relationships (e.g. Greenhouse) is unfavorable. When the relationship between exogenous feature series and the response series is unclear, using the history of the response to predict its future values produces valid and equally efficient intervals.

(2) We now show the tables. Table <ref> is for ridge, <ref> for lasso, <ref> for RF, and <ref> for NN. To ensure the regression algorithms have better prediction performance, we double the number of $\alpha$ chosen between $0.0001$ and 10 for ridge and lasso. We also double the number of trees and uplift the restriction on the maximum depth of each tree for RF. Doing so dramatically reduces average interval widths without affecting validity. For CP methods, results are averaged over 10 trials. We remark that we consider the performance on both univariate and multivariate versions together so that there is only one set of bold numbers (width and Winkler score) for each dataset.

(3) We use |keras| to construct RNN with the following parameter specifications:

* We add two hidden LSTM layers, followed by a dense output layer. Each LSTM layer has 100 hidden neurons, so the output from the first hidden layer is fed into the second hidden layer. We use the Tanh activation function for these two hidden layers and the Relu activation function for the dense layer. The optimizer is Adam with a fixed learning rate of $5\times 10^{-4}$ under the MSE loss. The batch size is fixed at 100 with a maximum of 100 epochs. We use early stopping if there is no improvement in training error after 10 epochs.

We summarize the results in grouped boxplots as we did before.

* Figure <ref> shows result for the Greenhouse dataset. We see that these two regression algorithms using different CP methods show similar average coverage. In terms of width, |EnPI| using NN produces shorter intervals on the multivariate version and |EnPI| using RNN produces shorter intervals on the univariate version.

* Figure <ref> shows result for the Appliances dataset. We see that RNN and NN behave very similarly in terms of coverage on the univariate version, but RNN produces more valid intervals than NN on the multivariate version. In terms of width, |EnPI| using NN produces shorter intervals on the multivariate version and |EnPI| produces almost the same intervals on the univariate version, regardless of using NN or RNN.

* Figure <ref> shows result for the Beijing air dataset. As before, coverage by both regression methods is very similar, but RNN produces more valid intervals than NN on the multivariate version. In terms of width, |EnPI| using RNN produces shorter intervals than NN regardless of both multivariate and univariate versions.

* Figure <ref> shows result for the Solar dataset. Both regression algorithms behave very similarly in terms of coverage, but |EnPI| using NN may produce shorter intervals on the multivariate version.

In general, we think |EnPI| (or CP methods in general) using RNN produces very competitive results as it using NN does. The intervals can even be more valid and efficient for certain datasets. We think these results further show that |EnPI| under complex neural networks can produce prediction intervals for these networks in a principled way, thus quantifying the uncertainty when using these methods. Lastly, we note that when the neural network depends on the order of the training data indices (e.g. RNN), ensemble-based methods such as |EnPI| that subsample training data at random may need revision so as to ensure the prediction performance of these neural networks.

Greenhouse data: coverage (left two) and width (right two) boxplots over 10 trials.

Appliances data: coverage (left two) and width (right two) boxplots over 10 trials.

Beijing air data: coverage (left two) and width (right two) boxplots over 10 trials.

Ridge: the smallest widths and Winkler scores among all methods are in bold.

3cMultivariate     3cUnivariate

Dataset     Method     Coverage     Width     Winkler Score     Coverage     Width     Winkler Score

4*[l]Greenhouse     ARIMA     0.78    122.14     5.75e+04

Ensemble     0.88     55.35     1.59e+04     0.91     139.72     4.15e+04

ICP     0.89     55.39     1.64e+04    0.90     147.77     4.75e+04

WeightedICP     0.89     55.65     1.62e+04     0.90     166.69     5.07e+04

4*[l]Appliances     ARIMA     0.93     236.25     2.71e+06

Ensemble     0.88     345.30     3.95e+06     0.90     147.10     2.56e+06

ICP     0.88     340.46     3.92e+06     0.90     150.55     2.57e+06

WeightedICP     0.89     345.40     3.96e+06    0.90     151.72     2.57e+06

4*Beijing air     ARIMA     0.95     77.55     7.38e+05

Ensemble     0.90     102.39     1.03e+06      0.90     54.97     6.90e+05

ICP     0.90     103.36     1.05e+06     0.90     57.01     6.97e+05

WeightedICP     0.90     101.97     1.03e+06      0.90     56.91     6.98e+05

4*Solar     ARIMA     0.87     145.18     2.19e+06

Ensemble     0.91     299.09     2.41e+06     0.91     192.00     2.04e+06

ICP     0.91     293.54     2.33e+06    0.90     189.96     1.99e+06

WeightedICP     0.91     292.32     2.33e+06     0.90     193.01     1.98e+06

Lasso: the smallest widths and Winkler scores among all methods are in bold.

3cMultivariate     3cUnivariate

Dataset     Method     Coverage     Width     Winkler Score     Coverage     Width     Winkler Score

4*[l]Greenhouse     ARIMA     0.78     122.14     5.75e+04

Ensemble     0.89     54.66     1.56e+04     0.90     135.43     4.06e+04

ICP     0.89     54.18     1.63e+04    0.89     150.57     4.76e+04

WeightedICP     0.89     56.52     1.65e+04      0.90     158.21     4.76e+04

4*[l]Appliances     ARIMA     0.93     236.25     2.71e+06

Ensemble     0.87     347.47     3.96e+06      0.90     146.33     2.56e+06

ICP     0.89     329.06     3.89e+06    0.90     148.12     2.56e+06

WeightedICP     0.89     322.42     3.92e+06      0.90     149.39     2.56e+06

4*Beijing air     ARIMA     0.95     77.55     7.38e+05

Ensemble     0.90     102.22     1.03e+06     0.90     55.02     6.89e+05

ICP     0.90     102.23     1.04e+06     0.90     56.72     6.95e+05

WeightedICP     0.90     102.37     1.03e+06     0.90     56.61     6.93e+05

4*Solar     ARIMA     0.87     145.18     2.19e+06

Ensemble     0.91     299.39     2.41e+06     0.91     192.26     2.04e+06

ICP     0.91     293.73     2.33e+06     0.90     193.69     1.98e+06

WeightedICP     0.91     293.87     2.34e+06     0.90     191.54     2.00e+06

Random Forest: the smallest widths and Winkler scores among all methods are in bold.

3cMultivariate     3cUnivariate

Dataset     Method     Coverage     Width     Winkler Score     Coverage     Width     Winkler Score

4*[l]Greenhouse     ARIMA     0.78     122.14     5.75e+04

Ensemble     0.88     65.75     2.20e+04     0.91     118.70     3.83e+04

ICP     0.89     76.79     2.70e+04     0.90     153.68     5.20e+04

WeightedICP     0.90     79.72     2.73e+04     0.91     173.49     5.61e+04

4*[l]Appliances     ARIMA     0.93     236.25     2.71e+06

Ensemble     0.82     646.01     7.19e+06     0.89     205.04     2.81e+06

ICP     0.87     767.71     8.24e+06     0.90     271.47     3.58e+06

WeightedICP     0.88     804.99     8.40e+06     0.90     278.31     3.65e+06

4*Beijing air     ARIMA     0.95     77.55     7.38e+05

Ensemble     0.89     92.18     9.99e+05    0.90     63.36     7.99e+05

ICP     0.90     122.93     1.33e+06     0.90     89.04     1.11e+06

WeightedICP     0.90     126.14     1.41e+06     0.90     91.79     1.17e+06

4*Solar     ARIMA     0.87     144.79     2.19e+06

Ensemble     0.90     198.62     1.96e+06     0.91     193.29     1.89e+06

ICP     0.90     237.34     2.28e+06     0.90     237.54     2.37e+06

WeightedICP     0.90     248.94     2.35e+06    0.91     237.71     2.34e+06

Neural Network: the smallest widths and Winkler scores among all methods are in bold.

3cMultivariate     3cUnivariate

Dataset     Method     Coverage     Width     Winkler Score     Coverage     Width     Winkler Score

4*[l]Greenhouse     ARIMA     0.78     122.14     5.75e+04

Ensemble     0.88     51.18     1.57e+04     0.83     131.12     5.35e+04

ICP     0.88     55.09     1.67e+04     0.84     133.31     5.03e+04

WeightedICP     0.88     56.72     1.67e+04     0.86     139.63     5.02e+04

4*[l]Appliances     ARIMA     0.93     236.25     2.71e+06

Ensemble     0.88     303.76     3.75e+06     0.89     169.01     3.09e+06

ICP     0.89     331.67     4.12e+06     0.89     185.74     3.16e+06

WeightedICP     0.89     325.41     4.12e+06     0.89     186.03     3.14e+06

4*Beijing air     ARIMA     0.95     77.55     7.38e+05

Ensemble     0.89     103.78     1.14e+06     0.89     91.67     1.00e+06

ICP     0.89     113.75     1.26e+06     0.90     100.65     1.07e+06

WeightedICP     0.89     116.05     1.28e+06     0.90     100.14     1.07e+06

4*Solar     ARIMA     0.87     145.18     2.19e+06

Ensemble     0.89     194.99     2.02e+06     0.89     206.98     2.40e+06

ICP     0.89     209.03     2.05e+06     0.89     221.71     2.42e+06

WeightedICP     0.90     210.21     2.05e+06    0.89     222.02     2.41e+06

Greenhouse data: coverage (left two) and width (right two) boxplots over 10 trials.

Appliances data: coverage (left two) and width (right two) boxplots over 10 trials.

Beijing air data: coverage (left two) and width (right two) boxplots over 10 trials.

Solar data: coverage (left two) and width (right two) boxplots over 10 trials.