This is a foundational article about an algorithmic inferential modeling approach that supports an agile consumer-centric demand forecasting and inventory planning process without conventional normal (Gaussian) modeling assumptions. Conventional modeling assumptions may no longer work effectively in today’s disrupted supply chain environment, where unusual events and outliers tend to dominate demand history. Demand forecasting is becoming a greater challenge for demand and supply planning decision support systems designed with standard normal (Gaussian) modeling assumptions.

Regular and intermittent demand forecasting procedures are used for demand and inventory planning processes in practice, when 

• modeling product mix based on future patterns of demand at the item and store level
• selecting the most appropriate inventory control policy based on type of stocked item
• setting safety stock levels for SKUs at multiple locations

Intermittent demand (also known as sporadic demand) comes about when a product experiences several periods of zero demand. Often in these situations, when demand occurs it is small, and sometimes highly variable in size. In inventory control, the level of safety stock depends on the service level you select, on the replenishment lead time as well as the reliability of the forecast.

Step 1. Data Exploration as an Essential Quality Check in the Forecasting Process

Consumer demand-driven historical data, in the retail industry for example, are characterized to a large extent by trends (consumer demographics, business cycles) and seasonal patterns (consumer habits: economics). Before making any specific modeling assumptions, demand forecasters and planners should first examine data quality in the context of the new environment along with an exploratory data analysis (EDA) examination of the changing data characteristics and quality in demand history. if you want to achieve greater agility in forecasting for the entire supply chain, this preparatory step can be time consuming but is an essential undertaking, 

Looking for insight into data quality in this spreadsheet, you see that an expected seasonal peak for December 2016 seems to appear in the following month (Jan 2017). Interchanging the December with the January value appears to have a significant impact on the underlying seasonal variation, including uncertainty and other factors (compare ‘before’ in column 2 with ‘after’ in column 3 in the table below).

Also, May 2017 appears unusually weak, but we would call on domain experts to advise the forecaster on that unusual event. In any case, from a preliminary examination into the quality of the data, we see that consumer habit (an economic factor) may constitute about two-thirds of the total variation in the demand history (Above result obtained with a Two-way ANOVA w/o replication algorithm).

First Takeaway: “Bad data will beat a good forecaster every time” (Paraphrasing W. Edward Deming)

• Embrace Change & Chance by improving data quality through exploratory data analysis(EDA) as a preliminary step essential in creating agility in the demand forecasting and planning process.

For intermittent demand forecasting, on the other hand, it is necessary to examine the nature of the interdemand intervals and its relation to the distribution of non-zero demand sizes over a specified lead-time. The widely used Croston methods will under closer scrutiny point to flawed assumptions about the independence of zero intervals and demand volumes. Assumptions about Interval sizes and nonzero demand volumes should be reconciled in practice with real data, not just with simulated data.

Characterizing intermittent data with SIB models differs fundamentally from Croston-based methods in that a dependence of interval sizes on demand can be made explicit and can be validated with real data. The SIB modeling approach does not assume that intervals and demand sizes are independent or generated by conventional data-generating models.

The evidence of the dependence on interval durations for demand size can be explored by examining the frequency distribution of ‘lag time’ durations. I define a ”Lag-time” Zero­­ Interval LZI as the zero-interval duration preceding a nonzero demand size. In my previous LinkedIn article on SIB models for intermittent and regular demand forecasting, I focused on measurement error models for forecasting and evaluating intermittent demand volumes based on a dependence of an interval duration distribution. 

A Structured Inference Base (SIB) Model for Lead Time Demand

Step 2. Leadtime Demand Forecasting Needs a Specified Lead-time (Time Horizon) For lead-time demand forecasting with a fixed horizon, a location scale measurement error model can be created for assessing the effectiveness and accuracy of the demand forecasting process. The Profile Forecast Error (FPE) dataused in modeling the accuracy and performance of lead-time demand forecasts can be represented as the output of a measurement model: FPE = β + σ ɛ, in which β and σ are unknown parameters,  and the input ɛ is a measurement error with a known or assumed non-Gaussian distribution.

Keeping in mind the pervasive presence of outliers and unusual values in today’s real-world forecasting environment, I will henceforth shy away from the conventional normal distribution assumptions.  Rather, for error distribution assumptions, I will be referring to a flexible family of distributions, known as theExponential family. This is a rich family of distributions particularly suited for SIB modeling; it contains many familiar distributions including the normal (Gaussian) distribution, as well as distributions with thicker tails and skewness, so much more appropriate in today’s disruptive forecasting environment. These are some of the reasons I regard this article as foundational. However, by following a step-by-step development using a real-world data example and nothing more than arithmetic, some algebra and the logarithm, I hope that you can follow the process.

The SIB modeling approach is algorithmic and data-driven, in contrast to conventional data-generating models with normality assumptions. The measurement model for Forecast Profile Error (or Miss) = β + σ ɛ is known as a location scale measurement model because of its structure. The Forecast Profile Error (FPE) model shows that the FPE data result from a translation and scaling of an input measurement error ɛ. For the three forecasting approaches (judgment, method and model), used as examples, I display the FPE data below.

The forecast profile errors in the spreadsheet are calculated with the formula where a(i) are the components in the actual alphabet profile (AAP) and f(i) are the components in the forecast alphabet profile (FAP). These have been explained before in several previous articles available on my LinkedIn Profile and the my Delphus website blog.

The sums of the rows can be interpreted as a measure of information or ignorance about the forecast profile error. The closer to zero the better and the sign indicates over or underforecasting. In this example, the techniques are overforecasting. The units are known as ‘nats’ (for natural logarithms).

