fable.intermittent: benchmarking probabilistic forecasting methods for intermittent time series

1. Introduction

• Asymmetric forecast-error costs make point forecasts insufficient; inventory control can require the upper tail of the demand distribution.The optimal forecast may be a quantile rather than the mean.Source: §1 · p. 2
• Probabilistic methods span statistical and machine-learning approaches, but few have public implementations, and existing implementations are scattered across frameworks.This makes systematic comparisons on the same data laborious.Source: §1 · p. 2
• fable.intermittent implements literature-based probabilistic methods within fable and fits and evaluates them on collections of series through one pipeline.The package uses the fable framework to provide a straightforward common workflow.Source: §1 · p. 2
• TWEES is exponential smoothing with a Tweedie predictive distribution; its separate recursions model the predictive mean and demand-occurrence probability.In the reported experiments, TWEES is among the most accurate models.Source: §1 · pp. 2–3
• tweedieDistr evaluates the Tweedie density up to 20 times faster than tweedie while retaining the same numerical accuracy.Repeated density evaluation is computationally demanding during model estimation, especially across large collections of series.Source: §1 · p. 3

2.1. Models

• Implemented models cover bootstrapping, Bayesian methods, exponential smoothing, ARMA-based approaches, and static empirical or parametric distributions.All implemented models are local; some exponential-smoothing and Bayesian methods model occurrence separately from demand size.Source: Table 1 · p. 3, §2.1 · p. 4

2.2. Datasets

• The four datasets cover automotive parts, pasta sales, Royal Air Force spare parts, and Walmart product sales, with differing intermittency and demand-size variability.Pasta and RAF have demand spikes (high CV2), and RAF has the most intermittent demand; Table 2 summarizes series count, length, zero share, ADI, and CV2.Source: §2.2 · pp. 4–5, Table 2 · p. 5
• The datasets are stored in tsibble format, with observations, timestamps, series identifiers, and, where available, exogenous covariates.The data are formatted for ingestion by the fable pipeline.Source: §2.2 · pp. 4, 6

2.3. Easy benchmarking with fable.intermittent

• The common fable syntax represents probabilistic forecasts as distributional objects, enabling the same evaluation pipeline to compare 11 models on collections of series.The example pipeline fits, forecasts, and evaluates models using fabletools functions.Source: §2.3 · pp. 6–7

3. tweedieDistr: a fast implementation of the Tweedie distribution

• For 1 < ρ < 2, the Tweedie distribution has a compound Poisson-Gamma form: Gamma sums model positive demand and a Poisson component supplies the zero mass.This structure supports intermittent demand; the Tweedie family is characterized by a power mean-variance relationship.Source: §3.1 · pp. 7–8
• At fixed mean and power, changing dispersion changes both the probability of zero and the distribution of positive values.Figure 2 illustrates the comparison for mean 2 and power 1.2.Source: §3.1 · p. 8, Figure 2 · p. 8
• The density evaluation can become a bottleneck when fitting Tweedie distributions across many series or repeatedly using the distribution as a time-series likelihood.For positive values, evaluation requires computing infinite sums.Source: §3.1 · p. 8
• dtweedieDistr averages 15, 300, and 600 times faster for density, cumulative, and inverse-cumulative functions, respectively.The speedups are attributed partly to C++ implementations, algorithm choices, and vectorisation over mean, dispersion, and power.Source: §3.1 · p. 9
• The density function matches tweedie’s accuracy, the cumulative function is 2% more accurate on average, and the inverse cumulative function can error near ρ = 2.tweedieDistr currently handles only ρ ∈ (1, 2), while tweedie supports ρ > 0.Source: §3.1 · p. 9

4. The TWEES model

• TWEES models observations with a Tweedie distribution whose mean follows a simple exponential smoothing recursion.The mean update uses learnable smoothing parameters and the training-data average.Source: §4 · p. 9
• A second damped smoothing process models occurrence separately, allowing occurrence and demand size to evolve independently.It is applied to the indicator that demand is positive.Source: §4 · p. 9
• The model sets Tweedie dispersion so its zero mass matches the forecast probability of no occurrence.This links the occurrence forecast to the distribution’s zero probability.Source: §4 · pp. 9–10
• The model jointly estimates its parameters by maximizing the Tweedie log-likelihood; longer-horizon forecasts simulate both smoothing processes forward.One-step-ahead forecasts use observed values directly.Source: §4 · p. 11
• TWEES restricts ρ to (1.2, 1.8): the upper bound limits density-evaluation terms, and the lower bound supports numerical stability.Before fitting, each series is rescaled by its median positive value.Source: §4 · p. 11

5.1. Setup

• Evaluation used two expanding-window splits at each dataset’s native horizon, averaging scores across windows and time series and using an MCS at α = 0.05 to assess significance.For methods requiring sample paths, 10^5 paths were drawn; quantile scores used levels 0.5, 0.75, 0.835, 0.975, and 0.995, scaled by the naive forecast’s MAE.Source: §5.1 · pp. 11–12

