Forecasting risk analysis for supply chains with intermittent demand

: This paper focuses on the forecasting risk analysis in supply chains with intermittent demand, which is typical for the inventory management of the ‘slow-moving items’, such as service parts or high-priced capital goods. The adopted demand model is based on the Generalised Beta-Binomial Distribution (GBBD), which is capable of incorporating the additive distortions in the demand historical records as parameters. For this setting, there are proposed explicit expressions for forecasting risk and the prediction function, which minimises the error impact on the risk. The efficiency of the proposed approach is confirmed by computer simulation and is illustrated by an application example for forecasting of the intermittent demand values for car spare parts


Introduction
Recently, performance measuring in manufacturing and logistics has been highlighted by a number of researchers.Owing to growth of enterprise collaboration and new challenges in managing supply chains, risk evaluation has become a very important issue.This paper focuses on minimising the forecasting risk for inventory management in supply chains with intermittent demand in the case of reporting errors in the past demand records.
For large-scale supply chains, accurate demand forecasting is motivated by strong financial incentives, since there is a high need to balance the customer service levels against investments over a large assortment of Stock-Keeping Units (SKU) (Harris, 1997;Hopp and Spearman, 1995;Silver et al., 1998).As follows from related research (Snyder, 2002), a large proportion of SKUs is usually made up of 'slow-moving items' with irregular (intermittent) demand patterns, that is, random sequences with a large proportion of zeros and great variability among the remaining non-zero integer values (Willemain et al., 2004).A typical case study of NSK Corporation (Ann Arbor, MI, a precision parts manufacturer for the automotive industry) reported by Smart (2002) shows that two-thirds of the 7000 SKUs demonstrate highly intermittent demand behaviour, with half of the zeros in the history data.Similar demand patterns are common for all the service parts businesses in the aerospace, automotive and high-tech electronics industries.The intermittent demands are also usual for high-priced capital goods, such as heavy machinery, where the inventory replenishment lead times can be up to several months, and the costs may reach several hundred thousand dollars.
In contrast to the standard marginal stock level setting, the inventory management of the slow-moving items requires a probability distribution of the SKU demand.It is used to generate the demand forecasts, which are key inputs to the inventory control models allowing to correct time and size of replenishment orders (reorder points and order quantities).Besides, the demand distribution is particularly important for estimating the customer service level (e.g. a 95% or 99% likelihood of not stocking out over a lead time).In such situation, the conventional statistical forecasting methods (exponential smoothing and moving averages) are limited to forecasting of the average demand per period only.They are not able to produce credible results for the complete set of all possible demand values, since they implicitly assumed the normal, bell-shaped distribution, which ignores both the integer nature of the historical data and special role of zeros.And with the intermittent data, the distribution shape can be essentially non-symmetric and even non-unimodal, with several maximums.So often, the conventional methods may produce misleading inputs to inventory control models, with rather costly consequences (Grange, 1998;Hopp et al., 1997).For this reason,it is important to apply an adequate demand probability distribution that ensure accurate reorder points forecast.
The first heuristic technique for the intermittent demand forecasting was developed by Croston (1972) who applied separate exponential smoothing to both the non-zero demand sizes and inter-arrival times between successive demands.Later, the problem was addressed by a number of authors (Willemain et al., 1994;Johnston and Boylan, 1996;Syntetos and Boylan, 2001) who concentrated on improving the forecast accuracy and proposed several modifications of the Croston's method, including EWMA-smoothing and log-transformation of the data (both the demand values and inter-arrival times) in order to avoid the positive-space constraint.However, the corresponding model variables were still deterministic and defined in a continuous space.Consequently, all of these methods do not produce forecast distributions and associated prediction intervals (Hyndman and Shenstone 2003).
The first stochastic demand models were based on the Poisson or NBD (negative-binomial) distributions, and they originated from research on stochastic inter-purchase times and proved to be very accurate in the fitting of the aggregated data describing frequently purchased goods (Agrawal and Smith, 1996;Dunn et al., 1983;Grange, 1998;Gupta, 1991;Wagner and Taudes, 1987).However, their basic assumption on exponential (or gamma-exponential) distribution of the inter-arrival times does not take into account the multimodality, which becomes a vital issue for personalisation of the marketing decisions.To overcome this problem, Telang et al. (2004) proposed recently a hierarchical probabilistic model of user's repeat visits that incorporates Weibull, Laplace and double-exponential mixed distribution to account for users' schedule.This model was successfully applied to forecasting of visits to massively popular internet websites, but it seems to be hardly applied to the intermittent demand modelling, which possesses the above-mentioned specificities.In this paper, the multimodality is resolved within the Bayesian framework (Azoury and Miller, 1984), using the Generalised Beta-Binomial Demand Distribution (GBBD) developed in our previous paper (Dolgui et al., 2005).
However, even when the selected stochastic model of the demand is adequate, it still can be misleading in cases of reporting errors in the demand records (Jacobs and Wagner, 1989).A number of researchers investigated the related problem of unobserved lost sales and proposed efficient identification algorithms for relevant stochastic demand models (Agrawal and Smith, 1996;Nahmias, 1994;Tripathi et al., 1994).However, no research has been reported on the performance of the Beta-Binomial Demand Model (BBM) in such settings.So, this paper focuses on the analysis of the Beta-Binomial Distribution-(BBD-) based forecasting risk and also proposes consistent estimators for the demand values, which are robust in respect to the distorted historical data.This paper also contributes in some neighbouring areas, such as implication of errors in survey data (Gaba and Winkler, 1992), where the data distortions may have an essential effect on inferences, such as too-high point estimates and too-narrow estimation intervals.
The remainder of the paper is organised as follows.Section 2 defines the general inventory optimisation problem and justifies applications of the GBBD.In Section 3, the BBM is described in detail and its main characteristics are presented.Section 4 focuses on minimisation of the forecasting risk and contains the main theoretical deliverables, such as the explicit expressions for the forecasting risk and the prediction function, which minimises the error impact on the risk.Section 5 contains simulation results.Section 6 describes an application example.Section 7 summarises the main contributions of this paper.

