Intermittent demand forecasting with tsintermittent

Nikolaos Kourentzes

Introduction

Intermittent demand time series record items that are requested infrequently, so that many periods have zero demand. Spare parts are the typical example, but the same pattern appears in many slow-moving items. Conventional forecasting methods, such as exponential smoothing, do not account for the zeros explicitly, and Croston (1972) proposed an alternative that has since been extended in several ways.

The tsintermittent package provides Croston’s method and its variants, the TSB method, simple exponential smoothing, Croston-style moving averages, the iMAPA temporal aggregation approach, and tools to classify intermittent series. A central feature of the package is that the parameters and initial values of all methods can be optimised for each time series, using cost functions designed for intermittent demand (Kourentzes, 2014). The package also provides empirical predictive distributions over the lead time.

This vignette explains the logic behind these choices and shows how to use the main arguments. The example series ts.data1 is a monthly intermittent demand series with 84 observations.

plot(ts.data1, ylab = "Demand")

Forecasts as demand rates

Croston’s method separates an intermittent series into two series: the non-zero demand sizes \(z_t\) and the intervals between them \(x_t\). Each is smoothed with single exponential smoothing, but only in periods where demand occurs:

\[\hat{z}_t = \hat{z}_{t-1} + \alpha_z (z_t - \hat{z}_{t-1}), \qquad \hat{x}_t = \hat{x}_{t-1} + \alpha_x (x_t - \hat{x}_{t-1}),\]

and the forecast is their ratio, \(\hat{y}_t = \hat{z}_t / \hat{x}_t\). In periods without demand the estimates stay unchanged. The forecast is the same for all future periods.

It is important to interpret this forecast correctly. A forecast of 0.5 does not say that half a unit will be demanded next period. It says that, on average, one unit is demanded every two periods: it is a demand rate. The forecast is therefore not trying to track the demand in each period, which is zero in most periods and positive in a few, but the average demand per period. This view drives both how the package optimises the methods and how it provides predictive distributions.

Methods

Croston and its variants are provided by crost(). Croston’s original method is biased, and two corrections are available through the argument type: the Syntetos-Boylan approximation ("sba"), which multiplies the forecast by \((1 - \alpha_x/2)\), and the Shale-Boylan-Johnston correction ("sbj"), which multiplies it by \((1 - \alpha_x/(2 - \alpha_x))\).

fit <- crost(ts.data1, type = "sba", outplot = TRUE)

fit$weights   # Smoothing parameters for demand and intervals
#> [1] 0.29283632 0.01067471
fit$initial   # Initial values for demand and intervals
#> [1] 5.937983 9.123422

TSB, provided by tsb(), replaces the intervals with the probability of demand \(d_t\), which is 1 in periods with demand and 0 otherwise. The demand probability is updated in every period, while the demand size is updated only when demand occurs, and the forecast is \(\hat{y}_t = \hat{d}_t \hat{z}_t\). Because the probability is updated every period, TSB reacts to a series that stops being demanded, which Croston’s method cannot do. The first smoothing parameter is for the demand size and the second for the demand probability, and the optimiser keeps the latter no larger than the former, following Teunter et al. (2011).

fit.tsb <- tsb(ts.data1, outplot = TRUE)

fit.tsb$weights
#> [1] 3.262541e-01 1.939481e-09

Simple exponential smoothing, provided by sexsm(), is a common benchmark. It makes no special treatment of zero demand.

Croston-style moving averages, provided by crost.ma(), replace exponential smoothing with moving averages of the demand sizes and intervals. The orders of the moving averages play the role of the smoothing parameters.

iMAPA, provided by imapa(), models the series at multiple temporal aggregation levels and combines the resulting forecasts. It is described in its own section below.

sexsm(ts.data1)$alpha
#> [1] 0.01837332
crost.ma(ts.data1)$order
#> [1] 8 7

Optimising the methods

What is optimised

By default all methods optimise both their smoothing parameters and their initial values. The relevant arguments are:

Optimising the initial values matters for intermittent series. Such series are often short and have few non-zero observations, and the smoothing parameters are typically low, so the initial values can affect the forecasts for a long time. Optimising them avoids relying on ad-hoc values.

# Optimise everything (default)
crost(ts.data1, type = "sba")$frc.out[1]
#> [1] 0.1908284

# Single parameter for demand and intervals
crost(ts.data1, type = "sba", nop = 1)$weights
#> [1] 0.1530906 0.1530906

# Fixed parameters, optimised initial values
crost(ts.data1, type = "sba", w = c(0.1, 0.05))$initial
#> [1] 0.3922488 1.2072842

# Everything fixed
crost(ts.data1, type = "sba", w = c(0.1, 0.05), init = c(5, 3), init.opt = FALSE)$frc.out[1]
#> [1] 0.5478768

Cost functions

Conventional exponential smoothing is optimised by minimising the squared or absolute one-step-ahead errors. For intermittent demand this is problematic. Most periods have zero demand, so an optimiser driven by period-by-period errors is pulled towards forecasts close to zero (Teunter and Duncan, 2009; Wallström and Segerstedt, 2010). More fundamentally, such errors compare a demand rate with the realised demand of each period, which are not of the same nature.

The package instead offers cost functions that compare like with like. The mean absolute rate (MAR) and mean squared rate (MSR) compare the forecast demand rate with the observed demand rate up to each period, measured by the cumulative mean of the series (Kourentzes, 2014):

\[r_t = \hat{y}_t - \frac{1}{t}\sum_{j=1}^{t} y_j, \qquad \text{MAR} = \sum_{t} |r_t|, \qquad \text{MSR} = \sum_{t} r_t^2.\]

The cumulative mean is unreliable in the first periods, where it is based on few observations, so it is winsorised over the first 30% of the sample. Because the comparison is on rates, the timing of individual demands does not affect the cost, while both the size and the frequency of demand do.

Two further cost functions measure the error over the forecast horizon \(h\). The cumulative error from origin \(t\) is

\[e_t^{(h)} = \sum_{j=1}^{h} y_{t+j} - \sum_{j=1}^{h} \hat{y}_{t+j|t},\]

and the mean squared cumulative error (MSCE) and mean absolute cumulative error (MACE) are the mean of \((e_t^{(h)})^2\) and of \(|e_t^{(h)}|\) across all available origins. These are natural choices when the forecasts are used to decide stock over a lead time, as discussed in the next section, and they use the same cumulative errors that define the predictive distribution. All methods in the package produce a constant forecast over the horizon, so the cumulative forecast from origin \(t\) is simply \(h\) times the forecast made at \(t\).

Unlike the other cost functions, MSCE and MACE depend on the forecast horizon: changing h changes the optimised parameters. This is intended, but it is worth understanding what it implies. Svetunkov, Kourentzes, and Killick (2024) show that cost functions based on multiple steps ahead errors shrink the parameters of exponential smoothing models towards zero. The reason is that the variance of the multi-step errors grows with the smoothing parameters, so minimising those errors pulls the parameters down, more strongly the longer the horizon. Their analysis is for linear state space models, but the methods in this package are built on exponential smoothing and the same mechanism applies: a noisier demand rate forecast inflates the cumulative error over \(h\) periods, so the optimiser favours smaller smoothing parameters as \(h\) grows.

Some shrinkage is useful, as it makes the forecasts less reactive to noise, which suits intermittent demand. However, Svetunkov, Kourentzes, and Killick (2024) also show that the shrinkage can become excessive when the sample is short and the horizon is long. Intermittent demand series are often short, so MSCE and MACE should be used with lead times that are small relative to the length of the series.

The conventional mean squared error (MSE) and mean absolute error (MAE) are also available for comparison. MAE is not recommended for intermittent demand, as it is particularly prone to favouring forecasts close to zero.

The argument cost selects the cost function, with MAR as the default:

sapply(c("mar", "msr", "mse", "mae"),
       function(x){crost(ts.data1, type = "sba", cost = x)$weights})
