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

# fable.intermittent: benchmarking probabilistic forecasting methods for intermittent time series
Source: [https://arxiv.org/html/2609.28607](https://arxiv.org/html/2609.28607)
Stefano DamatoEmail:[stefano\.damato@supsi\.ch](mailto:[email protected])Corresponding author:Corresponding authorAffiliation:SUPSI, Istituto Dalle Molle di Studi sull’Intelligenza Artificiale \(IDSIA\), Lugano, SwitzerlandLorenzo ZambonEmail:[lorenzo\.zambon@supsi\.ch](mailto:[email protected])Affiliation:SUPSI, Istituto Dalle Molle di Studi sull’Intelligenza Artificiale \(IDSIA\), Lugano, SwitzerlandGiorgio CoraniEmail:[giorgio\.corani@supsi\.ch](mailto:[email protected])Affiliation:SUPSI, Istituto Dalle Molle di Studi sull’Intelligenza Artificiale \(IDSIA\), Lugano, SwitzerlandDario AzzimontiEmail:[dario\.azzimonti@supsi\.ch](mailto:[email protected])Affiliation:SUPSI, Istituto Dalle Molle di Studi sull’Intelligenza Artificiale \(IDSIA\), Lugano, Switzerland

###### Abstract

Intermittent time series are common in spare\-parts demand and retail sales\. Since the cost of forecast errors is typically asymmetric, decisions such as inventory control require the full predictive distribution rather than a point forecast\. Many probabilistic forecasting methods have been proposed; their implementations, however, are scattered across different software frameworks, making it difficult to compare them systematically\. We introducefable\.intermittent, anRpackage that implements several probabilistic forecasting methods for intermittent series within thefableframework\. The package allows several models to be fitted and evaluated on a collection of time series through a single, simple forecasting pipeline\. We also introduce TWEES, a new exponential smoothing model with a Tweedie predictive distribution\. Fitting TWEES requires repeated evaluation of the computationally demanding Tweedie density\. We also release theRpackagetweedieDistr, whose implementation of the Tweedie distribution is substantially faster than the existing one while preserving the same numerical accuracy\. We evaluate the methods implemented infable\.intermittenton four datasets, also released in the package\.

###### Keywords:

Intermittent demand; Probabilistic forecasting; Tweedie distribution; Exponential smoothing; Forecast reconciliation

\\nopreprintlinetrue

## 1Introduction

Intermittent time series are non\-negative valued time series characterized by positive values interspersed with zeros\([Boylan and Syntetos, 2021](https://arxiv.org/html/2609.28607#bib.bib7)\)\. They commonly arise, for instance, in supply\-chain applications\([Syntetos et al\., 2016](https://arxiv.org/html/2609.28607#bib.bib48)\)and retail\([Fildes et al\., 2022](https://arxiv.org/html/2609.28607#bib.bib19)\)\. Most methods for intermittent time series\([Croston, 1972](https://arxiv.org/html/2609.28607#bib.bib11);[Syntetos and Boylan, 2005](https://arxiv.org/html/2609.28607#bib.bib49);[Teunter et al\., 2011](https://arxiv.org/html/2609.28607#bib.bib51)\)provide point forecasts, which are insufficient to support decision making\. The costs of forecasting errors are usually asymmetric, so that the optimal forecast is a quantile rather than the mean of the predictive distribution\([Kolassa, 2016](https://arxiv.org/html/2609.28607#bib.bib29)\); inventory control, for instance, requires the upper tail of the predictive distribution of the demand\([Prak and Teunter, 2019](https://arxiv.org/html/2609.28607#bib.bib38)\)\.

Probabilistic models for intermittent time series span several methodologies, from statistical to machine learning approaches\([Lang et al\., 2026](https://arxiv.org/html/2609.28607#bib.bib31)\): static empirical and parametric distributions\([Boylan and Babai, 2022](https://arxiv.org/html/2609.28607#bib.bib8);[Kolassa, 2016](https://arxiv.org/html/2609.28607#bib.bib29)\), bootstrap methods\([Willemain et al\., 2004](https://arxiv.org/html/2609.28607#bib.bib55);[Zhou and Viswanathan, 2011](https://arxiv.org/html/2609.28607#bib.bib58)\), Bayesian dynamic models\([Harvey and Fernandes, 1989](https://arxiv.org/html/2609.28607#bib.bib22);[Babai et al\., 2021](https://arxiv.org/html/2609.28607#bib.bib5)\), exponential smoothing models for count data\([Snyder et al\., 2012](https://arxiv.org/html/2609.28607#bib.bib44);[Svetunkov and Boylan, 2023](https://arxiv.org/html/2609.28607#bib.bib47)\), ARMA\-based models\([Sbrana, 2025](https://arxiv.org/html/2609.28607#bib.bib41);[Sbrana and Babai, 2026](https://arxiv.org/html/2609.28607#bib.bib42)\), Gaussian processes for intermittent data\([Damato et al\., 2025](https://arxiv.org/html/2609.28607#bib.bib12)\)\. Probabilistic forecasts can also be produced by a global deep\-learning model such as DeepAR\([Salinas et al\., 2020](https://arxiv.org/html/2609.28607#bib.bib40)\), trained with a distribution head suitable for intermittent data\.

Only a few methods, however, have publicly available implementations\. InR\([R Core Team, 2026](https://arxiv.org/html/2609.28607#bib.bib39)\), thesmoothpackage\([Svetunkov, 2026](https://arxiv.org/html/2609.28607#bib.bib46)\)implements the intermittent exponential smoothing \(iETS\) model\([Svetunkov and Boylan, 2023](https://arxiv.org/html/2609.28607#bib.bib47)\); whilefableCount\([Almeida and Vieira, 2024](https://arxiv.org/html/2609.28607#bib.bib2)\)makes available the models GLARMA and INGARCH for time series of counts\. InPython,GluonTS\([Alexandrov et al\., 2020](https://arxiv.org/html/2609.28607#bib.bib1)\)provides global deep\-learning models\. Other models, e\.g\.[Damato et al\. \(2025\)](https://arxiv.org/html/2609.28607#bib.bib12)and[Sbrana and Babai \(2026\)](https://arxiv.org/html/2609.28607#bib.bib42), are only released in project repositories\. Available implementations, moreover, are spread across different frameworks, making it laborious to compare several probabilistic methods for intermittent time series on the same data\.

We fill this gap withfable\.intermittent, anRpackage which extends thefableframework\([O’Hara\-Wild et al\., 2026a](https://arxiv.org/html/2609.28607#bib.bib35)\)by implementing several probabilistic forecasting methods from the literature\. The package allows these methods to be fitted and evaluated on collections of time series with a straightforward pipeline\.

fable\.intermittentalso provides TWEES, a novel probabilistic model for intermittent time series\. It is an exponential smoothing model with a Tweedie predictive distribution\. The Tweedie distribution is suitable for intermittent data as it combines a point mass at zero with a continuous mixture of Gamma distributions over the positive real line\. Recent studies\([Damato et al\., 2025](https://arxiv.org/html/2609.28607#bib.bib12);[Damato et al\., 2026](https://arxiv.org/html/2609.28607#bib.bib13)\)show that it improves the estimation of the highest quantiles compared with distributions traditionally used for forecasting of intermittent time series, such as the negative binomial and hurdle\-shifted negative binomial\. Following the design of[Snyder et al\. \(2012\)](https://arxiv.org/html/2609.28607#bib.bib44), we use two different exponential smoothing recursions to model the mean of the predictive distribution and the probability of demand occurrence\. In our experiments, TWEES is among the most accurate models\.

The evaluation of the Tweedie density, however, is computationally demanding\([Dunn and Smyth, 2005](https://arxiv.org/html/2609.28607#bib.bib15)\)\. This is problematic when the Tweedie distribution is the predictive distribution of a time series model, as the density must be evaluated repeatedly during estimation, especially when the model is fitted on a large collection of time series\. We address this limitation withtweedieDistr, anRpackage providing a novel implementation of the Tweedie distribution\. For density evaluation, our implementation is up to 20 times faster than that of thetweediepackage\([Dunn, 2026](https://arxiv.org/html/2609.28607#bib.bib14)\), while retaining the same numerical accuracy\. Bothfable\.intermittentandtweedieDistrare released under the LGPL\-3\.0 License111[https://gnu\.org](https://gnu.org/)and available on CRAN222[https://cran\.r\-project\.org/web/packages/fable\.intermittent/index\.html](https://cran.r-project.org/web/packages/fable.intermittent/index.html), [https://cran\.r\-project\.org/web/packages/tweedieDistr/index\.html](https://cran.r-project.org/web/packages/tweedieDistr/index.html)\., the official archive ofRpackages\.

The paper is organized as follows\. Sec\.[2](https://arxiv.org/html/2609.28607#S2)presents the models implemented infable\.intermittent, the datasets available within the package, and an example of the forecasting pipeline\. Sec\.[3](https://arxiv.org/html/2609.28607#S3)describes our implementation of the Tweedie distribution, available intweedieDistr, and its speedup over the implementation provided bytweedie\. Sec\.[4](https://arxiv.org/html/2609.28607#S4)introduces the TWEES model\. Sec\.[5](https://arxiv.org/html/2609.28607#S5)evaluates the models offable\.intermittenton several datasets\. Sec\.[6](https://arxiv.org/html/2609.28607#S6)shows how the resulting forecasts can be probabilistically reconciled usingfable\.bayesRecon\([Azzimonti et al\., 2026](https://arxiv.org/html/2609.28607#bib.bib4)\)\. Finally, Sec\.[7](https://arxiv.org/html/2609.28607#S7)presents our conclusions\.

## 2Literature models and datasets provided byfable\.intermittent

Table 1:Methods from the literature implemented in the package, grouped by methodological family\.### 2\.1Models

The package implements the probabilistic models for intermittent time series listed in Tab\.[1](https://arxiv.org/html/2609.28607#S2.T1)\. The last column indicates which models rely on*Croston’s decomposition*\([Croston, 1972](https://arxiv.org/html/2609.28607#bib.bib11)\), which decomposes the series into an occurrence process, indicating whether demand is positive, and a demand size process, containing the positive demand values\. According to the model classification by[Januschowski et al\. \(2020\)](https://arxiv.org/html/2609.28607#bib.bib26), all the implemented models are local\.

Bootstrapping methods for intermittent demand, reviewed by[Hasni et al\. \(2019\)](https://arxiv.org/html/2609.28607#bib.bib23), resample past demand values to estimate the demand size process\. In particular, we implement WSS\([Willemain et al\., 2004](https://arxiv.org/html/2609.28607#bib.bib55)\), which models the occurrence process as a two\-state Markov Chain, and VZ\([Zhou and Viswanathan, 2011](https://arxiv.org/html/2609.28607#bib.bib58)\), which bootstraps demand intervals\. The Bayesian methods of[Harvey and Fernandes \(1989\)](https://arxiv.org/html/2609.28607#bib.bib22)use a prior for the parameters of the predictive distribution, and update it as the time series evolves; GAMPOISB uses a Gamma prior for the Poisson parameter, while BETANBB uses a Beta prior for the probability parameter of a negative binomial model, which count parameter is learned via maximum likelihood\.[Snyder et al\. \(2012\)](https://arxiv.org/html/2609.28607#bib.bib44)couple exponential smoothing with predictive distributions suitable for count data: in HSPES, two processes smooth the rate of a Poisson on the demand size and a probability of occurrence to get a hurdle\-shifted Poisson forecast distribution, while in NEGBINES a single process determines the mean of a negative binomial distribution\.

ARMA\-based models use a Gaussian likelihood, but the mean is tailored to properties of intermittent demand:[Sbrana \(2025\)](https://arxiv.org/html/2609.28607#bib.bib41)additionally models the occurrence with a Markov chain, while[Sbrana and Babai \(2026\)](https://arxiv.org/html/2609.28607#bib.bib42)constrain the mean of the forecast distribution to be positive\. Finally, static distributions treat time series data as i\.i\.d\. observations; they are considered a solid benchmark for intermittent time series\([Spiliotis et al\., 2021](https://arxiv.org/html/2609.28607#bib.bib45)\)\. We implement EMPSD, the simple empirical distribution of the data\([Boylan and Babai, 2022](https://arxiv.org/html/2609.28607#bib.bib8)\), and PARAMSD that fits a parametric distribution\([Kolassa, 2016](https://arxiv.org/html/2609.28607#bib.bib29)\)\. Our implementation fits five candidate parametric distributions \(Poisson, negative binomial, their hurdle\-shifted versions and a discretized Tweedie distribution\) and selects among them by minimizing the Bayes Information Criterion\([Schwarz, 1978](https://arxiv.org/html/2609.28607#bib.bib43), BIC,\)\. See[A](https://arxiv.org/html/2609.28607#A1)for more detailed descriptions\.

Some of the methods above share the same routines, which are run several times for each time series fit\. Our package improves the computational efficiency of these methods by implementing the repeated routines, such as exponential smoothing and ARMA recursions and Bayesian updates, inC\+\+via the packageRcpp\([Eddelbuettel and Francois, 2011](https://arxiv.org/html/2609.28607#bib.bib17)\)\.

### 2\.2Datasets

Table 2:Overview of the datasets included in the package\.NNandLLdenote respectively the number of time series in the data set and their length\. The last three columns report the overall proportion of zero observations, the median across series of the average inter\-demand interval \(ADI\), and the median squared coefficient of variation of the non\-zero demand sizes \(CV2\\mathrm\{CV\}^\{2\}\)\.The package provides four datasets of intermittent time series \(Tab\.[2](https://arxiv.org/html/2609.28607#S2.T2)\) in a format ready to be ingested by thefablepipeline\.*Auto*\([Türkmen et al\., 2021](https://arxiv.org/html/2609.28607#bib.bib52)\)contains monthly demand for automotive spare parts; its series are short and only mildly intermittent\.*Pasta*\([Mancuso et al\., 2021](https://arxiv.org/html/2609.28607#bib.bib33)\)contains daily sales of 118 pasta products, including promotion covariates\.*RAF*\([Syntetos and Boylan, 2005](https://arxiv.org/html/2609.28607#bib.bib49)\)contains monthly demand for spare parts of the Royal Air Force\.TinyM5is a subset from the data set of the M5 competition\([Makridakis et al\., 2022](https://arxiv.org/html/2609.28607#bib.bib32)\)and contains daily sales of different products at Walmart; we use the selection of[Joachimiak \(2022\)](https://arxiv.org/html/2609.28607#bib.bib27)\. Both*Pasta*and*RAF*are characterized by demand spikes \(highCV2\\mathrm\{CV\}^\{2\}\), which challenge the tails of the predictive distributions\.*RAF*has also the most intermittent demand\.

Figure 1:Classification of the series of the four data sets according to average inter\-demand interval \(ADI, log scale\) and the squared coefficient of variation of the non\-zero demand sizes \(CV2, linearly\-adjusted log scale\)\. Dashed lines mark the cut\-offs of[Syntetos et al\. \(2005\)](https://arxiv.org/html/2609.28607#bib.bib50)\(ADI=1\.32\\mathrm\{ADI\}=1\.32,CV2=0\.49\\mathrm\{CV\}^\{2\}=0\.49\)\. A small number of series withCV2=0\\mathrm\{CV\}^\{2\}=0\(constant non\-zero demand size\) appear in the plot exactly on zero because of the linearly adjusted log scale\.Fig\.[1](https://arxiv.org/html/2609.28607#S2.F1)shows the distribution of the average inter\-demand interval \(ADI\) and the squared coefficient of variation of non\-zero demand sizes \(CV2\) across the time series of the different datasets, following the classification of[Syntetos et al\. \(2005\)](https://arxiv.org/html/2609.28607#bib.bib50)\.

The data sets are released in thetsibble\([Wang et al\., 2020](https://arxiv.org/html/2609.28607#bib.bib53)\)format, a tidy data structure to store time series\. In each data set, the columnvaluestores the actual observations, and the columnindexsaves their timestamps\. One or more additional columns serve as time series identifiers; others may provide exogenous covariates\.

### 2\.3Easy benchmarking withfable\.intermittent

The packagefable\.intermittentextendsfable\([O’Hara\-Wild et al\., 2026a](https://arxiv.org/html/2609.28607#bib.bib35)\), the tidy time series modelling framework of thetidyverts333[https://tsibble\.tidyverts\.org/](https://tsibble.tidyverts.org/)ecosystem\. All models share a common syntax and can be straightforwardly fitted to collections of time series stored astsibbleobjects\. Their probabilistic forecasts are represented asdistributionalobjects\([O’Hara\-Wild et al\., 2026c](https://arxiv.org/html/2609.28607#bib.bib37)\)\. This common representation allows the models to be evaluated through the same pipeline, based onfabletools\([O’Hara\-Wild et al\., 2026b](https://arxiv.org/html/2609.28607#bib.bib36)\)functions, and compared with other models available infableor its extensions\. Listings[2](https://arxiv.org/html/2609.28607#S2.F2)and[3](https://arxiv.org/html/2609.28607#S2.F3)show the basic pipeline to obtain a tsibble,res, with the results of Sec\.[5\.2](https://arxiv.org/html/2609.28607#S5.SS2)\. In a few lines, we train and evaluate 11 models, the ten models of Tab\.[1](https://arxiv.org/html/2609.28607#S2.T1)and the new model TWEES introduced in Sec\.[4](https://arxiv.org/html/2609.28607#S4)\.

library\(fable\)

library\(fable\.intermittent\)

library\(dplyr\)

h<\-28

train<\-pasta\|\>filter\(index<=max\(index\)\-h\)

test<\-pasta\|\>filter\(index\>max\(index\)\-h\)

head\(train,1\)

indexbrandproductvaluepromotion

<date\><chr\><chr\><dbl\><int\>

12014\-01\-02B1170

Figure 2:Load thepastadata set and create a train/test split\. Lines 12\-16 show the shape of the tsibbletrainwhere the columnvaluecontains the values of the time series\.fit<\-train\|\>

model\(

wss=WSS\(value\),

vz=VZ\(value\),

betanbb=BETANBB\(value\),

gampoisb=GAMPOISB\(value\),

hspes=HSPES\(value\),

negbines=NEGBINES\(value\),

twees=TWEES\(value\),

empsd=EMPSD\(value\),

paramsd=PARAMSD\(value\),

marwal=MARWAL\(value\),

nnarma=NNARMA\(value\)

\)

fc<\-fit\|\>

forecast\(h=h\)

res<\-fc\|\>

accuracy\(data,measures=experiment\_measures\)\|\>

group\_by\(\.model\)\|\>

summarise\(across\(where\(is\.numeric\),\\\(x\)mean\(x,na\.rm=TRUE\)\)\)

Figure 3:The pipeline for fitting, forecasting and evaluation on the Pasta dataset\. Here,experiment\_measuresis a list offabletoolsforecast evaluation measures;\.modelis the column where model names are stored, over which we group the results and compute the mean\.

## 3tweedieDistr: a fast implementation of the Tweedie distribution

### 3\.1The Tweedie distribution

Figure 4:Two Tweedie distributions with the same mean and power parameter, but different dispersion\.The Tweedie is a family of exponential dispersion distributions\([Jørgensen, 1987](https://arxiv.org/html/2609.28607#bib.bib28)\)characterized by a power mean\-variance relationship\. IfY∼Tw⁡\(μ,ϕ,ρ\)Y\\sim\\mathrm\{Tw\}\(\\mu,\\phi,\\rho\), then

Var⁡\(Y\)=ϕ​μρ,\\mathrm\{Var\}\(Y\)=\\phi\\mu^\{\\rho\},\(1\)whereμ\>0\\mu\>0is the mean,ϕ\>0\\phi\>0is the dispersion parameter, andρ\>0\\rho\>0is the power parameter\. As in[Damato et al\. \(2025\)](https://arxiv.org/html/2609.28607#bib.bib12), we limit the power toρ∈\(1,2\)\\rho\\in\(1,2\): for this choice of the power parameter, the Tweedie distribution can be parameterized as a Compound Poisson Gamma model:

Y=∑i=1NXi,with\\displaystyle Y=\\sum\_\{i=1\}^\{N\}X\_\{i\},\\quad\\text\{with\}\(2\)Xi​∼i\.i\.d\.​Gamma​\(α,β\),N∼Pois⁡\(λ\),\\displaystyle X\_\{i\}\\overset\{i\.i\.d\.\}\{\\sim\}\\mathrm\{Gamma\}\(\\alpha,\\beta\),\\quad N\\sim\\mathrm\{Pois\}\(\\lambda\),whereλ\\lambda,α\\alpha, andβ\\betaare functions of the original parameters\([Dunn and Smyth, 2005](https://arxiv.org/html/2609.28607#bib.bib15)\)\. From this parametrisation it can be seen why the Tweedie is suitable for intermittent demand: the sum of Gamma random variables models the positive mass and allows for a flexible right tail behaviour, while the Poisson distribution provides an atomic mass in zero computed from Eq\. \([2](https://arxiv.org/html/2609.28607#S3.E2)\) as

P⁡\(Y=0\)=P⁡\(N=0\)=e−λ=e−μ2−ρϕ⁡\(2−ρ\)P\(Y=0\)=P\(N=0\)=e^\{\-\\lambda\}=e^\{\-\\frac\{\\mu^\{2\-\\rho\}\}\{\\phi\(2\-\\rho\)\}\}\(3\)Fig\.[4](https://arxiv.org/html/2609.28607#S3.F4)shows two Tweedie distributions with the same mean \(μ=2\\mu=2\) and power \(ρ=1\.2\\rho=1\.2\), but different dispersionϕ\\phi; this results in different probabilities of zero and different distributions of positive values\.

### 3\.2A novel, fast implementation intweedieDistr

Fory\>0y\>0, the evaluation of the Tweedie density is computationally cumbersome because it requires computing infinite sums\. TheRpackagetweedie\([Dunn, 2026](https://arxiv.org/html/2609.28607#bib.bib14)\)uses evaluation strategies based on series expansion\([Dunn and Smyth, 2005](https://arxiv.org/html/2609.28607#bib.bib15)\)and Fourier inversion of the characteristic function\([Dunn and Smyth, 2008](https://arxiv.org/html/2609.28607#bib.bib16)\)\. The running times of these implementations, however, can become a bottleneck when fitting a Tweedie distribution to each time series of a large collection, or when the Tweedie distribution is used as the likelihood of a time series model and must therefore be evaluated repeatedly\.

Our packagetweedieDistrprovides a substantially faster implementation of the Tweedie distribution\. Similarly to the package of[Dunn \(2026\)](https://arxiv.org/html/2609.28607#bib.bib14),tweedieDistrfollows theRstatistical convention implementing functions for the evaluation of the density function \(dtweedie\(\)\), the cumulative distribution function \(ptweedie\(\)\) and its inverse \(qtweedie\(\)\), and for random sample generation \(rtweedie\(\)\)\.

In Fig\.[5](https://arxiv.org/html/2609.28607#S3.F5)we show the speedup on the three core functionsdtweedie\(\),ptweedie\(\), andqtweedie\(\)provided bytweedieDistrovertweediefor different values of power, dispersion, and sample size\. See[C](https://arxiv.org/html/2609.28607#A3)for more details about our experimental setup\.

The implementations oftweedieDistrare on average 15, 300, and 600 times faster than their counterpart intweedie\. The speedup can be attributed, in part, to the use ofC\+\+viaRcpp/RcppArmadillo\([Eddelbuettel and Francois, 2011](https://arxiv.org/html/2609.28607#bib.bib17)\)in our implementations\. Moreover, a substantial advantage is provided by different choices in terms of algorithms, detailed in[B](https://arxiv.org/html/2609.28607#A2)\. Finally,tweedieDistrvectorises over mean, dispersion and power, whereastweedierequires anR\-level loop for the power parameter, as it only treats it as a scalar value\.

![Refer to caption](https://arxiv.org/html/2609.28607v1/tweedie_speedup_heatmap_stacked_raw_named.png)Figure 5:Improvement of the running times of the core density \(dtweedie\(\)\), cumulative \(ptweedie\(\)\) and quantile \(qtweedie\(\)\) functions for different choices ofϕ\\phiandρ\\rho\. For additional details on the setup of this evaluation, see[C](https://arxiv.org/html/2609.28607#A3)\.The functiondtweedie\(\)fromtweedieDistrachieves the same accuracy as the one intweedie\. Instead, the cumulative function evaluation \(ptweedie\(\)\) is 2% more accurate intweedieDistrthan intweedieon average, when compared to Monte Carlo evaluations\. Moreover, thetweedieimplementation of the inverse cumulative distribution \(qtweedie\(\)\) sometimesS: returns and error whenρ\\rhoapproaches 2\. We observed such failures in the upper row of the third panel in Fig\.[5](https://arxiv.org/html/2609.28607#S3.F5); they are more frequent in cells closer to the top left corner, that is, whenϕ\\phigrows larger\.D: nella fig 3 però abbiamo dei numeri, vuol dire che i numeri in quella cella sono calcolati su un numero di esperimenti piu’ piccolo?S: siHowever, the scope of thetweediepackage is wider, as it implements the functions above for anyρ\>0\\rho\>0, whiletweedieDistrcurrently handles onlyρ∈\(1,2\)\\rho\\in\(1,2\)\.

## 4The TWEES model

TWEES is an exponential smoothing model with a Tweedie predictive distribution\. A smoothed level controls the mean of a Tweedie likelihood, following a simple exponential smoothing recursion as in[Snyder et al\. \(2012\)](https://arxiv.org/html/2609.28607#bib.bib44)\. Thus, denoting the data of time series asy1,…,yTy\_\{1\},\\dots,y\_\{T\}, TWEES modelsYi∼Tw⁡\(μi,ϕi,ρ\)Y\_\{i\}\\sim\\mathrm\{Tw\}\(\\mu\_\{i\},\\phi\_\{i\},\\rho\), whit the recursion

μi=αμ​yi−1\+θμ​μ¯\+\(1−αμ−θμ\)​μi−1\\mu\_\{i\}=\\alpha\_\{\\mu\}y\_\{i\-1\}\+\\theta\_\{\\mu\}\\bar\{\\mu\}\+\(1\-\\alpha\_\{\\mu\}\-\\theta\_\{\\mu\}\)\\mu\_\{i\-1\}whereαμ\\alpha\_\{\\mu\}andθμ\\theta\_\{\\mu\}are learnable smoothing parameters, andμ¯\\bar\{\\mu\}is the average of the training data\. This model is also similar to the TweedieGP model of[Damato et al\. \(2025\)](https://arxiv.org/html/2609.28607#bib.bib12), where the mean of the Tweedie distribution is modeled with a Gaussian process\. Recall that the variance ofYiY\_\{i\}is tied to its mean through Eq\. \([1](https://arxiv.org/html/2609.28607#S3.E1)\); here, however, the dispersionϕi\\phi\_\{i\}is not kept fixed, but derived from a separate occurrence process, described next\.

To let occurrence and demand size evolve independently, we run a second damped exponential smoothing on the occurrence indicatorot=𝟙yt\>0o\_\{t\}=\\mathbbm\{1\}\_\{y\_\{t\}\>0\}, as in the hurdle\-shifted Poisson model of[Snyder et al\. \(2012\)](https://arxiv.org/html/2609.28607#bib.bib44):

πi=απ​oi−1\+θπ​π¯\+\(1−απ−θπ\)​πi−1,\\pi\_\{i\}=\\alpha\_\{\\pi\}o\_\{i\-1\}\+\\theta\_\{\\pi\}\\bar\{\\pi\}\+\(1\-\\alpha\_\{\\pi\}\-\\theta\_\{\\pi\}\)\\pi\_\{i\-1\},whereαπ\\alpha\_\{\\pi\}andθπ\\theta\_\{\\pi\}are smoothing parameters too, andπ¯\\bar\{\\pi\}is the average of the occurrence process\. Then, we constrainϕi\\phi\_\{i\}so that the mass at zero of the Tweedie matches the occurrence forecastπi\\pi\_\{i\}: since the zero\-mass of a Tweedie variable ise−λie^\{\-\\lambda\_\{i\}\}withλ\\lambdaas in Eq\. \([3](https://arxiv.org/html/2609.28607#S3.E3)\), settinge−λi=1−πie^\{\-\\lambda\_\{i\}\}=1\-\\pi\_\{i\}and solving forϕi\\phi\_\{i\}gives

ϕi=μi2−ρ\(2−ρ\)​\(−log⁡\(1−πi\)\)\.\\phi\_\{i\}=\\frac\{\\mu\_\{i\}^\{2\-\\rho\}\}\{\(2\-\\rho\)\\left\(\-\\log\(1\-\\pi\_\{i\}\)\\right\)\}\.
The parameters\(ρ,μ0,αμ,θμ,π0,απ,θπ\)\(\\rho,\\mu\_\{0\},\\alpha\_\{\\mu\},\\theta\_\{\\mu\},\\pi\_\{0\},\\alpha\_\{\\pi\},\\theta\_\{\\pi\}\)are estimated jointly by maximising the Tweedie log\-likelihood implied by the recursions above, subject toθη≥0,αη≥0\\theta\_\{\\eta\}\\geq 0,\\alpha\_\{\\eta\}\\geq 0, and the usual smoothing constraintsαη\+θη<1\\alpha\_\{\\eta\}\+\\theta\_\{\\eta\}<1forη∈\{μ,π\}\\eta\\in\\\{\\mu,\\pi\\\}\.

One\-step\-ahead forecasts are obtained directly from the fitted model, given the observed values; longer\-horizon forecasts are generated by simulating both exponential smoothing processes forward in time\.

As in[Damato et al\. \(2025\)](https://arxiv.org/html/2609.28607#bib.bib12), the series is rescaled by its median positive value before fitting; this operation is allowed as the distribution is absolutely continuous fory\>0y\>0\. In TWEES we restrict the power parameter toρ∈\(1\.2,1\.8\)\\rho\\in\(1\.2,1\.8\): the upper bound limits the number of terms required to evaluate the density, while the lower bound ensures numerical stability\.

## 5Experiments

### 5\.1Setup

We evaluate the models introduced in Sec\.[2\.1](https://arxiv.org/html/2609.28607#S2.SS1)and[4](https://arxiv.org/html/2609.28607#S4)on the datasets introduced in Sec\.[2\.2](https://arxiv.org/html/2609.28607#S2.SS2)\. We usefable\.intermittent\(v\. 0\.3\.0\) andtweedieDistr\(v\. 0\.2\.0\); our code is available in this anonymous GitHub repository444[https://anonymous\.4open\.science/r/benchmarking\_intermittent\-DD78](https://anonymous.4open.science/r/benchmarking_intermittent-DD78); note that in this repository, the name of our packages is no longer anonymous\.\. All models are run with default options; for those requiring sample paths,10510^\{5\}are drawn\.

For each dataset, we perform forecast evaluation using an expanding window approach where we use the firstTTobservations as a training set and the followinghhobservations as a holdout test set\. In particular, we use two windows settingT=L−j​hT=L\-jhforj=1,2j=1,2\. The value ofhhis set to the dataset’s native forecast horizon:h=6h=6for Auto,h=12h=12for RAF, andh=28h=28for Pasta and TinyM5\.

Writing the training set byy1,…,yTy\_\{1\},\\allowbreak\\dots,\\allowbreak y\_\{T\}and the holdout test set byyT\+1,…,yT\+hy\_\{T\+1\},\\dots,\\allowbreak y\_\{T\+h\}, we use the Root Mean Squared Scaled Error\([Hyndman and Koehler, 2006](https://arxiv.org/html/2609.28607#bib.bib25)\)to evaluate point forecasts:

RMSSE=1h​∑t=T\+1T\+h\(yt−y^t\)21T−1​∑t=2T\(yt−yt−1\)2\.\\mathrm\{RMSSE\}=\\sqrt\{\\dfrac\{\\dfrac\{1\}\{h\}\\sum\_\{t=T\+1\}^\{T\+h\}\\left\(y\_\{t\}\-\\hat\{y\}\_\{t\}\\right\)^\{2\}\}\{\\dfrac\{1\}\{T\-1\}\\sum\_\{t=2\}^\{T\}\\left\(y\_\{t\}\-y\_\{t\-1\}\\right\)^\{2\}\}\}\.wherey^t\\hat\{y\}\_\{t\}denotes the mean of the forecast distribution at timett\. We assess the predictive distribution using the quantile score\([Gneiting and Raftery, 2007](https://arxiv.org/html/2609.28607#bib.bib20)\)

QSτ​\(y,y^τ\)=\{τ⁡\(y−y^τ\)if​y≥y^τ,\(1−τ\)​\(y^τ−y\)if​y<y^τ\.\\mathrm\{QS\}\_\{\\tau\}\(y,\\hat\{y\}\_\{\\tau\}\)=\\begin\{cases\}\\tau\\,\(y\-\\hat\{y\}\_\{\\tau\}\)&\\text\{if \}y\\geq\\hat\{y\}\_\{\\tau\},\\\\ \(1\-\\tau\)\\,\(\\hat\{y\}\_\{\\tau\}\-y\)&\\text\{if \}y<\\hat\{y\}\_\{\\tau\}\.\\end\{cases\}Following[Spiliotis et al\. \(2021\)](https://arxiv.org/html/2609.28607#bib.bib45), we evaluate the quantile levelsτ=0\.5,0\.75,0\.835,0\.975,0\.995\\tau=0\.5,0\.75,\\allowbreak 0\.835,\\allowbreak 0\.975,0\.995, and we scale the quantile scores by the Mean Absolute Error \(MAE\) of the naive forecast\. We omit lower\-tail levels \(τ=0\.005,0\.025,0\.165,0\.25\\tau=0\.005,0\.025,\\allowbreak 0\.165,\\allowbreak 0\.25\), as in the data sets considered their quantile forecast is zero for most time series and models\. The resulting scaled quantile scoresQSτ\\mathrm\{sQS\}\_\{\\tau\}is defined as follows:

sQSτ=∑t=T\+1T\+hQSτ​\(yt,y^t,τ\)∑t=2T\|yt−yt−1\|,\\mathrm\{sQS\}\_\{\\tau\}=\\frac\{\\sum\_\{t=T\+1\}^\{T\+h\}\\mathrm\{QS\}\_\{\\tau\}\(y\_\{t\},\\hat\{y\}\_\{t,\\tau\}\)\}\{\\sum\_\{t=2\}^\{T\}\|y\_\{t\}\-y\_\{t\-1\}\|\},wherey^t,τ\\hat\{y\}\_\{t,\\tau\}is the quantile forecast of levelτ\\tauat timestamptt\. We note thatsQS0\.5\\mathrm\{sQS\}\_\{0\.5\}corresponds to the Mean Absolute Scaled Error\([Hyndman and Koehler, 2006](https://arxiv.org/html/2609.28607#bib.bib25)\)\.

The scoring rules we introduced are evaluated over the test set of each window: subsequently, we compute a score for each time series by taking the mean over the two windows\. Since the scoring rules are scale\-independent, we evaluate each method on each dataset by taking the mean across different time series; this practice has been shown to be more robust than the use of ranks with MCB test\([Koning et al\., 2005](https://arxiv.org/html/2609.28607#bib.bib30)\), especially when dealing with high quantile levels\([Corani et al\., 2026](https://arxiv.org/html/2609.28607#bib.bib10)\)\. To test the significance of our findings, we use, for each scoring rule, the Model Confidence Set \(MCS,[Hansen et al\., 2011](https://arxiv.org/html/2609.28607#bib.bib21); confidenceα=0\.05\\alpha=0\.05\), a procedure that identifies a set of statistically equivalent best methods by iteratively eliminating the worst one\.

### 5\.2Results

Table 3:Auto data set: RMSSE andsQSτ\\mathrm\{sQS\}\_\{\\tau\}\(τ=0\.5,0\.75,0\.835,0\.975,0\.995\\tau=0\.5,0\.75,0\.835,0\.975,0\.995\)\.Table 4:Pasta data set: RMSSE andsQSτ\\mathrm\{sQS\}\_\{\\tau\}\(τ=0\.5,0\.75,0\.835,0\.975,0\.995\\tau=0\.5,0\.75,0\.835,0\.975,0\.995\)\.Table 5:RAF data set: RMSSE andsQSτ\\mathrm\{sQS\}\_\{\\tau\}\(τ=0\.5,0\.75,0\.835,0\.975,0\.995\\tau=0\.5,0\.75,0\.835,0\.975,0\.995\)\.Table 6:TinyM5 data set: RMSSE andsQSτ\\mathrm\{sQS\}\_\{\\tau\}\(τ=0\.5,0\.75,0\.835,0\.975,0\.995\\tau=0\.5,0\.75,0\.835,0\.975,0\.995\)\.We report the average scores of each method in Tab\.[3](https://arxiv.org/html/2609.28607#S5.T3)\(Auto\), Tab\.[4](https://arxiv.org/html/2609.28607#S5.T4)\(Pasta\), Tab\.[5](https://arxiv.org/html/2609.28607#S5.T5)\(RAF\), and Tab\.[6](https://arxiv.org/html/2609.28607#S5.T6)\(TinyM5\)\. For each metric, we report in bold the methods with the lowest score \(the best\), we underline the scores of those methods falling into the MCS, and we report the rank of each method in grey\.

Across all four datasets, no single method dominates uniformly, but two patterns can be observed\. First: static, simple benchmarks like EMPSD and PARAMSD are hard to beat on RMSSE and on central quantiles\. As already noted by[Boylan and Babai \(2022\)](https://arxiv.org/html/2609.28607#bib.bib8), static methods are particularly effective when the historical data is limited \(Auto, Tab\.[3](https://arxiv.org/html/2609.28607#S5.T3)\) or very sparse \(RAF, Tab\.[5](https://arxiv.org/html/2609.28607#S5.T5)\)\. The difference between the two methods is often minimal, except on the Auto dataset, where the parametric model is better than the empirical distribution on all metrics\.

Second, we observe that two methods, TWEES and NEGBINES, have better scores and almost never fall outside the top five\. For instance, on the TinyM5 dataset, Tab\.[6](https://arxiv.org/html/2609.28607#S5.T6), the two methods are always part of the MCS\. TWEES, in particular, ranks in the top three on RMSSE andsQS0\.835\\mathrm\{sQS\}\_\{0\.835\}across all datasets\.

Our results show that exponential smoothing models generally perform better\. For TWEES, in particular, the two\-process approach for occurrence and demand size is key: using a single process that only controls the mean, as in NEGBINES, led to poorer results in preliminary experiments\. The negative binomial method, however, remains very competitive, and almost always ranks better than TWEES on Pasta, Tab\.[4](https://arxiv.org/html/2609.28607#S5.T4)\. HSPES is the worst of the three ES methods: indeed, its modeling of the occurrence is effective, but the demand process is fit on comparatively few non\-zero observations and the light right tail of its Poisson likelihood struggles to cover the highest quantiles\.

A related problem affects the GAMPOISB model: the Gamma distribution driving the mean of the Poisson often has low variance, constraining the model to behave similarly to a Poisson model, thus affecting its performance on high quantile levels\. The other Bayesian model, BETANBB, is frequently effective: its results are particularly good on the extreme quantiles for Auto, and in the central ones for all data sets, confirming the suitability of the negative binomial distribution on intermittent demand\.

Bootstrapping methods sometimes rank in top positions, but they lack consistency across different metrics\. VZ, for instance, performs well on high quantiles only when time series are long \(Pasta, TinyM5\)\. Its mean forecasts are accurate only on the Auto dataset, the one with the least sparse demand\. Its main limitation is that, using Croston’s decomposition, the number of values to sample from is equal to the number of non\-zero values in the time series\. Similarly, the performance of WSS is unstable, probably because bootstrap demand values, after the jittering step, are then rounded to positive integers: this latter step biases the quantile forecasts of the demand size upwards\. Instead in RAF, where the lowest quantiles are only driven by the occurrence Markov Chain, WSS has good scores\.

ARMA\-based models are probably the class with the worst performance\. NNARMA and MARWAL sometimes perform well on RMSSE, with the latter ranking second on RAF\. Yet on probabilistic scores the two methods are at the bottom of the ranking, particularly at the central quantiles where they rank 8th to 11th on every dataset\. This gap between a competitive mean forecast and poor\-performing quantile forecasts can be traced back to the use of a Gaussian forecast distribution, whose symmetric tails fail to model frequent zeros and sporadic, possibly large, demand spikes\.

## 6A worked example of intermittent time series reconciliation

Figure 6:Hierarchy of the Pasta dataset\.The integration offable\.intermittentin thefableframework automatically gives access to time series reconciliation tools\. Within thefableframework, in particular,fable\.bayesReconprovides methods for the reconciliation of probabilistic forecasts of smooth, intermittent and mixed time series hierarchies\. Here, we show how to produce coherent probabilistic forecasts\([Athanasopoulos et al\., 2024](https://arxiv.org/html/2609.28607#bib.bib3)\)for a hierarchy of time series\. We use the pasta dataset, which has a natural hierarchical structure, with 118 bottom\-level time series \(daily sales of pasta products\) aggregated into 4 time series representing the sales of each brand, further aggregated into total sales\. The hierarchy is represented in Fig\.[6](https://arxiv.org/html/2609.28607#S6.F6)\. The code snippet in Listing[7](https://arxiv.org/html/2609.28607#S6.F7)shows the code for loading the data and creating the hierarchical structure and the train/test split\. The experimental setup is the same as the one used in the previous section\.

library\(fable\)

library\(fable\.intermittent\)

library\(dplyr\)

pasta\_hts<\-pasta\|\>

aggregate\_key\(brand/product,value=sum\(value\)\)

n\_test<\-28

train<\-pasta\_hts\|\>filter\(index<=max\(index\)\-n\_test\)

test<\-pasta\_hts\|\>filter\(index\>max\(index\)\-n\_test\)

Figure 7:Loading the data and computing the aggregated series withaggregate\_key\.The bottom time series are mostly non\-smooth, as shown in Fig\.[1](https://arxiv.org/html/2609.28607#S2.F1); therefore we fit a TWEES model\. On the other hand, the aggregated time series are smooth, so we fit an ETS model\([Hyndman and Athanasopoulos, 2021](https://arxiv.org/html/2609.28607#bib.bib24), Ch\. 8\)\. We reconcile these forecasts with the BUIS\([Zambon et al\., 2024a](https://arxiv.org/html/2609.28607#bib.bib56)\), MixCond and TDcond\([Zambon et al\., 2024b](https://arxiv.org/html/2609.28607#bib.bib57)\)algorithms, implemented by the functionsbayesRecon\_BUIS\(\),bayesRecon\_mixCond\(\)andbayesRecon\_TDcond\(\)of the packagefable\.bayesRecon\. As a comparison, we fit an ETS model to all the series of the hierarchy, which we then reconcile with MinT\([Wickramasuriya et al\., 2019](https://arxiv.org/html/2609.28607#bib.bib54)\)\.

fit<\-train\_hts\|\>

model\(ets=ETS\(value\),twees=TWEES\(value\)\)\|\>

mutate\(base\_mix=if\_else\(is\_aggregated\(product\),ets,twees\)\)

fit\_rec<\-fit\|\>reconcile\(BUIS=bayesRecon\_BUIS\(base\_mix\),

MixCond=bayesRecon\_MixCond\(base\_mix\),

TDcond=bayesRecon\_TDcond\(base\_mix\),

MinT=min\_trace\(ets\)\)

fc<\-fit\_rec\|\>forecast\(h=n\_test\)

Figure 8:Fit and reconcile mixed forecasts\.Table 7:RMSSE andsQSτ\\mathrm\{sQS\}\_\{\\tau\}computed on the hierarchy of the Pasta dataset \(Fig\.[6](https://arxiv.org/html/2609.28607#S6.F6)\)\. Base forecasts \(grey background\) are incoherent; mixed are computed with TWEES on the bottom series and ETS on the aggregated series\.Listing[8](https://arxiv.org/html/2609.28607#S6.F8)shows that we can easily fit all models and reconcile their forecasts with very few lines of code and an intuitive syntax\. Tab\.[7](https://arxiv.org/html/2609.28607#S6.T7)shows theRMSSE\\mathrm\{RMSSE\}andsQSτ\\mathrm\{sQS\}\_\{\\tau\}\(τ∈\{0\.5,0\.75,0\.835,0\.975,0\.995\}\\tau\\in\\\{0\.5,0\.75,0\.835,0\.975,0\.995\\\}\) for the base forecasts \(ETS and mixed\) and the reconciled forecasts \(BUIS, MixCond, TDcond and MinT\)\. Among the base forecasts, mixed is consistently better than ETS, as TWEES is more effective than ETS on the intermittent bottom series\. Moreover, reconciliation via BUIS and MixCond improves the base forecasts on all the scores except for the two highest quantiles, while TDcond and MinT are less effective\. Indeed, BUIS and MixCond are tailored for hierarchies with count\-valued bottom forecasts; conversely, MinT relies on Gaussian base forecasts, which often are not effective on intermittent time series\.

## 7Conclusion

We presentfable\.intermittent, anRpackage that implements a broad set of probabilistic methods for intermittent time series within thefableframework\. Its purpose is to bring together, under one interface, methods that were previously scattered across different packages and programming languages, or had no publicly available implementation at all\. We also introduce TWEES, a new exponential smoothing model with a Tweedie predictive distribution\. Finally, we provide an efficient implementation of the Tweedie distribution in the separateRpackagetweedieDistr\. Since the package adopts the tidy syntax of thetidyvertsecosystem, every method is fitted, forecast, evaluated and visualized through the same concise pipeline, and can be applied to large collections of series with a single call\. All the models are compatible with the rest of the ecosystem: they can be used together with other fable forecasting models, and their forecasts reconciled across a hierarchy using functions available in the same ecosystem\.

Several directions remain open for future work\. The static distributions performed strongly and could be refined\. For example, the models might account for seasonality or give more weight to recent observations\. PARAMSD, which currently selects among candidate distributions with a log\-likelihood criterion, could use a quantile score, since decisions typically rely on the upper tail of the predictive distribution\. Another direction is to extend some of the models to incorporate exogenous variables, such as promotions, which strongly affect demand\. Finally, we stress that our implementations are open\-source and ready to be extended by users; we encourage contributions: the guide on the GitHub repository of the package555[https://github\.com/StefanoDamato/fable\.intermittent/commit/ee95d66947a49740dd0d0bac50ade369c6a572e3](https://github.com/StefanoDamato/fable.intermittent/commit/ee95d66947a49740dd0d0bac50ade369c6a572e3)details how to build new models and integrate them in the package\.

## References

- Alexandrov et al\. \(2020\)Alexandrov, A\., Benidis, K\., Bohlke\-Schneider, M\., Flunkert, V\., Gasthaus, J\., Januschowski, T\., Maddix, D\.C\., Rangapuram, S\., Salinas, D\., Schulz, J\., Stella, L\., Türkmen, A\.C\., Wang, Y\., 2020\.GluonTS: Probabilistic and Neural Time Series Modeling in Python\.Journal of Machine Learning Research 21, 1–6\.
- Almeida and Vieira \(2024\)Almeida, G\., Vieira, M\., 2024\.fableCount: INGARCH and GLARMA Models for Count Time Series in Fable Framework\.R package version 0\.1\.0\.
- Athanasopoulos et al\. \(2024\)Athanasopoulos, G\., Hyndman, R\.J\., Kourentzes, N\., Panagiotelis, A\., 2024\.Forecast reconciliation: A review\.International Journal of Forecasting 40, 430–456\.
- Azzimonti et al\. \(2026\)Azzimonti, D\., Damato, S\., Zambon, L\., Carrara, C\., Corani, G\., 2026\.fable\.bayesRecon: Bayesian Reconciliation in the ’fable’ Framework\.R package version 0\.2\.0\.
- Babai et al\. \(2021\)Babai, M\., Chen, H\., Syntetos, A\., Lengu, D\., 2021\.A compound\-Poisson Bayesian approach for spare parts inventory forecasting\.International Journal of Production Economics 232, 107954\.
- Box and Jenkins \(1970\)Box, G\.E\.P\., Jenkins, G\.M\., 1970\.Time Series Analysis: Forecasting and Control\.Holden\-Day, San Francisco\.
- Boylan and Syntetos \(2021\)Boylan, J\., Syntetos, A\., 2021\.Intermittent Demand Forecasting: Context, methods and applications\.1 ed\., Wiley\.
- Boylan and Babai \(2022\)Boylan, J\.E\., Babai, M\.Z\., 2022\.Estimating the cumulative distribution function of lead\-time demand using bootstrapping with and without replacement\.International Journal of Production Economics 252, 108586\.
- Brent \(1971\)Brent, R\.P\., 1971\.An algorithm with guaranteed convergence for finding a zero of a function\.The Computer Journal 14, 422–425\.
- Corani et al\. \(2026\)Corani, G\., Damato, S\., Azzimonti, D\., Zambon, L\., 2026\.Model selection with proper scoring rules on data sets of time series: prefer the mean scaled score\.[arXiv:2606\.24715](http://arxiv.org/abs/2606.24715)\.
- Croston \(1972\)Croston, J\.D\., 1972\.Forecasting and stock control for intermittent demands\.Operational Research Quarterly \(1970\-1977\) 23, 289\.
- Damato et al\. \(2025\)Damato, S\., Azzimonti, D\., Corani, G\., 2025\.Forecasting intermittent time series with Gaussian Processes and Tweedie likelihood\.International Journal of Forecasting , S0169207025000974\.
- Damato et al\. \(2026\)Damato, S\., Rubattu, N\., Azzimonti, D\., Corani, G\., 2026\.Intermittent time series forecasting: local vs global models\.[arXiv:2601\.14031](http://arxiv.org/abs/2601.14031)\.
- Dunn \(2026\)Dunn, P\.K\., 2026\.Tweedie: Evaluation of Tweedie Exponential Family Models\.R package version 3\.1\.20\.
- Dunn and Smyth \(2005\)Dunn, P\.K\., Smyth, G\.K\., 2005\.Series evaluation of Tweedie exponential dispersion model densities\.Statistics and Computing 15, 267–280\.
- Dunn and Smyth \(2008\)Dunn, P\.K\., Smyth, G\.K\., 2008\.Evaluation of Tweedie exponential dispersion model densities by Fourier inversion\.Statistics and Computing 18, 73–86\.
- Eddelbuettel and Francois \(2011\)Eddelbuettel, D\., Francois, R\., 2011\.Rcpp: Seamless R and C\+\+ Integration\.Journal of Statistical Software 40, 1–18\.
- Feng \(2021\)Feng, C\.X\., 2021\.A comparison of zero\-inflated and hurdle models for modeling zero\-inflated count data\.Journal of Statistical distributions and applications 8, 8\.
- Fildes et al\. \(2022\)Fildes, R\., Ma, S\., Kolassa, S\., 2022\.Retail forecasting: Research and practice\.International Journal of Forecasting 38, 1283–1318\.Special Issue: M5 competition\.
- Gneiting and Raftery \(2007\)Gneiting, T\., Raftery, A\.E\., 2007\.Strictly proper scoring rules, prediction, and estimation\.Journal of the American Statistical Association 102, 359–378\.
- Hansen et al\. \(2011\)Hansen, P\.R\., Lunde, A\., Nason, J\.M\., 2011\.The model confidence set\.Econometrica 79, 453–497\.
- Harvey and Fernandes \(1989\)Harvey, A\.C\., Fernandes, C\., 1989\.Time series models for count or qualitative observations\.Journal of Business & Economic Statistics 7, 407–417\.
- Hasni et al\. \(2019\)Hasni, M\., Babai, M\., Aguir, M\., Jemai, Z\., 2019\.An investigation on bootstrapping forecasting methods for intermittent demands\.International Journal of Production Economics 209, 20–29\.
- Hyndman and Athanasopoulos \(2021\)Hyndman, R\.J\., Athanasopoulos, G\., 2021\.Forecasting: principles and practice\.3rd ed\., OTexts, Melbourne, Australia\.URL:[https://OTexts\.com/fpp3](https://otexts.com/fpp3)\. accessed on October 31, 2025\.
- Hyndman and Koehler \(2006\)Hyndman, R\.J\., Koehler, A\.B\., 2006\.Another look at measures of forecast accuracy\.International Journal of Forecasting 22, 679–688\.
- Januschowski et al\. \(2020\)Januschowski, T\., Gasthaus, J\., Wang, Y\., Salinas, D\., Flunkert, V\., Bohlke\-Schneider, M\., Callot, L\., 2020\.Criteria for classifying forecasting methods\.International Journal of Forecasting 36, 167–177\.
- Joachimiak \(2022\)Joachimiak, K\., 2022\.m5: ’M5 Forecasting’ Challenges Data\.R package version 0\.1\.1\.
- Jørgensen \(1987\)Jørgensen, B\., 1987\.Exponential dispersion models\.Journal of the Royal Statistical Society Series B\-Methodological 49, 127–145\.
- Kolassa \(2016\)Kolassa, S\., 2016\.Evaluating predictive count data distributions in retail sales forecasting\.International Journal of Forecasting 32, 788–803\.
- Koning et al\. \(2005\)Koning, A\.J\., Franses, P\.H\., Hibon, M\., Stekler, H\., 2005\.The M3 competition: Statistical tests of the results\.International Journal of Forecasting 21, 397–409\.
- Lang et al\. \(2026\)Lang, Y\., Mao, W\., He, J\., Deng, X\., Liu, W\., Zuo, M\., 2026\.A comprehensive review on intermittent time series forecasting from the perspective of machine learning\.Neurocomputing , 134783\.
- Makridakis et al\. \(2022\)Makridakis, S\., Spiliotis, E\., Assimakopoulos, V\., 2022\.M5 accuracy competition: Results, findings, and conclusions\.International Journal of Forecasting 38, 1346–1364\.
- Mancuso et al\. \(2021\)Mancuso, P\., Piccialli, V\., Sudoso, A\.M\., 2021\.A machine learning approach for forecasting hierarchical time series\.Expert Systems with Applications 182, 115102\.
- Mersmann \(2024\)Mersmann, O\., 2024\.microbenchmark: Accurate Timing Functions\.R package version 1\.5\.0\.
- O’Hara\-Wild et al\. \(2026a\)O’Hara\-Wild, M\., Hyndman, R\., Wang, E\., 2026a\.fable: Forecasting Models for Tidy Time Series\.R package version 0\.5\.0\.
- O’Hara\-Wild et al\. \(2026b\)O’Hara\-Wild, M\., Hyndman, R\., Wang, E\., 2026b\.fabletools: Core Tools for Packages in the ’fable’ Framework\.R package version 0\.8\.0\.
- O’Hara\-Wild et al\. \(2026c\)O’Hara\-Wild, M\., Kay, M\., Hayes, A\., Hyndman, R\., 2026c\.distributional: Vectorised Probability Distributions\.R package version 0\.7\.1\.
- Prak and Teunter \(2019\)Prak, D\., Teunter, R\., 2019\.A general method for addressing forecasting uncertainty in inventory models\.International Journal of Forecasting 35, 224–238\.Special Section: Supply Chain Forecasting\.
- R Core Team \(2026\)R Core Team, 2026\.R: A Language and Environment for Statistical Computing\.R Foundation for Statistical Computing\. Vienna, Austria\.
- Salinas et al\. \(2020\)Salinas, D\., Flunkert, V\., Gasthaus, J\., Januschowski, T\., 2020\.DeepAR: Probabilistic forecasting with autoregressive recurrent networks\.International Journal of Forecasting 36, 1181–1191\.
- Sbrana \(2025\)Sbrana, G\., 2025\.Markov Walk and Walmart sales prediction\.Journal of the Operational Research Society , 1–12\.
- Sbrana and Babai \(2026\)Sbrana, G\., Babai, M\.Z\., 2026\.Non\-negative autoregressive moving average models for intermittent demand: Forecast accuracy and inventory implications\.European Journal of Operational Research , S0377221726005412\.
- Schwarz \(1978\)Schwarz, G\., 1978\.Estimating the dimension of a model\.The Annals of Statistics 6, 461–464\.
- Snyder et al\. \(2012\)Snyder, R\.D\., Ord, J\.K\., Beaumont, A\., 2012\.Forecasting the intermittent demand for slow\-moving inventories: A modelling approach\.International Journal of Forecasting 28, 485–496\.
- Spiliotis et al\. \(2021\)Spiliotis, E\., Makridakis, S\., Kaltsounis, A\., Assimakopoulos, V\., 2021\.Product sales probabilistic forecasting: An empirical evaluation using the M5 competition data\.International Journal of Production Economics 240, 108237\.
- Svetunkov \(2026\)Svetunkov, I\., 2026\.smooth: Forecasting Using State Space Models\.URL:[https://CRAN\.R\-project\.org/package=smooth](https://cran.r-project.org/package=smooth), doi:[10\.32614/CRAN\.package\.smooth](http://dx.doi.org/10.32614/CRAN.package.smooth)\. r package version 4\.5\.0\.
- Svetunkov and Boylan \(2023\)Svetunkov, I\., Boylan, J\.E\., 2023\.iETS: State space model for intermittent demand forecasting\.International Journal of Production Economics 265, 109013\.
- Syntetos et al\. \(2016\)Syntetos, A\.A\., Babai, Z\., Boylan, J\.E\., Kolassa, S\., Nikolopoulos, K\., 2016\.Supply chain forecasting: Theory, practice, their gap and the future\.European Journal of Operational Research 252, 1–26\.
- Syntetos and Boylan \(2005\)Syntetos, A\.A\., Boylan, J\.E\., 2005\.The accuracy of intermittent demand estimates\.International Journal of Forecasting 21, 303–314\.
- Syntetos et al\. \(2005\)Syntetos, A\.A\., Boylan, J\.E\., Croston, J\.D\., 2005\.On the categorization of demand patterns\.Journal of the Operational Research Society 56, 495–503\.
- Teunter et al\. \(2011\)Teunter, R\.H\., Syntetos, A\.A\., Babai, M\.Z\., 2011\.Intermittent demand: Linking forecasting to inventory obsolescence\.European Journal of Operational Research 214, 606–615\.
- Türkmen et al\. \(2021\)Türkmen, A\.C\., Januschowski, T\., Wang, Y\., Cemgil, A\.T\., 2021\.Forecasting intermittent and sparse time series: A unified probabilistic framework via deep renewal processes\.PLOS ONE 16, e0259764\.
- Wang et al\. \(2020\)Wang, E\., Cook, D\., Hyndman, R\.J\., 2020\.A new tidy data structure to support exploration and modeling of temporal data\.Journal of Computational and Graphical Statistics 29, 466–478\.
- Wickramasuriya et al\. \(2019\)Wickramasuriya, S\.L\., Athanasopoulos, G\., Hyndman, R\.J\., 2019\.Optimal forecast reconciliation for hierarchical and grouped time series through trace minimization\.Journal of the American Statistical Association 114, 804–819\.
- Willemain et al\. \(2004\)Willemain, T\.R\., Smart, C\.N\., Schwarz, H\.F\., 2004\.A new approach to forecasting intermittent demand for service parts inventories\.International Journal of Forecasting 20, 375–387\.
- Zambon et al\. \(2024a\)Zambon, L\., Azzimonti, D\., Corani, G\., 2024a\.Efficient probabilistic reconciliation of forecasts for real\-valued and count time series\.Statistics and Computing 34, 21\.
- Zambon et al\. \(2024b\)Zambon, L\., Azzimonti, D\., Rubattu, N\., Corani, G\., 2024b\.Probabilistic reconciliation of mixed\-type hierarchical time series, in: Proceedings of the Fortieth Conference on Uncertainty in Artificial Intelligence, PMLR\. p\. 4078–4095\.
- Zhou and Viswanathan \(2011\)Zhou, C\., Viswanathan, S\., 2011\.Comparison of a new bootstrapping method with parametric approaches for safety stock determination in service parts inventory systems\.International Journal of Production Economics 133, 481–485\.

## Appendix AMethods from the literature

We describe here the methods from the literature implemented infable\.intermittent\. We denote byy1,…,yTy\_\{1\},\\dots,y\_\{T\}the values of a time series with training set of lengthTT\. We usehhto denote the length of the forecast horizon\. Forecasts are generated for the variablesYT\+1,…,YT\+hY\_\{T\+1\},\\dots,Y\_\{T\+h\}, where we use capital letters to denote a random variable\.

Several of the methods rely on*Croston’s decomposition*\([Croston, 1972](https://arxiv.org/html/2609.28607#bib.bib11)\), which we review here\. Lett1<t2<⋯<tnt\_\{1\}<t\_\{2\}<\\dots<t\_\{n\}be thennperiods at which a positive value is observed, i\.e\.yti\>0y\_\{t\_\{i\}\}\>0\. The*demand sizes*are the positive values themselves,di=ytid\_\{i\}=y\_\{t\_\{i\}\}fori=1,…,ni=1,\\dots,n; the*demand intervals*ℓi=ti−ti−1\\ell\_\{i\}=t\_\{i\}\-t\_\{i\-1\}\(withℓ0=0\\ell\_\{0\}=0\) are the number of periods between consecutive positive demands\. Alternatively, the*occurrence*process is used, that is the binary time seriesot=𝟙\{yt\>0\}o\_\{t\}=\\mathbbm\{1\}\_\{\\\{y\_\{t\}\>0\\\}\}fort=1,…,Tt=1,\\dots,T\.

### A\.1Bootstrapping methods

#### WSS\([Willemain et al\., 2004](https://arxiv.org/html/2609.28607#bib.bib55)\)

This bootstrapping method uses a Markov chain with states\{0,1\}\\\{0,1\\\}to model the occurrence\. Denoting\{Oi\}\\\{O\_\{i\}\\\}the Bernoulli variables of the Markov chain, the transition probabilities are estimated using the training set as

ℙ⁡\(Oi\+1=j\|Oi=k\)=∑i=1T−1𝟙\{oi=k,oi\+1=j\}∑i=1T𝟙\{oi=k\}\.\\mathbb\{P\}\(O\_\{i\+1\}=j\|O\_\{i\}=k\)=\\frac\{\\sum\_\{i=1\}^\{T\-1\}\\mathbbm\{1\}\_\{\\\{o\_\{i\}=k,o\_\{i\+1\}=j\\\}\}\}\{\\sum\_\{i=1\}^\{T\}\\mathbbm\{1\}\_\{\\\{o\_\{i\}=k\\\}\}\}\.The positive part of the forecast distribution is obtained sampling demand sizesdid\_\{i\}and adding a Gaussian noise with variancedi\\sqrt\{d\_\{i\}\}; the demand size is then rounded and negative values are set to 0\.

Thus, the forecast distribution is a mixture of a mass in zero and a bootstrap distribution on positive integers\. Forecasts are generated via autoregressive sampling\.

#### VZ\([Zhou and Viswanathan, 2011](https://arxiv.org/html/2609.28607#bib.bib58)\)

This method is based on Croston’s decomposition\. Predictions are generated with independent bootstrap samples of demand intervals and demand size: the former determine when positive values will appear, the latter determine their magnitude\. Training times are almost immediate\.

### A\.2Bayesian models

#### BETANBB\([Harvey and Fernandes, 1989](https://arxiv.org/html/2609.28607#bib.bib22)\)

A Bayesian model using a negative binomial forecast distribution, that isYi∼NB⁡\(r,pi\)Y\_\{i\}\\sim\\mathrm\{NB\}\(r,p\_\{i\}\)\. The count parameterrris estimated globally across the time series, while the probability parameterpip\_\{i\}follows a Beta distribution,pi∼Beta⁡\(ai,bi\)p\_\{i\}\\sim\\mathrm\{Beta\}\(a\_\{i\},b\_\{i\}\), whose parameters are updated following the equations

ai\+1=ω⁡\(ai\+v\)\+\(1−ω\),bi\+1=ω⁡\(bi\+yi\)a\_\{i\+1\}=\\omega\(a\_\{i\}\+v\)\+\(1\-\\omega\),\\quad b\_\{i\+1\}=\\omega\(b\_\{i\}\+y\_\{i\}\)\(4\)withω∈\(0,1\]\\omega\\in\(0,1\]andv\>0v\>0\. The parametersrr,vv,ω\\omega,a1a\_\{1\},b1b\_\{1\}are estimated to minimize the negative log\-likelihood of the observations\. Probabilistic forecasts are generated via autoregressive sampling\.

#### GAMPOISB\([Harvey and Fernandes, 1989](https://arxiv.org/html/2609.28607#bib.bib22)\)

In this model, a hierarchical Gamma\-Poisson structure is used:Yi∼Pois⁡\(λi\)Y\_\{i\}\\sim\\mathrm\{Pois\}\(\\lambda\_\{i\}\)withλi∼Γ⁡\(ai,bi\)\\lambda\_\{i\}~\\sim\\Gamma\(a\_\{i\},b\_\{i\}\)\. Similarly to Eq\. \([4](https://arxiv.org/html/2609.28607#A1.E4)\), the parameters of the distribution vary as

ai\+1=ω​ai\+yi,bi\+1=ω​bi\+1a\_\{i\+1\}=\\omega a\_\{i\}\+y\_\{i\},\\quad b\_\{i\+1\}=\\omega b\_\{i\}\+1withω∈\(0,1\]\\omega\\in\(0,1\]\. Initial parameters, asa1a\_\{1\}andb1b\_\{1\}, and the updating coefficientω\\omegaare learned via an optimizer to minimize the negative log\-likelihood, which takes the form of a negative binomial distribution\. One step\-ahead forecasts are computed in closed\-form, while for longer horizons autoregressive sampling is needed\.

### A\.3Exponential smoothing models

#### HSPES\([Snyder et al\., 2012](https://arxiv.org/html/2609.28607#bib.bib44)\)

In this model, demand size and occurrence are forecast separately\. To both, a simple exponential smoothing\([Hyndman and Athanasopoulos, 2021](https://arxiv.org/html/2609.28607#bib.bib24), Ch\. 8\)is applied, such that respectively

πi=απ​oi−1\+θπ​π¯\+\(1−απ−θπ\)​πi−1\\pi\_\{i\}=\\alpha\_\{\\pi\}o\_\{i\-1\}\+\\theta\_\{\\pi\}\\bar\{\\pi\}\+\(1\-\\alpha\_\{\\pi\}\-\\theta\_\{\\pi\}\)\\pi\_\{i\-1\}and

λi=αλ​\(di−1−1\)\+θλ​λ¯\+\(1−αλ−θλ\)​λi−1\\lambda\_\{i\}=\\alpha\_\{\\lambda\}\(d\_\{i\-1\}\-1\)\+\\theta\_\{\\lambda\}\\bar\{\\lambda\}\+\(1\-\\alpha\_\{\\lambda\}\-\\theta\_\{\\lambda\}\)\\lambda\_\{i\-1\}where for bothη∈\{π,λ\}\\eta\\in\\\{\\pi,\\lambda\\\},αη≥0,θη≥0\\alpha\_\{\\eta\}\\geq 0,\\theta\_\{\\eta\}\\geq 0, andαη\+θη<1\\alpha\_\{\\eta\}\+\\theta\_\{\\eta\}<1\. The exponential smoothing is calledundampedordampeddepending on whetherθη\\theta\_\{\\eta\}is equal to 0 or not\. Note that the shifted demand size time seriesdi−1d\_\{i\}\-1is potentially shorter thany1,…,yTy\_\{1\},\\dots,y\_\{T\}and is only updated whenyi\>0y\_\{i\}\>0\.

The two processes are combined into a hurdle\-shifted Poisson distribution \(HSP\) such thatYi∼Oi​\(1\+Zi\)Y\_\{i\}\\sim O\_\{i\}\(1\+Z\_\{i\}\)whereOi∼Ber⁡\(πi\)O\_\{i\}\\sim\\mathrm\{Ber\}\(\\pi\_\{i\}\)andZi∼Pois⁡\(λi\)Z\_\{i\}\\sim\\mathrm\{Pois\}\(\\lambda\_\{i\}\)\. Thusℙ⁡\(Yi=0\)=1−πi\\mathbb\{P\}\(Y\_\{i\}=0\)=1\-\\pi\_\{i\}andℙ⁡\(Yi=k\)=πi​e−λi​λik−1\(k−1\)\!\\mathbb\{P\}\(Y\_\{i\}=k\)=\\pi\_\{i\}e^\{\-\\lambda\_\{i\}\}\\frac\{\\lambda\_\{i\}^\{k\-1\}\}\{\(k\-1\)\!\}\.

The parameters of the exponential smoothing process,αη\\alpha\_\{\\eta\}andθη\\theta\_\{\\eta\}, and the initial valuesη0\\eta\_\{0\}forη∈\{π,λ\}\\eta\\in\\\{\\pi,\\lambda\\\}are learned to minimise the negative log\-likelihood, whileπ¯\\bar\{\\pi\}andλ¯\\bar\{\\lambda\}are set as the sample mean of the occurrence and the shifted demand size respectively\.

#### NEGBINES\([Snyder et al\., 2012](https://arxiv.org/html/2609.28607#bib.bib44)\)

In this model, a damped or undamped simple exponential smoothing process determines the mean parameter of a negative binomial distribution, which is related to the parametrisation asμ=r​1−pp\\mu=r\\frac\{1\-p\}\{p\}\. Thus,Yi∼NB⁡\(μi​p1−p,p\)Y\_\{i\}\\sim\\mathrm\{NB\}\(\\mu\_\{i\}\\frac\{p\}\{1\-p\},p\), where the mean is computed with the recursion

μi=α​yi−1\+θ​μ¯\+\(1−α−θ\)​μi−1\.\\mu\_\{i\}=\\alpha y\_\{i\-1\}\+\\theta\\bar\{\\mu\}\+\(1\-\\alpha\-\\theta\)\\mu\_\{i\-1\}\.The parametersp∈\(0,1\)p\\in\(0,1\)andμ1\\mu\_\{1\}of the negative binomial distribution, andα≥0\\alpha\\geq 0andθ≥0\\theta\\geq 0such thatα\+θ<1\\alpha\+\\theta<1are determined using an optimiser\.μ¯\\bar\{\\mu\}is the mean of the training set of the time series\.

The forecast distribution is available in closed form only forh=1h=1; for longer horizons, it is obtained via autoregressive sampling\. Being the negative binomial overdispersed, this model can expand its right tail more than HSPES; however, it is necessarily unimodal\.

### A\.4Static methods

#### EMPSD\([Boylan and Babai, 2022](https://arxiv.org/html/2609.28607#bib.bib8)\)

The empirical distribution is a widely used and simple probabilistic benchmark for intermittent time series\([Spiliotis et al\., 2021](https://arxiv.org/html/2609.28607#bib.bib45)\)\. The predictive distribution isℙ\(Yi=k\)=1T∑i=1T𝟙\{yi=k\}\\mathbb\{P\}\(Y\_\{i\}=k\)=\\frac\{1\}\{T\}\\sum\_\{i=1\}^\{T\}\\mathbbm\{1\}\_\{\\\{y\_\{i\}=k\\\}\}\.

#### PARAMSD\([Kolassa, 2016](https://arxiv.org/html/2609.28607#bib.bib29)\)

Fitting static distributions has been shown to be a strong benchmark:[Boylan and Babai \(2022\)](https://arxiv.org/html/2609.28607#bib.bib8)show that, for slow\-moving items, a well\-chosen parametric distribution estimates the lead\-time demand CDF as accurately as bootstrap resampling, particularly when the demand history is short\. This model treats the series as i\.i\.d\. draws from a count distribution: the predictive distribution is identical at every horizon,YT\+j∼FΘ^Y\_\{T\+j\}\\sim F\_\{\\hat\{\\Theta\}\}for alljj, whereΘ^\\hat\{\\Theta\}is a set of parameters learned via maximum likelihood\. Five candidate distributions are fitted to the historyy1,…,yTy\_\{1\},\\dots,y\_\{T\}: Poisson \(Pois⁡\(λ\)\\mathrm\{Pois\}\(\\lambda\)\), negative binomial \(NB⁡\(r,p\)\\mathrm\{NB\}\(r,p\)\), their hurdle\-shifted counterparts\([Feng, 2021](https://arxiv.org/html/2609.28607#bib.bib18)\), and Tweedie\.

Hurdle shifted models use a Croston\-type decomposition into occurrenceot=𝟙\{yt\>0\}o\_\{t\}=\\mathbbm\{1\}\_\{\\\{y\_\{t\}\>0\\\}\}and positive demand sizes, placing an explicit massπ^0=1−1T​∑tot\\hat\{\\pi\}\_\{0\}=1\-\\tfrac\{1\}\{T\}\\sum\_\{t\}o\_\{t\}at zero and fitting the count distribution to the shifted sizes\{yt−1∣yt\>0\}\\\{y\_\{t\}\-1\\mid y\_\{t\}\>0\\\}:

ℙ⁡\(Y=0\)=π0,ℙ⁡\(Y=k\)=\(1−π0\)​pΘ∖\{π\}​\(k−1\),k≥1\.\\mathbb\{P\}\(Y=0\)=\\pi\_\{0\},\\qquad\\mathbb\{P\}\(Y=k\)=\(1\-\\pi\_\{0\}\)p\_\{\\Theta\\setminus\\\{\\pi\\\}\}\(k\-1\),\\quad k\\geq 1\.The choice of the distribution is driven by an information criterion \(either the defaultBIC=−2​ℓ​\(Θ^\)\+\|Θ\|​log⁡T\\mathrm\{BIC\}=\-2\\ell\(\\hat\{\\Theta\}\)\+\|\\Theta\|\\log T, orAIC=−2​ℓ​\(Θ^\)\+2​\|Θ\|\\mathrm\{AIC\}=\-2\\ell\(\\hat\{\\Theta\}\)\+2\|\\Theta\|, where\|Θ\|\|\\Theta\|is the number of parameters\), and the minimiser is returned\.

To allow for a fair comparison with other distributions, we do not evaluate the log\-density of the Tweedie distribution on count data, but discretise it by integrating it over the interval\[y−1/2,y\+1/2\)\[y\-1/2,y\+1/2\)for anyy∈ℕy\\in\\mathbb\{N\}\.

### A\.5ARMA methods

#### MARWAL\([Sbrana, 2025](https://arxiv.org/html/2609.28607#bib.bib41)\)

The model works as the interplay of an ARMA\(1,1\) process\([Box and Jenkins, 1970](https://arxiv.org/html/2609.28607#bib.bib6)\),ziz\_\{i\}, and a Markov chainχi\\chi\_\{i\}on the occurrence switching between 0 and 1, such that

yi=χi​ω\+ziy\_\{i\}=\\chi\_\{i\}\\omega\+z\_\{i\}wherezi=λ​zi−1\+νi\+θ​νi−1z\_\{i\}=\\lambda z\_\{i\-1\}\+\\nu\_\{i\}\+\\theta\\nu\_\{i\-1\}andω\\omegais the mean of non\-zero values of the training set\. The model achieves fast training times, as the parametersλ\\lambdaandθ\\thetaare estimated in closed form\.

Furthermore, a validation loop is performed to determine the starting valuet0∈\{1,…,T\}t\_\{0\}\\in\\\{1,\\dots,T\\\}of the training setyt0,…,yTy\_\{t\_\{0\}\},\\dots,y\_\{T\}\.

The mean forecasty^T\+h\\hat\{y\}\_\{T\+h\}and its varianceσ^T\+h2\\hat\{\\sigma\}^\{2\}\_\{T\+h\}is computed in closed form for any horizon, and a Gaussian distribution is assumed for the forecasts such that

YT\+h∼𝒩⁡\(y^T\+h,σ^T\+h2\)\.Y\_\{T\+h\}\\sim\\mathcal\{N\}\\left\(\\hat\{y\}\_\{T\+h\},\\hat\{\\sigma\}^\{2\}\_\{T\+h\}\\right\)\.Negative samples or quantiles, however, are set to zero by design\.

#### NNARMA\([Sbrana and Babai, 2026](https://arxiv.org/html/2609.28607#bib.bib42)\)

This model drives the demand series through a Gaussian ARMA\(1,1\) process, whose forecast distribution is then constrained to remain non\-negative\. Writingνi\\nu\_\{i\}for the one step\-ahead innovation, the model is

yi−ω=ϕ⁡\(yi−1−ω\)\+νi\+θ​νi−1,y\_\{i\}\-\\omega=\\phi\(y\_\{i\-1\}\-\\omega\)\+\\nu\_\{i\}\+\\theta\\nu\_\{i\-1\},withϕ,ω\>0\\phi,\\omega\>0,θ≤0\\theta\\leq 0andϕ\+θ\>0\\phi\+\\theta\>0\. The three parameters are estimated by least squares on the innovationsνi\\nu\_\{i\}\. The mean and variance of the forecast distribution are computed in closed form for any horizon via the standard ARMA\(1,1\) recursion\.

Input time series undergo a multiplicative deseasonalisation, and seasonal factors, computed as the average ratio between the time series and its moving average, are then applied to the forecasts\. Similarly to the Markov Walk model, the distribution has its negative part set to zero\. However, we underline that this leads to an ill\-posed distribution, as collapsing the negative mass to zero would also modify its mean\.

## Appendix BImplementative differences betweentweedieDistrandtweedie

#### dtweedie\(\)

IntweedieDistr, thanks to the fastC\+\+implementation, the density functiondtweedie\(\)can be evaluated via its series expansion\([Dunn and Smyth, 2005](https://arxiv.org/html/2609.28607#bib.bib15)\)for any choice of parametersϕ\\phiandρ∈\(1,2\)\\rho\\in\(1,2\)\(the meanμ\\mudoes not enter the series expansion\)\. On the other hand,tweedieuses the series expansion only in a few cases\. In the other cases, an interpolation strategy with saddlepoint approximation is used on tabulated data embedded in the package\.

#### ptweedie\(\)

For the evaluation of the cumulative distribution function on positive values, thetweediepackage uses the Fourier inversion strategy described by[Dunn and Smyth \(2008\)](https://arxiv.org/html/2609.28607#bib.bib16)\. For thetweedieDistrimplementation, on the other hand, we use the Compound\-Poisson parametrisation to write the cumulative distribution function as

ℙ⁡\(Y≤y\)=e−λ\+∑n=1\+∞Tn​\(y\)​with\\displaystyle\\mathbb\{P\}\(Y\\leq y\)=e^\{\-\\lambda\}\+\\sum\_\{n=1\}^\{\+\\infty\}T\_\{n\}\(y\)\\ \\text\{with\}Tn​\(y\)=pPois​\(n∣λ\)​FGamma​\(y∣n​α,β\),\\displaystyle T\_\{n\}\(y\)=p\_\{\\mathrm\{Pois\}\}\(n\\mid\\lambda\)F\_\{\\mathrm\{Gamma\}\}\(y\\mid n\\alpha,\\beta\),fory\>0y\>0, wherepPoisp\_\{\\mathrm\{Pois\}\}denotes the mass function of a Poisson distribution andFGammaF\_\{\\mathrm\{Gamma\}\}the cumulative distribution function of a Gamma distribution\. We then identify the index of the largest term

M=arg⁡maxn≥1​Tn​\(y\)M=\\arg\\max\_\{n\\geq 1\}T\_\{n\}\(y\)and we approximate the series∑n=1\+∞Tn​\(y\)\\sum\_\{n=1\}^\{\+\\infty\}T\_\{n\}\(y\)as

∑n∈ΛTn​\(y\)​with​Λ:=\{n∈ℕ:Tn​\(y\)TM​\(y\)≥e−37\}\.\\sum\_\{n\\in\\Lambda\}T\_\{n\}\(y\)\\ \\text\{with\}\\ \\Lambda:=\\left\\\{n\\in\\mathbb\{N\}:\\frac\{T\_\{n\}\(y\)\}\{T\_\{M\}\(y\)\}\\geq e^\{\-37\}\\right\\\}\.The threshold ofe−37≃8⋅10−17e^\{\-37\}\\simeq 8\\cdot 10^\{\-17\}has been chosen as it guarantees double precision in a 64\-bit floating point arithmetic\([Dunn and Smyth, 2005](https://arxiv.org/html/2609.28607#bib.bib15)\)\.

#### qtweedie\(\)

The quantile function also differs among the two packages\. The major reason for the speedup is that both packages exploit their implementation of cumulative distribution function, which is faster intweedieDistr\. However, thetweediepackage relies only on Brent’s method\([Brent, 1971](https://arxiv.org/html/2609.28607#bib.bib9)\), whiletweedieDistruses aC\+\+implementation of Newton\-Raphson algorithm, where the bisection method kicks in where convergence is not reached in few steps \(∼4%\\sim 4\\%of the times\)\.

#### rtweedie\(\)

Sampling from a Tweedie distribution forρ∈\(1,2\)\\rho\\in\(1,2\)easily leverages the expression in Eq\. \([2](https://arxiv.org/html/2609.28607#S3.E2)\): first, drawNNas aPois⁡\(λ\)\\mathrm\{Pois\}\(\\lambda\), then drawyyfromGamma⁡\(n​α,β\)\\mathrm\{Gamma\}\(n\\alpha,\\beta\)\. Both packages use this strategy\.

## Appendix CBenchmarking speed and accuracy of a Tweedie implementation

We benchmarktweedieDistr’sdtweedie\(\),ptweedie\(\)andqtweedie\(\)against their namesake functions in the referencetweediepackage, on identical inputs, over a grid of dispersionϕ∈\{0\.2,0\.6,…,7\.0\}\\phi\\in\\\{0\.2,0\.6,\\dots,7\.0\\\}\(18 values\), powerρ∈\{1\.1,1\.2,…,1\.9\}\\rho\\in\\\{1\.1,1\.2,\\dots,1\.9\\\}\(9 values\)\. For each\(ϕ,ρ\)\(\\phi,\\rho\)pair, the experiment is repeated 200 times to average out sampling variability; this number is up to 54% smaller for the combinations of parameters where the quantile function gives errors\. We do not iterate over different values ofμ\\mu, always set to 1, as it does not enter the series expansion of the density[Dunn and Smyth \(2005\)](https://arxiv.org/html/2609.28607#bib.bib15)\.

For the density and cumulative distribution function, for each repetition we draw a probability of a zero observationp0∼Unif⁡\(0\.2,0\.8\)p\_\{0\}\\sim\\mathrm\{Unif\}\(0\.2,0\.8\), then generate the test input of sizenn,x1,…,xnx\_\{1\},\\dots,x\_\{n\}as independent draws from a zero\-inflated Gamma distribution, that is,Xi​∼iid​O⋅ZX\_\{i\}\\overset\{\\mathrm\{iid\}\}\{\\sim\}O\\cdot Z, withO∼Ber⁡\(p0\)O\\sim\\mathrm\{Ber\}\(p\_\{0\}\)andZ∼Gamma⁡\(1,1\)Z\\sim\\mathrm\{Gamma\}\(1,1\); this is supposed to mimic the properties of Tweedie\-generated data\. For the quantile function, we generate the inputq1,…,qnq\_\{1\},\\dots,q\_\{n\}fromQi​∼iid​Unif​\(0,1\)Q\_\{i\}\\overset\{\\mathrm\{iid\}\}\{\\sim\}\\mathrm\{Unif\}\(0,1\)\. Both packages are evaluated on this same input\. Wall\-clock time is measured withmicrobenchmark\([Mersmann, 2024](https://arxiv.org/html/2609.28607#bib.bib34)\), running 10 replications per comparison; the reported speedup is the ratio of the reference package’s mean time totweedieDistr’s mean time\.

To account for the dependence of the computational speedup on the sample size, we run these experiments forn=2,10,50,200,500,2000n=2,10,50,200,500,2000; the results shown in Fig\.[5](https://arxiv.org/html/2609.28607#S3.F5)are generated usingn=200n=200fordtweedie\(\)andptweedie\(\), andn=10n=10forqtweedie\(\)\. We observe that the gap in the speedup increases with the sample size for the cumulative function and its inverse, with the new implementations being up to 1000 times faster than the original one\. The difference in speed between the two versions, on the contrary, decreases withnnfor the density; however, even for the largest sample size we use, 2000, our novel implementation remains from 1\.5 to 20 times faster than thetweediepackage implementation\.

The correctness offable\.intermittentis verified by an automated test suite that checks its output againsttweediewith a tolerance of10−810^\{\-8\}\. In the experiments grid described above, the differences in the density evaluation between the two packages were below the machine epsilon used in the experiments\. For the cumulative probability function and the quantile function, the median absolute difference is again smaller than the machine epsilon; but the maximum difference is respectively in the order of5×10−25\\times 10^\{\-2\}and7×10−47\\times 10^\{\-4\}\. In the few cases where there is a disagreement among the two packages, Monte Carlo simulations confirmedfable\.intermittentto be more accurate\.