Inventory optimisation problem
A basic inventory optimisation problem can be stated as follows (Grange, 1998).Denote by m the number of SKU and let D l and Q l be the demand rate and stock guidance ('order-up-to' quantity) of the lth SKU over a fixed period of time (e.g. one week).Define the total fill rate as ( ) where E{D l } is the expected demand for the lth SKU, P{D l ≤ Q l } is the individual service level of this SKU, and Then the problem is to minimise the total inventory cost ensuring the target total fill rate F: In order to solve problem (2) that is rather common for the combinatorial optimisation, one must know the probability distribution of the demand Dl for each SKU.Thus, this paper focuses on statistical aspects of the inventory optimisation, namely consistent estimation of the distribution parameters and minimisation of the forecasting risk with respect to the demand data distortions.
As was proved by Grange (1998), an adequate demand model may be based on the BBD.An essential advantage of the BBM versus the classical Binomial Model (BM) is that the BBM takes into account possible correlation of the demand within a time period and between the periods, while the BM assumes that all requests are independent.The BBM also has some advantages over the Poisson-based models, which assumes only binary demand values and perform poorly for many industrial data sets (Willemain et al., 2004).In this paper, we apply a Generalisation of the Beta-Binomial Model (GBBM), proposed in our earlier research (Dolgui et al., 2005).

