Time series

In most machine learning you can shuffle the rows. In time series you cannot: order is the signal, and it is also what makes ordinary validation lie to you. Split a series at random and information from the future leaks into training, so the model scores beautifully in testing and falls apart the first time it forecasts something real. This section covers the structure worth modelling — trend, seasonality, autocorrelation — validation that respects time, and the naive baseline that a surprising number of forecasting projects never beat and never check against.

The structure#

• Trend: the long-run direction.
• Seasonality: repeating patterns at fixed periods, daily, weekly, yearly.
• Autocorrelation: today resembles yesterday. This is what makes forecasting possible and also what makes naive validation wrong.

Autocorrelation is the one to sit with, because both halves of that sentence follow from the same fact. If today resembles yesterday, then yesterday is a useful feature, which is why forecasting works at all. And if today resembles yesterday, then a validation row whose neighbours are in training is not a test of anything: the model can answer it by remembering the days on either side, which is a skill it will never be asked for again.

The mistake that invalidates everything#

The size of that effect is worth measuring rather than believing. Below is a series whose level drifts — a random walk, which is what most business series look like once seasonality is removed. The model is deliberately trivial: predict a day using the nearest day it saw in training. That is roughly what any flexible model does on a series with no other features, so it is a fair stand-in.

import numpy as np

rng = np.random.default_rng(0)
n = 500
level = np.cumsum(rng.normal(0, 1, n))          # a drifting level
y = level + rng.normal(0, 0.5, n)

def nearest_day(train_idx, test_idx):
    """Predict each test day with the value of the closest day in training."""
    out = np.empty(len(test_idx))
    for j, i in enumerate(test_idx):
        k = train_idx[np.argmin(np.abs(train_idx - i))]
        out[j] = y[k]
    return out

def rmse(a, b):
    return float(np.sqrt(np.mean((a - b) ** 2)))

perm = rng.permutation(n)
rand_test, rand_train = np.sort(perm[:100]), np.sort(perm[100:])
print("random split", round(rmse(y[rand_test], nearest_day(rand_train, rand_test)), 2))

time_train, time_test = np.arange(400), np.arange(400, 500)
print("time split  ", round(rmse(y[time_test], nearest_day(time_train, time_test)), 2))

# random split 1.33
# time split   5.47

from sklearn.model_selection import TimeSeriesSplit, cross_val_score

# Wrong
cross_val_score(model, X, y, cv=5)

# Right: each fold trains on the past and validates on the future
cross_val_score(model, X, y, cv=TimeSeriesSplit(n_splits=5))

fold 1  train [-------]  test [--]
fold 2  train [----------]  test [--]
fold 3  train [-------------]  test [--]
fold 4  train [----------------]  test [--]

Each fold trains only on what came before its test window, and the training
set grows. This is forward chaining, and it is what TimeSeriesSplit does.

Features#

df["lag_1"] = df.value.shift(1)
df["lag_7"] = df.value.shift(7)
df["rolling_7"] = df.value.shift(1).rolling(7).mean()   # shift FIRST
df["dow"] = df.date.dt.dayofweek
df["is_holiday"] = df.date.isin(holidays)

That caption is easy to nod at and easy to get wrong at half past six on a Friday, so here are the two columns side by side. The series has one spike, and on the spike day the leaky feature has already noticed.

import numpy as np

values = np.array([10.0, 12, 11, 40, 13, 12, 11])

def rolling_including_today(v, w=3):
    return [round(float(v[i - w + 1 : i + 1].mean()), 1) if i >= w - 1 else None
            for i in range(len(v))]

def rolling_shifted(v, w=3):
    return [round(float(v[i - w : i].mean()), 1) if i >= w else None
            for i in range(len(v))]

print("value     ", [float(x) for x in values])
print("leaky     ", rolling_including_today(values))
print("shifted   ", rolling_shifted(values))

# value      [10.0, 12.0, 11.0, 40.0, 13.0, 12.0, 11.0]
# leaky      [None, None, 11.0, 21.0, 21.3, 21.7, 12.0]
# shifted    [None, None, None, 11.0, 21.0, 21.3, 21.7]

Which lags to build is a question with a cheap empirical answer: compute the correlation between the series and itself shifted by k, for k from one to about three seasons, and take the peaks. Daily retail data will show a spike at 7 and a smaller one at 14 and 21. Hourly traffic shows 24 and 168. Build lags at those, plus rolling means at a few windows to give the model a smoothed level as well as individual days.