5.2. Results

• On Auto, the parametric static model beats the empirical distribution on all metrics, although their difference is often minimal elsewhere.Static methods are particularly effective when historical data is limited, as on Auto, or demand is very sparse, as on RAF.Source: §5.2 · pp. 12–13, Table 3 · p. 12, Table 5 · p. 13
• Across four datasets, no method dominates; static benchmarks are hard to beat on RMSSE and central quantiles, while TWEES ranks in the top three on RMSSE and sQS0.835.TWEES is also almost always among the top five; on TinyM5, TWEES and NEGBINES are always in the MCS.Source: §5.2 · pp. 12–13, Table 3 · p. 12, Table 4 · p. 13
• Exponential smoothing models generally perform better; TWEES’s two-process occurrence-and-size design outperformed a mean-only approach in preliminary experiments, while NEGBINES usually ranks better on Pasta.The negative binomial method remains competitive overall.Source: §5.2 · p. 13, Table 4 · p. 13
• HSPES models occurrence effectively, but its demand-size process uses few non-zero observations and a light-tailed Poisson likelihood that struggles with high quantiles.This limits HSPES among the three exponential-smoothing methods.Source: §5.2 · pp. 13–14
• GAMPOISB’s low-variance Gamma component can make it behave like a Poisson model at high quantiles, whereas BETANBB performs well on Auto’s extreme quantiles and central quantiles across datasets.The authors note the suitability of the negative binomial distribution for intermittent demand.Source: Table 3 · p. 12, §5.2 · p. 14
• Bootstrap methods lack consistency: VZ’s high-quantile performance is strongest on longer series, while WSS’s rounding can bias demand-size quantiles upward.VZ’s mean forecasts are accurate only on Auto, the least sparse dataset, and its Croston decomposition limits resampling to the observed non-zero values; WSS performs well on RAF’s lowest quantiles, which are driven by its occurrence Markov Chain.Source: §5.2 · p. 14
• ARMA-based methods can perform well on RMSSE but rank 8th to 11th on central quantiles across datasets, with Gaussian symmetric tails poorly matching zeros and demand spikes.Their probabilistic scores are particularly weak despite sometimes competitive mean forecasts.Source: §5.2 · p. 14

6. A worked example of intermittent time series reconciliation

• The hierarchy contains 118 bottom-level pasta-product series, aggregated to four brands and then total sales.The bottom series receive TWEES forecasts and the aggregated series ETS forecasts; the comparison fits ETS to all series and reconciles with MinT.Source: §6 · pp. 15–16, Table 7 · p. 16
• The example uses BUIS, MixCond, and TDcond to reconcile mixed forecasts, with forecasts produced 28 steps ahead.The fable integration provides access to reconciliation tools, and the example presents fitting and reconciliation through a concise pipeline.Source: §6 · p. 16
• Mixed forecasts are consistently better than ETS; BUIS and MixCond improve base forecasts on all scores except the two highest quantiles, while TDcond and MinT are less effective.BUIS and MixCond are tailored to count-valued bottom forecasts; MinT relies on Gaussian base forecasts, which are often ineffective for intermittent series.Source: Table 7 · p. 16, §6 · pp. 16–17

7. Conclusion

• fable.intermittent gathers probabilistic forecasting methods in the fable ecosystem, where they can be fitted, forecast, evaluated, and visualized through one pipeline.The models can also be combined with other fable models and reconciled across hierarchies using the same ecosystem.Source: §7 · p. 17
• The paper introduces TWEES, an exponential smoothing model with a Tweedie predictive distribution, and provides a faster Tweedie implementation in tweedieDistr.The Tweedie distribution implementation is released as a separate R package.Source: §7 · p. 17
• Future directions include refining static distributions, selecting PARAMSD models with a quantile score, and extending some models to use exogenous variables such as promotions.The authors note that promotions strongly affect demand.Source: §7 · p. 17

Appendix A. Methods from the literature

• The appendix defines the training series and forecast horizon, then reviews Croston’s decomposition into positive demand sizes and intervals between demands.It also defines the binary occurrence process as 1 when demand is positive and 0 otherwise.Source: Appendix A · p. 22
• WSS models binary occurrence with a two-state Markov chain and forms a zero-mass/positive-demand mixture; Croston-based bootstrapping separately resamples intervals and demand sizes.WSS adds Gaussian noise to sampled demand sizes, rounds them, and sets negative values to zero; forecasts are generated by autoregressive sampling.Source: Appendix A · p. 22, VZ (Zhou and Viswanathan, 2011) · p. 22

VZ (Zhou and Viswanathan, 2011).

• The Bayesian negative-binomial model uses a globally estimated count parameter and time-varying Beta-distributed probabilities, with parameters fitted by negative log-likelihood.Its forecasts use autoregressive sampling.Source: VZ (Zhou and Viswanathan, 2011) · p. 22

