Who

This talk presents results and software developed with the help of the whole group

Dario Azzimonti (Maintainer bayesRecon and fable.bayesRecon)

Stefano Damato (Maintainer fable.intermittent)

Lorenzo Zambon

Giorgio Corani

Special thanks also to Chiara Carrara and Anubhab Biswas.

What are intermittent time series?

Intermittent time series: contain mostly integer-values with many zero values.

Many more formal definitions (e.g. ADI\(>1.32\), Syntetos et al. (2005)).

Examples of intermittent time series

Often arise when modelling sales of SKU in demand forecasting
(Notable examples: M5, VN2, RAF or Auto parts datasets)

Today’s running example: Pasta dataset (Mancuso et al. 2021)

print(pasta_tsibble, n=3)

# A tsibble: 215,350 x 4 [1D]
# Key:       Brand, SKU [118]
  date       Brand   SKU value
  <date>     <chr> <int> <dbl>
1 2014-01-02 B1        1     7
2 2014-01-03 B1        1     5
3 2014-01-04 B1        1     9
# ℹ 215,347 more rows

Dataset in tsibble format will be avaialble on the next version of fable.intermittent.

How do we forecast intermittent time series?

• Point forecasts are not reasonable: decision making requires uncertainty

How do we forecast intermittent time series?

• Point forecasts are not reasonable: decision making requires uncertainty

• We need to return probabilistic forecasts.

How do we forecast intermittent time series?

• Point forecasts are not reasonable: decision making requires uncertainty

• We need probabilistic forecasts and custom models for intermittent time series

Package fable.intermittent

10 probabilistic intermittent forecasting models

TWEES = TWeedie Exponential Smoothing

ETS state-space model, Gaussian observation noise replaced by a Tweedie distribution

\[y_t \sim \text{Tweedie}(\mu_t, \phi, \rho), \quad \text{where}\]

• \(\mu_t\) follows a (optionally damped) exponential smoothing recursion, as in classical ETS.

• Dispersion \(\phi\) and power \(\rho\) (smoothing parameters) estimated jointly by maximum likelihood.

Tweedie distribution

• Parameters: mean \(\mu>0\), dispersion \(\phi>0\), power \(\rho \in(1,2)\)

• Point mass at zero and continuous density for \(y>0\)

• Density has no closed form but it is computable (Dunn and Smyth (2005), Damato et al. (2025))

Tweedie implementation

Fast implementation: dtweedie/ptweedie/qtweedie/rtweedie compiled in C++