• Lags: the value k steps ago. Cheap, and the single most useful family of features on almost any series.
• Rolling statistics: mean, standard deviation, min and max over a window ending before today. The standard deviation is the one people forget, and it is what lets a model know the series has become volatile.
• Differences: value minus the value a season ago. This removes the level, which matters because a tree model cannot extrapolate beyond the range of levels it was trained on.
• Calendar: day of week, day of month, month, holiday, days to the next holiday, and any known future event such as a planned promotion. These are the only features you are certain to have at forecast time.

Seasonality#

Seasonality is the part of the series that repeats on a known period, and the useful thing about it is that it is knowable in advance. Separate it from the trend before you model anything, because what looks like a model failing on Saturdays is usually a weekly pattern nobody gave it.

import numpy as np

rng = np.random.default_rng(2)
n = 364
t = np.arange(n)
dow = t % 7
weekly = np.array([0, 2, 3, 4, 8, 20, 15], dtype=float)      # the truth, unobserved
series = 100 + 0.05 * t + weekly[dow] + rng.normal(0, 2, n)

trend = np.convolve(series, np.ones(7) / 7, mode="same")     # a centred 7-day mean
resid = series - trend
profile = np.array([resid[dow == d].mean() for d in range(7)])

print(profile.round(2))
print((weekly - weekly.mean()).round(2))

# [-6.46 -4.87 -4.75 -3.73  1.5  13.54  8.63]
# [-7.43 -5.43 -4.43 -3.43  0.57 12.57  7.57]

Two ways to hand that to a model. Dummy variables, one per day of week or per month, which are exact and cost a column each — fine for seven days, awkward for 365 days of the year. Or Fourier terms, a handful of sines and cosines at the seasonal frequency, which give a smooth periodic shape in four or six columns and are the standard answer for yearly seasonality on daily data.

import numpy as np

def fourier_terms(t, period, k):
    """2k columns capturing a smooth repeating pattern of the given period."""
    out = []
    for i in range(1, k + 1):
        out.append(np.sin(2 * np.pi * i * t / period))
        out.append(np.cos(2 * np.pi * i * t / period))
    return np.column_stack(out)

t = np.arange(730)
X_year = fourier_terms(t, period=365.25, k=3)       # 6 columns for yearly shape
X_week = fourier_terms(t, period=7, k=2)            # 4 columns for weekly shape
print(X_year.shape, X_week.shape)                   # (730, 6) (730, 4)

Beat the naive baseline first#

Predicting that tomorrow equals today, or equals the same day last week, is often surprisingly strong. Many forecasting projects produce models that never beat it, and nobody checked. Compute it before anything else.

import numpy as np

def naive(history, horizon):
    """Tomorrow equals today. The hardest baseline on a random walk."""
    return np.full(horizon, history[-1])

def seasonal_naive(history, horizon, season=7):
    """Each future day equals the same weekday last week."""
    last = history[-season:]
    return np.array([last[i % season] for i in range(horizon)])

def drift(history, horizon):
    """Extend the straight line from the first to the last observation."""
    slope = (history[-1] - history[0]) / (len(history) - 1)
    return history[-1] + slope * np.arange(1, horizon + 1)

h = np.array([100.0, 104, 99, 103, 108, 130, 120, 102, 106, 101, 105, 110, 132, 122])
print(naive(h, 3).round(1))             # [122. 122. 122.]
print(seasonal_naive(h, 3).round(1))    # [102. 106. 101.]
print(drift(h, 3).round(1))             # [123.7 125.4 127.1]

You know how to validate a forecast honestly. This is what to put inside it. The classical methods, exponential smoothing and ARIMA, are fast, interpretable and genuinely hard to beat on many business series, so they are a baseline rather than a formality. The alternative is to reshape the problem into ordinary supervised learning with lag features and fit a tree, which handles many related series and external variables easily. The section also covers the horizon, since accuracy one day ahead says little about thirty, and giving an interval rather than a single number nobody should trust.

What you have                                  Start with
one series, a few hundred points               exponential smoothing or ARIMA
one series, external drivers you know ahead    gradient boosting on lag features
thousands of related series                    one global gradient boosting model
thousands of series, years of history, a team  a deep global model, and prove it
demand that is zero most days                  an intermittent-demand method
a series you cannot explain at all             the seasonal naive baseline, honestly

Classical models#

