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)), demand, color='lightgrey', label='Demand')
plt.plot(CFE, label='CFE')
plt.plot(PIS, label='PIS')
plt.axhline(0, lw=1, color='black')

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)),n), demand, color='lightgrey')
plt.plot(demand_rate, label='Demand Rate')

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)), ts, color='lightgrey')
plt.xlabel('Time index')

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
    ts : (N,) array_like
        1-D input array
    alpha : float
        Smoothing factor, `0 < alpha < 1`, default = 0.1
    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
            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)


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")

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), ts, color='lightgrey')

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)), 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}')

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.title('Cumulative Forecasting Error')

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.


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')

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), ts, color='lightgrey')


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.


Accuracy metrics for different values of alpha.

0.00 6.619574 117.737804 33.752689 219.360000 222.720000 20.640000 47 -5921.640000
0.05 7.251768 106.405893 10.639970 85.954346 112.280559 18.878666 47 -3341.478878
0.10 7.605402 107.205868 5.956724 39.335916 76.681052 14.980889 47 -2124.001490
0.15 7.779682 109.055526 4.648616 19.791792 57.684920 1.897217 47 -1471.296153
0.20 7.892810 110.959703 4.765412 10.188768 49.825913 -5.575011 43 -1080.649959
0.25 7.999936 112.866346 5.763382 4.767469 45.481206 -9.864628 41 -824.777258
0.30 8.083558 114.841793 7.393092 1.290773 42.830600 -12.350334 40 -644.604293
0.35 8.147493 116.928055 9.512085 -1.169582 40.743894 -13.791315 36 -510.876608
0.40 8.191765 119.127456 12.033105 -2.997149 38.893421 -14.600775 36 -408.532700
0.45 8.215478 121.421458 14.901910 -4.345081 37.250558 -15.004828 33 -329.528871
0.50 8.218292 123.792666 18.085363 -5.277029 35.785136 -15.637853 33 -269.266190
0.55 8.227726 126.239831 21.565192 -5.830786 34.466570 -16.121885 29 -224.748918
0.60 8.282826 128.785547 25.336396 -6.046069 33.264806 -16.197522 28 -193.635146
0.65 8.333973 131.480332 29.410686 -5.976651 32.236909 -15.898489 26 -173.748405
0.70 8.414663 134.407264 33.825407 -5.697197 32.437841 -15.269306 27 -162.811487
0.75 8.485696 137.690303 38.657840 -5.310121 32.643048 -14.370995 26 -158.269206
0.80 8.548356 141.508910 44.045148 -4.955429 32.827849 -13.312039 26 -157.124730
0.85 8.622289 146.123005 50.212600 -4.826025 32.954745 -12.952348 23 -155.738630
0.90 8.734745 151.916982 57.518128 -5.192089 32.963934 -12.682085 24 -149.535388
0.95 8.858563 159.482068 66.532620 -6.441564 32.758937 -12.551001 21 -132.516604
1.00 9.081844 169.782741 78.202335 -9.153333 32.180000 -12.620000 20 -96.326667


On Intermittent Demand Model Optimisation and Selection, 2014