Error Metrics for Intermittent Demand (CFE, PIS, MSR)
Intermittent demand is present in many retail settings and when forecasting stock requirements, businesses must strike a balance between the cost of goods and loss of sales due to a lack of stock. When training a forecasting model, using common accuracy measures, such as the MAE or MSE, do not always translate to ideal real-world outcomes. Due to the intermittence, forecasts at or near zero may reduce the error but would result in shelves being empty.
We will explore a new range of error metrics targeted towards intermittent demand, such as the Cumulative Forecasting Error (CFE), Periods in Stock (PIS) and Mean Squared Rate (MSR). Using Croston's method to create a forecast, we will see how each measure reacts as we adjust the smoothing factor, alpha.
Cumulative Forecasting Error (CFE)
The CFE is the sum of the difference between a forecast and the observed demand.
$$CFE = \sum_{i=1}^t (y_t - \hat{y_t})$$For example, in the table below we have 7 days of sales with a given forecast of 2 units per day. The error is calculated each day, then the cumulative sum is taken for the period.
| Day (t) | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|
| Demand ($y$) | 2 | 0 | 0 | 5 | 0 | 3 | 5 |
| Forecast ($\hat{y}$) | 2 | 2 | 2 | 2 | 2 | 2 | 2 |
| Error ($y - \hat{y}$) | 0 | -2 | -2 | 3 | -2 | 1 | 3 |
| CFE ($\sum Error$) | 0 | -2 | -4 | -1 | -3 | -2 | 1 |
At the end of this 7-day period, we have a CFE of 1.
A positive CFE means there was an under-estimation of the required stock during the period. A negative number means an over-estimation. You could have wild swings in the CFE but end up back at 0 by the end of the period. You would look at this zero error and think - wow.. this is a great model! But in actual fact there may have been periods where you were vastly under-stocked or over-stocked. One way to balance this would be to track CFEmin and CFEmax, which we will investigate later.
Number of Shortages
Along with tracking the CFEmin and CFEmax, we can also look at the Number of Shortages (NOS). The NOS is how many days that we were under-stocked throughout the period, which can be calculated by counting the number of positive CFE occurrences.
In the example above, there was only one day with a positive CFE, so the NOS would be equal to 1.
Periods in Stock (PIS)
Instead of looking at the CFE along with its min and max values, we can calculate the Periods in Stock. This is the accumulation of the CFE, which will give insight into what is happening throughout the period. It the negative cumulative CFE, given by
$$PIS = -\sum_{i=1}^t CFE$$Applying this to the example from before, we get
| Day (t) | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|
| CFE | 0 | -2 | -4 | -1 | -3 | -2 | 1 |
| PIS | 0 | 2 | 6 | 7 | 10 | 12 | 11 |
A positive PIS over a given period is an indication of being over-stocked. In this case, we have a PIS of 11.
You can think of a fictitous store preparing stock each day. In the example above, they get two units ready for sale on day one and manage to sell both, resulting in a PIS of 0. Perfect. On the second day they prepare another two units but this time there are no sales, so the PIS grows by two. On the third day, they again prepare another two units and again there are no sales. So they have have a total of four units unsold that day, which brings the PIS up to six.
To visualise the CFE and and PIS, we can do a quick plot using Numpy and Matplotlib.
import numpy as np
import matplotlib.pyplot as plt
# Example
demand = np.array([2,0,0,5,0,3,5])
forecast = np.array([2] * 7)
CFE = np.cumsum(demand - forecast)
PIS = -np.cumsum(CFE)
# Plot the results
fig, ax = plt.subplots(figsize=(9,5))
plt.bar(np.arange(7), demand, color='lightgrey', label='Demand')
plt.plot(CFE, label='CFE')
plt.plot(PIS, label='PIS')
plt.axhline(0, lw=1, color='black')
plt.legend()
plt.show()
While the CFE ends up back near zero, the PIS has grown due to being over-stocked for most of the period.
Mean Squared Error and the Demand Rate
The Mean Squared Error is often used for cost functions, but when it comes to intermittent data it can prove to be problematic. Due to all the zeros in the data, the error might be minimised by simply casting a forecast of zero, or one close to zero. On paper it looks great, we've reduced the error. But in the real world it would translate to empty shelves. And you can't sell what you don't have.
When making a forecast with Croston's method or one of it's adaptations, we are left with a single value which represents the demand rate. For example, a forecast of 0.5 doesn't mean we predict a sale of half a unit, but perhaps 5 units over the next 10 days.
So instead of using the actual demand, we can create a new series for the demand rate.
# Example of demand rate
demand = np.array([5,0,2,1,0,0,0,3,0,1,3,0,4])
n = len(demand)
demand_rate = np.cumsum(demand) / np.arange(1,n+1)
fig, ax = plt.subplots(figsize=(9,4))
plt.bar(np.arange(0,n), demand, color='lightgrey')
plt.plot(demand_rate, label='Demand Rate')
plt.legend()
plt.show()
Note that the value for the demand rate at t=0 will be equivalent to the initial demand. This can be problematic, so when calculating the MSR the first few values of the demand rate/forecast can often be trimmed. This is implemented below by ignoring the first 10% of values.
Parameter Selection for Croston's Method
To explore each of these error metrics, we will use an example time series from the M5 dataset. I've chosen row index 7 and the final 50 days.
import pandas as pd
# Read in M5 dataset, select the final 100 days of row index 5
sales = pd.read_csv('data/sales_train_validation.csv')
ts = sales.iloc[7,-50:].values.astype('int')
# Plot sales series
plt.figure(figsize=(12,4))
plt.bar(np.arange(50), ts, color='lightgrey')
plt.xlabel('Time index')
plt.ylabel('Sales')
plt.show()
Croston's method will be performed on this series to produce a forecast, the code has been taken from the Croston's Method tutorial.
def croston(ts, alpha=0.1):
"""
Perform Croston's method on a time series, ts, and return
a forecast
Parameters
----------
ts : (N,) array_like
1-D input array
alpha : float
Smoothing factor, `0 < alpha < 1`, default = 0.1
Returns
-------
forecast : (N+1,) ndarray
1-D array of forecasted values
"""
# Initialise arrays for demand, z, and period, p. Starting
# demand is first non-zero demand value, starting period is
# mean of all demand intervals
ts_trim = np.trim_zeros(ts, 'f')
n = len(ts_trim)
z = np.zeros(n)
p = np.zeros(n)
p_idx = np.flatnonzero(ts)
p_diff = np.diff(p_idx, prepend=-1)
p[0] = np.mean(p_diff)
z[0] = ts[p_idx[0]]
q = 1
for i in range(1,n):
if ts_trim[i] > 0:
z[i] = alpha*ts_trim[i] + (1-alpha)*z[i-1]
p[i] = alpha*q + (1-alpha)*p[i-1]
q = 1
else:
z[i] = z[i-1]
p[i] = p[i-1]
q += 1
f = z / p
nan_arr = [np.nan] * (len(ts)-n+1)
return np.concatenate((nan_arr, f))
The error() function below will take the original time series and the forecast and return a dictionary of the calculated error metrics.
from scipy.spatial.distance import sqeuclidean
def error(ts, forecast):
"""
Return a dictionary of error metrics
"""
cfe = np.cumsum(ts - forecast)
cfe_max_idx = np.abs(cfe).argmax()
pis = -np.cumsum(cfe)
n = len(ts)
d_rate = np.cumsum(ts) / np.arange(1,n+1) # Demand rate
trim = int(n*0.1) # Number to trim from start of demand rate, ie 10%
err = {
'MAE': np.sum(np.abs(ts - forecast))/n,
'MSE': sqeuclidean(ts, forecast)/n,
'MSR': sqeuclidean(d_rate[trim:], forecast[trim:])/(n-trim),
'CFE': cfe[-1],
'CFE_MAX': np.max(cfe),
'CFE_MIN': np.min(cfe),
'NOS': len(cfe[cfe>0]),
'PIS': pis[-1],
}
return err
Iterating Through Alpha
We will test different values for alpha, creating a forecast and then returning the dictionary of error metrics each time. We'll use values between 0 and 1 in increments of 0.05. The results will then be stored in a dataframe, df. The full table of results can be seen at the bottom of the tutorial in the appendix.
import pandas as pd
alphas = np.arange(0,1.05,0.05)
df = pd.DataFrame([])
for a in alphas:
forecast = croston(ts, alpha=a)
idx = np.argmax(np.isfinite(forecast)) # Index of first non-NaN value
err = error(ts[idx:], forecast[idx:-1])
err = pd.DataFrame(err, index=[a])
df = pd.concat([df, err], axis=0)
MAE, MSE and MSR
The first metrics we will plot are the Mean Absolute Error, Mean Squared Error and Mean Squared Rate. The plots will show the error as a function of alpha.
metrics = ['MAE', 'MSE', 'MSR']
fig, axs = plt.subplots(nrows=3, figsize=(9,8), sharex=True)
for i, m in enumerate(metrics):
df[m].plot(ax=axs[i], legend=m)
axs[i].legend(loc="upper left")
fig.supxlabel('Alpha')
fig.supylabel('Error')
plt.show()
We have three different error curves. The MAE shows a minimum when $\alpha=0$, the MSE has its minimum at $\alpha=0.05$ and the MSR at $\alpha=0.15$.
Let's plot the forecasts using the optimal level found to see if we can figure out what is going on here.
fig, ax = plt.subplots(figsize=(9,5))
for m in metrics:
alpha_min = df[m].abs().idxmin()
forecast = croston(ts, alpha_min)
l = f'alpha={alpha_min:.2f}, selected by {m}'
plt.plot(forecast, label=l)
plt.bar(np.arange(len(ts)), ts, color='lightgrey')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Demand')
plt.show()
It's been established that using the MAE or MSE as the error metric for optimising alpha can result in a bias towards a zero forecast. We can see that reducing the MAE results in a forecast closer to zero, which may lead to be under-stocked. We will see if this is the case later when looking at the CFE and PIS.
When trying to minimise the MSR, we are fitting a forecast to the demand rate. We can plot do a plot of the forecast against the demand rate to see how it fit.
# Construct demand rate array
d_rate = np.cumsum(ts) / np.arange(1,len(ts)+1)
trim = int(len(ts)*0.1) # Trim % from the start of d_rate
# Plot results
alpha_min = df['MSR'].abs().idxmin()
fig, ax = plt.subplots(figsize=(9,5))
plt.bar(np.arange(len(ts)), ts, color='lightgrey')
plt.plot(np.arange(trim,len(ts)), d_rate[trim:],
label='Demand Rate')
plt.plot(croston(ts, alpha_min), linestyle='--',
label=f'Forecast, alpha={alpha_min:.2f}')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Demand')
plt.show()
CFE, CFEmin and CFEmax
We'll now look at the Cumulative Forecasting Error, to gauge whether we were over- or under-stocked at the end of the period.
fig, axs = plt.subplots(figsize=(9,5))
plt.plot(df['CFE'], label='CFE')
plt.plot(df['CFE_MIN'], linestyle='--', label='CFE_MIN')
plt.plot(df['CFE_MAX'], linestyle='--', label='CFE_MAX')
plt.axhline(linewidth=0.5, color='black')
plt.legend()
plt.title('Cumulative Forecasting Error')
plt.xlabel('Alpha')
plt.ylabel('Error')
plt.show()
Recall that a CFE of zero means that at the end of the period the sum of the stock forecasted was equal to the actual demand. In the above plot, the CFE reaches an absolute minimum around $\alpha=0.35$. When we tried to minimise the MAE, we got an optimal value of $\alpha=0.0$. Looking at the CFE, this would have left us severley under-stocked as that is when the CFE is at it's maximum there.
A CFE of zero does not always indicate the ideal level of stock, as there can be instances of over- or under-stocking during that period, which is why the CFEmin and CFEmax are also tracked. The CFEmax was always positive, which means that regardless of the value of alpha, there were periods of being under-stocked.
PIS and NOS
We'll now have a look at the final metrics, PIS and NOS.
metrics = ['PIS', 'NOS']
fig, axs = plt.subplots(nrows=2, figsize=(9,8))
for i, m in enumerate(metrics):
axs[i].plot(df[m], label=m)
axs[i].axhline(linewidth=0.5, color='black')
axs[i].legend()
plt.xlabel('Alpha')
fig.supylabel('Error')
plt.show()
The PIS had it's absolute minimum at $\alpha=1$. When the PIS is negative, it indicates being under-stocked. The NOS fell as alpha increased, and also had its minimum at $\alpha=1$. But since the NOS was greater than 1, it shows that we still recorded stock shortages.
We can plot the forecast using the optimal value of alpha found by reducing the CFE and PIS.
metrics = ['CFE', 'PIS']
fig, ax = plt.subplots(figsize=(9,5))
for m in metrics:
alpha_min = df[m].abs().idxmin()
forecast = croston(ts, alpha_min)
l = f'alpha={alpha_min:.2f}, selected by {m}'
plt.plot(forecast, label=l)
plt.bar(np.arange(len(ts)), ts, color='lightgrey')
plt.legend()
plt.xlabel('Time')
plt.ylabel('Demand')
plt.show()
Conclusion
While the MAE and MSE suggested using low value for $\alpha$, when we looked at the CFE and PIS we discovered that this would leave us severley under-stocked. As the value of $\alpha$ increased, the CFE improved and reached zero at $\alpha=0.35$. However, looking at the PIS we could see we would still have been under-stocked. Minimising the PIS suggests setting $\alpha=1$, but even at this level we would still have experienced stock-shortages, as indicated by the NOS.
As mentioned in the introduction, forecasting intermittent demand is nuanced. There are many trade-offs in selecting an error metric. But I hope this has provided more insight and a few more tools for making a better decision.
Appendix
Accuracy metrics for different values of alpha.
df
References
On Intermittent Demand Model Optimisation and Selection, 2014 https://kourentzes.com/forecasting/2014/06/11/on-intermittent-demand-model-optimisation-and-selection/