The Prediction Interval Coverage Audit: Published Bayesian Prediction Intervals Exhibit Systematic Undercoverage in Time Series Forecasting

Spike and Tyke

Abstract. Bayesian prediction intervals for time series forecasting carry an implicit promise: a nominal 95% interval should contain the realized value 95% of the time. We audited 120 published forecasting papers that report Bayesian prediction intervals, recomputing empirical coverage on held-out data using original code and data where available ($n=47$) and calibrated simulation otherwise ($n=73$). The median actual coverage for nominal 95% intervals was 81.2% (IQR: 74.6%--87.3%), representing a 13.8 percentage point gap. Undercoverage worsened with forecast horizon, dropping 2.3 percentage points per doubling of horizon length. ARIMA-family models exhibited less undercoverage than neural forecasters (median coverage 86.1% vs 74.4%). Decomposition analysis identified posterior predictive distributions that are too narrow---driven primarily by model misspecification rather than computational approximation error---as the dominant mechanism. These findings suggest that published Bayesian prediction intervals should be interpreted as approximate 80% intervals rather than 95% intervals, and that recalibration procedures should become standard practice.

1. Introduction

A 95% Bayesian prediction interval communicates that the forecaster assigns 95% posterior predictive probability to the realized outcome falling within the stated bounds. When many such intervals are collected, the fraction containing the true value---the empirical coverage---should converge to the nominal level. Departures from this guarantee indicate either computational error in posterior inference or, more fundamentally, that the model is misspecified in ways that compress the predictive distribution.

1.1 Motivation

The trustworthiness of prediction intervals underpins downstream decisions in energy markets [1], epidemiology [2], and supply chain management [3]. If published intervals systematically undercover, decision-makers who rely on them face unquantified risk. Despite the centrality of coverage to interval quality, no large-scale audit of published Bayesian prediction intervals exists. Simulation studies have examined coverage under known data-generating processes [4], but the gap between simulation conditions and applied practice remains unmeasured.

1.2 Scope and Contributions

This audit targets Bayesian prediction intervals specifically, excluding frequentist confidence-based intervals (which have separate coverage semantics) and point forecasts with ad hoc uncertainty bands. We examine 120 papers published between 2015 and 2025, spanning macroeconomic, energy, epidemiological, and financial forecasting domains. The analysis quantifies (a) the magnitude of undercoverage, (b) its dependence on forecast horizon, model class, and domain, and (c) the relative contribution of model misspecification versus computational approximation to the observed gap.

2. Related Work

Gneiting and Raftery [4] established the framework for evaluating probabilistic forecasts through proper scoring rules, emphasizing calibration and sharpness. Their work provides the theoretical foundation: a well-calibrated forecast system produces prediction intervals whose empirical coverage matches nominal levels across all quantile levels, not merely at the 95th percentile.

Hyndman and Athanasopoulos [5] documented that prediction intervals from exponential smoothing and ARIMA models in the forecast R package undercover at long horizons due to accumulated parameter estimation uncertainty that the plug-in approach ignores. Their proposed bootstrap correction improves coverage but remains imperfect because it cannot address structural misspecification.

Kuleshov et al. [6] demonstrated that neural network uncertainty estimates are poorly calibrated and proposed post-hoc recalibration via isotonic regression on a held-out calibration set. Their method, while effective, treats the model as a black box and does not diagnose the source of miscalibration.

Bracher et al. [7] analyzed prediction interval coverage in the context of COVID-19 forecasting hubs, finding that ensemble models achieved better coverage (89%) than individual models (median 78%) at nominal 95% levels. This domain-specific finding aligns with our broader cross-domain results.

Talts et al. [8] developed simulation-based calibration (SBC) as a diagnostic for Bayesian computation, checking that posterior samples are consistent with the prior predictive distribution. SBC detects computational errors but does not address model misspecification, which our analysis identifies as the primary driver of undercoverage.

3. Methodology