Simple implementation in fable.intermittent

  # Fit ETS and TWEES models on training time series
  fit <- train |>
    model(ets = ETS(value), twees = TWEES(value))
  
  # Forecast one month (28 days) ahead
  fc <- fit |> forecast(h = 28)

  # Plot forecast and time series of brand B2 and SKU 19
  fc |> filter(Brand == "B2", SKU == 19) |>
    autoplot(train |> filter(Brand == "B2", SKU == 19, level = 90) +
    autolayer(test |> filter(Brand == "B2", SKU == 19), value)

Benchmark on Pasta dataset

Experimental setup: 118 time series, training on 1797 days

Models: all 10 available on fable.intermittent plus ETS

# Fit all models on the training data
fit_inter <- train |>
  model(
    ets = ETS(value),
    empdistr = EMPDISTR(value),
    staticdistr = STATICDISTR(value),
    betanbb = BETANBB(value),
    gampois = GAMPOISB(value),
    wss = WSS(value),
    vz = VZ(value),
    hspes = HSPES(value),
    negbines = NEGBINES(value),
    twees = TWEES(value),
    marwal = MARWAL(value)
  )

Benchmark on Pasta dataset

# Forecast the next 6 months
fc_inter <- fit_inter |> 
  forecast(h = 28)

acc_inter <- fc_inter |>
  accuracy(test_ts, 
  measures = list(
      SPL95=scaled_pinball_loss,
      CRPS=CRPS),
  .train = list(train_ts$value), 
  level=c(95))

Pasta dataset: a hierarchy

118 SKU divided between 4 brands summing up to a total

pasta_hts <- pasta_tsibble |>
  aggregate_key(Brand / SKU, value = sum(value))

How to reconcile?

We cannot use Gaussian reconciliation (MinT): bottom are intermittent.

We need mixed type reconciliation.

Reconciliation via conditioning: bayesRecon

• Base forecasts are incoherent with independent predictive distributions for each series.

• Reconciliation via conditioning takes the joint incoherent distribution and conditions on the hierarchy constraints (\(\boldsymbol u = \boldsymbol{Ab}\)), keeping only coherent points.

• Works for discrete, continuous and mixed-type predictive distributions

Implemented in R package bayesRecon: no forecast model fit, the user provides forecasts.

fable.bayesRecon

Implements reconciliation via conditioning into fable

• Gaussian with unknown covariance (bayesRecon_t)
  Bayesian approach to account for the uncertainty of the covariance matrix. Student’s t reconciled distribution (Carrara et al. 2026)

• Discrete forecasts (bayesRecon_BUIS)
  Bottom Up Importance Sampling: hierarchical constraints incorporated as sampling constraints (Zambon, Azzimonti, and Corani 2024)

• Mixed type forecasts (bayesRecon_MixCond, bayesRecon_TDcond)
  Gaussian forecasts reconciled with MinT, discrete forecasts reconciled with bottom-up (MixCond) or top-down (TDcond) updates.(Zambon, Azzimonti, Rubattu, et al. 2024)

Reconciling the Pasta hierarchy

Hierarchy: 118 disaggregated time series, 4 mid-level aggregates, one top level.

# Fit base forecasting models
# base_mix -> bottom TWEES, rest ETS
fit <- train_hts |>
  model(ets = ETS(value), twees = TWEES(value)) |>
  mutate(base_mix = if_else(is_aggregated(SKU), ets, twees))

# reconcile forecasting models
# mixCond on base_mix
# MinT    on ETS
fit_rec <- fit |> reconcile(mixCond = bayesRecon_MixCond(base_mix),
                            MinT    = min_trace(ets))

# Forecast 28 steps ahead
n_test <- 28
fc <- fit_rec |> forecast(h = n_test)

Results of reconciled forecasts

Intermittent time series in fable framework

fable.intermittent and fable.bayesRecon

• 10 probabilistic methods for intermittent time series implemented in fable environment

• Easy benchmarking and new methods development

• Tweedie distribution and TWEES model competitive with state of the art

• All works seamlessly with reconciliation functions in fable.bayesRecon

Help needed! Add/propose new methods by posting an issue on Github

https://github.com/StefanoDamato/fable.intermittent/issues

Thanks for your attention

Models implemented

Tweedie distribution

Member of the exponential dispersion family, parametrized by mean \(\mu>0\), dispersion \(\phi>0\), power \(p\in(1,2)\): \[\mathbb{E}[X]=\mu, \qquad \text{Var}(X) = \phi\,\mu^{p}\]

For power \(p \in (1,2)\) Tweedie has a compound Poisson-Gamma representation with

• Point mass at zero: \[\mathbb{P}(X=0) = \exp\!\left(-\frac{\mu^{2-p}}{\phi(2-p)}\right)\]

• For \(X>0\): density has no closed form but computable (Dunn and Smyth 2005)

Captures both “no demand” periods and a heavy right tail of spikes.

TWEES details

TWEES(formula, damped = TRUE, scaling = TRUE)

• damped = TRUE (default): damped trend in the smoothing recursion;

• scaling = TRUE (default): series rescaled by its max before fitting, for numerical stability

• \(h=1\) forecast: exact Tweedie distribution; \(h>1\): obtained by simulating the recursion forward (no closed form, as in ETS prediction intervals)

• Brings the Tweedie-likelihood idea (introduced in Damato et al. (2025) with Gaussian Processes) into a lightweight exponential smoothing model

fable.intermittent::ptweedie vs tweedie::ptweedie

fable.intermittent::dtweedie vs tweedie::dtweedie