Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
99,132 characters · 25 sections · 32 citation commands
\graphicspath{{figures/}}
Conformal prediction VovkGammermanShafer2005, AngelopoulosBates2023 constructs distribution-free prediction intervals with exact finite-sample coverage guarantees under exchangeability: when the joint distribution of calibration and test observations is invariant to permutations, the probability that $Y_{T+h}$ falls outside the interval is exactly $\alpha$, for every finite sample size. This exact guarantee does not extend to time series. Serial dependence, volatility clustering, and gradual distributional drift all violate exchangeability, so classical split conformal prediction may deliver coverage substantially different from the nominal $1-\alpha$ level. The practitioner's natural response is rolling-origin calibration: estimate the conformal quantile from the $m$ most recent pseudo-out-of-sample errors, so that the calibration sample tracks the current error distribution rather than averaging over a possibly stale history. This approach is intuitively compelling and widely used XuXie2021, XuXie2023, but the window length $m$ is typically chosen by informal cross-validation or left as a tuning parameter without theoretical guidance, and there is no existing result quantifying how close coverage is to $1-\alpha$ as a function of $m$, $T$, and the degree of non-stationarity.
The present paper addresses this gap by proposing and analysing rolling-origin conformal prediction as a method for time-series forecasting. We formalise the construction, characterise its coverage properties, and derive principled guidance on the window length $m$. Rather than seeking an exact finite-sample guarantee --- which cannot hold without exchangeability --- we establish precisely how much coverage deviates from $1-\alpha$. Under local stationarity of the score distribution Dahlhaus1997, uniform $\alpha$-mixing Rio1993, DoukhanMassartRio1994, quantile regularity, and a near-stationarity condition $m = o(T)$, we prove a four-term decomposition of the coverage deviation $\Delta_{T,m,h} := \lvert \mathbb{P}(Y_{T+h} \in \widehat{C}_{T+h|T}(1-\alpha)) - (1-\alpha) \rvert$ and derive the minimax-optimal window rule $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ with coverage-error rate $O(T^{-\beta/(2\beta+1)})$ (recovering $T^{2/3}$ and $O(T^{-1/3})$ at $\beta=1$). Each term in the decomposition has a direct empirical counterpart and a clear operational meaning. The deviation bound is two-sided: it controls both undercoverage ($\mathbb{P} < 1-\alpha$) and overcoverage ($\mathbb{P} > 1-\alpha$), reflecting that neither direction of error is structurally excluded in the time-series setting.
\paragraph{Relation to existing work.} GibbsCandes2021 propose adaptive conformal inference, which updates the quantile online via a fixed gradient step; the update rule tracks distribution shift but has no explicit bound on the instantaneous coverage error or guidance on the step size. BarberCandesRamdas2023 establish the general theory of conformal prediction beyond exchangeability, providing marginal validity under weighted exchangeability conditions; their framework does not address the rolling calibration setting or the window-length choice. The rolling-calibration approach is studied by XuXie2021, XuXie2023 in an online sequential setting, but without a finite-sample decomposition or an optimal-window result. Our contributions are: (i) a complete formulation of rolling-origin conformal prediction for non-exchangeable time series; (ii) a four-term coverage-error decomposition under H\"{o}lder-$\beta$ drift with minimax-optimal window $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ (Theorems (ref) and (ref)); (iii) Bahadur representations under both $\alpha$-mixing (Proposition (ref)) and physical dependence (Proposition (ref)); (iv) an oracle inequality for the implemented data-driven window selector, with the required concentration condition proved in Appendix (ref) (Theorem (ref)); and (v) empirical validation on six real series and a stratified M4 sample confirming the theoretical scaling and demonstrating consistent Winkler-score gains over full-history calibration.
The local-stationarity framework follows Dahlhaus1997 and Vogt2012. The mixing theory relies on Rio1993 and DoukhanMassartRio1994. The Bahadur representation uses the Bernstein inequality of MerlvedePeligradRio2009. Mixing properties of the GARCH and tvARCH processes relevant to the empirical application are established in FryzlewiczSubbaRao2011.
The remainder of the paper is organised as follows. Section (ref) presents the method in full: Algorithm (ref) (base) and Algorithm (ref) (volatility-scaled), practical guidance on window selection and score choice, and the oracle decomposition used in the theory. Section (ref) states the assumptions, proves the main coverage-error theorem and minimax lower bound, establishes Bahadur representations under $\alpha$-mixing and physical dependence, derives the optimal window rule, and provides extensions to polynomial mixing and volatility-scaled scores. Section (ref) presents the empirical analysis. Section (ref) concludes. Appendix (ref) contains the self-contained proof of the Bahadur representation (Proposition (ref)). Appendix (ref) proves the uniform concentration condition required for the oracle inequality (Theorem (ref)).
Let $(Y_t)_{t=1}^{T}$ be a univariate time series observed up to time $T$. A forecasting model produces an $h$-step-ahead point prediction $\hat{Y}_{t+h|t}$ using only information available through time $t$. The rolling-origin conformal score is the absolute forecast error \[ \hat{S}_{t,h} := \lvert Y_{t+h} - \hat{Y}_{t+h|t} \rvert. \] Rolling-origin conformal prediction constructs a prediction interval at the target origin $T$ by treating the $m$ most recent scores as a local calibration sample and taking their empirical $(1-\alpha)$-quantile as the half-width. The complete procedure is stated in Algorithm (ref).
The empirical calibration distribution and quantile are formally \[ \widehat{F}_{T,m,h}(x) := \frac{1}{m}\sum_{j=1}^{m} \mathbf{1}_{\{\hat{S}_{T-j,h} \leq x\}}, \qquad \hat{q}_{T,m,h} := \inf\bigl\{x : \widehat{F}_{T,m,h}(x) \geq 1-\alpha\bigr\}, \] and the rolling-origin conformal interval is
Coverage analysis reduces to studying $\mathbb{P}(\hat{S}_{T,h} \leq \hat{q}_{T,m,h})$, since $Y_{T+h} \in \widehat{C}_{T+h|T}(1-\alpha) \iff \hat{S}_{T,h} \leq \hat{q}_{T,m,h}$.
When nonstationarity is primarily driven by time-varying scale (e.g.\ GARCH volatility clustering), normalising the scores by an estimated conditional volatility reduces the effective drift rate $L_{F,h}$ and can substantially narrow the intervals. Algorithm (ref) implements this variant; it differs from Algorithm (ref) only in Steps 1 and 3--5.
The theoretical effect of volatility scaling is characterised in Corollary (ref): replacing $L_{F,h}$ with $L^{sc}_{F,h} \ll L_{F,h}$ in the coverage bound reduces the drift term and shifts the optimal window toward larger $m$.
\paragraph{Window selection.} Theorem (ref) establishes that the theoretically optimal window is $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ (equal to $T^{2/3}$ at $\beta=1$), but this result depends on unknown constants ($L_{F,h}$, $A_h(\infty)$) that are not directly estimable from the data. In practice we recommend Winkler cross-validation: evaluate the Winkler score on a held-out validation fold of the pseudo-out-of-sample errors for each candidate $m$ drawn from a fine grid in $[0.1, 4.0] \times T^{2\beta/(2\beta+1)}$ (using $\beta = 1$ as default when $\beta$ is unknown), and select the minimiser. This is a consistent estimator of $m^{\star}$ in the sense that the cross-validation selected $\hat{m}$ achieves the same asymptotic coverage-deviation rate as $m^{\star}$, since the Winkler score is a proper scoring rule whose population minimiser coincides with the optimal $m^{\star}$.
\paragraph{Score function.} Algorithm (ref) uses the absolute error $\hat{S}_{t,h} = \lvert Y_{t+h} - \hat{Y}_{t+h|t} \rvert$ as the conformal score. This is natural for symmetric, unimodal conditional distributions. Other score functions are compatible with the framework: the signed error $Y_{t+h} - \hat{Y}_{t+h|t}$ produces one-sided intervals; a normalised residual $\hat{S}_{t,h}/\hat{\sigma}_{t,h}$ reduces to the VS-ROCP of Algorithm (ref); quantile-regression residuals $\max(\tau(Y_{t+h}-q_{t,h}(\tau)), (\tau-1)(Y_{t+h}-q_{t,h}(\tau)))$ produce asymmetric intervals. The theoretical analysis of Section (ref) applies to any score satisfying Assumptions (ref)--(ref).
\paragraph{Forecasting model.} The method is model-agnostic: $\mathcal{M}$ in Algorithm (ref) can be any procedure that produces a point prediction $\hat{Y}_{t+h|t}$ using only past data. The theoretical cost of using an estimated rather than oracle model appears in Term (IV) of Theorem (ref) via the pair $(r_T, \eta_T)$ in Assumption (ref). For AR($p$) models, $r_T = O(T^{-1/2})$ and $\eta_T = 0$ under standard mixing conditions. For ARMA-GARCH, the same holds under the conditions of FryzlewiczSubbaRao2011. For nonlinear or machine-learning models, Term (IV) must be bounded case by case, but the remaining three terms of the decomposition are unaffected by the model choice.
\paragraph{Relationship to rolling-origin evaluation.} Rolling-origin evaluation is the standard protocol for assessing time-series forecast accuracy TashmanFildes2000, and the pseudo-out-of-sample scores $\hat{S}_{t,h}$ in Step 1 of Algorithm (ref) are exactly the quantities computed in any rolling-origin evaluation exercise. Rolling-origin conformal prediction therefore adds no computational overhead to a forecasting pipeline that already performs rolling-origin evaluation: the conformal interval is obtained by taking the empirical quantile of the scores that are already being recorded, restricted to the most recent $m^{\star}$ of them.
\paragraph{Oracle decomposition.} To separate forecasting error from calibration error, fix an oracle predictor $f^{\star}_{t,h}$ (e.g.\ the true conditional mean, not computable in practice) and define the oracle score $S^{o}_{t,h} := \lvert Y_{t+h} - f^{\star}_{t,h} \rvert$ with marginal distribution $F_{t,h}(x) := \mathbb{P}(S^{o}_{t,h} \leq x)$. The corresponding oracle empirical distribution and quantile are \[ \widehat{G}_{T,m,h}(x) := \frac{1}{m}\sum_{j=1}^{m} \mathbf{1}_{\{S^{o}_{T-j,h} \leq x\}}, \qquad \hat{q}^{o}_{T,m,h} := \inf\bigl\{x : \widehat{G}_{T,m,h}(x) \geq 1-\alpha\bigr\}, \] and the average calibration-window distribution and its $(1-\alpha)$-quantile are \[ \bar{F}_{T,m,h}(x) := \frac{1}{m}\sum_{j=1}^{m} F_{T-j,h}(x), \qquad q^{\circ}_{T,m,h} := \inf\bigl\{x : \bar{F}_{T,m,h}(x) \geq 1-\alpha\bigr\}. \] The quantity $q^{\circ}_{T,m,h}$ is the population quantile of the time-averaged calibration-window distribution. It differs from the current quantile of $F_{T,h}$ by an amount proportional to the local drift rate $L_{F,h}$ and the window length $m$, which is the source of the bias term in Theorem (ref).
Before stating the main result it is worth pausing to understand what the theory needs to control, because the four assumptions correspond directly to the four terms that appear in Theorem (ref). Each assumption is minimal in the sense that dropping it would make the corresponding term uncontrollable.
The first condition concerns how the distribution of forecast errors changes over time. In a perfectly stationary series the calibration scores from two years ago are just as informative as those from last week. Real series are not stationary: business cycles shift macroeconomic volatility, GARCH dynamics cluster financial volatility, and structural change alters the unconditional error distribution entirely. We capture this through a H\"{o}lder-$\beta$ smoothness condition on the score distributions $F_{t,h}$, following Dahlhaus1997.
The parameter $\beta$ controls how smoothly the error distribution evolves. When $\beta = 1$ (Lipschitz drift), a gap of $k$ time steps produces a distributional shift of order $k/T$ --- linear in time. Larger $\beta$ describes smoother evolution, smaller $\beta$ more erratic change. The constant $L_{F,h}$ is the drift rate: a series undergoing rapid structural change has large $L_{F,h}$, while a slowly drifting macro series has small $L_{F,h}$. Both parameters are unobservable, which is precisely why the window $m$ cannot be set optimally without data-driven selection; the main theorem makes their role in the optimal window explicit.
The second condition controls the dependence structure within the calibration window. Forecast errors from adjacent time periods are typically correlated: an AR model that overpredicts today will likely overpredict tomorrow. If this dependence is too strong, the $m$ calibration scores effectively contain less than $m$ independent pieces of information, and the empirical quantile is more variable than the $m^{-1/2}$ rate one might naively expect.
The quantity $A_h(m)$ is the effective dependence burden of the calibration window: a summable mixing sequence has $A_h(\infty) < \infty$ (the window's information content grows at the standard $m^{-1/2}$ rate), while a polynomially mixing sequence with $\alpha_h(k) \asymp k^{-a}$ has $A_h(m) \asymp m^{1-a}$ (the effective rate degrades to $m^{-a/2}$, slower for strongly dependent processes). In either case, $A_h(m)$ enters the bound as a multiplicative inflation of the quantile-noise term. All ARMA and GARCH processes with appropriate parameter restrictions are summably mixing FryzlewiczSubbaRao2011.
Bounding the quantile-noise term also requires that the score distribution does not have a flat density near the quantile of interest. If the density were zero at the $(1-\alpha)$-quantile, small deviations in the empirical CDF would translate into large deviations in the quantile --- the empirical quantile would be extremely sensitive to estimation noise.
The upper bound $\bar{f}_{h}$ and lower bound $\underline{f}_{h}$ on the density appear explicitly in every term of Theorem (ref) as the ratio $\bar{f}_{h}/\underline{f}_{h}$. A distribution with a very flat density (small $\underline{f}_{h}$) near its quantile produces large estimation noise for a given calibration window --- intuitively, many observations fall near the boundary between covering and not covering, so small sample fluctuations matter a great deal. Conversely, a very peaked density (large $\bar{f}_{h}$) means small drifts in the distribution translate into larger shifts in the quantile, amplifying the drift-bias term. The condition is satisfied by all continuous distributions with bounded, bounded-away-from-zero densities in a neighbourhood of the relevant quantile, including normal, Student-$t$, and most parametric forecast error distributions.
The fourth and final condition concerns the gap between the oracle forecasts used in the analysis and the estimated forecasts used in practice. The theory is developed with oracle predictions because it isolates the calibration mechanism cleanly. The practitioner uses estimated predictions, and the difference introduces a fourth source of coverage error.
This condition is deliberately stated in a model-free way: the pair $(r_T, \eta_T)$ encapsulates the entire estimation cost of the forecasting model, regardless of its internal structure. For AR($p$) models, $r_T = O(T^{-1/2})$ and $\eta_T = 0$ under standard mixing conditions, so Term (IV) in the bound decays at rate $T^{-1/2}$ --- faster than the dominant $T^{-\beta/(2\beta+1)}$ rate and thus asymptotically negligible. For ARMA-GARCH, the same holds under the conditions of FryzlewiczSubbaRao2011. For machine-learning models, the pair must be supplied separately, but the three calibration terms in the bound remain unaffected: the window-length choice is determined by the first three terms alone, regardless of how accurately the model predicts.
With the four assumptions in place, the strategy for bounding the coverage error is to linearise the empirical quantile $\hat{q}^{o}_{T,m,h}$ around a deterministic target $q^{\circ}_{T,m,h}$. This is the content of the Bahadur representation: it says that the difference $\hat{q}^{o}_{T,m,h} - q^{\circ}_{T,m,h}$ is, to leading order, just a scaled version of the centred empirical CDF evaluated at $q^{\circ}_{T,m,h}$, plus a remainder that is negligible relative to $m^{-1/2}$.
Why does this linearisation matter? Because once we have it, bounding $\lvert \hat{q}^{o}_{T,m,h} - q^{\circ}_{T,m,h} \rvert$ reduces to bounding $\lvert \widehat{G}_{T,m,h}(q^{\circ}_{T,m,h}) - (1-\alpha) \rvert$, which is a sum of centred indicator random variables --- an object we know how to control using Rio's covariance inequality under mixing. Without the linearisation, quantile deviations are genuinely nonlinear and much harder to handle.
Establishing the Bahadur representation rigorously for a locally non-stationary, $\alpha$-mixing sequence is the primary technical step in the paper. The difficulty is that standard proofs assume stationarity; here the calibration scores $S^o_{T-1,h}, \ldots, S^o_{T-m,h}$ have different marginal distributions, so one cannot apply classical results directly. The proof, given in Appendix (ref), uses a stationary approximation to reduce to the distribution $F_{T,h}$, then controls the approximation error via Assumption (ref), and finally applies a bracketing-entropy maximal inequality to bound the remainder uniformly.
The condition $m/T \to 0$ is satisfied by the optimal window $m^{\star} \asymp T^{2\beta/(2\beta+1)}$, since $m^{\star}/T = T^{-1/(2\beta+1)} \to 0$. The proof uses a bracketing-entropy argument for the VC class of half-lines combined with the Merlev\`{e}de--Peligrad--Rio Bernstein inequality MerlvedePeligradRio2009; see Appendix (ref) for details.
The Bahadur representation converts a problem about quantile deviations into a problem about empirical CDFs, which in turn can be controlled by the four assumptions. The main theorem assembles this into a single closed-form bound that separates the four distinct sources of coverage error. Each source corresponds to a question a practitioner might naturally ask:
Terms (I) and (III) pull in opposite directions as $m$ changes: larger $m$ means less quantile noise but more drift bias. Their balance determines the optimal window, derived in the next subsection. The theorem makes this trade-off precise.
Theorem (ref) tells us how large the coverage error can be for any choice of $m$. The natural next question is: what value of $m$ makes the bound as small as possible? The answer is not $m = T$ (use all the data), nor $m = 1$ (use only the most recent observation), but a balance between the two --- and the location of that balance depends on $\beta$ and $T$ in a precise way.
Under the short-memory regime $A_{h}(\infty) := \sum_{k=1}^{\infty}\alpha_h(k) < \infty$, the Bahadur remainder $\bar{f}_{h} B_{m,h} = O(m^{-3/4}(\log m)^{3/4})$ is negligible relative to Term (II), which is $O(m^{-1/2})$. Terms (I) and (IV) are controlled by the forecasting model and do not depend on the window choice. Under H\"{o}lder-$\beta$ drift the window-dependent part of (ref) is
Minimising $R_h(m;\beta)$ over $m > 0$: setting $\partial R_h/\partial m = 0$ gives $\tfrac{1}{2}\Gamma_h m^{-3/2} = \beta L_{F,h} m^{\beta-1}/T^\beta$, hence $m^{\beta+1/2} = \Gamma_h T^\beta/(2\beta L_{F,h})$, yielding
Substituting back gives the optimised coverage-error rate
At $\beta = 1$ (Lipschitz drift) these reduce to $m^{\star} \asymp T^{2/3}$ and rate $O(T^{-1/3})$, the values previously reported in the literature. For $\beta > 1$ (smoother drift) both the window and the rate are larger; for $0 < \beta < 1$ (rougher drift) the window is smaller and the rate is slower. Section (ref) establishes that $T^{-\beta/(2\beta+1)}$ is not merely an upper bound but the minimax-optimal rate --- no conformal procedure can achieve better coverage accuracy on the class $\mathcal{F}(L,\beta)$.
The series-specific constant $C_h := (\Gamma_h/(L_{F,h}\beta))^{2/(2\beta+1)}$ governs the magnitude of $m^\star$: series with strong GARCH clustering have large $A_{h}(\infty)$ and hence large $\Gamma_h$, making $m^\star$ smaller; series with slow drift have small $L_{F,h}$, making $m^\star$ larger.
The optimal window result says that rolling-origin conformal prediction with $m^\star \asymp T^{2\beta/(2\beta+1)}$ achieves coverage error $O(T^{-\beta/(2\beta+1)})$. A natural and important question is whether this rate could be improved by a cleverer procedure --- one that, say, uses a variable window or a non-uniform weighting of the calibration scores. The answer is no.
The rate $T^{-\beta/(2\beta+1)}$ is a fundamental statistical limitation of the problem: it is the minimax rate for this class of processes, in the same sense that $n^{-1/2}$ is the minimax rate for estimating a mean from $n$ i.i.d.\ observations. No conformal procedure --- regardless of how it constructs its calibration set, weights its observations, or adapts its quantile estimate --- can achieve a uniformly better coverage-error rate on $\mathcal{F}(L,\beta)$. The rolling-origin scheme with Winkler-optimal $m^\star$ is therefore not just practically convenient; it is statistically efficient in the strongest possible sense.
The proof follows the classical two-point Le Cam strategy Tsybakov2009: construct two hypotheses in $\mathcal{F}(L,\beta)$ that are close enough to be nearly indistinguishable yet far enough apart that any procedure covering one must miscover the other. The signal strength $\Delta$ measures the separation; optimising $\Delta$ subject to the model-class and indistinguishability constraints gives the rate. The key technical ingredient that makes the construction work for general $\beta$ is using a H\"{o}lder-$\beta$ bump function to encode the perturbation, rather than the more common Lipschitz bump.
Theorem (ref) establishes an upper bound on the coverage deviation. We now show that the rate $T^{-\beta/(2\beta+1)}$ is sharp: no conformal procedure can achieve better accuracy uniformly over the class $\mathcal{F}(L,\beta)$.
Define the minimax coverage risk as \[ R_T^{\star}(\beta) := \inf_{C_T}\,\sup_{P \in \mathcal{F}(L,\beta)} \lvert \mathbb{P}_P(Y_{T+h} \in C_T) - (1-\alpha) \rvert, \] where $\mathcal{F}(L,\beta)$ is the class of all locally stationary processes satisfying Assumption (ref) with parameters $L$ and $\beta$.
When the mixing coefficients decay polynomially, $\alpha_h(k) \asymp k^{-a}$ with $0 < a \leq 1$, the process has long memory in the sense relevant to our problem: the cumulative dependence measure satisfies $A_h(m) \asymp m^{1-a}$, which grows with $m$ rather than remaining bounded. The dominant stochastic term in (ref) is then $(\bar{f}_{h}/\underline{f}_{h})\sqrt{A_h(m)/m} \asymp m^{-a/2}$ rather than $m^{-1/2}$, and the window-dependent bound becomes \[ R_h(m;\beta,a) \asymp \Gamma_h^{(a)}\,m^{-a/2} + L_{F,h}\left(\frac{m}{T}\right)^{\!\beta}, \qquad \Gamma_h^{(a)} := \frac{\bar{f}_{h}}{\underline{f}_{h}}\sqrt{C_{a}\,A_{h}(\infty)^{(a)}}, \] where $C_a > 0$ is a constant and $A_{h}(\infty)^{(a)}$ is the appropriate long-memory analogue of $A_h(\infty)$. Balancing $m^{-a/2}$ against $(m/T)^\beta$ gives $m^{a/2+\beta} \asymp T^\beta$, hence the optimal window
The formula (ref) applies for $0 < a < 1$, where $A_h(m) \asymp m^{1-a}$ grows with $m$. Within this regime, smaller $a$ corresponds to slower-decaying mixing ($\alpha_h(k) \asymp k^{-a}$ decays more slowly for smaller $a$), hence stronger long-range dependence and a larger optimal window. Conversely, as $a \nearrow 1$ from below, the long-memory effect weakens and $m^{\star} \asymp T^{2\beta/(2\beta+a)}$ approaches the boundary $T^{2\beta/(2\beta+1)}$. At $\beta=1$ (Lipschitz drift), the formula reduces to $T^{2/(2+a)}$, the standard long-memory optimal window.
One should not take $a \to \infty$ inside (ref): for $a > 1$ the mixing coefficients are summable ($\sum_k \alpha_h(k) < \infty$), $A_h(m)$ is bounded, and the long-memory formula no longer applies. In that regime the stochastic term reverts to $m^{-1/2}$ and the optimal window is the short-memory rule $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ from Section (ref). The polynomial mixing formula interpolates smoothly between slow mixing ($a$ near 0) and the boundary of summability ($a = 1$), but does not extend beyond it. At the boundary $a = 1$, $A_h(m) \asymp \log m$ and the dominant term is $\sqrt{\log(m)/m}$, giving $m^{\star} \asymp T^{2/3}(\log T)^{1/3}$, a logarithmic correction to the short-memory rule. Table (ref) summarises the window rules across regimes.
When the primary source of nonstationarity is time-varying scale rather than a drifting conditional distribution, normalising the conformal scores can substantially reduce the drift term (III) in (ref) while leaving terms (I) and (II) unaffected.
\paragraph{Setup.} Suppose the oracle forecast error factors as $Y_{t+h} - f^{\star}_{t,h} = \sigma_{t,h}\,\varepsilon_{t,h}$, where $\sigma_{t,h} > 0$ is a time-varying scale parameter and $\varepsilon_{t,h}$ is a mean-zero innovation. If the distribution of $\varepsilon_{t,h}$ is locally stationary with drift constant $L^{sc}_{F,h}$, then the oracle score $S^o_{t,h} = \sigma_{t,h}\lvert \varepsilon_{t,h} \rvert$ has drift constant $L_{F,h} \asymp \sigma_{t,h} L^{sc}_{F,h} + \dot{\sigma}_{t,h}$, which may be much larger than $L^{sc}_{F,h}$ when $\sigma_{t,h}$ is itself varying rapidly. The volatility-scaled oracle score \[ S^{sc,o}_{t,h} := \frac{S^o_{t,h}}{\sigma_{t,h}} = \lvert \varepsilon_{t,h} \rvert \] has drift constant $L^{sc}_{F,h}$ directly, removing the contribution of $\dot{\sigma}_{t,h}$ from the bias term. In practice one substitutes the estimated scale $\hat{\sigma}_{t,h}$ (from GARCH, stochastic volatility, or a realised-volatility estimator) to obtain $\hat{S}^{sc}_{t,h} := \hat{S}_{t,h}/\hat{\sigma}_{t,h}$. The scaling introduces an additional estimation error of order $\lvert \hat{\sigma}_{t,h} - \sigma_{t,h} \rvert/\sigma_{t,h}$, which enters through Assumption (ref) and is typically $o(1)$ for consistent volatility estimators.
The Bahadur representation (Proposition (ref)) is proved under $\alpha$-mixing. An alternative, and in some respects more natural, framework for nonlinear time series is the physical dependence framework of Wu2005PNAS, which characterises dependence through the stability of a functional representation $Y_{t,T} = G(t/T, \varepsilon_t, \varepsilon_{t-1}, \ldots)$ under perturbations of the innovation sequence.
Physical dependence is neither implied by nor implies $\alpha$-mixing in general, but holds for ARMA, GARCH, and many nonlinear time-series models Wu2005PNAS. For practitioners working with threshold autoregressive, bilinear, or neural-network models --- all of which are naturally described by functional representations but may not satisfy standard mixing conditions --- this framework is the more appropriate one.
The following proposition shows that the Bahadur representation of Proposition (ref) remains valid under this alternative framework, with the functional CLT rate $O(\sqrt{\log m/m})$ replacing the $O(m^{-1/2})$ rate of the i.i.d.\ case.
Theorems (ref) and (ref) together establish that the optimal window is $m^\star \asymp T^{2\beta/(2\beta+1)}$ and that no procedure can do better. But this is a population-level statement: it assumes knowledge of $L_{F,h}$, $\beta$, and $\Gamma_h$, none of which are observable. The practitioner must therefore select $m$ from the data alone.
The natural question is whether data-driven selection destroys the optimality guarantee. It does not. After the calibration step, the candidate intervals $\{\hat{q}_m : m \in \mathcal{M}_T\}$ are fixed functions of the observed calibration data, and the validation fold estimates their future predictive performance given those realised intervals. Theorem (ref) of Appendix (ref) shows that the Winkler criterion concentrates uniformly around this conditional risk at a rate faster than $R_T^\star(\beta)$, so the data-driven selector $\hat{m}$ performs as well as if $m^\star$ were known.
Let $\mathcal{M}_T = \{m_1, \ldots, m_K\}$ be a finite candidate grid with $K = K_T$ satisfying $\log K_T = o(T^{1/(2\beta+1)})$ (e.g.\ the 30-point grid in $[0.1, 4.0]\times T^{2\beta/(2\beta+1)}$ used in Section (ref)). Let $\mathcal{F}_{\mathrm{cal}}$ be the $\sigma$-field generated by all data used in calibration, define the conditional validation risk $\mathcal{R}_T(m) := \mathbb{E}[\mathcal{W}_T(m) \mid \mathcal{F}_{\mathrm{cal}}]$, and the Winkler cross-validation criterion \[ \mathcal{W}_T(m) := \frac{1}{n_{\mathrm{val}}} \sum_{t \in \mathcal{T}_{\mathrm{val}}} W_\alpha\!\bigl(Y_{t+h},\, \hat{Y}_{t+h|t} - \hat{q}_{m,\alpha},\, \hat{Y}_{t+h|t} + \hat{q}_{m,\alpha}\bigr), \] where $\mathcal{T}_{\mathrm{val}}$ is a held-out validation fold of size $n_{\mathrm{val}} \asymp T$ and $W_\alpha(y, l, u) := (u-l) + (2/\alpha)[\max(l-y,0) + \max(y-u,0)]$ is the Winkler score Winkler1972. Define the data-driven window selector $\hat{m} := \operatorname{arg\,min}_{m \in \mathcal{M}_T}\mathcal{W}_T(m)$.
We evaluate rolling-origin conformal prediction on two complementary datasets. The real-data application uses six series representing distinct forms of nonstationarity: US CPI log-inflation and unemployment (FRED monthly, $T \approx 709$); S&P 500 and VIX daily log-returns ($T \approx 6200$--$6400$); daily WIG20 log-returns ($T \approx 6400$); and European electricity load ($T \approx 9000$). The M4 scaling analysis uses a stratified sample of 93 series from the M4 competition benchmark MakridakisEtal2020 across all five frequency classes (Yearly, Quarterly, Monthly, Weekly, Daily), chosen to span a T range of approximately 20 to 9000 and allow a regression test of the $T^{2/3}$ prediction.
Three calibration schemes are compared throughout. Full history uses all available pseudo-out-of-sample errors; this is the natural baseline reflecting standard practice. Rolling uses the Winkler-optimal $m^{\star}$ chosen from a 30-point grid in $[0.10, 4.0] \times T^{2/3}$ (using $\beta=1$ as the empirical benchmark). Volatility-scaled rolling normalises scores by GARCH conditional volatility before computing the calibration quantile, implementing Corollary (ref). Adaptive conformal inference GibbsCandes2021 was evaluated but excluded from reported results: the fixed-step update rule $q_{t+1} = q_t + \gamma(\alpha - \mathbf{1}[\text{miss}])$ exhibited persistent bimodal coverage on fat-tailed series (coverage $\to 1$ or $\to 0$ depending on the window), a known structural limitation under heavy tails and regime shifts that cannot be resolved by adjusting $\gamma$.
Forecasts are produced by linear AR($p$) with lag selected by BIC and ARMA(1,1)--GARCH(1,1). Intervals are evaluated at horizons $h \in \{1, 5, 22\}$ using empirical coverage, mean interval half-width, rolling local coverage over 50-observation windows, and the Winkler interval score Winkler1972.
Rolling-origin calibration outperforms full-history calibration in 31 of 36 comparisons (86%) by Winkler score. Among the 31 wins, the median improvement is 12.3% (range 0.5%--16%). The five cases in which full-history wins are all at the shortest horizon $h = 1$ for macro series, where the rolling window buys little adaptivity and the quantile-noise cost is comparatively high --- consistent with Remark (ref) below.
Table (ref) reports the dataset-level results at $h=1$ that aggregate to these comparisons. Rolling-origin calibration achieves the lowest Winkler score on every dataset; the gap is most pronounced for financial series (S&P 500/VIX/WIG20), where the Winkler-optimal window collapses to the grid boundary $m = 17$, and is smallest for the macro and electricity series, where the optimal window is large relative to $T^{2/3}$ and the rolling/full-history schemes nearly coincide.
Rolling coverage tracks the 90% target with high precision at short and medium horizons (Figures (ref) and (ref)): at $h = 1$ and $h = 5$, every series lies within $\pm 2\%$ of the nominal level (mean absolute deviation 0.008 and 0.012, respectively). At $h = 22$ performance is more mixed: 50% of series lie within $\pm 2\%$ and 67% within $\pm 5\%$, with the unemployment series at $h = 22$ providing the main failure case (coverage 0.81 for AR, 0.82 for ARMA-GARCH, versus the nominal 0.90). At this horizon the 22-step-ahead error distribution is dominated by model misspecification --- unemployment is highly persistent and its long-horizon error distribution is wide and non-stationary --- and the rolling quantile cannot adapt quickly enough. This illustrates a boundary of Assumption (ref): the local-stationarity condition requires the score distribution to drift at rate $O(T^{-1})$, which fails when structural shifts accumulate faster than the sample grows.
For the three financial series (S&P 500, VIX, WIG20), the Winkler-optimal window falls below the lower boundary of the 30-point evaluation grid ($m < 17$, ratio $< 0.05$), meaning the Winkler score continues to decline at the smallest evaluated $m$. This is not a numerical artefact: GARCH volatility clustering concentrates information about the current error scale in the very most recent observations, making very short windows ($m \approx 5$--$15$) theoretically preferable. The ratio $C_h = (\Gamma_h/L_{F,h})^{2/3}$ is small for these series because $\Gamma_h$ --- which grows with the mixing constant $A_h(\infty)$ --- is large relative to the drift rate $L_{F,h}$: strong dependence penalises large windows heavily while the score distribution drifts slowly. These series are excluded from the $m^{\star}$ ratio analysis in Table (ref) as the grid boundary prevents identification of the true optimum.
Figure (ref) shows empirical coverage as a function of calibration window $m$ for each dataset. At the short horizon $h=1$, coverage is near 90% across all window lengths for the macro and WIG20 series, while the financial and electricity series show a mild hump --- slightly overcovering at short windows and converging to the target as $m$ grows. The dashed vertical line marks $T^{2/3}$: in every panel the Winkler minimum (Figure (ref)) lies close to this benchmark, consistent with the $m^\star \asymp T^{2/3}$ prediction at $\beta=1$.
Volatility-scaled and plain rolling scores produce identical results for the AR model on all financial series, as the AR model does not produce GARCH-conditional volatility forecasts and the scaled scheme falls back to the rolling quantile. For the ARMA-GARCH model on financial series, volatility scaling produces marginally narrower intervals with equivalent coverage, consistent with Corollary (ref). The electricity ARMA-GARCH volatility-scaled results are excluded because GARCH conditional volatility at the scale of $\sim 70{,}000$ MW causes numerical overflow in the multi-step variance forecast.
Figure (ref) summarises coverage at the Winkler-optimal $m^\star$ for each scheme and dataset. Rolling and volatility-scaled rolling are nearly identical across all four datasets at $h=1$, lying within 1 percentage point of the 90% target. Full history is slightly closer to target for the electricity and WIG20 series at this horizon, consistent with the theoretical prediction that full history outperforms rolling when drift is slow relative to calibration-window noise.
Figure (ref) shows rolling local coverage (mean $\pm$ one standard deviation computed over 50-observation windows). The error bars quantify conditional stability: a wide bar indicates that coverage fluctuates substantially over time even if marginal coverage is near 90%. Rolling and volatility-scaled rolling show narrower error bars than full history for the financial and electricity datasets, reflecting the method's adaptivity to local changes in the error distribution.
We test the prediction $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ empirically. Under $\beta = 1$ (Lipschitz drift) this specialises to $m^{\star} \asymp T^{2/3}$, which we use as the benchmark since $\beta$ is not separately identified from the empirical window scaling. We regress $\log m^{\star}$ on $\log T$ across 93 M4 series drawn from all five frequency classes (Yearly through Daily), spanning $T \approx 20$ to $9{,}000$ and $T^{2/3} \approx 8$ to $440$. The Winkler-optimal $m^{\star}$ is identified from a 30-point grid in $[0.10, 4.0] \times T^{2/3}$, yielding a continuous ratio distribution with only 7 of 93 series at a grid boundary (compared to 60 of 62 under the original 7-point grid). The log-$T$, log-$m^{\star}$ correlation across all 93 series is $r = 0.925$, confirming that sample size is the dominant predictor of the optimal window.
The pooled regression $\log m^\star_i = \beta_0 + \beta_1 \log T_i + \varepsilon_i$ on the $n = 93$ series gives $\hat\beta_1 = 0.776$ (heteroskedasticity-robust HC1 SE $= 0.033$, $95\%$ CI $[0.710, 0.841]$, $R^2 = 0.855$), which does not contain $2/3$. This specification is misspecified because $\log m^\star$ also depends on the series-specific constant $\log C_h$, which is correlated with $\log T$ across the M4 panel: Daily and Hourly series tend to have both large $T$ and small $C_h$ (strong GARCH clustering inflates $A_{h}(\infty)$ and shrinks $m^\star$), while Yearly and Quarterly series have small $T$ and large $C_h$. This induces a positive correlation between $\log T_i$ and the omitted regressor $\log C_{h,i}$, biasing $\hat\beta_1$ upward.
Adding frequency fixed effects $\delta_{f(i)}$ for the five M4 classes absorbs the cross-frequency mean of $\log C_h$ (using frequency as a proxy for the structural regime that determines $\Gamma_h/L_{F,h}$). The augmented specification $\log m^\star_i = \beta_0 + \beta_1 \log T_i + \delta_{f(i)} + \varepsilon_i$ gives $\hat\beta_1 = 0.614$ (HC1 SE $= 0.097$, $95\%$ CI $[0.424, 0.805]$, within-$R^2 = 0.861$, $n = 93$, $\mathrm{df} = 87$), which contains $2/3$. Clustering standard errors at the frequency level (5 clusters) widens the CI to $[0.388, 0.840]$ and does not overturn the conclusion. The bias direction (0.776 $\to$ 0.614 after controlling for $C_h$) matches the sign predicted by the omitted-variable formula.
Table (ref) reports the full set of regression specifications: the pooled OLS, the FE-augmented regression, and the within-frequency slopes. The Daily class — the largest by $T$ and the most homogeneous in dependence structure — yields a within-frequency slope of $0.670$ that contains $2/3$ at the centre of its $95\%$ CI. Within-frequency slopes for smaller classes (Yearly, Quarterly, Weekly) are noisier owing to small sample sizes ($n \leq 20$ per class) and limited $T$ range within class.
Table (ref) reports the implied series-specific constants $\hat{C}_f = \exp(\hat{a}_f)$. The wide spread of $\hat{C}_f$ across frequencies (from $0.05$ for Quarterly to $34.14$ for Weekly) confirms that $C_h$ varies substantially across the M4 panel, justifying the FE specification used to recover the structural slope $\beta_1 = 2/3$.
Figure (ref) (left panel) shows the log-log scatter of $m^\star$ against $T^{2/3}$ for all 93 series, coloured by frequency class. The clustering by frequency is the visual signature of the series-specific constant $C_h$: Daily series (purple) sit systematically below the $m = T^{2/3}$ reference line, reflecting small $C_h$ from GARCH clustering, while Yearly series (red) sit above it, reflecting larger $C_h$ from slower drift. The right panel is a forest plot of the slope estimates with 95% confidence intervals; the pooled-FE estimate and the Daily within-frequency estimate both span the theoretical $2/3$ value.
The within-frequency slopes are heterogeneous, ranging from 0.237 (Weekly) to 1.483 (Quarterly), with only the Daily class yielding a 95% CI $[0.422, 0.917]$ that individually contains $2/3$. This within-class heterogeneity reflects genuine variation in $C_h$ across series: the $T^{2\beta/(2\beta+1)}$ rule is a cross-frequency average regularity, and the series-specific constant encodes the ratio $\Gamma_h/L_{F,h}$ from (ref) directly. Mean coverage across all 93 series ranges from 0.900 (Daily) to 0.915 (Quarterly), with no frequency class showing systematic deviation beyond 1.5 percentage points.
Figure (ref) shows coverage and mean interval half-width as a function of forecast horizon $h \in \{1, 5, 22\}$, averaged across all datasets. Coverage remains close to 90% at $h=1$ and $h=5$ for both schemes, with rolling tracking the target more closely than full history at $h=22$. The right panel shows that rolling intervals are narrower than full-history intervals at longer horizons, confirming that the Winkler score gain documented in the real-data section reflects genuine interval efficiency, not overcoverage.
Table (ref) confirms the visual pattern numerically: rolling-origin calibration delivers Winkler-score improvements of $14\%$, $13\%$, and $15\%$ at $h \in \{1, 5, 22\}$ respectively, with coverage remaining within $1\%$ of nominal at the shorter horizons. The $h = 22$ row shows the same coverage degradation discussed above — rolling at $0.871$ vs nominal $0.90$ — which the analysis attributes to model misspecification at long horizons rather than calibration-window choice.
We have proposed and analysed rolling-origin conformal prediction for time-series forecasting, establishing its theoretical properties under local nonstationarity and weak dependence. The method constructs prediction intervals by calibrating against the $m$ most recent pseudo-out-of-sample forecast errors, adapting automatically to the serial dependence, volatility clustering, and distributional drift that invalidate classical conformal guarantees. The exact finite-sample coverage guarantee of classical conformal prediction rests irreducibly on exchangeability and cannot hold in the time-series setting; the contribution of this paper is to establish what replaces it: a precise characterisation of the coverage deviation as a function of the window length $m$, the sample size $T$, and the structure of the dependence and drift, together with the optimal window rule that minimises the deviation.
The main result (Theorem (ref)) bounds the coverage deviation $\lvert \mathbb{P}(Y_{T+h} \in \widehat{C}_{T+h|T}(1-\alpha)) - (1-\alpha) \rvert$ by a sum of four interpretable terms, under the general H\"{o}lder-$\beta$ drift model. The optimal calibration window $m^{\star} \asymp T^{2\beta/(2\beta+1)}$ yields coverage-error rate $O(T^{-\beta/(2\beta+1)})$; at $\beta = 1$ these reduce to $T^{2/3}$ and $O(T^{-1/3})$. Theorem (ref) establishes that this rate is minimax-optimal: no conformal procedure can achieve better accuracy uniformly over $\mathcal{F}(L,\beta)$. The Bahadur representation is established under both $\alpha$-mixing (Proposition (ref)) and physical dependence (Proposition (ref)), covering ARMA, GARCH, and a broad class of nonlinear time-series models. Theorem (ref) and Appendix (ref) together provide an oracle inequality for the implemented Winkler cross-validation window selector: conditional on the calibration sample, $\hat{m}$ achieves near-minimal validation loss, and in expectation over calibration samples its average loss is within $o(T^{-\beta/(2\beta+1)})$ of the best candidate window.
The empirical analysis on six real series and 93 M4 competition series confirms the theory's main predictions. Rolling-origin calibration outperforms full-history calibration in 86% of comparisons (median Winkler improvement 12.3%), achieves coverage within $\pm 2\%$ of the 90% target on every series at horizons $h \in \{1, 5\}$, and the cross-frequency regression recovers slope $0.614$ (95% CI $[0.424, 0.805]$) consistent with the theoretical $2\beta/(2\beta+1)$ at $\beta = 1$.
One direction remains open. Multivariate extensions to HPD prediction regions for vector outputs require reworking of the oracle-quantile theory in higher dimensions, where the Bahadur representation and the density-bound argument of Theorem (ref) do not directly apply.
All code and data-fetching scripts required to reproduce the empirical results are publicly available at \url{https://github.com/profsms/conform-ROE-repl}. Running the script 01_fetch_data.jl downloads and caches every dataset used in the paper; subsequent analysis scripts are deterministic given a fixed random seed.
Macroeconomic series (CPI inflation, unemployment) are retrieved from the Federal Reserve Economic Data (FRED) database at \url{https://fred.stlouisfed.org} via the FRED REST API (a free API key, exported as the FRED_API_KEY environment variable, is required; the script falls back to a synthetic series if the key is absent). Financial return series (S&P 500, VIX, WIG20) are fetched from the Yahoo Finance API. The electricity load series is the ENTSO-E hourly load dataset. The M4 competition series are sampled from the M4 competition repository at \url{https://github.com/Mcompetitions/M4-methods}.
The author is grateful to an anonymous colleague whose critical reading of an earlier draft identified an overclaim in Theorem (ref) and prompted the formal reduction lemma in the proof of Theorem (ref).
The AI model Sonnet 4.7 and Opus 4.7 were used during the preparation of Appendix (ref) for exploratory generation of proof-strategy ideas and feedback on the conditional-versus-unconditional formulation of the uniform concentration result. All mathematical content, proofs, assumptions, and final formulations were developed, verified, and written by the author, who takes full responsibility for the correctness of the manuscript.
The author received no specific funding for this work.
The authors declare no conflict of interest.