Exponential smoothing weights recent observations more heavily and handles trend and seasonality. ARIMA models a series through its own past values and past errors. Both are strong, fast and interpretable baselines, and on many business series they are hard to beat.

Simple exponential smoothing is one line, and writing it out shows exactly what the single parameter does. The level is a running blend of what just happened and what you already believed.

import numpy as np

def ses(series, alpha):
    """level_t = alpha * observation + (1 - alpha) * level_{t-1}"""
    level = series[0]
    out = [level]
    for x in series[1:]:
        level = alpha * x + (1 - alpha) * level
        out.append(level)
    return np.array(out)

s = np.array([10.0, 12, 11, 40, 13, 12, 11])     # a series with one spike
print(ses(s, 0.3).round(2))
print(ses(s, 0.8).round(2))

# [10.   10.6  10.72 19.5  17.55 15.89 14.42]   alpha 0.3: slow, ignores the spike
# [10.   11.6  11.12 34.22 17.24 13.05 11.41]   alpha 0.8: fast, chases the spike

The fuller versions add a trend term and a seasonal term with their own smoothing parameters, which is Holt-Winters. ARIMA arrives at similar forecasts from a different direction: differencing the series until it is stationary, then modelling what is left as a combination of its own past values and its own past errors. The practical point is not which is theoretically nicer but that both fit in milliseconds, both produce prediction intervals for free from their error model, and both are hard to beat on a single business series with a few hundred points.

Treat it as supervised learning#

Build lag and rolling features, add calendar features, and fit a gradient-boosted tree. This handles many related series at once and incorporates external variables easily, which is why it is common in practice.

date        sales     ->   lag_1  lag_7  roll_7  dow  holiday  promo   y
2026-03-08   120            118    115    117.4    6     0        0     120
2026-03-09   118            120    119    117.9    0     0        1     118
2026-03-10   131            118    122    118.6    1     0        1     131

One row per date, every feature known before that date. Now it is a table,
and everything you know about tabular modelling applies.

Gradient boosting on that table is the honest default for applied work, for four reasons. It takes external variables — price, promotion, weather, a competitor's opening — which classical models handle awkwardly. It fits thousands of series in one model, so a product with six weeks of history borrows the shape of everything similar. It handles non-linear effects and interactions without being told they exist. And it is fast enough to backtest properly, which matters more than any single point of accuracy, because a model you can only evaluate once is a model you cannot trust.

import numpy as np
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 300
t = np.arange(n).astype(float)
y = 50 + 0.4 * t + 6 * np.sin(2 * np.pi * t / 7) + rng.normal(0, 3, n)

X = np.column_stack([t % 7, np.sin(2 * np.pi * t / 7), np.cos(2 * np.pi * t / 7)])
train, test = slice(0, 250), slice(250, 300)

def rmse(a, b):
    return round(float(np.sqrt(np.mean((a - b) ** 2))), 2)

gbm = GradientBoostingRegressor(random_state=0).fit(X[train], y[train])
print("GBM on the level    ", rmse(y[test], gbm.predict(X[test])))

trend = LinearRegression().fit(t[train].reshape(-1, 1), y[train])
resid = y - trend.predict(t.reshape(-1, 1))
gbm2 = GradientBoostingRegressor(random_state=0).fit(X[train], resid[train])
pred = trend.predict(t[test].reshape(-1, 1)) + gbm2.predict(X[test])
print("GBM on the residual ", rmse(y[test], pred))
print("the trend alone     ", rmse(y[test], trend.predict(t[test].reshape(-1, 1))))

# GBM on the level     59.81
# GBM on the residual   3.11
# the trend alone       4.93

When a deep model earns its place#

Not on one series with three years of monthly data. The forecasting competitions that made global deep models credible were run on tens of thousands of related series with years of history each, and that is the setting where the extra capacity has something to learn from. Below roughly that scale, a deep model is a slower, less interpretable way of arriving at a worse number.

• Many related series, so the model can learn shape from thousands of examples rather than one. This is the main condition, and it is the one people skip.
• Inputs that do not fit in a table: raw text about a product, an image of a store shelf, a long history where the relevant pattern is a hundred steps back.
• A need for the whole horizon at once. Sequence models output all thirty days in one pass, which avoids the compounding errors of feeding forecasts back in.
• Enough traffic that a point of accuracy is worth a retraining pipeline, a GPU and someone who understands it at three in the morning.