#>             mar       msr          mse       mae
#> [1,] 0.29283632 0.1456198 6.890404e-01 0.9999934
#> [2,] 0.01067471 0.1140640 5.751319e-08 0.9271863
crost(ts.data1, h = 3, type = "sba", cost = "msce")$weights
#> [1] 0.8041929 0.4680887
crost(ts.data1, h = 3, type = "sba", cost = "mace")$weights
#> [1] 0.9993288 0.9526520

Predictive distributions over the lead time

Why over the lead time

Since the forecasts are demand rates, a predictive distribution for the demand in a single period is not meaningful: the forecast does not say when demand will occur, only how much is demanded on average per period. What is meaningful, and what matters for stock control, is the total demand over the lead time. Over a lead time of \(h\) periods the forecasts sum to the expected cumulative demand, which can be compared directly with the cumulative realised demand.

The package therefore provides the predictive distribution of the cumulative demand over the forecast horizon, and when quantiles are requested h is understood as the lead time. More precisely, h should be the period the stock decision must cover. With continuous review this is the lead time. With periodic review it is the lead time plus the review period, often called the protection interval. For example, with monthly data, a lead time of two months, and monthly reviews, h = 3.

How it is calculated

The distribution is empirical. For every in-sample origin \(t\) the cumulative errors \(e_t^{(h)}\) over the next \(h\) periods are calculated, as in the definition above, using the forecasts made at origin \(t\). The quantiles of these errors are added to the cumulative forecast from the end of the sample:

\[Q_\tau = \sum_{j=1}^{h} \hat{y}_{n+j|n} + q_\tau\left(e^{(h)}\right),\]

where \(q_\tau\) is the empirical quantile at probability \(\tau\). Since demand cannot be negative, quantiles are truncated at zero.

Working with the cumulative errors directly matters. The forecast errors at different steps ahead are correlated, even when the model is correct, so the variance of the error over the lead time is not the one-step-ahead variance multiplied by the lead time, nor the sum of the variances at each step. Both shortcuts ignore the covariances between the errors and typically understate the uncertainty. The empirical distribution of the cumulative errors contains these covariances by construction, without assuming a model for the demand (Saoud, Kourentzes, and Boylan, 2022). The package takes this one step further and reads the quantiles directly from that distribution, so no distribution needs to be assumed for the errors either.

Any probability can be requested, but quantiles below the median are of little use here. They give a negative safety stock, that is, stock below the expected demand over the lead time, and they describe the lower tail of the cumulative demand, which is not what stock control is concerned with.

The two parts of \(Q_\tau\) map directly onto stock control. \(Q_\tau\) is the stock needed to cover the demand over the lead time with probability \(\tau\): the order-up-to level under periodic review, or the reorder point under continuous review. The error quantile \(q_\tau(e^{(h)})\) is the safety stock, the stock held on top of the expected demand over the lead time. Both are reported.

Quantiles are calculated only when the argument probs is given, for example probs = c(0.9, 0.95). The output then contains the cumulative forecast, the quantiles, the safety stock, and the cumulative errors used to calculate them.

fit <- crost(ts.data1, h = 3, type = "sba", probs = c(0.9, 0.95))
fit$lt$frc        # Cumulative forecast over the lead time
#> [1] 0.5724852
fit$lt$quantile   # Quantiles of the cumulative demand over the lead time
#>      90%      95% 
#> 1.418639 2.028401
fit$lt$ss         # Safety stock for each probability
#>       90%       95% 
#> 0.8461539 1.4559154
summary(fit)
#> Croston's method (SBA)
#> 
#> Series: 84 periods, 9 with demand.
#> Mean inter-demand interval (p): 8.556, CV^2 of non-zero demand: 0.7318, PKa class: SBA
#> 
#> Optimised with cost function MAR.
#> Smoothing parameters: demand 0.2928, interval 0.01067 (optimised)
#> Initial values: demand 5.938, interval 9.123 (optimised)
#> In-sample MAR: 2.561
#> 
#> Forecast for h = 3: 0.1908 per period
#> Demand over the lead time of 3 periods (81 cumulative errors):
#>                        90%    95%
#> Cumulative forecast 0.5725 0.5725
#> Quantile            1.4186 2.0284
#> Safety stock        0.8462 1.4559

