1 Introduction↩︎

Prediction intervals for time-series forecasts are central to applications such as retail and finance. Conformal prediction (CP) provides distribution-free coverage under exchangeability [1], [2]. Extensions address temporal dependence [3]–[5], online adaptation [6], [7], weighted quantiles [8], [9], time-series datasets [10], [11], and distribution shifts [12].

Many applications require uncertainty for aggregated targets such as annual sums and year-over-year growth rates. This is challenging because aggregation induces dependence and nonlinear transformations, so combining pointwise intervals is invalid, and CP guarantees do not directly extend to multi-step dependent settings.

We propose SA-MSCP, a simulation-enhanced multi-step split conformal method [5] for aggregated forecasting. The method uses expanding-window cross-validation to collect residuals across forecast origins, then simulates future paths with a block bootstrap over residual sequences to better preserve dependence before constructing prediction intervals from empirical quantiles of aggregated trajectories.

Related work includes copula-based calibration for joint coverage [11], conformalized interval arithmetic for sums under group exchangeability [13], and hierarchical reconciliation approaches [14]. In contrast, SA-MSCP focuses on temporal aggregation within a univariate series.

Our contributions are: (i) a practical MSCP-based procedure for aggregated targets that improves empirical coverage without finite-sample guarantees, and (ii) an empirical evaluation on M4 and a proprietary dataset of 2,000 series showing consistent coverage gains over a simulated-path baseline with significance assessed via Wilcoxon tests.

2 Methodology↩︎

Split conformal prediction (SCP) constructs prediction intervals from calibration residuals and achieves marginal coverage under exchangeability [2], [15]. In time series, temporal dependence violates this assumption, motivating multi-step conformal approaches that approximate residual distributions via rolling or expanding windows.

SA-MSCP follows this framework. Residuals \(\hat{\varepsilon}_{k,h} = y_{k+h} - \hat{y}_{k+h}\) are collected via expanding-window cross-validation across forecast origins \(k\) and horizons \(h\). In our implementation, the cross-validation uses an expanding window with an initial calibration window of \(10\) observations; at each step the window grows by one observation and the next \(h\) points are used for validation. We center each horizon-wise residual column, extract all valid consecutive residual blocks of length \(b\), and simulate \(S\) future paths by resampling these blocks with replacement and stitching them together along the horizon. This block bootstrap preserves local residual dependence better than horizon-wise independent resampling.

In this paper, we apply the method to monthly time series with annual aggregates, but the same approach can be used with other temporal resolutions and aggregation units. Each simulated path yields monthly forecasts, which are aggregated to annual totals \(\widehat{Y}_{s} = \sum_{m=1}^{12} \widehat{y}_{s,m}\) and transformed into growth rates using consecutive annual values. Prediction intervals are then obtained from empirical quantiles of the simulated samples at the desired coverage levels for both pointwise and aggregated targets.

We use block size \(b=12\) for M4, matching its 36-month forecast horizon and preserving annual dependence patterns, and \(b=3\) for the proprietary data. The horizon is only 12 months for the proprietary data, so using \(b=12\) would leave too few distinct bootstrap combinations; \(b=3\) instead preserves within-quarter dependence while retaining sufficient resampling diversity. This simulation-based calibration approximates the distribution of aggregated quantities under temporal dependence, enabling uncertainty quantification beyond pointwise horizons.

We compute empirical coverage rates and interval widths at multiple confidence levels (e.g., 90%, 95%, 99%) across series and forecast scenarios (original, aggregated, growthrate). We assess statistical significance of differences versus a simulated-path baseline using Wilcoxon signed-rank tests.

3 Experiments↩︎

We fit Auto-ARIMA models using the ARIMA implementation in the fable package [16]. Future trajectories are simulated using a modified version of the package’s generate() function with block-bootstrap innovations drawn from cross-validated residual blocks. For each time series, we generate \(S = 10{,}000\) simulated paths.

We evaluate SA-MSCP on M4 monthly sales data for aggregated annual totals (Aggregated Sales) and year-over-year growth rates (Y-o-Y Sales Growth), which are the primary targets of interest, and also report original monthly forecasts (Raw Sales) for reference. We compare against a baseline simulated-path approach (Baseline) [17] that uses conditional simulation without cross-validated residual calibration.

A direct split-conformal baseline built on aggregated nonconformity scores is less suitable in our setting because we start from monthly series and evaluate annual aggregates: once the data are aggregated to the yearly level, only a small number of calibration points remain for each series, making direct SCP at the aggregated level statistically inefficient and unstable.