The bar to clear is a tuned gradient boosting model on the same features, backtested over the same rolling origins, at the same horizons. That comparison is easy to arrange and frequently not arranged, which is why published deep forecasting results and applied experience disagree as often as they do.

Horizon matters#

Forecasting tomorrow is easier than forecasting next month, and errors grow with the horizon. Evaluate at every horizon you will actually use rather than reporting a single averaged number that hides how quickly accuracy decays.

There are two ways to get a multi-step forecast and the choice has consequences. The recursive way forecasts one step, appends the prediction, and forecasts again; it needs one model and it compounds its own errors, so by step thirty it may be forecasting from a history it invented. The direct way trains a separate model per horizon, each predicting h steps ahead from features available now; it costs h models and does not compound, and it is usually the better answer once the horizon is more than a few steps.

                 Recursive                   Direct
models            1                          one per horizon
error behaviour   compounds with each step   independent per step
features          can use lag_1              only lags at least h old
consistency       forecasts flow smoothly    steps can disagree with each other
use when          short horizons             the horizon is the thing you care about

Intervals, not points#

Coverage is a measurement, not an assumption. Take every backtest fold, count how often the actual value fell inside the stated interval, and compare that with the number on the label. An eighty percent interval that contains sixty percent of outcomes is not a cautious forecast, it is a wrong one, and the planner who sized their safety stock from it will find out in a way you do not want to be present for.

import numpy as np

def coverage(actual, low, high):
    return float(np.mean((actual >= low) & (actual <= high)))

rng = np.random.default_rng(1)
actual = rng.normal(100, 10, 500)              # the truth: standard deviation 10

honest = coverage(actual, 100 - 1.2816 * 10, 100 + 1.2816 * 10)
optimistic = coverage(actual, 100 - 1.2816 * 6, 100 + 1.2816 * 6)
print(round(honest, 3), round(optimistic, 3))

# 0.852 0.616
# The second interval is labelled 80% and covers 62%. Nothing in the point
# forecast changed; only the claimed uncertainty was wrong.

Intervals from a model's own error assumptions are usually too narrow, because they account for the noise and not for the model being wrong about the level. Two fixes work in practice. Fit quantile models directly, one per quantile, so the interval is learned rather than derived. Or take the residuals from the backtest, use their empirical quantiles as the interval width, and widen with the horizon — the crude version of conformal prediction, and considerably better than a formula that assumes the model is correct.

Forecasting is where a careless evaluation fools you most easily, because the ordinary way of splitting data lets the model see the future. This section sets up the discipline instead: the seasonal naive baseline you have to beat before anything clever is worth trying, backtesting that respects the order of time, a metric that can be compared across series of different sizes, lag features built without leakage, and quantile models for when a range matters more than a single number.

Beat this first#

def seasonal_naive(history, horizon, season=7):
    """Forecast each future day as the same weekday last week."""
    last_season = history[-season:]
    return [last_season[i % season] for i in range(horizon)]

For daily data with a weekly pattern, this is surprisingly hard to beat. Any model that does not beat it clearly is not worth its complexity.

Clearly is doing work in that sentence. A model that beats seasonal naive by three percent on one test window has not beaten it, because the fold-to-fold spread is usually larger than three percent. Run several origins, look at the spread, and treat a win smaller than the spread as no result. This is also the honest way to end a project early: if nothing beats the baseline after a fortnight, the baseline is the deliverable, and a seasonal naive forecast in production costs nothing to run and never has an outage.

Which baseline to use depends on the series, and picking the flattering one is a way of cheating that is easy to do by accident. Use seasonal naive on anything with a weekly or yearly shape. Use plain naive, tomorrow equals today, on anything that looks like a random walk, which includes most prices and most financial series and which seasonal naive loses to badly. If you are unsure, compute all of them; they are three lines each and the most expensive one to skip is the one you would have lost to.

Backtest like you will forecast#

train ------------------|test 14 days|
train ----------------------------|test 14 days|
train --------------------------------------|test 14 days|

Each fold trains only on the past, forecasts the next 14 days, then the origin moves forward.

One split tells you how the model did in one fortnight, which may have been unusual. Several rolling origins show how it does across conditions, and match how the model will be used: retrained and forecasting forward repeatedly.

import numpy as np

def rolling_origin(series, forecaster, score, horizon=14, folds=5,
                   step=14, min_train=120):
    """Each fold trains on everything up to its origin and scores the next
    horizon days. The origin moves forward by step each time."""
    out = []
    for f in range(folds):
        end = min_train + f * step
        if end + horizon > len(series):
            break
        history, actual = series[:end], series[end : end + horizon]
        out.append((end, score(actual, forecaster(history, horizon), history)))
    return out

