EconBase
← Back to paper

Feature-based intermittent demand forecast combinations: bias, accuracy and inventory implications

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

61,929 characters · 15 sections · 106 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Feature-based intermittent demand forecast combinations: accuracy and inventory implications

abstractIntermittent demand forecasting is a ubiquitous and challenging problem in production systems and supply chain management. In recent years, there has been a growing focus on developing forecasting approaches for intermittent demand from academic and practical perspectives. However, limited attention has been given to forecast combination methods, which have achieved competitive performance in forecasting fast-moving time series. The current study aims to examine the empirical outcomes of some existing forecast combination methods and propose a generalized feature-based framework for intermittent demand forecasting. The proposed framework has been shown to improve the accuracy of point and quantile forecasts based on two real data sets. Further, some analysis of features, forecasting pools and computational efficiency is also provided. The findings indicate the intelligibility and flexibility of the proposed approach in intermittent demand forecasting and offer insights regarding inventory decisions.
keywordsIntermittent demand forecasting; Forecast combinations; Time series features; Diversity; Empirical evaluation

\setcounter{page}{1}

Introduction

Intermittent demand with several periods of zero demand is ubiquitous in practice. Over half of inventory consists of spare parts, in which demand is typically intermittent nikolopoulos2021we. Given the high purchase and shortage costs associated with intermittent demand applications, accurate forecasts could be coupled with improved inventory management in the field of manufacturing jiang2020intermittent, aerospace wang2016select, retailing sillanpaa2018forecasting and so on balugani2019periodic,babai2019new.

What makes intermittent demand challenging to forecast is that there are two sources of uncertainty: the sporadic demand occurrence, and the demand arrival timing. Seminal work on intermittent demand forecasting by croston1972forecasting proposed to forecast the sizes of demand and the inter-demand intervals separately. Then some scholars followed this idea and put forward some developments. For example, Syntetos-Boylan Approximation (SBA) proposed by syntetos2005accuracy delivered approximately unbiased estimates and constituted the benchmark in subsequently proposed methodologies for intermittent demand forecasting.

syntetos2005categorization proposed a categorization of demand patterns to facilitate the selection of croston1972forecasting's method and SBA syntetos2005accuracy. A classification rule was expressed in terms of the average inter-demand interval and the squared coefficient of variation of demand sizes syntetos2005categorization. kostenko2006note developed the SBC categorization scheme syntetos2005categorization and suggested a simple and more accurate rule, which has been widely used in the research of intermittent demand petropoulos2015forecast, spiliotis2021product.

However, croston1972forecasting's method and SBA update demand sizes and intervals, which leads to inapplicability in periods of zero demand when considering inventory obsolescence. To overcome this shortcoming, teunter2011intermittent proposed a new method called Teunter-Syntetos-Babai (TSB) to update the demand probability instead of the demand interval. TSB has been proved to have good empirical performance for the demands within linear and sudden obsolescence babai2014intermittent.

The aforementioned forecasting methods for intermittent demand are all parametric methods, which estimate the parameters of a specific distribution. Instead, non-parametric intermittent demand methods directly estimate empirical distribution based on past data, with no need for any assumption of a standard probability distribution. The bootstrapping methods, and the overlapping and non-overlapping aggregation methods dominate the research field of non-parametric intermittent demand forecasting willemain2004new, hasni2019performance, hasni2019spare, boylan2021intermittent,boylan2016performance.

In particular, temporal aggregation is a promising approach to intermittent demand forecasting, in which a lower-frequency time series can be aggregated to a higher-frequency time series. Latent characteristics of the demand, such as trend and seasonality, appear at higher levels of aggregation. nikolopoulos2011aggregate first introduced temporal aggregation to intermittent demand forecasting and proposed the Aggregate-Disaggregate Intermittent Demand Approach (ADIDA). To tackle the challenge of determining the optimal aggregation level, petropoulos2015forecast considered combinations of forecasts from multiple temporal aggregation levels simultaneously. This approach is called the Intermittent Multiple Aggregation Prediction Algorithm (IMAPA). The overall results of their work suggested that combinations of forecasts from different frequencies led to improved forecasting performance.

Recently, some attention has been paid to applying machine learning approaches to improve forecasting accuracy for intermittent demand, such as neural networks lolli2017single, support vector machines kaya2018intermittent, jiang2020intermittent, and so on.

Despite that intermittent demand forecasting has obtained some research achievements in recent decades nikolopoulos2011aggregate,petropoulos2015forecast,kourentzes2021elucidate, there is still much scope for improvements nikolopoulos2021we. For example, limited attention has been given to combination schemes for intermittent demand forecasting. The literature indicates that forecast combination can improve forecast accuracy in modeling fast-moving time series bates1969combination,de2000review, petropoulos2022forecasting,li2022improving. In this study, we aim to examine whether the forecast combination improves intermittent demand forecasts. The main contributions of our work are: (1) providing a discussion and comparison of forecast combination methods in the context of intermittent demand forecasting, (2) developing a feature-based combination framework for intermittent demand, which can determine optimal combination weights evaluated by the given error measure, and (3) improving the accuracy of both point and quantile forecasts to support real inventory decisions.