The empirical quantiles are only as reliable as the number of errors behind them. A series of 24 periods with a lead time of 5 provides fewer than 20 overlapping cumulative errors, and an empirical 95% quantile from so few is close to the maximum observed error. The package warns when the sample is too small for the requested probabilities. This is the price of the empirical approach: it needs more data than approximations built on the one-step-ahead errors (Saoud, Kourentzes, and Boylan, 2022).

Plotting the forecasts

When quantiles are available there are two ways to view the forecasts.

The per-period view shows the series, the in-sample and out-of-sample demand rate, and each lead-time quantile divided by the lead time, spread over the forecast horizon. This is a convenient way to place the quantile next to the forecast, but it is not a per-period quantile.

The lead-time view shows the rolling cumulative demand over windows of \(h\) periods, the corresponding cumulative forecasts made at the start of each window, and the cumulative forecast and quantiles for the next \(h\) periods. This is the view in which the distribution is defined.

Either view can be drawn when the model is fitted, through outplot. outplot = TRUE draws the per-period view, and outplot = "period" or outplot = "leadtime" selects a view explicitly. The same views are available afterwards with plot(), which also accepts graphical arguments such as main or ylim:

fit <- crost(ts.data1, h = 3, type = "sba", probs = c(0.9, 0.95), outplot = "leadtime")

plot(fit)                     # Per-period view

plot(fit, view = "leadtime")  # Lead-time view

plot(fit, view = "leadtime", main = "Demand over a 3-month lead time")

Reusing fitted models

A fitted model can produce forecasts and lead time quantiles for any horizon without optimising it again. forecast() gives the point forecasts, the expected demand per period, for any h. It deliberately gives no prediction intervals, since a per-period interval would be misleading for demand rates. The quantiles of the demand over the lead time, for any lead time and probabilities, are given by lt.dist(), which uses the series and in-sample forecasts stored in the fitted model:

fit <- crost(ts.data1, type = "sba")
forecast(fit, h = 6)
#>            Jan       Feb       Mar       Apr       May       Jun
#> 1987 0.1908284 0.1908284 0.1908284 0.1908284 0.1908284 0.1908284
lt.dist(fit, h = 3, probs = c(0.9, 0.95))$quantile
#>      90%      95% 
#> 1.418639 2.028401

The parameters of a model optimised with MSCE or MACE depend on the horizon used when fitting it, and they are not reoptimised for a new h.

A fitted model can also be applied to new data, for example when new observations become available, with the input model of the forecasting functions, or model.fit for iMAPA. The smoothing parameters and the initial values of the model are used without optimisation. This is appropriate when the new data extend the series the model was fitted to, as both start at the same period. To use the smoothing parameters of a model on a different series, set init.opt = TRUE, so that only the initial values are optimised.

fit <- crost(window(ts.data1, end = c(1985, 12)), type = "sba")
crost(ts.data1, model = fit)$frc.out[1]
#> [1] 0.2052833
forecast(fit, h = 3, newdata = ts.data1)
#>            Jan       Feb       Mar
#> 1987 0.2052833 0.2052833 0.2052833
crost(ts.data2, model = fit, init.opt = TRUE)$initial
#> [1] 3.580031 9.250680

Classifying series

The function idclass() classifies series according to their average inter-demand interval \(p\) and the squared coefficient of variation of their non-zero demands \(CV^2\). Several schemes are available through type: the Syntetos-Boylan-Croston classification ("SBC"), the Kostenko-Hyndman classification ("KH" and its approximation "KHa"), and the Petropoulos-Kourentzes classification ("PK" and its approximation "PKa"). The latter extends the Kostenko-Hyndman classification for series without any zero demand, \(p = 1\), which are assigned to SES instead of Croston or SBA (Petropoulos and Kourentzes, 2015). Such series are rare at the original frequency, but common once a series is temporally aggregated, which is why this classification is used by iMAPA. The exact schemes optimise SBA for each series, so they are slower.