A metric that means something#

import numpy as np

def mase(actual, forecast, history, season=7):
    """Error relative to the in-sample seasonal naive error. Below 1 beats the baseline."""
    naive_err = np.mean(np.abs(history[season:] - history[:-season]))
    return np.mean(np.abs(actual - forecast)) / naive_err

Percentage errors explode when actual values are near zero, and plain error cannot be compared across a high-volume and a low-volume product. Scaling by the naive error fixes both, and says directly whether you beat the baseline.

Put the two together and the value of several origins becomes visible immediately. The same forecaster, on the same series, scores between 0.74 and 1.39 depending only on which fortnight it was asked about.

import numpy as np

def seasonal_naive(history, horizon, season=7):
    last = history[-season:]
    return np.array([last[i % season] for i in range(horizon)])

def flat(history, horizon):
    return np.full(horizon, history[-1])

rng = np.random.default_rng(0)
t = np.arange(260)
series = 100 + 0.1 * t + 12 * np.sin(2 * np.pi * t / 7) + rng.normal(0, 4, 260)

for end, m in rolling_origin(series, seasonal_naive, mase):
    print("train to", end, "MASE", round(m, 3))

scores = [m for _, m in rolling_origin(series, seasonal_naive, mase)]
flat_scores = [m for _, m in rolling_origin(series, flat, mase)]
print("seasonal naive mean", round(float(np.mean(scores)), 3))
print("flat naive mean    ", round(float(np.mean(flat_scores)), 3))

# train to 120 MASE 0.828
# train to 134 MASE 1.197
# train to 148 MASE 0.74
# train to 162 MASE 1.388
# train to 176 MASE 0.798
# seasonal naive mean 0.99
# flat naive mean     1.854

Lags without leakage#

When forecasting 14 days ahead, a feature like 'yesterday's sales' is not available for day 14 at forecast time. Either build features only from data at least as old as the horizon, or forecast one step at a time and feed predictions back in, knowing errors compound.

import numpy as np

def features_for(series, i, horizon):
    """Features for predicting day i, using only data available h days earlier.
    Every lag is at least horizon old, so this works for the direct approach."""
    cut = i - horizon                       # the last day you will actually know
    if cut < 28:
        return None
    past = series[:cut + 1]
    return {
        "lag_h": past[-1],                  # h days before the target
        "lag_h_7": past[-8],                # and a week before that
        "roll_7": past[-7:].mean(),
        "roll_28": past[-28:].mean(),
        "vol_28": past[-28:].std(),
        "dow": i % 7,                       # known for any future date
    }

series = np.arange(100).astype(float)
print(features_for(series, i=60, horizon=14)["lag_h"])    # 46.0, not 59.0

Intervals from quantile models#

from sklearn.ensemble import HistGradientBoostingRegressor

models = {q: HistGradientBoostingRegressor(loss="quantile", quantile=q).fit(X_train, y_train)
          for q in (0.1, 0.5, 0.9)}
low, mid, high = (models[q].predict(X_test) for q in (0.1, 0.5, 0.9))
coverage = np.mean((y_test >= low) & (y_test <= high))   # should be near 0.8

An inventory planner needs to know how bad it could plausibly be, not just the middle guess. Check that an 80% interval really contains about 80% of actual values in the backtest.

Three quantile models are three separate fits, and nothing forces them to agree. On a hard row the 0.1 model can predict above the 0.9 model, which is called quantile crossing and produces an interval with a negative width. Nobody notices until a downstream system divides by it. Sort the three predictions row by row before using them; it is one line and it is always correct, because a set of quantiles that are out of order was wrong in a way sorting cannot make worse.

import numpy as np

low = np.array([12.0, 30.0, 41.0])
mid = np.array([18.0, 28.0, 45.0])        # the middle model disagrees on row 2
high = np.array([25.0, 27.0, 52.0])       # and the upper one is lowest there

stacked = np.sort(np.column_stack([low, mid, high]), axis=1)
low, mid, high = stacked[:, 0], stacked[:, 1], stacked[:, 2]
print(low, mid, high)
# [12. 27. 41.] [18. 28. 45.] [25. 30. 52.]

print(bool((high >= low).all()))      # True, now guaranteed