The rest of the paper is organized as follows. Section (ref) reviews a series of forecast combination methods discussed in this work. Section (ref) proposes a generalized forecast combination framework for intermittent demand. In Section (ref), we apply our framework to two real datasets and present results based on point forecasts and quantile forecasts. Section (ref) concludes the paper.

A review of forecast combinations

Combining forecasts from different methods or models has achieved satisfactory results in practice. wang2022forecast provided an up-to-date review of forecast combinations including combining point forecasts and combining probabilistic forecasts. The Simple Average (SA) has been proved to be a hard-to-beat forecast combination method clemen1989combining, stock2004combination, lichtendahl2020some, which simply combines forecasts with an equal weight of $1/M$ ($M$ is the number of forecasting methods to be combined). clemen1989combining reviewed over two hundred articles and concluded that SA should be used as a benchmark when proposing more complex weighting schemes. palm1992combine emphasized that SA could reduce the variance of forecasts and avoid the uncertainty of weight estimation. The phenomenon that SA outperforms more complicated combination methods is referred to “forecast combination puzzle" stock2004combination, smith2009simple,claeskens2016forecast.

Because SA is sensitive to extreme values, some attention has been paid to other more robust combination schemes, including the median and trimmed means stock2004combination, lichtendahl2020some, PETROPOULOS2020110. jose2008simple studied two mean-based methods, trimmed and Winsorized means, and verified their improved combined forecasts. The simple combination schemes based on the mean and median are easy to calculate and avoid parameter estimation error. However, there is still no consensus on which of the mean and the median of individual forecasts performs better.

In the field of intermittent demand forecasting, forecast combination methods have been largely overlooked. To the best of our knowledge, only SA has been applied to improving intermittent demand forecasting petropoulos2015forecast. Recently, the organizers of the M5 competition makridakis2021m5 used SA as the combination benchmark, such as the average of exponential smoothing (ES) and ARIMA. The M5 competition focused on sales forecasts involving a mass of intermittent time series. As shown in M5 results, combinations performed better or equally well with the individual methods that they consist of makridakis2022m5.

To further exploit the value of forecast combinations, a handful of research has focused on finding optimal weights for combining different forecasting models over the past half-century. The seminal work by bates1969combination proposed the idea of weighted forecast combinations. newbold1974experience continued this stream of research and investigated more forecasting models and multiple forecast horizons. In their work, a weighted combination can be expressed as a linear function such that

equation*[equation* omitted — 155 chars of source]

where $\mathbf{\hat{y}}_{T+1}$ is the column vector of forecasts at time $T+1$ generated from $M$ forecasting models, and ${\mathbf{w}}$ is the column vector of weights.

granger1984improved investigated some regressive approaches to obtain linear combinations. They demonstrated that the method with a constant term and unrestricted weights performed better. The combination weights can be estimated by Ordinary Least Squares (OLS). Linear combination has a long and successful history in forecasting. However, the issue related to determining the best set of forecasting models to combine is also worthy of attention. Lasso-based methods can do this trick by producing the selection and shrinkage toward zero tibshirani1996regression. diebold2019machine proposed a variant of Lasso, partially-egalitarian LASSO (peLASSO), which set the weights of some forecasting methods to zero and shrunk the survivors toward equality. They provided an empirical assessment to forecast Eurozone GDP growth and found that peLASSO outperformed SA and the median diebold2019machine.