set.seed(1)
dataset <- cbind(simID(20, 60, idi = 1.05, cv2 = 0.1),
                 simID(20, 60, idi = 1.20, cv2 = 0.1),
                 simID(20, 60, idi = 2.00, cv2 = 0.8))
colnames(dataset) <- paste0("ts", 1:60)
cls <- idclass(dataset, type = "PKa", outplot = "summary")

cls$summary
#>         Series
#> Croston     33
#> SBA         26
#> SES          1
cls.sbc <- idclass(dataset, type = "SBC", outplot = "detail")

Temporal aggregation and iMAPA

Temporal aggregation sums an intermittent series into non-overlapping buckets of several periods, for example from monthly to quarterly demand. At higher aggregation levels there are fewer zeros, the intervals between demands shrink, and the variability of the demand sizes changes, so the series becomes easier to forecast. A single aggregation level must be chosen, though, and the best level differs between series and is not known in advance.

iMAPA avoids this choice by using many aggregation levels at once (Petropoulos and Kourentzes, 2015). For each aggregation level from minimumAL to maximumAL it:

  1. aggregates the series, dropping the first observations where needed so that all buckets are complete;
  2. classifies the aggregated series with the PKa scheme and fits Croston, SBA, or SES accordingly, so that the method can change with the aggregation level;
  3. forecasts the aggregated series and divides the forecast equally over the periods of each bucket, returning it to the original frequency.

The forecasts from all levels are then combined, with the mean or the median as set by comb. Aggregation levels with fewer than four observations, or fewer than two non-zero demands, are skipped. Because the methods selected at the higher levels see a series with less intermittency, this combines model selection with the benefits of aggregation without requiring the user to choose either. The maximum aggregation level still matters: very high levels leave few observations and over-smooth the series. By default it is half the length of the series.

The argument outplot gives four views of the results, numbered as in earlier versions of the package, plus the lead time view described above. All are also available with plot():

fit <- imapa(ts.data1, outplot = 1)   # Series and combined forecast

plot(fit, view = "levels")            # Forecasts from each aggregation level

plot(fit, view = "selection")         # Summary of the methods selected

plot(fit, view = "selection.detail")  # Classification per level

summary(fit)                          # Characteristics and method per level
#> iMAPA
#> 
#> Aggregation levels:
#>      AL  n     p     cv2 model        a1        a2 initial.z initial.x
#> AL1   1 84 8.556 0.73183   SBA 2.928e-01 1.067e-02     5.938     9.123
#> AL2   2 42 4.875 0.91349   SBA 1.673e-01 1.979e-01     4.606     5.352
#> AL3   3 28 3.250 0.91349   SBA 1.848e-01 1.951e-01     5.361     4.354
#> AL4   4 21 2.857 0.78316   SBA 1.805e-01 2.056e-01     5.844     3.670
#> AL5   5 16 2.143 0.14000   SBA 1.402e-01 1.510e-01     1.259     1.939
#> AL6   6 14 1.857 1.06574   SBA 1.985e-01 2.015e-01     6.695     2.998
#> AL7   7 12 1.833 0.86782   SBA 2.275e-01 2.369e-01     7.124     2.510
#> AL8   8 10 1.667 0.24000   SBA 2.074e-01 2.712e-01     1.143     1.002
#> AL9   9  9 1.500 0.24000   SBA 2.651e-01 2.059e-01     1.214     1.089
#> AL10 10  8 1.333 0.09600   SBA 1.103e-01 2.608e-09     1.178     1.178
#> AL11 11  7 1.750 0.31276   SBA 2.112e-01 1.636e-01     1.880     1.727
#> AL12 12  7 1.167 0.86782   SBA 1.994e-01 2.240e-01     7.507     2.000
#> AL13 13  6 1.200 0.21605   SBA 2.382e-05 2.328e-01     2.220     1.682
#> AL14 14  6 1.000 0.86782   SES 2.387e-01        NA     4.059        NA
#> AL15 15  5 1.000 0.06173   SES 3.943e-01        NA     1.500        NA
#> AL16 16  5 1.000 0.37500   SES 2.291e-01        NA     1.500        NA
#> AL17 17  4 1.000 0.31276   SES 3.135e-01        NA     1.500        NA
#> AL18 18  4 1.000 0.31276   SES 3.135e-01        NA     1.500        NA
#> 
#> Forecasts combined with the mean.
#> Forecast for h = 10: 0.1532 per period