What Can Be Learned About the Measurement Process given the Forecast Profile Errors and the Observed Data?

Step 3. Setting Up the Model with Real Data

We can analyze the performance of a SIB model for lead-time forecasts as follows. In practice, we have multiple measurements of observed forecast profile errors (over a time horizon m = 12 in the spreadsheet example):

and where ɛ = {ɛ1, ɛ2, ɛ3, . . . ɛ12} are now 12 realizations of measurement errors from an assumed distribution in the Exponential family.

Step 4. A Critical Data Reduction Step

What information can we uncover about the forecasting process? Like a detective, we can explore a SIB model and find that, based on the observed data, there is a clue revealed now about the unknown, but realized measurement errors ɛ. This is evidence that will guide us to the next important SIB modeling step: It points to a decomposition of the measurement error distribution into two components: (1) a marginal distribution for the observed components and (2) a conditional distribution (based on the observed components) for the remaining unknown measurement error distribution, which depends on the parameters β and σ. So, what are these observed components of the error distribution that we uncover?

The insight or essential information is gleaned from the structure of the model and the recorded forecast profile error data. If we now select a suitable location measure m(.), and a scale measure s(.), we can make a calculation that yields important observables about the measurement process for each forecasting technique used. The SIB model shows, with some elementary algebraic manipulations, that the observables can be expressed in equations like this (using a 12-month lead-time):

That is, letting d = (d1, d2, … , d12) represent the left hand-side and right-hand side equations,we can reduce the SIB model to only two equations with two unknown variables m(ɛ) and s(ɛ) that represent the remaining unknown information in the measurement error model

Step 5. Conditioning on What You Know to be True

We do not need to go any further with details at this point. The conditional distribution (given the known d = (d1, d2, … , d12) for the variables m(ɛ) and s(ɛ) can be derived from the second equation using an assumed distribution for ɛ.

Using the selected location measure m(.), and a scale measure s(.), we can make a calculation that yields important observables about the measurement process for each method or model. If I select m(.) = L-Skill score as the location measure, the calculated L-Skill scores for the three approaches are obtained by summing the values on each row.

The L-Skill scores can range over the whole real line, the smaller the score in absolute value the better.Like the forecast profile accuracy measure D(a|f), the L-Skill score is zero if the actual and forecast alphabet profiles are identical. But, unlike the accuracy measure D(a|f), it can be shown that the L-Skill score will always be zero for constant level profiles, like MAVG-12.

The s(.) = Profile Accuracy D(a|f) divergences for the scale measure are found to be

We note that the 12 left-hand side equations named d(FPE) are equal to the right-hand side equations  d(ɛ). In the spreadsheet example for the 12 Lead-time demand values for Part #0174, d = d(FPE) = d(ɛ) is shown in the table. Then d = (d1, d2, … , d12), where we can now calculate

Step 6. Deriving Posterior Distributions for the Unknown Parameters β and σ Along with Posterior Prediction Limits and Likelihood Analyses

The derivations and formulae can be found in U of Toronto Professor D.A.S. Fraser’s 1979 book Inference and Linear Models, Chapter 2, and in his peer-reviewed journal articles dealing with statistical inference and likelihood methods. These are not mainstream results in modern statistical literature, but that does not diminish their value in practice.

Statistical inference refers to the theory, methods, and practice of forming judgments about the parameters of a population and the reliability of statistical relationships

These algorithms can be readily implemented in a modern computing environment, which was not the case more than four decades ago when I was first exposed to them. With normal (Gaussian) error distribution assumptions, there are closed form solutions (i.e. solvable in a mathematical analysis) that have a semblance to more familiar Bayesian inference solutions.

Step 7. Application to Inventory Planning

The SIB inferential analysis will yield a posterior distribution (conditional on the observed d) for the unknown parameters β and σ from which we can derive unique confidence bounds for β/σ (L-Skill score) and σ (Forecast Profile Accuracy). These confidence bounds will give us the service levels we require to set desired level of safety stock.

1.      For the ETS(A,A,M) model and data example, the reduced error distribution for location measure m(ɛ) and scale measure s(ɛ) is conditional on observed d = (d1, d2, … , d12):

·       Location component:            m(ɛ) = [m(FPE) – β]/ σ = [0.001 – β]/ σ

·       Scale component:                  s(ɛ) = s(FPE)/ σ = 0.044/ σ

2.      Define “Safety Factor” SF = √12 m(ɛ) /s(ɛ) = √12 {0.001– β 0}/ 0.044, where β 0 = max β under a selected contour boundary

Then,   β 0 = 0.001+ SF * 0.044/ √12 Is the desired level of safety stock for the service level you select.

Final Takeaway:

A data-driven non-Gaussian SIB modeling approach for lead-time demand forecasting is based on

• No independence assumptions made on demand size and demand interval distributions
• Non-normal (non-Gaussian) distributions assumed throughout the inference process.
• Using bootstrap resampling (MCMC) from a SIB model if the underlying empirical distributions are to be used
• Deriving posterior lead-time demand distributions with bootstrap sampling techniques
• Determining upper tail percentiles of a lead-time demand distribution from family of exponential distributions for lead-time demands or from empirical distributions

Hans Levenbach, PhD is Owner/CEO of Delphus, Inc and Executive Director, CPDF Professional Development Training and Certification Programs. Dr. Hans is the author of a new book (Change&Chance Embraced) on Demand Forecasting in the Supply Chain and created and conducts hands-on Professional Development Workshops on Demand Forecasting and Planning for multi-national supply chain companies worldwide. Hans is a Past President, Treasurer and former member of the Board of Directors of the International Institute of Forecasters.