Across all targets and levels, SA-MSCP delivers higher coverage than Baseline with a corresponding increase in interval width, but both methods miss the target coverage level. The coverage and width values are presented in Tables 1 and 2, respectively. In absolute terms, coverage gains range from roughly 6.9 to 17.5 percentage points (Column Coverage Delta in Table 3). Relative to Baseline, SA-MSCP improves coverage while widening intervals. The coverage cost summarizes the relative width increase per coverage gain: costs are lowest at the 10% level, increase at 5%, and are highest at 1% (Column Coverage Cost in Table 3). Raw Sales exhibits the lowest cost across levels, whereas Y-o-Y Sales Growth is the most costly—especially at 1%. These wider intervals are expected under compounded uncertainty from temporal aggregation.

3.0.0.1 Real-world monthly sales.

We additionally evaluate on 2,000 actual monthly sales series. As with M4, SA-MSCP attains higher coverage with wider intervals at higher confidence. For Raw Sales, SA-MSCP achieves 92.7%, 94.5%, and 96.1% coverage at the 10%, 5%, and 1% levels, respectively. This performance is better than on M4, but it still undercovers except at the 10% miscoverage level. Aggregated Sales and Y-o-Y Sales Growth exhibit slightly lower coverages but the same qualitative pattern. In Table 4, Aggregated Sales and Y-o-Y Sales Growth have identical coverage values because the proprietary test horizon is only 12 months, so annual aggregation produces a single forecast-year total. The corresponding growth rate is then obtained by dividing each simulated annual sum by the same observed annual sum from the previous year. This is a monotone transformation, so if an interval for aggregated sales covers the annual sum, then the corresponding growth-rate interval also covers the associated growth rate; conversely, if the aggregated-sales interval misses the annual sum, then the growth-rate interval also misses the associated growth rate.

Relative to Baseline, SA-MSCP improves coverage by 4.4–14.0 percentage points (Column Coverage Delta in Table 6). The coverage cost is smallest for Aggregated Sales at 10% and increases as nominal levels tighten; the largest cost occurs at the 1% level for Y-o-Y Sales Growth, which is lower in magnitude than on M4 but shows the same qualitative trend.

Overall, SA-MSCP achieves materially higher coverage than Baseline across both M4 and proprietary real-world sales data, at the expense of wider intervals. On M4, the width increase per unit coverage gain is smallest for Aggregated at 10% and largest for Growthrate at 1%. On the proprietary dataset, the corresponding extremes are likewise smallest for Aggregated at 10% and largest at 1%, consistent with aggregation effects observed in both settings.

In the absence of exchangeability, conformal prediction methods do not enjoy finite-sample coverage guarantees, and nominal levels should therefore be interpreted as target rather than guaranteed coverage. Consistent with this theoretical limitation, our empirical results show that while SA-MSCP substantially improves coverage relative to simulated-path baselines, achieved coverage remains below nominal levels as shown in Tables 1 and 4. This behavior reflects the intrinsic difficulty of uncertainty quantification in non-stationary and dependent time-series settings, rather than a failure of the proposed method.

In practice, users requiring closer alignment between nominal and achieved coverage may apply post-hoc calibration strategies, such as increasing the nominal coverage level, or applying online adaptation schemes introduced in [6]. Investigating principled calibration schemes that account for temporal dependence is an important direction for future work.

3.0.0.2 Significance testing.

Wilcoxon signed-rank tests comparing SA-MSCP vs Baseline yield \(p\)-values \(\ll 0.001\) for coverage and width across Aggregated and Growthrate targets at 10%, 5%, and 1% levels. Differences are statistically significant.

4 Conclusion↩︎

We propose SA-MSCP, a simulation-enhanced multi-step split conformal method that combines expanding-window cross-validation, block-bootstrap future-path simulation, and calibration on aggregated trajectories. SA-MSCP improves empirical coverage for monthly, aggregated annual, and growth-rate targets, at the cost of wider intervals, reflecting a fundamental trade-off in conformal prediction under dependence and aggregation. Residuals across forecast origins and horizons are not strictly exchangeable; expanding-window cross-validation partially accounts for temporal dependence by pooling residuals over time, while the block bootstrap preserves local cross-horizon dependence in simulated paths. Sequential CP methods (e.g., SPCI, PID) target pointwise horizons and long-run coverage, making extension to aggregated quantities nontrivial. Copula-based calibration could model cross-horizon dependence for sums and growth rates, but requires multiple time series, whereas our setting focuses on a single series [11].

Overall, SA-MSCP provides a practical and scalable approach to uncertainty quantification for aggregated forecasting tasks.

References↩︎