Beta-binomial demand model
Let us assume that the demand data are arranged in an integer matrix {s following the generalised beta distribution where i denotes the item number (within a group of similar SKUs), j defines the observation period, n is the model parameter, ( , ) B α β is the complete beta function; α, β and π 0 , π 1 are the shape and range distribution parameters respectively (0 ≤ π 0 < π 1 ≤ 1).This distribution takes into account the intermittent data nature by implicitly assuming that the observation time-segment may be uniformly divided in n time segments (an hour, a day, etc.), within which not more than one order can be received.It should be noted that such formulation follows the concept of the hierarchical (or mixed) statistical distributions where the probability of a random event is also treated as a random variable.
Under such assumptions, the generalised BBM may be expressed as the weighted sum of the shifted beta-binomial PDFs: where P r is the probability that the demand for the time-segment is equal to r, and the weights are expressed as To obtain this expression, it is necessary to apply the linear transformation π π Δ = − So, the marginal distribution of r may be presented as ( ) Finally, the obtained sum of the integrals may be expressed via the beta-functions to yield (4), ( 5).
For computing convenience, the ratios of the beta-functions ( , )/ B α β can be further simplified and Equation ( 4) can be rewritten as: Using the same approach with the dummy random variable θ, the demand mean E(r) can be expressed via the conditional expectation as 0 and after substitution of E(θ), this transformation yields the following average demand value: where the terms /( ) nα α β + and /( ) + can be interpreted as the average non-zero/zero probabilities for the standard BBD-model.Similarly, the demand variance V(r) can be expressed via the conditional mean and variance as ( ) and after computation of the conditional components and substitution of E(θ) and V(θ) one can obtain: ( ) ( ) where the first term An essential feature of the GBBM is that it is capable of describing additive binary distortions in the demand records.As proved in our previous work (Pashkevich and Dolgui, 2005), to take into account distortions, the GBBM should be parameterised in the following way: π 0 = ε 0 and π 1 = 1−ε 1 , where ε 0 is the probability that a false order was erroneously recorded, and ε 1 is the probability that an order was wrongly omitted in the record.

Minimising forecasting risk
Traditionally, the demand probabilities are forecasted using the BBM Bayes predictor (Collet, 2002), which due to the conjugate property of the classical BBD is reduced up to ( ) ( ) /( ), p s s n where , α β are the estimates of the BBM parameters, and s is the last observed demand value.This predictor minimises the mean square forecast error in the case of the consistent estimation of α, β.However, in the case of the data distortions, the classical predictor looses its optimality.Hence, a new predictor was developed, which minimises the forecasting risk in the case of the reporting errors.
It was proved (Kharin and Pashkevich, 2003) that for the BBM with the reporting errors, the mean square optimal forecast obeys the following probability distribution ( ) and the predictor function is the weighted sum where the weights ω si are defined as follows 0 ( , ) To evaluate the related forecast risk, let us note that in the non-distorted data case the risk is computed from the expression 2 0 ( )( 1)( ), R n αβ α β α β α β = + + + + + which includes the distribution shape parameters α, β only.Then, if the data are corrupted by the additive distortions with levels ε 0 , ε 1 , it is proved that the forecasting risk of the classical Bayes predictor is increased in accordance with the expression Since the distortions also influence the model parameter identification and α, β-estimators may produce non-consistent result, the forecasting risk is further increased.It was proved (Pashkevich and Kharin, 2004) that in this case the forecasting risk can be minimised by using the proposed prediction function which in the case of relatively small error levels is approximated by the expression which is the product of the classical Bayes predictor and the correction factor depending on the error levels ε 0 , ε 1 .Application of the proposed robust predictor and its influence on the minimisation of the forecasting risk is illustrated in the following section.

Simulation study
To evaluate efficiency of the proposed forecasting and risk evaluation procedures, they were applied to a number of randomly generated data sets that originated from the beta-binomial demand distribution with parameters α = 0.5, β = 9.5, n = 10.Then, these data were contaminated by imposing the additive binary distortions with levels ε 0 , ε 1 ∈[0, 0.05].Totally, 10,000 random samples were generated and used for both the model parameter estimation and demand value forecasting.Two basic cases were investigated based on the assumptions of the known and unknown distortion levels.
For the case of the known error levels, simulation results are shown in Figure 1, where R denotes the forecasting risk and ε is the distortion level, assuming that ε 0 = ε 1 .
As follows from the figure, the distortions influence is non-negligible and they may cause essential increase of the forecasting risk.For instance, in the error-free case, the risk is equal to 0.046 while for ε = 0.05 the classical technique produces results with the risk about 0.071.In contrast, the proposed robust predictor yields 1.3 times better result, with the risk value of 0.054.It is comparable with the minimum possible bound defined by the Bayesian theory (only by 0.008 higher than the lower bound), that confirms the proposed method efficiency.In the case of the unknown distortion levels, the experimental data sets were first used for the identification of the error levels.Then, the estimates of ε 0 , ε 1 were passed on to the forecasting routine.The obtained results are also optimistic, while even theoretically they cannot be obviously as good as in the previous case because of the additional source of the uncertainty related to the stochastic distortions.As follows from the computer experiments, for ε = 0.05 the proposed estimation method (given in our previous paper, Dolgui et al., 2004) in combination with the prediction technique, proposed here, is able to reduce the forecasting risk down to 0.55.Hence, the developed forecasting method demonstrates a robust behaviour in respect to the data distortions, which is attractive for industrial applications.

Application example
To validate the efficiency of the developed forecasting method in real world settings, the demand data were used for two service parts reported by Snyder (2002).The review period was assumed to be one month, and values for 36 months were available.
The distribution parameters α, β and distortion levels ε 0 = ε 1 were estimated using the demand data for the first 24 months (for several values of n).
Since the demand values possess stable seasonality, the Bayes forecast for the third year (months 25-36) was generated using the corresponding month demand of the second year (months 13-24).Two types of the predictors were used: the classical Bayes and the developed ones.The performance of each forecast was evaluated by comparing the square root of the mean square forecast error, that is, the forecasting risk (R).
Tables 1 and 2 present the actual demand for the months 25-36 and the forecast results for different values of the parameter n.It should be noted that the initial n-value is defined as the maximum demand in the past samples (months 1-24), and then it is slightly increased to consider the case of possible larger demand in the future.As follows from the tables, the proposed technique ensures essentially lower forecast error for the part 2, and a little lower error for part 1 when compared to the classical approach.As an example, for part 2, n = 8, the classical model leads to the error 1.72, while the proposed model ensures a value of 1.30 that results in 33% risk reduction.For the considered demand data, the risk sensitivity was also analysed with respect to the parameter n.As follows from relevant computations, the classical demand model is rather non-sensitive to selecting of n-value, while the generalised model (incorporating reporting errors) should be properly tuned to achieve the low value of the forecast risk.
In particular, it is prudent to define n close to the maximum demand in the learning sample, whereas unjustified increases of n may reduce the model efficiency.However, this complication in model tuning is compensated by reduction of the forecast risk.

Conclusion
Forecasting intermittent demand is a challenging problem of inventory management.
Current practice in such modelling favours the exponential smoothing of the demand values or applying exponential smoothing separately to the intervals between non-zero demands and their sizes (Croston's method).Alternatives recommended in inventory control literature are based on the bootstrap or Poisson distribution, which assumes exponential inter-arrival times.This paper describes another demand model, based on the generalised BBD, which is capable of incorporating the additive distortions as parameters.In combination with the Bayesian technique, this model is also able to describe the intermittent demand patterns with essentially non-exponential (poly-modal) distribution of the inter-arrival times.This paper contributes to the forecasting risk analysis in the presence of demand record errors, which are modelled as binary additive distortions that contaminate the actual demand data described by the generalised BBD.For this setting, explicit expressions were derived for the forecasting risk and the prediction function, which minimise the error impact.The efficiency of the proposed approach is confirmed by computer simulation and is illustrated by an application example for forecasting the intermittent demand values for car spare parts.
m}, where each observed value s ij is presented as a sum of n binary random variables Bernoulli trials with the random probabilities p j w rl computed via the Binomial theorem.
be treated as the weighted variance of the standard BBD with the scale factor equal to the square range of p, and the remainder ones are the weighted variances of the binomial distribution with parameters π 0 , π 1 .
the ascending/descending factorials of order k.Related proof is based on computing of the mean square error of the Bayes estimator taking into account the probability distributions (4), (11).

Figure 1
Figure 1 Forecasting risk for classical and proposed predictor functions (theoretical curves and experimental 95% confidence intervals)