GAMPOISB (Harvey and Fernandes, 1989).

• GAMPOISB uses a hierarchical Gamma-Poisson structure with parameters that vary over time; one-step forecasts are closed-form, while longer horizons require autoregressive sampling.Its parameters and updating coefficient are learned by minimizing negative log-likelihood.Source: GAMPOISB (Harvey and Fernandes, 1989) · p. 23

Appendix A.3. Exponential smoothing models

• HSPES smooths demand occurrence and positive demand size separately, then combines them into a hurdle-shifted Poisson forecast distribution.The smoothing parameters and initial values are learned by minimizing negative log-likelihood.Source: GAMPOISB (Harvey and Fernandes, 1989) · p. 23

NEGBINES (Snyder et al., 2012).

• NEGBINES uses damped or undamped exponential smoothing for a negative-binomial mean; its forecasts are closed-form only one step ahead and necessarily unimodal.For longer horizons, forecasts use autoregressive sampling; the negative-binomial distribution can have a wider right tail than HSPES.Source: NEGBINES (Snyder et al., 2012) · p. 23

EMPSD (Boylan and Babai, 2022).

• EMPSD treats observations as i.i.d. draws from a count distribution, giving the same predictive distribution at every forecast horizon.A well-chosen parametric distribution has been reported as a strong benchmark for slow-moving items, particularly with short demand histories.Source: PARAMSD (Kolassa, 2016) · p. 24

PARAMSD (Kolassa, 2016).

• PARAMSD fits five candidate distributions: Poisson, negative binomial, their hurdle-shifted counterparts, and Tweedie.For count data, the Tweedie density is discretized by integration over [y − 1/2, y + 1/2).Source: PARAMSD (Kolassa, 2016) · p. 24
• PARAMSD selects the fitted distribution with the lowest BIC by default, or AIC when specified.Both criteria use the fitted log-likelihood and parameter count.Source: PARAMSD (Kolassa, 2016) · p. 24

MARWAL (Sbrana, 2025).

• MARWAL combines an ARMA(1,1) process with a Markov chain switching occurrence between zero and one; its ARMA parameters are estimated in closed form.The process uses ω, the mean of non-zero training values, and a validation loop selects the starting value t0.Source: MARWAL (Sbrana, 2025) · p. 24
• MARWAL computes forecasts in closed form for any horizon and assumes Gaussian forecast distributions, setting negative samples or quantiles to zero.Source: MARWAL (Sbrana, 2025) · p. 24

NNARMA (Sbrana and Babai, 2026).

• NNARMA uses a Gaussian ARMA(1,1) process with positive-mean parameter constraints; setting negative forecast mass to zero makes the distribution ill-posed because it changes the mean.Source: NNARMA (Sbrana and Babai, 2026) · p. 25
• NNARMA estimates its three parameters by least squares on innovations and computes forecast means and variances in closed form for any horizon.Its inputs are multiplicatively deseasonalised, then forecasts receive seasonal factors based on average ratios to the moving average.Source: NNARMA (Sbrana and Babai, 2026) · p. 25

• tweedieDistr evaluates the density through its series expansion for all ϕ and ρ in (1, 2), whereas tweedie uses that expansion only in some cases and otherwise interpolates tabulated data with a saddlepoint approximation.Source: Appendix B · p. 25
• For positive y, tweedieDistr expresses the cumulative distribution function using a Compound-Poisson parametrisation with Poisson mass and Gamma cumulative-distribution components.Source: Appendix B · p. 25
• tweedieDistr computes quantiles with a C++ Newton-Raphson algorithm, falling back to bisection when convergence is not reached in a few steps; tweedie uses Brent’s method.The fallback is used about 4% of the time.Source: Appendix B · p. 26
• Both packages sample by drawing a Poisson count and then a Gamma value using the Compound-Poisson representation.Source: Appendix B · p. 26

• The benchmark compared identical inputs over 18 dispersion and 9 power values, repeating each parameter-pair experiment 200 times, except where quantile-function errors reduced repetitions by up to 54%.Source: Appendix C · p. 26
• Timing used 10 replications across sample sizes n = 2, 10, 50, 200, 500, and 2000; the reported speedup is tweedie’s mean time divided by tweedieDistr’s.The benchmark generated density and cumulative-function results at n = 200 and quantile results at n = 10.Source: Appendix C · p. 26
• Up to 1000 times faster for cumulative distribution and quantile evaluation; density evaluation remains 1.5 to 20 times faster at n = 2000.Source: Appendix C · p. 26
• Density differences were below machine epsilon; median absolute differences for cumulative probabilities and quantiles were also below machine epsilon, while maximum differences were about 5 × 10^-2 and 7 × 10^-4, respectively.Monte Carlo simulations favored fable.intermittent in the few cases of disagreement.Source: Appendix C · p. 27