Then report coverage per horizon, not pooled. Intervals are usually about right one day ahead and far too narrow at thirty, because the model's uncertainty about its own level grows with the horizon and a fixed-width interval does not. A single pooled coverage of 79 percent can hide 95 percent at day one and 55 percent at day thirty, and day thirty is the one the planner is sizing stock against.

In practice you rarely forecast one thing. You forecast every product in every store, including the ones launched last week with almost no history. This section covers the choice that follows: a separate model per series, which is simple but learns nothing across them, or one global model that shares what it learns between them. It also covers hierarchies that have to add up, and demand that is zero on most days.

One model per series, or one for all#

With ten thousand products, fitting a separate classical model to each is simple but learns nothing across products. A single global model, such as gradient boosting with product features and lags, learns shared patterns and helps products with short histories. Global models have done very well in forecasting competitions for this reason.

                      Local: one model per series    Global: one model, all series
short-history series  fits badly or not at all       borrows shape from similar ones
per-series quirks     captured exactly               captured only if a feature says so
training time         10,000 small fits              one fit
deployment            10,000 artefacts to version    one artefact
a new product         no model until it has history  predicted from day one
diagnosis             read one model                 read feature importances, once
scales to             hundreds of series             millions

The global model needs the series identity supplied as features rather than as an identifier, because an id is meaningless to a tree. Give it what actually distinguishes the series: category, price band, store size, how long it has existed, its average volume, its typical weekday shape. Then a product launched last week is predicted from products that resemble it, which is the whole point and is exactly what a local model cannot do.

Hierarchies must add up#

If store forecasts do not sum to the region forecast, planners get two different answers. Reconciliation methods adjust forecasts at all levels so they add up, often improving accuracy at the same time. Libraries in the statsforecast family implement them.

The two simple approaches are worth knowing by hand, because they encode different beliefs about where the signal is. Bottom-up forecasts each leaf and adds, which is right when the leaves are stable and noisy in independent ways. Top-down forecasts the total, which is smoother and easier to predict, then splits it by historical shares. Neither is generally better, and the honest answer is to check which wins at every level of your own hierarchy.

import numpy as np

def bottom_up(leaf_forecasts):
    return float(np.sum(leaf_forecasts)), np.asarray(leaf_forecasts, dtype=float)

def top_down(total_forecast, historical_shares):
    shares = np.asarray(historical_shares, dtype=float)
    shares = shares / shares.sum()
    return float(total_forecast), total_forecast * shares

leaves = np.array([100.0, 60.0, 40.0])     # per-store forecasts
total, per_store = bottom_up(leaves)
print(total, per_store)                    # 200.0 [100.  60.  40.]

total, per_store = top_down(210.0, leaves) # the total model says 210
print(total, per_store.round(1))           # 210.0 [105.  63.  42.]

Mostly zeros#

Many products sell nothing on most days. Standard models forecast a small fraction every day, which is never what happens. Methods for intermittent demand, such as Croston's, forecast the size and timing of sales separately. Evaluate them on the decision they support, such as stock-outs, not only on error.

Croston's method is two exponential smoothers: one over the sizes of the non-zero sales, one over the gaps between them. The forecast rate is the smoothed size divided by the smoothed gap. What it buys you is not a better average — the average is easy — but the two numbers a stock decision actually needs: how much when it sells, and how long you wait.

import numpy as np

def croston(series, alpha=0.1):
    """Returns the demand rate per period, from smoothed size and smoothed gap."""
    size = interval = None
    gap = 0
    for x in series:
        gap += 1
        if x > 0:
            size = x if size is None else alpha * x + (1 - alpha) * size
            interval = gap if interval is None else alpha * gap + (1 - alpha) * interval
            gap = 0
    return size / interval, size, interval

d = np.array([0, 0, 0, 5, 0, 0, 0, 0, 7, 0, 0, 6, 0, 0, 0, 0, 0, 4.0])
rate, size, interval = croston(d)
print(round(float(d.mean()), 3))     # 1.222  a flat model forecasts this every day
print(round(rate, 3))                # 1.229  the same rate, arrived at differently
print(round(size, 2), round(interval, 2))   # 5.15 4.19

Which is also why MASE and RMSE are poor judges here. A model that forecasts 1.22 every day has a low absolute error and orders the wrong amount at the wrong time. Score these on the decision: how often the item was out of stock, and how much capital sat in a warehouse. If those two numbers are not available, forecasting intermittent demand accurately is not something anyone can confirm you have done.