The aforementioned weighted forecast combinations need to generate multiple forecasts in the training period, which multiplies the computation time. Especially for highly intermittent demand, the covariance matrix of forecast errors is often singular (can not calculate the inversion in bates1969combination's methods), because the obtained errors may have many zero values. Similarly, for regressive approaches, the standardizing process can not be implemented when the true values for training are always zero. Therefore, bates1969combination's weighted combination methods and regression-based methods are not applicable for highly intermittent data.

Recent studies indicate that using all time series in the dataset to estimate the combination weights show outstanding performance in forecasting fast-moving time series montero2021principles,talagala2022fformpp,wang2022uncertainty. One mainstream is feature-based forecast combinations. For example, monteromanso2020fforma: developed FFORMA (Feature-based FORecast Model Averaging), which used 42 features to estimate the optimal combination weights based on a meta-learning algorithm.

Different feature-based combination approaches applied different time series features to improve forecasting performance wang2009rule, petropoulos2014horses, kang2020gratis, li2022bayesian. However, a significant characteristic of intermittent demand is that there exist a large number of zeros and irregular patterns, which makes the feature sets used in previous literature inapplicable for intermittent demand. theodorou2021exploring proposed a methodological approach for feature extraction and selection to explore the representativeness of M5 dataset. On the basis of the FFORMA framework monteromanso2020fforma:, kang2022forecast used the diversity of forecasting models as the only feature. The diversity has proved to be a novel type of efficient feature, which can not only improve the forecast accuracy but also reduce the computational complexity.

The potential of time series features and the diversity of forecasts have not been investigated when producing forecast combinations for intermittent demand. In our work, we extract a set of time series features selected for intermittent demand and calculate the diversity based on a pool of intermittent demand forecasting methods. To this end, a forecast combination framework for intermittent demand can be constructed by mapping the two types of features to the combination weights based on eXtreme Gradient Boosting (XGBoost) chen2016xgboost. The proposed framework can be applied to both point and quantile forecast combinations.

Forecast combination for intermittent demand

Time series features for intermittent demand

Several studies have investigated the features of intermittent demand Kourentzes2016tsintermittent,Hara2021feasts,theodorou2021exploring. First, we consider the two most popular attributes to divide intermittent demand in the SBC classification scheme syntetos2005categorization,kostenko2006note. Then we review the 42 features selected for exploring the feature spaces of M5 competition data theodorou2021exploring. To ensure the interpretability and compute as few features as possible, we remove the features based on complex statistical methods, such as STL decomposition and Fourier transform, and take out Boolean variables with minimal information. The reserved features in theodorou2021exploring are used in our work.

Therefore, we consider nine explainable time series features for intermittent demand forecasting, which are listed in (ref). They imply the intermittency, volatility, regularity and obsolescence of intermittent demand. Given a time series $\left\{ {{y}_{t}},t=1,2,\cdots ,T \right\}$, we describe the nine features as follows.

itemize${F}_{1}$, ${F}_{2}$: The two features are average Inter-Demand Interval (IDI) and squared Coefficient of Variation (CV$^2$) to measure intermittency and demand size volatility in the SBC classification scheme syntetos2005categorization,kostenko2006note. • ${F}_{3}$: Entropy-based measures have been applied to quantify the regularity and unpredictability of time-series data kang2017visualising, theodorou2021exploring. We use the approximate entropy in this paper. A relatively small value of ${F}_{3}$ indicates that the demand series includes more regularity and is more forecastable. • ${F}_{4}$, ${F}_{5}$: The two features describe the ratios of some specific values in a given time series. ${F}_{4}$ measures the percentage of zero values. ${F}_{5}$ denotes the percentage of values lying outside $[\mu_y - \sigma_y, \mu_y + \sigma_y]$, where $\mu_y$ and $\sigma_y$ are the mean and standard deviation of time series $\{y_t\}$, respectively. • ${F}_{6}$: This feature provides the coefficient of a linear least squares regression, which measures the linear time trend of the variances of component chunks for the target series. For monthly data in the following experiments, we set the chunk length $L=12$. Moreover for daily data $L=10$, consistent with theodorou2021exploring. • ${F}_{7}$: ${F}_{7}$ first calculates the consecutive changes of the demand, i.e., the first difference of the demand series. Then the mean absolute value of the consecutive changes is taken.

The last two features focus on the presence of recent demand to capture the obsolescence, which is a challenging problem in the field of intermittent demand babai2014intermittent,babai2019new.

itemize${F}_{8}$: ${F}_{8}$ calculates the sum of squares of the last chunk out of $K$ chunks expressed as a ratio with the sum of squares over the whole series. We set $K = 4$ for the Royal Air Force (RAF) dataset and $K = 10$ for M5 competition data, so that the length of the last chunk of each series is longer than the forecasting horizon. • ${F}_{9}$: ${F}_{9}$ computes the percentage of consecutive zero values at the end of the series, i.e., the number of consecutive zero values at the end over the length of the time series.

Based on the nine time series features tailored for intermittent demand, each target time series can be represented using a nine-dimensional vector. The feature vector can be used as the input to train the forecast combination model in the proposed framework.

table[table omitted — 1,799 chars of source]

Diversity for intermittent demand

In this paper, we extend kang2022forecast's work and use the diversity of the pool of methods to develop a forecasting combination method for intermittent demand. The scaled diversity between any two forecasting methods is defined as:

equation[equation omitted — 234 chars of source]

where $H$ is the forecast horizon, ${{\hat{y}}}_{ih}$ is the $h$-th step forecast generated from the $i$-th forecasting model, and $\left\{ {{y}_{t}},t=1,2,\cdots ,T \right\}$ is a series of observed values.

Assuming that the forecasting pool contains $M$ methods, we apply (ref) to each two of them. For each target time series, we can construct a diversity vector consisting of $M(M - 1)/2$ pairwise diversity measures. This vector can be viewed as the feature vector for the corresponding series in the proposed framework.

The main merits of applying forecast diversity to intermittent demand forecasting are twofold. The first aspect is simplicity in principle, as the calculation only depends on forecasting values, with no need to compute a separate set of features. The second is general applicability. The diversity can be obtained automatically from intermittent demand forecasts and comprehended quickly by forecasters without expertise. Choosing relevant features to match the actual inventory management problem may lead to feature selection bias, especially when forecasters’ information is inadequate. Therefore, in contrast to time series features, the diversity shows remarkably simplicity and interpretability in intermittent demand forecasting.

Evaluation metrics for intermittent forecasting

In previous studies, various forecasting evaluation metrics have been used for intermittent demand. The chosen metric of forecast errors may influence the ranked performance of the forecasting methods. silver1998inventory pointed out that no single metric was universally best. wallstrom2010evaluation discussed a series of forecasting error measurements, especially for intermittent demand and split them into two categories, traditional (accuracy) and bias error measurements. As traditional measures, wallstrom2010evaluation considered mean absolute deviation (MAD), mean square error (MSE), symmetric Mean Absolute Percentage Error (sMAPE). As bias error measures, they examined the Cumulated Forecast Error (CFE), Number Of Shortages (NOS), and Periods In Stock (PIS). kourentzes2014intermittent evaluated model selection results based on two accuracy metrics. The first is the Mean Absolute Scaled Error (MASE), which was suggested to be the standard measure for the data with different scales and zero values hyndman2006another. The second is the scaled Absolute Periods In Stock (sAPIS), which is a scale-independent variant of PIS.

However, kolassa2016evaluating explored traditional accuracy measures and argued that measures such as MAD, MASE and MAPE are unsuitable for intermittent demand. A flat zero forecast is frequently “best” for the measure of MAE when the demand is highly intermittent, because zero is the conditional median of the demand. Therefore, especially for intermittent demand, an MAE-minimizing method, which is the conditional median, prefers a lower forecast than the MSE-minimizing method, which is the expectation. In recent M5 competition with much intermittent demand, the accuracy of point forecasts is required to be evaluated using Root Mean Squared Scaled Error (RMSSE) makridakis2021m5.

kolassa2020best emphasized that different error measures reward different point forecasts, and different measures should not be applied to a single point forecast petropoulos2022forecasting. In our work, we use RMSSE makridakis2021m5 to measure the performance of point forecasts, which can be obtained as:

equation[equation omitted — 230 chars of source]

where $H$ is the forecasting horizon. ${{\hat{y}}}_{T+h}$ is the $h$-th step forecast generated from a series of observed values $\left\{ {{y}_{t}},t=1,2,\cdots ,T \right\}$, and ${{y}_{T+h}}$ is the true value.

Generalized forecast combination framework

The quality of forecast combination has been demonstrated to depend on the individual forecasts as well as the diversity between forecasts lemke2010meta, kourentzes2019another, kang2022forecast. Therefore, defining an appropriate forecasting pool is one of the most crucial steps in the forecast combination process. Firstly, we define a broad pool for intermittent demand forecasting. The pool includes traditional forecasting models, which are Naive, seasonal Naive (sNaive), Simple Exponential Smoothing (SES), Moving Averages (MA), AutoRegressive Integrated Moving Average (ARIMA), ExponenTial Smoothing (ETS), and intermittent demand forecasting methods, which are Croston’s method (CRO), optimized Croston’s method (optCro), SBA, TSB, ADIDA, IMAPA. The 12 forecasting methods in the pool are considered as statistical benchmarks in the M5 competition makridakis2021m5. In contrast to CRO with fixed smoothing parameters, the parameters in optCro are optimized to allow for more flexibility. Implementations for these methods exist in the {\fontseries{b}\selectfont forecast} Hyndman2020forecast and {\fontseries{b}\selectfont tsintermittent} Kourentzes2016tsintermittent packages in \proglang{R}. Then the pooling methods kourentzes2019another, lichtendahl2020some, diebold2019machine can be applied to reduce the number of forecasting methods and further improve the quality of the forecasting pool. We study the effect of three popular pooling algorithms in Section (ref).

In the proposed forecast combination framework, we build an XGBoost model to learn the relationship between features and combination weights. This approach transforms the combination problem into a classification problem by setting the best forecasting method as the target class for each time series. The two types of features in Section (ref) and Section (ref) are all valid inputs for the forecast combination model. We name the approach based on the nine time series features in Section (ref) as Feature-based Intermittent DEmand forecasting (FIDE). The diversity-based method is called DIVersity-based Intermittent DEmand forecasting (DIVIDE).

Then, given a forecast error metric, the optimization objectives for the FIDE and DIVIDE are

subequations\begin{align} &\underset{{w}_{F}}{\mathop{\arg \min }}\,\sum\limits_{n=1}^{N}{\sum\limits_{i=1}^{M}{w{{\left( {{F}_{n}} \right)}_{i}}}}\times error_{n,i}, and \\ &\underset{{w}_{D}}{\mathop{\arg \min }}\,\sum\limits_{n=1}^{N}{\sum\limits_{i=1}^{M}{w{{\left( {{D}_{n}} \right)}_{i}}}}\times error_{n,i}, \end{align}

respectively, where ${F}_{n}$ is the feature vector, ${D}_{n}$ is the diversity vector of the $n$-th time series, $N$ is the number of time series, and $M$ is the number of forecasting methods. $\text{error}_{n,i}$ is the forecast error of the $i$-th method for the $n$-th time series. RMSSE is used as the error measure of point forecasts in this paper. RMSSE focuses on the expectation, which is consistent with the candidate methods in the forecasting pool. Once the model has been trained, weights can be produced for a new series to generate the combined forecast. The process can be implemented based on the \proglang{R} package {\fontseries{b}\selectfont M4metalearning} by monteromanso2020fforma:.

Based on the forecasting pool for intermittent demand, we put forward a generalized forecast combination framework containing FIDE and DIVIDE. The flowchart of the proposed framework is presented in (ref). In the training phase, we generate forecasts based on the methods in the intermittent demand forecasting pool and calculate errors required in the objective function. In FIDE, we compute the features selected for intermittent demand and learn the relationship between the features and combination weights by (ref). In DIVIDE, the combination model can be obtained based on the diversity of different forecast methods (see (ref)), where the pairwise diversity values of the methods in the pool are used as time series features. Therefore, DIVIDE can be viewed as a special case of FIDE. In the forecasting phase, we calculate the features or the diversity for the new time series, and get the combination weights through the pre-trained XGBoost model. Finally, we utilize the optimal weights to average the forecasts from different methods in the pool and achieve the combined forecast results.

To evaluate the forecasting performance of the proposed framework, the time series need to be divided into three periods. Let $H$ be the forecast horizon and $T$ be the length of data. The first $T-H$ observations are used for training the forecast combination model. Then the $T-H$ observations are split into $T-H-H$ for training the forecasting methods in the pool and $H$ for testing. The final $H$ observations are used to evaluate the forecasting results.

The merits of the proposed framework include: (1) using a diverse forecasting pool, consisting of intermittent demand forecasting methods and traditional time series forecasting models, (2) considering a customizable objective function depending on actual inventory management requirements, (3) selecting intelligible time series features especially for intermittent demand, and (4) calculating the diversity with simple form only based on forecasting values.

figure[figure omitted — 393 chars of source]

Empirical evaluation

Real dataset

The proposed methods are applied to two real datasets. The first RAF dataset has been previously investigated in the literature kourentzes2021elucidate,petropoulos2015forecast,teunter2009forecasting. It contains 5000 monthly time series, with 84 observations each. Moreover, the second is M5 competition data, involving the unit sales of 3049 products sold by Walmart in the USA between 2011-01-29 and 2016-06-19 (1969 days). The dataset was organized in the form of hierarchical time series in M5 competition makridakis2021m5. We only consider the bottom level, i.e., 30,490 product-store unit sales in this paper.

In the following experiment, we examine three forecast horizons of 3, 6 and 12 months ahead for RAF dataset and 28-day-ahead forecasts for M5 data as required in M5 competition makridakis2021m5. The final observations of the horizon length are used to evaluate the forecasting performance. The seasonal periods are 12 and 7 for monthly and daily demand, respectively. Moreover, we preprocess the data before forecasting, removing the initial zero values and making the first non-zero demand as the initial value theodorou2021exploring. This is due to a lack of information that the initial zeros mean demands or sales are zero, or the product has not been in stores yet.

Following the SBC scheme kostenko2006note, the time series can be divided into four categories based on IDI and CV$^2$. (ref) describes the distributions of RAF and M5 datasets, respectively. The boundaries of different categories in (ref) are IDI$ = 4/3$, CV$^2 = 0.5$ kostenko2006note. The equal sign is placed on the less-than sign when classifying, e.g., if IDI$ > 4/3$ and CV$^2 <= 0.5$, the time series is “intermittent”. As shown in (ref), the RAF dataset exhibits high intermittence and contains 2729 intermittent, and 2271 lumpy series. While the M5 data has a wider-ranging distribution, including 22,206 intermittent, 5359 lumpy, 897 erratic, and 2028 smooth series.

figure[figure omitted — 937 chars of source]

Point forecasting

We compare our methods with individual models, SA, Median and FFORMA. Other combination methods reviewed in Section (ref), such as bates1969combination's original forecast combinations and regression-based methods granger1984improved,diebold2019machine, are omitted here, which are not applicable for highly intermittent data. We present the forecasting accuracy of different methods based RAF dataset and M5 competition data in (ref) and (ref) respectively.

As shown in (ref), the best individual method is ADIDA for all forecast horizons based on RMSSE. The simple combination methods (SA and Median) can not beat the best individual method. While the proposed methods based on intermittent demand features and the diversity consistently outperform others. FIDE shows obvious superiority compared with FFORMA using 42 time series features and offers the best forecasting results overall. Therefore, our chosen features are more appropriate to describe intermittent demand compared with the time series features in FFORMA designed for fast-moving data. Furthermore, the improved forecast accuracy of DIVIDE indicates that the diversity is a simple and efficient tool for intermittent demand forecasting combinations.

The forecasting results of M5 competition data in (ref) are organized by SBC classification scheme kostenko2006note. For each column (a subset of data), the combination model is optimized respectively. The last three rows show the RMSSE of the top three winning methods in the M5 competition for comparison, e.g., “M5-w1” denotes the first ranked method. The results of last three rows in (ref) are calculated based on corresponding M5 submissions. It should be noted that the M5 competition took the hierarchy of data into consideration and used weighted RMSSE to rank participants. The weights put more emphasis on the series that account for higher monetary sales makridakis2021m5. Therefore, comparing these methods in absolute terms is not entirely fair, as they were obtained from different application contexts and optimization objectives.

As shown in (ref), the best individual method is IMAPA for all classifications, which is inconsistent with the RAF data. The results emphasize the risk of choosing a single forecasting method and elicit the necessity of forecast combination. Our proposed methods achieve the competitive forecasting performance based on RMSSE when compared with the top three ranked methods in the M5 competition. Based on different classifications of M5 data, the performance of the proposed methods exhibits significant differences. The proposed FIDE and DIVIDE outperform FFORMA for the intermittent and lumpy data, which is consistent with the RAF dataset. The forecasting results based on the two datasets provide good evidence for the superiority of the proposed framework in intermittent demand forecasting. However, for the erratic and smooth data, our methods perform slightly worse than FFORMA. We acknowledge the limitations of the features used in the proposed framework, which are more applicable to intermittent demand.

table[table omitted — 1,367 chars of source]
table[table omitted — 1,737 chars of source]

The features and diversity in the proposed framework have been proven to improve the accuracy of intermittent demand forecasting. In the following experiments, we continue to provide a sensitivity analysis of multiple features, including diversity viewed as another form of features. We investigate the relationship between RMSSE and the number of features used in the proposed FIDE and DIVIDE, as shown in (ref). In FIDE, we set the feature number to be 3, 6 and 9(all). While in DIVIDE, the feature number varies across 5, 10, 20, 30, 40, 50 and 66(all). We use two ways to alter the number of features. One is to select features in the order of feature importance in XGBoost model, and the other is by random feature selection. The importance of each feature to the FIDE or DIVIDE is measured by the gain of features in the XGBoost model chen2016xgboost.

As shown in (ref), with the increase of the feature number, RMSSE of FIDE decreases when selecting features randomly (blue lines), but increases when selecting features in order of the importance (red lines). The findings emphasize the importance of choosing appropriate features as input in the proposed FIDE. For RAF dataset, the three most important features are Percent.zero.end ($F_9$), Ratio.last.chunk ($F_8$) and Linear.chunk.var ($F_6$). While for M5 competition data, the top three features are Percent.zero.end ($F_9$), IDI ($F_1$) and Ratio.last.chunk ($F_8$). Thus, the features to capture the recent demand are more critical for constructing the combination model in FIDE. However, the relationship between RMSSE and the feature number in DIVIDE is markedly different from that of FIDE. As the number of features increases, the overall trend of RMSSE is downward, though there is a non-significant increase when considering the importance of features for RAF dataset. Therefore, we recommend applying the whole diversity in the proposed framework.

figure[figure omitted — 733 chars of source]

To analyze the efficiency of the examined forecasting methods, we investigate the relationship between RMSSE and computation time based on RAF and M5 datasets. The results are computed indicatively for RAF dataset when $H=12$. As shown in (ref), the time consumption of our methods mainly comes from the individual forecasting methods, which is the limitation of forecast combinations. Simple combination schemes, such as SA and median, can save nearly half the time, but they perform worse than the best individual method for intermittent demand. The proposed framework generates forecasts in the training and forecasting periods, increasing the time consumption. However, the process of model training is in the off-line phase in real applications. Therefore, the increased computational time of training does not affect the forecasting efficiency. Compared with FFORMA, our methods based on features for intermittent demand and diversity are more computationally efficient, especially for the M5 dataset with large amounts of long time series. Based on these findings, decision makers can consider the trade-off between accuracy and computational cost in actual inventory management.

figure[figure omitted — 501 chars of source]

Reducing the number of forecasting methods used in the proposed framework can significantly save computational time. For instance, the computation time of our methods can be halved by removing the three most time-consuming methods in the forecasting pool. In the following experiment, we aim to investigate the potential of shrinking the forecasting pool to further improve the accuracy of the proposed framework. We call this process pooling, deriving from kourentzes2019another's research.

We study the effect of three pooling algorithms based on the RAF (indicatively for $h=12$) and M5 dataset, which are forecast islands proposed by kourentzes2019another, a screened method from lichtendahl2020some, and a Lasso-based method by diebold2019machine. The forecast islands kourentzes2019another remove some poorly performing models from the pool, which is shortened to “Islands". It conducts ${C}'=\left\{ 0,\Delta C \right\}$ for a series of ordered forecasts based on a criterion of forecasting performance $C$ and includes all forecasts until ${C}'\ge T$. $T=\text{Q3}+1.5\text{IQR}$ is related to the outlier of the boxplot, where Q3 is the third quartile and IQR is the interquartile range. The screened method lichtendahl2020some, shortened to “Screened", screens out forecasting models with highly correlated errors (correlation coefficient is over 0.95). The Lasso-based method diebold2019machine, “Lasso” for short, sets the regression coefficients of some forecasts to zero via a standard Lasso software, i.e., \proglang{R} package {\fontseries{b}\selectfont glmnet} Friedman2010Regularization, and the survivors form the final forecast pool.

(ref) presents the forecasting accuracy of original FIDE and DIVIDE and those considering the three pooling algorithms. In (ref), we find that none of the pooling methods significantly improves the forecasting performance. The proposed framework automatically reduces the weights of some methods to minimal values, which can be regarded as a generalized pooling method customized for each time series. The Islands remove the worst performing methods in the pool, but reduce the accuracy, especially for RAF dataset. For highly intermittent data, the poorly performing methods, such as Naive, sNaive, make important contributions to the forecast combination. Therefore, unless the pool of forecasting methods is too large which would render the computation process time-consuming, there is no need to add modeling complexity by implementing pooling approaches on top of our proposed framework.

table[table omitted — 801 chars of source]

Quantile forecasting

Based on the improving accuracy of point forecasts, the proposed combination methods are shown to be effective in providing robust forecasts to support decisions. However, in real supply chain management, estimating the right part of the demand distribution is also necessary for determining safety stock levels, which has been largely ignored in the research barrow2016distributions, spiliotis2021product. fildes2019retail reviewed retail demand forecasting and emphasized the connection of quantile, density, or volatility forecasting to the inventory control.

The intermittent demand forecasting methods (such as CRO, optCro, SBA, TSB, ADIDA and IMAPA) can not directly output quantile forecasts. We apply trapero2019empirical's empirical approach to estimate the desired quantiles. They recommended a kernel density estimation to model the forecast error distribution. We generate quantile forecasts by adjusting the point forecasts based on the respective quantiles calculated from the empirical distribution of residual errors, as follows:

equation[equation omitted — 118 chars of source]

where $T$ is the length of observations, ${{Q}_{T+h}}\left( u \right)$ is the probabilistic forecast for quantile $u$ at time $T+h$, ${\hat{y}}_{T+h}$ is the $h$-th step point forecast, ${\hat{q}_{|e}}\left( u \right)$ is the estimated $u$-th quantile of the residual errors. As shown in (ref), we assume that the demand pattern that occurred in the past will continue in the future. The obtained forecasts are based on in-sample approximations without requiring computing multiple forecasts. The approach has been verified to perform well for the RAF and M5 datasets spiliotis2021product, kourentzes2021elucidate.

The proposed framework can be extended to quantile forecast combinations by mapping the features to the errors of quantile forecasts. In FIDE and DIVIDE, we still compute the nine features in Section (ref) based on historical data and the diversity of different point forecasts as shown in Section (ref). The training and testing processes are consistent with Section (ref). We use the Scaled Pinball Loss (SPL) function to measure the precision of the quantile forecasts, which is required in the M5 competition makridakis2021m5. The SPL can be obtained as follows:

footnotesize\begin{equation} \begin{aligned} & SPL(u)= \\ & \frac{\sum\limits_{h=1}^{H}{u\left( {{y}_{T+h}}-{{Q}_{T+h}}\left( u \right) \right)\mathbf{1}\left\{ {{Q}_{T+h}}\left( u \right)\le {{y}_{T+h}} \right\}+\left( 1-u \right)\left( {{Q}_{T+h}}\left( u \right)-{{y}_{T+h}} \right)\mathbf{1}\left\{ {{Q}_{T+h}}\left( u \right)>{{y}_{T+h}} \right\}}}{H\cdot\frac{1}{T-1}\sum\nolimits_{t=2}^{T}{\left| {{y}_{t}}-{{y}_{t-1}} \right|}}, \\ \end{aligned} \end{equation}

where ${{y}_{T+h}}$ is the actual future value of the examined time series at point $T+h$, ${{Q}_{T+h}}\left( u \right)$ is the generated forecast for quantile $u$, $H$ is the forecasting horizon, $T$ is the length of the number of historical observations, and 1 is the indicator function (being 1 if true is within the postulated interval and 0 otherwise).

The following experiment based on the RAF and M5 datasets focuses on four quantiles, i.e. $u_1 = 0.750$, $u_2 = 0.835$, $u_3 = 0.975$, and $u_4 = 0.995$. $u_1$ and $u_2$ provide a good sense of the mid-right part of the distribution, while $u_3$ and $u_4$ provide information about its right tail, which is essential for the risk of extreme outcomes. We customize the objective function by assigning the error measure to SPL based on the correlated quantiles. Therefore, we obtain different combination weights for the four quantiles, respectively. The forecasting results in (ref) are computed based on 12-month-ahead forecasts for RAF dataset and 28-day-ahead forecasts for M5 data.

We can find in (ref) that the performance of individual methods changes considerably based on different quantiles. The finding indicates that each method is more appropriate for estimating different parts of the distribution of the series, which echoes the weakness of choosing a single method. The proposed combination methods exhibit steady performance across the four quantiles. For RAF dataset, FIDE and DIVIDE consistently outperform the rest in (ref). For M5 dataset, we add the top three ranked methods in the M5 competition for comparison, the results of which derive from spiliotis2021product. The proposed DIVIDE outperforms the first ranked method in M5 competition at quantiles 0.835, 0.975 and 0.995. While at quantile 0.750, our methods also provide competitive forecasting results. The improved performance of high quantile forecasts can contribute to practical inventory decisions for higher levels of service.

table[table omitted — 2,597 chars of source]

Conclusion

This paper focuses on forecast combinations for intermittent demand. We review a handful of forecasting methods, and investigate the performance of some existing forecast combination methods for intermittent demand. We introduce time series features and diversity to propose a generalized forecast combination framework, which can automatically determine the optimal combination weights. We conduct an empirical investigation based on real-life data to analyze the forecast accuracy and gain insights related to inventory decisions.

The results of point forecasts are measured by RMSSE, which focuses on the expectation. The proposed framework notably outperforms other combination methods and the best individual method, especially for the RAF dataset with highly intermittent series. Moreover, for M5 competition data, our methods achieve a competitive performance compared with the top three ranked methods in the M5 competition. In addition, the proposed framework can be regarded as a generalized pooling method customized for each time series by reducing the weights of some methods to minimal values. The empirical evaluation based on RAF and M5 datasets provides good evidence of the superiority and flexibility of the proposed framework. We acknowledge that our combination methods increase the computational time compared with individual methods. Decision makers should consider the trade-off between accuracy and computational cost in actual inventory management.

The proposed framework has also been applied to quantile forecast combinations, especially for high quantiles to estimate the right part of the demand distribution. We use SPL to measure the quantile forecasting performance and make it used in the optimization objective. The examined results show that our methods can provide accurate forecasts of both central tendency and high quantiles, which directly connect with the inventory decision.

The good performance of our proposed framework can be attributed to: (i) defining a forecasting pool suitable for intermittent demand, which consists of several intermittent demand forecasting methods and traditional time series forecasting models, (ii) applying diversity and time series features to determine the optimal combination weights automatically, and (iii) applying to both point and quantile forecasts to support inventory decisions. The diversity and the features selected for intermittent demand are all effective inputs of the proposed framework. Extracting the diversity independent of historical data makes it more flexible for intermittent demand forecasting, especially when the training set is limited in positive demands. In addition, the features in FIDE are all easily understood. The two features focusing on the presence of recent demand are proved more critical for constructing the forecast combination model. These advantages of the proposed methods lead to broad application prospects in intermittent demand forecasting.

However, we recognize the lack of a comprehensive evaluation of inventory performance in the current study. petropoulos2019inventory combined financial, operational, and service metrics to form a holistic measure for inventory control objectives. ducharme2021forecasting focused on stock-out events and proposed a novel metric called Next Time Under Safety Stock. The utility measures are essential to achieve a direct link between inventory holding costs and service levels in the production system. Such analysis needs to proceed based on restocking policies, which are not available for the RAF and M5 datasets without any background information of inventory. Future research should investigate the inventory performance of our proposed framework in the field of a specific inventory management problem. Another limitation of this paper is lacking an automatic procedure to choose features for modeling FIDE. Several scholars have investigated selecting features automatically from a large number of features lubba2019catch22,theodorou2021exploring. Although these approaches seem more general, they take over much computational time, and the selected features are often difficult to understand in the applications. Based on the results of our work, the nine features in FIDE are efficient and can be used as the benchmark pool of features for intermittent demand. In further research, we will study a standard procedure to select features automatically for the proposed framework, aiming to achieve both interpretability and computational efficiency.

Data Availability Statement

The RAF dataset has been used in previous literature teunter2009forecasting,petropoulos2015forecast,kourentzes2021elucidate and is available upon request. The M5 competition makridakis2021m5 data involves the unit sales of 3049 products between 2011-01-29 and 2016-06-19 (1969 days). The first 1941 observations for model training can be obtained from \url{https://github.com/Mcompetitions/M5-methods}; the final 28 observations is available upon request.

Acknowledgments

Yanfei Kang is supported by the National Natural Science Foundation of China (No. 72171011). Feng Li is supported by the Beijing Universities Advanced Disciplines Initiative (No. GJJ2019163) and the Emerging Interdisciplinary Project of CUFE. This research was supported by Alibaba Group through the Alibaba Innovative Research Program and the high-performance computing (HPC) resources at Beihang University.

Declaration of Interest Statement

No potential conflict of interest was reported by the authors.