3.1 Paper Selection and Data Collection

We searched Google Scholar, Scopus, and arXiv for papers published 2015--2025 containing the terms "Bayesian prediction interval" or "posterior predictive interval" combined with "time series" or "forecasting." From an initial pool of 843 papers, we applied inclusion criteria: (i) the paper reports at least one nominal 95% prediction interval on a time series forecasting task, (ii) either the code/data are publicly available or the methodology is described in sufficient detail for simulation replication, (iii) the forecasting target is a continuous variable. This yielded 120 papers.

For each paper, we recorded the model class, forecast horizon(s), domain, sample size, whether the posterior was obtained by MCMC, variational inference, or other approximation, and the reported prediction intervals.

3.2 Coverage Computation

For the 47 papers with available code and data, we performed direct replication. We held out the final 20% of each time series (or used the authors' stated test set) and computed empirical coverage:

$$\hat{C} = \frac{1}{T_{\text{test}}} \sum_{t=1}^{T_{\text{test}}} \mathbf{1}\left[ y_t \in \left( \hat{q}$${0.025}^{(t)}, \hat{q}{0.975}^{(t)} \right) \right]

where $\hat{q}_{\alpha}^{(t)}$ denotes the $\alpha$-quantile of the posterior predictive distribution at time $t$.

For the 73 papers without available code, we employed simulation calibration. We fit the reported model class to a synthetic time series generated from a data-generating process (DGP) calibrated to match the reported summary statistics (sample size, autocorrelation structure, noise level). Coverage was computed on held-out synthetic data. To validate this simulation approach, we applied it to the 47 papers with direct replication and confirmed that simulation-based coverage estimates fell within 3.1 percentage points of direct estimates (mean absolute difference: 1.7pp).

3.3 Horizon Dependence Model

To quantify how coverage degrades with forecast horizon $h$, we fit a mixed-effects logistic regression:

$$\text{logit}(\hat{C}$${i,h}) = \beta_0 + \beta_1 \log_2(h) + \beta_2 X_i + u_i + \varepsilon{i,h}

where $i$ indexes papers, $h$ indexes horizon, $X_i$ is a vector of paper-level covariates (model class, domain, sample size), $u_i \sim N(0, \sigma_u^2)$ is a random intercept, and $\varepsilon_{i,h}$ is residual error. The coefficient $\beta_1$ estimates the change in log-odds of coverage per doubling of horizon.

3.4 Decomposition of Undercoverage Sources

We decomposed undercoverage into three sources: (a) computational approximation error, (b) parameter estimation uncertainty ignored by plug-in methods, and (c) model misspecification.

For models with MCMC posteriors, computational approximation error was assessed via split-$\hat{R}$ diagnostics [9]. Parameter estimation uncertainty was evaluated by comparing plug-in prediction intervals (fixing parameters at posterior means) against full posterior predictive intervals. The residual gap---between full posterior predictive coverage and nominal coverage---was attributed to model misspecification:

$$\Delta_{\text{misspec}} = (1 - \hat{C}$${\text{full}}) - \Delta{\text{comp}} - \Delta_{\text{param}}

where $\Delta_{\text{comp}}$ is the coverage loss attributable to MCMC error and $\Delta_{\text{param}}$ is the gain from integrating over parameter uncertainty relative to plug-in.

3.5 Model Class Comparison

We categorized models into four classes: (i) ARIMA/SARIMA (n=38), (ii) state-space/structural time series (n=29), (iii) Gaussian processes (n=18), and (iv) neural forecasters including BNNs, DeepAR, and transformer-based models (n=35). Coverage was compared across classes using a Kruskal-Wallis test followed by Dunn's post-hoc comparisons with Bonferroni correction.

4. Results

4.1 Overall Coverage Gap

Nominal 95% prediction intervals achieved a median empirical coverage of 81.2% (IQR: 74.6%--87.3%) across all 120 papers. Only 14 papers (11.7%) achieved coverage above 93%, while 31 papers (25.8%) fell below 75%. The coverage distribution was left-skewed, with a long tail of severely undercovering models.

A one-sample Wilcoxon signed-rank test confirmed that the median coverage was significantly below 95% ($V = 142$, $p < 10^{-15}$). The 95% confidence interval for median coverage was [79.4%, 83.1%].

4.2 Horizon Dependence

Coverage decreased monotonically with forecast horizon. The mixed-effects logistic regression yielded $\hat{\beta}_1 = -0.134$ (SE = 0.018, $p < 10^{-10}$), corresponding to a 2.3 percentage point decrease in coverage per doubling of horizon on the probability scale (evaluated at mean coverage).

The degradation was not linear in horizon but approximately linear in $\log_2(h)$, consistent with the logistic model specification. At one-step-ahead, median coverage was 88.4%---still substantially below 95%---indicating that baseline miscalibration exists even before horizon effects accumulate.

4.3 Model Class Differences

The Kruskal-Wallis test rejected the null of equal coverage distributions across model classes ($\chi^2 = 34.7$, $df = 3$, $p < 10^{-7}$). Post-hoc Dunn tests revealed that the primary contrast was between ARIMA-family models and neural forecasters.

ARIMA models' superior coverage is partly explained by their well-understood distributional assumptions: when the true DGP is approximately linear and Gaussian, the ARIMA prediction interval formula is analytically derived and asymptotically correct. Neural forecasters, by contrast, rely on approximate Bayesian inference (variational dropout, MC dropout, or deep ensembles interpreted as approximate posteriors) whose coverage properties are not analytically tractable.

4.4 Decomposition of Undercoverage

Among the 47 directly replicated papers, we decomposed the average 14.2pp undercoverage gap into three components:

Computational approximation error (MCMC convergence issues, variational approximation bias) accounted for 2.1pp (15% of the gap). Models using variational inference showed larger computational contributions (3.4pp) than MCMC-based models (0.9pp).

Ignored parameter estimation uncertainty accounted for 3.3pp (23% of the gap). This was most pronounced in papers that used plug-in posterior means rather than full posterior integration for prediction.

Model misspecification accounted for 8.8pp (62% of the gap). The dominant form of misspecification was distributional: the assumed noise model (typically Gaussian) was too thin-tailed relative to the observed residual distribution. In 71% of undercovering models, the residual kurtosis exceeded the assumed distribution's kurtosis by a factor of 1.5 or more.

The misspecification contribution $\Delta_{\text{misspec}}$ can be further decomposed by noting that under the true DGP $P^*$, the Bayesian predictive distribution $P_M$ satisfies:

$$D_{\text{KL}}(P^* | P_M) = \underbrace{D_{\text{KL}}(P^* | P_{\mathcal{M}})}$${\text{model class gap}} + \underbrace{D{\text{KL}}(P_{\mathcal{M}} | P_M)}_{\text{within-class error}}

where $P_{\mathcal{M}}$ is the best model within the class. The model class gap dominated in 78% of cases, confirming that the issue is structural rather than parametric.

4.5 Domain Effects

Forecasting domain modulated coverage independent of model class. Energy forecasting papers achieved the highest median coverage (84.7%), while financial forecasting had the lowest (76.3%). This pattern is consistent with the relative stationarity of the underlying processes: energy demand exhibits strong, stable seasonal patterns that statistical models capture well, while financial returns are heavy-tailed and regime-switching.

4.6 Recalibration Potential

We applied conformal recalibration [10] to a subset of 30 models, splitting the test set into calibration (40%) and evaluation (60%) portions. Post-recalibration coverage improved to a median of 93.1% (IQR: 91.2%--95.4%), with interval width increasing by a median factor of 1.34. The cost of recalibration is thus wider intervals---intervals that honestly reflect the model's actual predictive uncertainty rather than its stated uncertainty.

The recalibrated interval half-width $w_{\text{recal}}$ relates to the original $w_{\text{orig}}$ via:

$$w_{\text{recal}} = w_{\text{orig}} \cdot \hat{q}_{1-\alpha/2}\left( \frac{|y_t - \hat{\mu}$$t|}{w{\text{orig}}} \right)

where the quantile is computed on the calibration set. This multiplicative correction averaged 1.34, meaning original intervals needed to be 34% wider to achieve nominal coverage.

5. Discussion

5.1 Interpretation of Findings

The central result---81% median coverage for nominal 95% intervals---means that published Bayesian prediction intervals are better understood as approximately 80% intervals. This reinterpretation has practical consequences: a decision-maker who treats a published 95% interval as a worst-case planning range will experience outcomes outside that range roughly one in five times, not one in twenty.

The horizon dependence of undercoverage reflects the compounding of model specification errors. At short horizons, even a misspecified model's one-step-ahead predictive distribution may be adequate because the signal-to-noise ratio is high. As horizon increases, the model must predict its own accumulated errors, and misspecification in the error structure becomes the dominant source of uncertainty. The $\log_2(h)$ scaling we observe is consistent with a multiplicative error accumulation model:

$$\text{Var}$${\text{true}}(y{t+h}) = \text{Var}{\text{model}}(y{t+h}) \cdot \left(1 + \delta \right)^{\log_2(h)}

where $\delta$ represents the per-doubling variance underestimation factor. Our estimate of $\delta \approx 0.12$ implies that the true predictive variance grows 12% faster per horizon doubling than the model anticipates.

5.2 Why Neural Forecasters Undercover More

Neural forecasters' worse coverage (74.4% vs 86.1% for ARIMA) arises from three interacting factors. First, approximate Bayesian inference methods (MC dropout, variational Bayes) systematically underestimate posterior width. Foong et al. [11] demonstrated that mean-field variational inference in neural networks produces predictive distributions whose variance is a lower bound on the true posterior predictive variance. Second, neural models are typically more flexible, so the gap between the training distribution and the test distribution matters more---overfitting to temporal patterns that do not generalize produces overconfident forecasts. Third, the lack of analytical prediction interval formulas means practitioners rely entirely on simulation from the approximate posterior, compounding approximation errors.

5.3 Recommendations for Practice

Based on these findings, we propose three recommendations. First, empirical coverage should be reported alongside prediction intervals as standard practice. The PIT (probability integral transform) histogram [4] provides a visual diagnostic, but a single coverage number at the nominal level is the minimum reporting standard. Second, conformal recalibration or similar post-hoc adjustment should be applied whenever a calibration set is available. The 34% median width increase we observed is the price of honesty---intervals that were 34% too narrow were being presented as well-calibrated. Third, model comparison should include coverage as a metric alongside point forecast accuracy. A model with 5% worse RMSE but 10pp better coverage may be preferable for decision-making under uncertainty.

5.4 Limitations

First, the simulation-based coverage estimates for the 73 papers without available code carry uncertainty that we have quantified (mean absolute error 1.7pp vs direct replication) but cannot eliminate. Alternative approaches such as contacting authors for code or restricting to fully reproducible papers would reduce this uncertainty at the cost of sample size and potential selection bias.

Second, our decomposition of undercoverage sources relies on assumptions about the independence of computational, parametric, and misspecification errors. In practice these interact: a poorly converged MCMC chain may mask misspecification, and misspecification may cause convergence difficulties. A joint sensitivity analysis framework, as developed by Kallioinen et al. [12], could provide a more rigorous decomposition.

Third, we measured coverage only at the 95% nominal level. Coverage miscalibration may differ at other quantile levels---in particular, tail coverage at 99% may be even worse. A full calibration audit across multiple quantile levels would require substantially more held-out data per paper.

Fourth, publication bias in our sample is possible: papers reporting prediction intervals may disproportionately represent those where intervals appeared well-calibrated. If so, our coverage estimates are optimistic, and the true median coverage in all Bayesian forecasting applications is lower than 81%.

Fifth, we treated all time points in the test set equally. Conditional coverage---whether intervals are well-calibrated in specific regimes such as high-volatility periods---may reveal even larger gaps. Christoffersen's [13] conditional coverage test applied to a subset of 20 papers rejected conditional coverage in 14 cases (70%), compared to 85% rejection of unconditional coverage.

6. Conclusion

Published Bayesian prediction intervals in time series forecasting systematically undercover, with a median actual coverage of 81% at nominal 95% levels. This gap is not a minor calibration issue but a substantive misrepresentation of predictive uncertainty. The problem intensifies with longer horizons and is more severe for neural forecasters than classical statistical models. Model misspecification---not computational error---is the primary driver. Recalibration is technically straightforward but requires reporting standards that distinguish nominal from empirical coverage. Until such standards are adopted, consumers of Bayesian prediction intervals should apply a mental discount: a published 95% interval is, on average, an 80% interval.

References

[1] Hong, T., Pinson, P., Fan, S., Zareipour, H., Troccoli, A., and Hyndman, R. J., 'Probabilistic Energy Forecasting: Global Energy Forecasting Competition 2014 and Beyond,' International Journal of Forecasting, 32(3), 896--913, 2016.

[2] Held, L., Meyer, S., and Bracher, J., 'Probabilistic Forecasting in Infectious Disease Epidemiology: The 13th Armitage Lecture,' Statistics in Medicine, 36(22), 3443--3460, 2017.

[3] Petropoulos, F., Apiletti, D., Assimakopoulos, V., et al., 'Forecasting: Theory and Practice,' International Journal of Forecasting, 38(3), 845--1130, 2022.

[4] Gneiting, T. and Raftery, A. E., 'Strictly Proper Scoring Rules, Prediction, and Estimation,' Journal of the American Statistical Association, 102(477), 359--378, 2007.

[5] Hyndman, R. J. and Athanasopoulos, G., Forecasting: Principles and Practice, 3rd ed., OTexts, Melbourne, 2021.

[6] Kuleshov, V., Fenner, N., and Ermon, S., 'Accurate Uncertainties for Deep Learning Using Calibrated Regression,' Proceedings of the 35th International Conference on Machine Learning, 2796--2804, 2018.

[7] Bracher, J., Ray, E. L., Gneiting, T., and Reich, N. G., 'Evaluating Epidemic Forecasts in an Interval Format,' PLOS Computational Biology, 17(2), e1008618, 2021.

[8] Talts, S., Betancourt, M., Simpson, D., Vehtari, A., and Gelman, A., 'Validating Bayesian Inference Algorithms with Simulation-Based Calibration,' arXiv:1804.06788, 2018.

[9] Vehtari, A., Gelman, A., Simpson, D., Carpenter, B., and Burkner, P.-C., 'Rank-Normalization, Folding, and Localization: An Improved $\hat{R}$ for Assessing Convergence of MCMC,' Bayesian Analysis, 16(2), 667--718, 2021.

[10] Vovk, V., Gammerman, A., and Shafer, G., Algorithmic Learning in a Random World, Springer, 2005.

[11] Foong, A. Y. K., Li, Y., Hernandez-Lobato, J. M., and Turner, R. E., 'On the Expressiveness of Approximate Inference in Bayesian Neural Networks,' Advances in Neural Information Processing Systems, 33, 15897--15908, 2020.

[12] Kallioinen, N., Paananen, T., Burkner, P.-C., and Vehtari, A., 'Detecting and Diagnosing Prior and Likelihood Sensitivity with Power-Scaling,' Statistics and Computing, 34, 57, 2024.

[13] Christoffersen, P. F., 'Evaluating Interval Forecasts,' International Economic Review, 39(4), 841--862, 1998.