The fitted models can be reused with model.fit, which skips the classification and optimisation:

imapa(ts.data1, model.fit = fit$model.fit)$frc.out
#>  [1] 0.1531747 0.1531747 0.1531747 0.1531747 0.1531747 0.1531747 0.1531747
#>  [8] 0.1531747 0.1531747 0.1531747

Lead time quantiles are available for iMAPA as well. The in-sample forecasts of each aggregation level are made at the start of each bucket rather than at every period, so the cumulative errors are calculated from forecasts made at every in-sample origin, using the fitted model of each level with fixed parameters. This takes longer than for the other methods, and paral can be used to speed it up. The cumulative cost functions MSCE and MACE are not offered for iMAPA, as temporal aggregation already targets the demand over multiple periods.

fit <- imapa(ts.data1, h = 3, maximumAL = 12, probs = c(0.9, 0.95))
fit$lt$quantile
#>      90%      95% 
#> 1.336003 1.951991
plot(fit, view = "leadtime")

Forecasting many series

data.frc() applies any of the methods to all columns of a data frame or matrix and collects the forecasts. Additional arguments are passed to the chosen method. With method = "auto", each series is forecast with Croston, SBA, or SES, according to its PKa classification.

out <- data.frc(dataset, method = "crost", type = "sba", h = 6)
out$frc.out[, 1:3]
#>        ts1      ts2      ts3
#> 1 34.53657 84.18071 11.61606
#> 2 34.53657 84.18071 11.61606
#> 3 34.53657 84.18071 11.61606
#> 4 34.53657 84.18071 11.61606
#> 5 34.53657 84.18071 11.61606
#> 6 34.53657 84.18071 11.61606

When probs is provided, the quantiles of the demand over the lead time of all series are collected in lt.out. Series that cannot be classified with method = "auto", because they have fewer than two non-zero values, are forecast with SES.

out <- data.frc(dataset, method = "auto", h = 3, probs = c(0.9, 0.95))
out$lt.out[, 1:3]
#>          ts1      ts2      ts3
#> 90% 126.2882 349.9759 43.30442
#> 95% 130.4752 366.7076 47.94016

References

Croston, J. D. (1972). Forecasting and stock control for intermittent demands. Operational Research Quarterly, 23(3), 289–303.

Kourentzes, N. (2014). On intermittent demand model optimisation and selection. International Journal of Production Economics, 156, 180–190. doi:10.1016/j.ijpe.2014.06.007

Petropoulos, F., and Kourentzes, N. (2015). Forecast combinations for intermittent demand. Journal of the Operational Research Society, 66(6), 914–924. doi:10.1057/jors.2014.62

Svetunkov, I., Kourentzes, N., and Killick, R. (2024). Multi-step estimators and shrinkage effect in time series models. Computational Statistics, 39(3), 1203–1239. doi:10.1007/s00180-023-01377-x

Saoud, P., Kourentzes, N., and Boylan, J. E. (2022). Approximations for the lead time variance: A forecasting and inventory evaluation. Omega, 110, 102614. doi:10.1016/j.omega.2022.102614

Teunter, R. H., and Duncan, L. (2009). Forecasting intermittent demand: a comparative study. Journal of the Operational Research Society, 60, 321–329.

Teunter, R. H., Syntetos, A. A., and Babai, M. Z. (2011). Intermittent demand: linking forecasting to inventory obsolescence. European Journal of Operational Research, 214(3), 606–615.

Wallström, P., and Segerstedt, A. (2010). Evaluation of forecasting error measurements and techniques for intermittent demand. International Journal of Production Economics, 128(2), 625–636.