Forecasting Hourly Store Traffic to Staff the Floor: Time Series from Math to Production 🛒

Statistics
Time Series
Forecasting
ARIMAX
Prophet
Gradient Boosting
MLOps
An end-to-end time series case study: forecast hourly customer arrivals at a grocery store, turn the forecast into a staffing plan with queueing theory, compare ARIMAX, Prophet and gradient boosting with rolling-origin validation, then save, ship, load and monitor the model.
Author

MLStation

Published

October 3, 2026

Every evening at 18:00, the manager of a mid-sized grocery store publishes tomorrow’s shift plan. Put too few people on the tills and queues grow, baskets are abandoned and loyal customers drift away. Put too many on and the store pays for idle hours. The whole decision rests on one number per hour: how many customers will walk in?

That makes it a time series forecasting problem with money attached. In this article we solve it end to end: from the business question, through the math of three forecasting models, to a saved model that a scheduling app loads and runs. 🚀

Tip🎯 What you will learn
  • How to frame a forecast around a business decision, not a metric
  • Why ordinary regression with a random train/test split gives wrong answers on time series
  • The math behind ARIMAX, Prophet and gradient boosting with lag features, kept short
  • How to validate forecasts honestly with rolling-origin evaluation, MASE and pinball loss
  • How to turn a forecast into staff numbers with the Erlang C queueing model and a cost-based quantile
  • How to save, ship, load and monitor the model, with an animated production console

All data in this article are simulated: two years of hourly arrivals for a fictional store, FreshCart, open 07:00–23:00. The simulation includes daily and weekly rhythm, seasons, weather, promotions, holidays, local events, paydays and a competitor opening. Because we generated the data, we also know the true arrival rate, which lets us measure how close each model gets to the best possible forecast.

🛒 1. The problem: a Friday with fixed staffing

FreshCart staffs the tills the same way every day. Watch what happens on a busy Friday with promotions running. Press Play, then try changing the number of cashiers.

Show the code: simulate two years of hourly arrivals 
import numpy as np, pandas as pd, json, math, time, hashlib, tempfile, pathlib
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from ms_charts import COLORS, SEQ, PLOTLY_CONFIG

rng = np.random.default_rng(2026)
OPEN, CLOSE = 7, 23                     # open 07:00-23:00 -> 16 hourly slots
H = CLOSE - OPEN
days = pd.date_range("2024-01-01", "2025-12-31", freq="D"); nd = len(days)
doy = days.dayofyear.values

# --- daily drivers, all known the evening before ---
temp = 19 + 11*np.sin(2*np.pi*(doy-110)/365.25) + rng.normal(0, 3.2, nd)   # forecast high, °C
rain = (rng.random(nd) < 0.22 + 0.10*np.cos(2*np.pi*(doy-120)/365.25)).astype(int)
promo = np.array([int(d.dayofweek in (3, 4, 5) and d.isocalendar().week % 5 in (0, 2)) for d in days])
hol = pd.to_datetime(["2024-01-01", "2024-05-27", "2024-07-04", "2024-09-02", "2024-11-28", "2024-12-25",
                      "2025-01-01", "2025-05-26", "2025-07-04", "2025-09-01", "2025-11-27", "2025-12-25"])
holiday = days.isin(hol).astype(int); pre_hol = days.isin(hol - pd.Timedelta(days=1)).astype(int)
event = np.zeros(nd, int); event[rng.choice(nd, 40, replace=False)] = 1    # game nights on the town calendar
payday = ((days.day == 1) | (days.day == 15)).astype(int)
competitor = (days >= pd.Timestamp("2025-03-15")).astype(int)            # a discounter opens nearby

hours = np.arange(OPEN, CLOSE)
wk_prof = .55*np.exp(-.5*((hours-12.5)/1.6)**2) + 1.0*np.exp(-.5*((hours-17.8)/1.4)**2) + .25  # weekday: lunch + after work
we_prof = 1.0*np.exp(-.5*((hours-11.5)/2.4)**2) + .55*np.exp(-.5*((hours-16)/2.0)**2) + .25   # weekend: late morning
dow_lvl = np.array([0.0, -.05, -.02, .03, .12, .22, .10])

rows, eta = [], 0.0
for d in range(nd):
    dow = days[d].dayofweek; prof = we_prof if dow >= 5 else wk_prof; tz = (temp[d]-20)/10
    for j, h in enumerate(hours):
        eta = 0.9*eta + rng.normal(0, 0.06)                               # unobserved surprises that persist for hours (AR(1), log scale)
        log_lam = (3.65 + np.log(prof[j]) + dow_lvl[dow] + 0.10*d/365.25 + 0.05*np.sin(2*np.pi*(doy[d]-330)/365.25)
                   - 0.06*tz**2 + 0.03*tz - 0.16*rain[d]*(1 + 0.6*(h >= 16)) + 0.20*promo[d]
                   + 0.45*pre_hol[d]*(h >= 14) - 0.9*holiday[d]
                   + 0.12*promo[d]*(dow >= 5) + 0.10*promo[d]*(h >= 17)
                   - 0.10*max(tz, 0)*(13 <= h <= 16) + 0.08*max(tz, 0)*(h >= 19)
                   + 0.55*event[d]*(16 <= h <= 19) - 0.25*event[d]*(h >= 20)
                   + 0.10*payday[d] - 0.08*competitor[d]*(dow < 5) + eta)
        lam = np.exp(log_lam)
        y = rng.poisson(rng.gamma(50.0, lam/50.0))                        # negative binomial counts
        rows.append((days[d] + pd.Timedelta(hours=int(h)), d, dow, h, temp[d], rain[d], promo[d], holiday[d],
                     pre_hol[d], event[d], payday[d], competitor[d], y, lam))
df = pd.DataFrame(rows, columns=["ts", "day", "dow", "hour", "temp", "rain", "promo", "holiday", "pre_holiday",
                                 "event", "payday", "competitor", "y", "lam"])
df["t"] = np.arange(len(df)); df["doy"] = df.ts.dt.dayofyear; df["weekend"] = (df.dow >= 5).astype(int)
print(f"{len(df):,} hourly rows · {df.ts.min():%Y-%m-%d} to {df.ts.max():%Y-%m-%d} · mean {df.y.mean():.1f} customers/hour · max {df.y.max()}")
11,696 hourly rows · 2024-01-01 to 2025-12-31 · mean 29.4 customers/hour · max 163
Show the code: pick a busy promo Friday 
fri = df[(df.dow == 4) & (df.promo == 1) & (df.ts.dt.year == 2025)].groupby("day").y.sum().idxmax()
friday = df[df.day == fri][["hour", "y"]].rename(columns={"y": "arrivals"})
ojs_define(friday=friday.to_dict(orient="records"),
           friday_label=f"{df[df.day == fri].ts.iloc[0]:%A %d %B %Y}")

With 3 cashiers the evening rush overwhelms the tills and dozens of customers walk out. With 7 the store survives the rush but pays for many idle hours in the morning. No fixed number works. The right answer changes every hour, so we must forecast every hour and staff to the forecast.

🧭 The decision we are supporting

Question Answer for FreshCart
⏰ When is the decision made? 18:00 the day before, when tomorrow’s shifts are published
🔭 Forecast horizon The 16 opening hours of the next day (1 to 29 hours ahead)
🧮 What do we forecast? Customer arrivals per hour, with a probability range, not just one number
👥 What do we decide? Cashiers and floor assistants per hour
💰 What does a wrong answer cost? Understaffing: lost sales and goodwill. Overstaffing: wages for idle time
📥 What is known at 18:00? History up to now, tomorrow’s calendar (holiday, promo, events, payday) and the weather forecast

The last row matters most. A feature is only allowed if it is known when the forecast is made. Tomorrow’s actual sales are not; tomorrow’s promotion calendar is.

flowchart TD
  A["🛒 Business problem<br/>staff each hour"] --> B["📊 Data<br/>hourly arrivals + drivers"]
  B --> C["🔍 Explore<br/>seasonality, effects, autocorrelation"]
  C --> D["🧮 Model<br/>ARIMAX · Prophet · Gradient boosting"]
  D --> E["✅ Validate<br/>rolling origin, MASE, pinball, coverage"]
  E --> F["🏆 Choose<br/>accuracy + calibration + operations"]
  F --> G["👥 Decide<br/>Erlang C + cost-based quantile"]
  G --> H["💾 Save<br/>versioned artifact + model card"]
  H --> I["🚚 Serve<br/>scheduler loads model, predicts"]
  I --> J["📡 Monitor<br/>errors, drift, retrain"]
  J -. "new data" .-> B

📊 2. Look at the data first

A forecaster’s first job is to understand the structure of the series. Use the range slider to zoom from two years down to a single week.

Show the code 
daily = df.groupby(df.ts.dt.date).agg(y=("y", "sum"), holiday=("holiday", "max")).reset_index()
daily["ts"] = pd.to_datetime(daily.ts); daily["ma7"] = daily.y.rolling(7, center=True).mean()
fig = go.Figure()
fig.add_trace(go.Scatter(x=daily.ts, y=daily.y, name="Daily customers", line=dict(color=COLORS["indigo"], width=1), opacity=.55, hovertemplate="%{x|%a %d %b %Y}<br>Daily customers: %{y:.0f}<extra></extra>"))
fig.add_trace(go.Scatter(x=daily.ts, y=daily.ma7, name="7-day average", line=dict(color=COLORS["rose"], width=2.5), hovertemplate="%{x|%a %d %b %Y}<br>7-day average: %{y:.1f}<extra></extra>"))
hd = daily[daily.holiday == 1]
fig.add_trace(go.Scatter(x=hd.ts, y=hd.y, mode="markers", name="Public holiday", marker=dict(color=COLORS["amber"], size=9, symbol="diamond"), hovertemplate="%{x|%a %d %b %Y}<br>Holiday traffic: %{y:.0f}<extra></extra>"))
fig.add_vline(x=pd.Timestamp("2025-03-15").timestamp()*1000, line=dict(color=COLORS["slate"], dash="dot"))
fig.add_annotation(x=pd.Timestamp("2025-03-15"), y=daily.y.max(), text="Competitor opens", showarrow=False, yshift=10, font=dict(color=COLORS["slate"]))
fig.update_layout(height=400, margin=dict(l=60, r=20, t=30, b=40), yaxis_title="Customers per day",
                  xaxis=dict(rangeslider=dict(visible=True), rangeselector=dict(buttons=[
                      dict(count=7, label="1w", step="day", stepmode="backward"),
                      dict(count=1, label="1m", step="month", stepmode="backward"),
                      dict(count=6, label="6m", step="month", stepmode="backward"), dict(step="all", label="All")])),
                  legend=dict(orientation="h", y=1.12))
fig.show(config=PLOTLY_CONFIG)
Figure 1: Daily customer totals (simulated). Notice the weekly sawtooth, holiday dips and the gentle yearly wave. Zoom in to see the weekday dip after a competitor opened in March 2025.

The hourly pattern is different on weekdays and weekends: a lunch bump and a big after-work peak on weekdays, a late-morning peak on Saturdays and Sundays.

Show the code 
hm = df[df.holiday == 0].pivot_table(index="dow", columns="hour", values="y", aggfunc="mean")
fig = go.Figure(go.Heatmap(z=hm.values, x=[f"{h}:00" for h in hm.columns], y=["Mon", "Tue", "Wed", "Thu", "Fri", "Sat", "Sun"],
                           colorscale=[[0, "#EEF2FF"], [.5, COLORS["indigo"]], [1, "#3730A3"]], colorbar=dict(title="Avg"),
                           hovertemplate="%{y} %{x}<br>%{z:.1f} customers<extra></extra>"))
fig.update_layout(height=330, margin=dict(l=60, r=20, t=20, b=40), yaxis=dict(autorange="reversed"))
fig.show(config=PLOTLY_CONFIG)
Figure 2: Average customers by hour and day of week. Weekdays peak around 18:00; weekends peak before noon.

How strongly do the drivers move traffic? Each bar compares average customers per hour with and without the driver.

Show the code 
base = df[(df.holiday == 0)]
def lift(col): g = base.groupby(col).y.mean(); return 100*(g.loc[1]/g.loc[0]-1)
eff = pd.DataFrame({"driver": ["🏷️ Promotion", "🌧️ Rain", "🏟️ Event day", "💵 Payday", "🏪 Competitor open"],
                    "lift": [lift("promo"), lift("rain"), lift("event"), lift("payday"), lift("competitor")]})
pre = 100*(df[df.pre_holiday == 1].y.mean()/base[base.pre_holiday == 0].y.mean()-1)
eff.loc[len(eff)] = ["🎄 Day before holiday", pre]
eff = eff.sort_values("lift")
fig = go.Figure(go.Bar(x=eff.lift, y=eff.driver, orientation="h", marker_color=[COLORS["rose"] if v < 0 else COLORS["teal"] for v in eff.lift],
                       text=[f"{v:+.0f}%" for v in eff.lift], textposition="outside", cliponaxis=False, hovertemplate="%{y}<br>Raw average lift: %{x:+.1f}%<extra></extra>"))
fig.update_layout(height=320, margin=dict(l=170, r=50, t=20, b=40), xaxis_title="Change in average customers per hour (%)",
                  xaxis=dict(range=[min(eff.lift.min(), 0) - 12, eff.lift.max() + 12]))
fig.show(config=PLOTLY_CONFIG)
Figure 3: Average effect of each driver on hourly traffic, as a percentage change (simulated data, raw averages before adjusting for other drivers).
Note📝 Raw averages can mislead

The competitor bar mixes two things: the competitor effect and the store’s own growth over time. Separating overlapping effects is exactly what a model does. Raw comparisons are for spotting candidates, not for measuring effects.

🔁 Autocorrelation: the fingerprint of a time series

A time series has memory when knowing the past helps us anticipate the present. Autocorrelation measures that memory.

Start with one concrete question:

If the store was unusually busy at 18:00 yesterday, is 18:00 today also likely to be unusually busy?

Let y_t mean the number of customers in the current opening hour t. Then y_{t-k} is the observation k opening-hours earlier. The autocorrelation at lag k is

\rho_k = \operatorname{Corr}(y_t,y_{t-k}) = \frac{\operatorname{Cov}(y_t,y_{t-k})} {\sqrt{\operatorname{Var}(y_t)\operatorname{Var}(y_{t-k})}}.

For a stable series, the two variances are approximately the same, giving the familiar simplified form

\rho_k \approx \frac{\operatorname{Cov}(y_t,y_{t-k})}{\operatorname{Var}(y_t)}.

The sample ACF replaces the unknown population quantities with what we can calculate from the observed data:

\hat\rho_k = \frac{\sum_{t=k+1}^{n}(y_t-\bar y)(y_{t-k}-\bar y)} {\sum_{t=1}^{n}(y_t-\bar y)^2}.

Read the symbols before the chart

  • t: the current time index.
  • k: the lag — how far backward we look.
  • y_t: customers observed at time t.
  • \bar y: average hourly customer count.
  • n: number of observations.
  • \rho_k: true autocorrelation at lag k; \hat\rho_k is its estimate from our data.
  • \rho_k\approx 1: strong positive memory; \rho_k\approx0: little linear memory; \rho_k<0: high values tend to be followed by low values at that lag.

FreshCart is open for 16 hourly slots per day. Therefore lag 16 compares an hour with the same store hour yesterday, while lag 112 = 16\times7 compares it with the same hour last week. Peaks there are the statistical fingerprint of daily and weekly repetition.

Show the code 
def acf(x, nlags):
    x = np.asarray(x, float) - np.mean(x); d = np.sum(x**2)
    return np.array([np.sum(x[k:]*x[:len(x)-k])/d for k in range(nlags+1)])
r = acf(df.y, 130)
cols = [COLORS["amber"] if k % 112 == 0 and k else (COLORS["rose"] if k % 16 == 0 and k else COLORS["indigo"]) for k in range(131)]
fig = go.Figure(go.Bar(x=np.arange(131), y=r, marker_color=cols, hovertemplate="lag %{x}<br>ρ = %{y:.2f}<extra></extra>"))
band = 1.96/np.sqrt(len(df))
for s in (band, -band): fig.add_hline(y=s, line=dict(dash="dash", color=COLORS["slate"]))
fig.add_annotation(x=16, y=r[16], text="1 day", yshift=14, showarrow=False); fig.add_annotation(x=112, y=r[112], text="1 week", yshift=14, showarrow=False)
fig.update_layout(height=320, margin=dict(l=60, r=20, t=20, b=40), xaxis_title="Lag (opening hours)", yaxis_title="Autocorrelation")
fig.show(config=PLOTLY_CONFIG)
Figure 4: Autocorrelation of hourly arrivals. Peaks at multiples of 16 (one day) and the tallest at 112 (one week). Dashed lines: approximate 95% band for white noise, ±1.96/√n.

⚠️ 3. Why ordinary regression is the wrong tool

The title needs one important qualification:

Linear regression itself is not forbidden for time series.
The problem is using an ordinary regression as if rows were independent, exchangeable observations and then validating it with a random split.

Suppose we begin with a perfectly reasonable model:

y_t = \beta_0+\beta_1\,\text{hour}_t+\beta_2\,\text{promo}_t+\cdots+\varepsilon_t.

Here:

  • y_t is customer traffic in hour t;
  • \beta_0 is the baseline level;
  • each \beta_j describes how a feature changes the expected traffic;
  • \varepsilon_t is everything the features did not explain — the residual or surprise.

If those surprises were independent and had roughly constant variance, ordinary regression could work well. FreshCart violates those conditions.

Ordinary tabular regression thinks

Row 101 and row 102 are simply two observations.

Once their feature columns are known, the leftover errors are unrelated.

A time series says

Row 102 happened after row 101.

A surprise at 17:00 can persist into 18:00, and the same hour yesterday or last week can still carry information.

Why this matters — four concrete failure modes

1. The residuals remember the past.
Ordinary least squares commonly uses the working assumption

\operatorname{Cov}(\varepsilon_t,\varepsilon_s)=0 \qquad(t\neq s),

meaning the unexplained part of one hour tells us nothing about another. But if an unrecorded local event makes 17:00 unusually busy, 18:00 may also be busy. Then

\operatorname{Cov}(\varepsilon_t,\varepsilon_{t-1})>0.

The coefficient estimate can still be useful, but the familiar OLS uncertainty formula

\operatorname{Var}(\hat{\boldsymbol\beta}) = \sigma^2(X^\top X)^{-1}

is no longer the correct variance formula. With correlated errors it becomes

\operatorname{Var}(\hat{\boldsymbol\beta}) = (X^\top X)^{-1}X^\top\Omega X(X^\top X)^{-1},

where \Omega records the covariance among residuals. The practical consequence is more important than memorizing the matrix: ordinary standard errors and prediction intervals can look more confident than they should, and we are throwing away useful temporal memory.

2. Random train/test splitting lets the future surround the test row.
In ordinary machine learning, random splitting is often excellent. In forecasting it changes the task. The model might train on Tuesday and Thursday while Wednesday is in the test set. That asks, “Can I reconstruct a missing point between known points?” Production asks, “Can I predict a point beyond everything I have seen?”

Production rule: if the forecast is created Tuesday at 18:00, every feature used for Wednesday must already be available Tuesday at 18:00.

This is why a centered moving average is dangerous: it can quietly include Wednesday information while pretending to be a Tuesday feature.

3. Customer arrivals are counts.
Hourly arrivals are non-negative integers. Busy hours also tend to have more variability than quiet hours. A Gaussian linear model on raw counts does not naturally enforce either fact. It can produce negative predictions and constant-width uncertainty bands. We can repair this with transformations, Poisson/count models, or other suitable likelihoods.

4. Time can change the data-generating process.
Traffic grows, seasons repeat, and structural breaks happen. When a competitor opens, the relationship before March may not be identical to the relationship afterward. A model must represent those changes or be retrained/adapted.

So the lesson is not “do not use regression.” It is:

\boxed{ \text{Respect time dependence} + \text{respect forecast availability} + \text{use an appropriate response model} + \text{validate forward in time}. }

The next diagnostic makes the first two problems visible.

Show the code: naive OLS, residual autocorrelation and the leakage test 
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.model_selection import KFold

Xn = pd.get_dummies(df[["hour", "dow"]].astype(str)).astype(float)
Xn[["promo", "rain", "holiday", "temp"]] = df[["promo", "rain", "holiday", "temp"]]
Xn["t"] = df.t/len(df); Xn["const"] = 1.0
bn = np.linalg.lstsq(Xn.values, df.y.values, rcond=None)[0]
res_naive = df.y.values - Xn.values @ bn
dw = np.sum(np.diff(res_naive)**2)/np.sum(res_naive**2)        # Durbin-Watson: 2 = no autocorrelation
r_res = acf(res_naive, 40)
neg_pred = int(((Xn.values @ bn) < 0).sum())

# leakage: same model, random 5-fold CV vs. honest "train on the past, test on the next 3 months"
feat_leak = ["hour", "dow", "doy", "temp", "rain", "promo", "holiday", "pre_holiday", "event", "payday", "t"]
gb = lambda: HistGradientBoostingRegressor(loss="poisson", max_iter=300, learning_rate=.05, random_state=0)
df["y_smooth3"] = df.y.rolling(3, center=True, min_periods=1).mean()   # a "smoothed traffic" column: quietly uses the next hour
def kfold_mae(cols):
    out = []
    for tr, te in KFold(5, shuffle=True, random_state=0).split(df):
        m = gb().fit(df.iloc[tr][cols], df.y.iloc[tr]); out.append(np.mean(np.abs(m.predict(df.iloc[te][cols]) - df.y.iloc[te])))
    return np.mean(out)
rand_mae = kfold_mae(feat_leak); leak_mae = kfold_mae(feat_leak + ["y_smooth3"])
cut = pd.Timestamp("2025-10-01")
trn, tst = df[df.ts < cut], df[df.ts >= cut]
m = gb().fit(trn[feat_leak], trn.y); time_mae = np.mean(np.abs(m.predict(tst[feat_leak]) - tst.y))
print(f"Durbin-Watson = {dw:.2f}  (2.0 = no autocorrelation)")
print(f"Lag-1 residual autocorrelation = {r_res[1]:.2f}")
print(f"Negative predictions from linear model on counts: {neg_pred}")
print(f"Random 5-fold MAE           = {rand_mae:.2f}")
print(f"  + smoothed-traffic column  = {leak_mae:.2f}")
print(f"Honest past -> future test  = {time_mae:.2f}")
Durbin-Watson = 0.81  (2.0 = no autocorrelation)
Lag-1 residual autocorrelation = 0.60
Negative predictions from linear model on counts: 54
Random 5-fold MAE           = 5.95
  + smoothed-traffic column  = 4.74
Honest past -> future test  = 6.59
Show the code 
fig = make_subplots(rows=1, cols=2, column_widths=[.6, .4], subplot_titles=("Residual ACF of naive regression", "MAE: random split vs. future test"))
fig.add_trace(go.Bar(x=np.arange(1, 41), y=r_res[1:], marker_color=COLORS["rose"], showlegend=False), 1, 1)
for s in (band, -band): fig.add_hline(y=s, line=dict(dash="dash", color=COLORS["slate"]), row=1, col=1)
fig.add_trace(go.Bar(x=["Random<br>+ leak", "Random", "Past →<br>future"], y=[leak_mae, rand_mae, time_mae],
                     marker_color=[COLORS["rose"], COLORS["amber"], COLORS["indigo"]], text=[f"{v:.2f}" for v in (leak_mae, rand_mae, time_mae)], textposition="outside", showlegend=False), 1, 2)
fig.update_xaxes(title_text="Lag (hours)", row=1, col=1); fig.update_xaxes(tickangle=0, row=1, col=2); fig.update_yaxes(title_text="Customers (MAE)", row=1, col=2)
fig.update_layout(height=360, margin=dict(l=60, r=20, t=50, b=50))
fig.show(config=PLOTLY_CONFIG)
Figure 5: Left: residuals of the naive regression are strongly autocorrelated (independent errors would stay inside the dashed band). Right: random splits, especially with a column that peeks at neighbouring hours, report a better error than the model achieves on genuinely future data.

How to read the diagnostic

The Durbin–Watson statistic

d= \frac{\sum_t(e_t-e_{t-1})^2}{\sum_t e_t^2} \approx2(1-\hat\rho_1)

compresses lag-1 residual dependence into one number. Here e_t is the residual at hour t, and \hat\rho_1 is residual autocorrelation one hour apart.

  • d\approx2: little lag-1 autocorrelation.
  • d<2: positive autocorrelation — positive errors tend to follow positive errors.
  • d>2: negative autocorrelation.

For FreshCart, d= 0.81, so the naive regression leaves strong temporal structure unexplained.

The right panel answers a different question: did our validation imitate production? Random splitting reports MAE 5.95, while the honest past→future test reports 6.59. The centered “smoothed traffic” feature looks even better under random CV (4.74) because it secretly contains future information.

That is data leakage: information unavailable at forecast time enters training or evaluation and makes the model look better than it can be after deployment.

Important✅ The fix is not “never use regression”

All three models below are regressions in a broad sense. What changes is that they (a) model the dependence over time, (b) respect the count nature of the data, and (c) are validated by training on the past and testing on the future.

🧮 4. Three ways to forecast

We give all three models the same information: the calendar, tomorrow’s weather forecast, promotions, events, paydays and the competitor flag. They differ in how they turn it into a forecast. Each model works on z_t = \log(1+y_t) or uses a count loss, which handles the “variance grows with the mean” problem.

Tip🧭 One notation rule for this section

Whenever you see a subscript t, read it as “at time t.” A hat, as in \hat y_t, means a model estimate or forecast. Bold symbols such as \mathbf{x}_t are vectors (several feature values together), while ordinary symbols such as y_t are single numbers.

Think of ARIMAX as regression + memory.

First explain the predictable part using tomorrow’s known features:

z_t=\mathbf{x}_t^\top\boldsymbol\beta+\eta_t, \qquad z_t=\log(1+y_t).

  • y_t: customers in hour t.
  • z_t=\log(1+y_t): a transformed count; the log compresses large peaks and guarantees positive forecasts after transforming back.
  • \mathbf{x}_t: the feature vector known for hour t — hour, weekday, weather forecast, promotion, holiday, event, etc.
  • \boldsymbol\beta: learned feature effects.
  • \eta_t: the remaining time-dependent surprise after those features are explained.

Now give the surprise memory. For this store, the useful version is deliberately simple:

\boxed{ \eta_t = \phi_1\eta_{t-1} + \phi_{16}\eta_{t-16} + \varepsilon_t }

Read it in English:

today’s unexplained traffic = some of the previous-hour surprise + some of the same-hour-yesterday surprise + a new random shock.

\phi_1 controls one-hour memory, \phi_{16} controls same-hour-yesterday memory, and \varepsilon_t is a new unpredictable innovation with standard deviation \sigma.

That is the main ARIMAX idea used in this article.

Textbooks often compress the same idea into

\phi(B)\Phi(B^s)(1-B)^d\eta_t = \theta(B)\Theta(B^s)\varepsilon_t.

Here B is the backshift operator: B\eta_t=\eta_{t-1}. s is a seasonal period, d is the number of differences, \phi and \Phi describe non-seasonal and seasonal autoregressive memory, while \theta and \Theta describe moving-average error terms.

You do not need this compact polynomial notation to follow our FreshCart model; it is included so the implementation connects to standard ARIMA literature.

Features \mathbf{x}_t. We include hour × weekend indicators, day of week, smooth yearly seasonality, trend, temperature, rain, promotion, event, payday, holiday and competitor status, plus a few business-motivated interactions.

How it is fitted. Conceptually we alternate between two jobs: estimate the feature effects, then estimate how the remaining residuals depend on recent residuals. Production libraries such as statsmodels can estimate the full model by maximum likelihood.

Forecast uncertainty. For future hour t+h, let \hat\eta_{t+h} be the predicted residual and let v_h be its variance. As the horizon h grows, more unknown future shocks can accumulate, so uncertainty usually widens. On the log scale,

\hat z_{t+h} = \mathbf{x}_{t+h}^{\top}\hat{\boldsymbol\beta} + \hat\eta_{t+h}.

Transforming back gives an approximate mean

E[1+y_{t+h}] = \exp\!\left(\hat z_{t+h}+\frac{v_h}{2}\right),

and a forecast quantile

Q_q(y_{t+h}) \approx \exp\!\left(\hat z_{t+h}+z_q\sqrt{v_h}\right)-1,

where q is the desired probability level and z_q is the corresponding standard-normal quantile.

Why use it here? ✅ Feature effects are interpretable; temporal memory is explicit; uncertainty is inexpensive to compute.
Limitation: ⚠️ complex interactions usually need to be designed by hand.

Prophet (Taylor and Letham, 2018) is a decomposable curve-fitting model:

z_t = g(t) + s(t) + h(t) + \mathbf{x}_t^\top\boldsymbol\beta + \varepsilon_t .

Read this as a sum of components:

  • g(t): long-term trend — the slowly changing baseline.
  • s(t): repeating seasonality such as weekly or yearly cycles.
  • h(t): holiday/event effects.
  • \mathbf{x}_t^\top\boldsymbol\beta: effects of extra regressors such as rain or promotions.
  • \varepsilon_t: unexplained noise.

The symbol t is time; it is not a fitted coefficient.

Trend g(t) is piecewise linear. A changepoint c_j is simply a candidate time where the slope is allowed to change. Prophet does not force a bend at every candidate; regularization keeps most changes small:

g(t) = \big(k + \mathbf{a}(t)^\top\boldsymbol\delta\big)\,t + \big(m + \mathbf{a}(t)^\top\boldsymbol\gamma\big), \qquad a_j(t) = \mathbb{1}[t \ge c_j],

with a sparse prior \delta_j \sim \operatorname{Laplace}(0,\tau) so the trend only bends when the data insist, as after the competitor opened.

Seasonality s(t) is a Fourier series for each period P (daily, weekly, yearly):

s_P(t) = \sum_{n=1}^{N}\Big[a_n \cos\!\Big(\tfrac{2\pi n t}{P}\Big) + b_n \sin\!\Big(\tfrac{2\pi n t}{P}\Big)\Big].

Holidays h(t) get their own effects. Prophet fits the maximum a posteriori estimate with L-BFGS. Our NumPy version uses a Gaussian (ridge) prior on the changepoints instead of a Laplace prior, which keeps it to one linear solve.

✅ Very easy for analysts, robust to gaps and outliers. ⚠️ Errors are treated as independent, so it cannot use “the last hour was unusually busy.”

Gradient-boosted trees (LightGBM, XGBoost, scikit-learn’s HistGradientBoostingRegressor) turn forecasting into supervised learning by adding time-safe lag features.

The beginner’s picture is:

small tree 1 fixes the biggest mistakes
       ↓
small tree 2 focuses on what remains
       ↓
small tree 3 corrects the remaining pattern
       ↓
many small corrections → final forecast

Mathematically the final score is an additive model:

\boxed{ F_M(\mathbf{x}) = F_0+\nu\sum_{m=1}^{M}f_m(\mathbf{x}) }

where F_0 is the starting prediction, f_m is tree m, M is the number of trees and \nu is the learning rate controlling how much each new tree is allowed to change the model.

For readers who want the optimization connection, boosting chooses each new tree to approximate the direction that most reduces the current loss:

g_t^{(m)} = -\frac{\partial \ell(y_t,F)}{\partial F}\Big|_{F_{m-1}}.

So g_t^{(m)} is not another observed feature; it is the negative loss gradient for training example t at boosting step m.

For counts we use the Poisson deviance, which models \log\lambda_t = F(\mathbf{x}_t):

\ell(y,F) = e^{F} - yF .

For staffing ranges we train extra models with the pinball (quantile) loss:

\ell_q(y,\hat y) = \max\big(q\,(y-\hat y),\,(q-1)(y-\hat y)\big).

Lag features must exist at 18:00 the day before: same hour yesterday (y_{t-16}), same hour last week (y_{t-112}) and yesterday’s mean. Using y_{t-1} would leak, because at 18:00 we do not yet know tomorrow’s 9:00 count.

✅ Learns interactions and non-linear effects automatically, and scales to thousands of stores in one model. ⚠️ Less interpretable, and its quantile intervals need calibration checks.

Family Examples When to consider it
🧠 Deep sequence models Temporal Fusion Transformer, N-HiTS Many related series, rich covariates, a team that can maintain deep models
🌐 Foundation models Chronos, TimesFM Zero-shot forecasts for new stores with little history
📦 Hierarchical reconciliation Store → region → chain, MinT Plans that must add up across levels

For one store with two years of data, the three models above are usually as accurate and far simpler to run. The textbook Forecasting: Principles and Practice covers the classical side in depth.

Show the code: the three models in NumPy and scikit-learn 
def fourier(x, period, K):
    return {f"{f}{k}": (np.sin if f == "s" else np.cos)(2*np.pi*k*np.asarray(x, float)/period) for k in range(1, K+1) for f in "sc"}

HW_COLS = [f"hw_{h}_{w}" for h in range(OPEN, CLOSE) for w in (0, 1)]
ARIMAX_FEATURES = HW_COLS + ["yr_s1", "yr_c1", "yr_s2", "yr_c2", "dow_fri", "dow_sat", "dow_sun", "trend", "temp", "temp2",
                             "rain", "rain_eve", "promo", "promo_we", "promo_eve", "heat_aft", "heat_eve", "event_pm",
                             "event_late", "payday", "holiday", "prehol_pm", "comp_wd"]

def arimax_features(d):
    """Feature matrix for the dynamic regression. Pure arithmetic, so it is easy to re-implement anywhere."""
    h, w, dow = d.hour.values, d.weekend.values, d.dow.values
    tz = (d.temp.values - 20)/10; hot = np.maximum(tz, 0)
    F = {c: ((h == int(c.split("_")[1])) & (w == int(c.split("_")[2]))).astype(float) for c in HW_COLS}
    for k, v in fourier(d.doy.values, 365.25, 2).items(): F["yr_" + k] = v
    F.update(dow_fri=(dow == 4)*1., dow_sat=(dow == 5)*1., dow_sun=(dow == 6)*1., trend=d.t.values/(365.25*H),
             temp=tz, temp2=tz**2, rain=d.rain.values*1., rain_eve=d.rain.values*(h >= 16), promo=d.promo.values*1.,
             promo_we=d.promo.values*w, promo_eve=d.promo.values*(h >= 17), heat_aft=hot*((h >= 13) & (h <= 16)),
             heat_eve=hot*(h >= 19), event_pm=d.event.values*((h >= 16) & (h <= 19)), event_late=d.event.values*(h >= 20),
             payday=d.payday.values*1., holiday=d.holiday.values*1., prehol_pm=d.pre_holiday.values*(h >= 14),
             comp_wd=d.competitor.values*(dow < 5))
    return np.column_stack([F[c] for c in ARIMAX_FEATURES]).astype(float)

class ARIMAX:
    """z_t = log(1+y_t) = x_t'beta + eta_t,  eta_t = phi1*eta_{t-1} + phi16*eta_{t-16} + eps_t  (two-step FGLS)."""
    def fit(self, d):
        X, z = arimax_features(d), np.log1p(d.y.values)
        beta = np.linalg.lstsq(X, z, rcond=None)[0]
        for _ in range(2):                                      # alternate: AR on residuals, then GLS refit of beta
            e = z - X @ beta
            A = np.column_stack([e[H-1:-1], e[:-H]]); phi = np.linalg.lstsq(A, e[H:], rcond=None)[0]
            Xf = X[H:] - phi[0]*X[H-1:-1] - phi[1]*X[:-H]; zf = z[H:] - phi[0]*z[H-1:-1] - phi[1]*z[:-H]
            beta = np.linalg.lstsq(Xf, zf, rcond=None)[0]
        e = z - X @ beta; A = np.column_stack([e[H-1:-1], e[:-H]]); phi = np.linalg.lstsq(A, e[H:], rcond=None)[0]
        self.beta, self.phi, self.sigma = beta, phi, float(np.std(e[H:] - A @ phi)); return self
    def forecast_day(self, prev_day, day, qs=(0.1, 0.8, 0.9)):
        """Forecast the 16 hours of `day` given the previous day's actuals (what is known at 18:00... plus closing)."""
        e_prev = np.log1p(prev_day.y.values) - arimax_features(prev_day) @ self.beta
        mu = arimax_features(day) @ self.beta
        eta, v, out_eta, out_v = e_prev[-1], 0.0, [], []
        for j in range(H):
            eta = self.phi[0]*eta + self.phi[1]*e_prev[j]; v = self.phi[0]**2*v + self.sigma**2
            out_eta.append(eta); out_v.append(v)
        m, sd = mu + np.array(out_eta), np.sqrt(out_v)
        zq = {0.1: -1.2816, 0.8: 0.8416, 0.9: 1.2816}
        return np.expm1(m + sd**2/2), {q: np.expm1(m + zq[q]*sd) for q in qs}

def prophet_features(d, cps):
    t = d.t.values/(365.25*H); hrs = d.hour.values.astype(float); w = d.weekend.values
    cols = [np.ones(len(d)), t] + [np.maximum(0, t - c) for c in cps]
    for v in fourier(hrs, 24, 6).values(): cols += [v, v*w]                    # daily shape, different on weekends
    for v in fourier(d.dow.values + hrs/24, 7, 3).values(): cols.append(v)      # weekly
    for v in fourier(d.doy.values, 365.25, 3).values(): cols.append(v)          # yearly
    X = arimax_features(d); keep = [ARIMAX_FEATURES.index(c) for c in ARIMAX_FEATURES if not c.startswith(("hw_", "yr_", "dow_", "trend"))]
    return np.column_stack(cols + [X[:, keep]])

class ProphetStyle:
    """g(t) piecewise-linear trend + Fourier seasonalities + holidays/regressors, MAP fit with a ridge prior on changepoints."""
    def fit(self, d, n_cp=10, tau=0.05):
        t = d.t.values/(365.25*H); self.cps = np.quantile(t[:int(.85*len(t))], np.linspace(.05, .95, n_cp))
        X, z = prophet_features(d, self.cps), np.log1p(d.y.values)
        P = np.eye(X.shape[1])*1e-6; P[2:2+n_cp, 2:2+n_cp] += np.eye(n_cp)*len(d)*tau**2   # prior shrinks trend changes
        self.w = np.linalg.solve(X.T @ X + P, X.T @ z); self.sigma = float(np.std(z - X @ self.w)); return self
    def predict(self, d, qs=(0.1, 0.8, 0.9)):
        m = prophet_features(d, self.cps) @ self.w; zq = {0.1: -1.2816, 0.8: 0.8416, 0.9: 1.2816}
        return np.expm1(m + self.sigma**2/2), {q: np.expm1(m + zq[q]*self.sigma) for q in qs}

# lag features that are known at 18:00 the day before
df["y_lag_day"] = df.y.shift(H); df["y_lag_week"] = df.y.shift(7*H)
df["y_mean_yday"] = df.groupby("day").y.transform("mean").shift(H)
GBM_FEATURES = ["hour", "dow", "weekend", "doy", "temp", "rain", "promo", "holiday", "pre_holiday", "event", "payday",
                "competitor", "y_lag_day", "y_lag_week", "y_mean_yday"]
def gbm(loss="poisson", q=None):
    kw = dict(loss=loss, max_iter=500, learning_rate=0.04, max_leaf_nodes=31, min_samples_leaf=40,
              l2_regularization=1.0, categorical_features=[0, 1], random_state=0)
    if q is not None: kw["quantile"] = q
    return HistGradientBoostingRegressor(**kw)
class GBMForecaster:
    def fit(self, d):
        d = d.dropna(subset=["y_lag_week", "y_mean_yday"])
        self.mean = gbm().fit(d[GBM_FEATURES], d.y)
        self.q = {q: gbm("quantile", q).fit(d[GBM_FEATURES], d.y) for q in (0.1, 0.8, 0.9)}; return self
    def predict(self, d):
        return self.mean.predict(d[GBM_FEATURES]), {q: m.predict(d[GBM_FEATURES]) for q, m in self.q.items()}
print(f"ARIMAX: {len(ARIMAX_FEATURES)} features · Gradient boosting: {len(GBM_FEATURES)} features")
ARIMAX: 55 features · Gradient boosting: 15 features

In a real project you would use maintained libraries. The NumPy versions above follow the same math, so this article runs with only NumPy, pandas and scikit-learn.

# ARIMAX with exact maximum likelihood (pip install statsmodels)
from statsmodels.tsa.statespace.sarimax import SARIMAX
fit = SARIMAX(np.log1p(train.y), exog=X_train, order=(1, 0, 0),
              seasonal_order=(1, 0, 0, 16)).fit(disp=False)
fc = fit.get_forecast(steps=16, exog=X_tomorrow)
mean, ci = np.expm1(fc.predicted_mean), np.expm1(fc.conf_int(alpha=0.2))

# Prophet (pip install prophet)
from prophet import Prophet
m = Prophet(changepoint_prior_scale=0.05, daily_seasonality=True, weekly_seasonality=True)
for c in ["temp", "rain", "promo", "event", "payday"]:
    m.add_regressor(c)
m.add_country_holidays("US")
m.fit(train.rename(columns={"ts": "ds", "y_log": "y"}))
fc = m.predict(future_df)            # yhat, yhat_lower, yhat_upper

# LightGBM (pip install lightgbm)
import lightgbm as lgb
mean_model = lgb.LGBMRegressor(objective="poisson", n_estimators=800, learning_rate=0.03)
p90_model  = lgb.LGBMRegressor(objective="quantile", alpha=0.9, n_estimators=800)
mean_model.fit(train[features], train.y)

✅ 5. Validate like you will forecast

Rolling-origin evaluation

The honest test copies production: stand at a date, train on everything before it, forecast the next period, move forward, repeat. This is rolling-origin (or time-series) cross-validation.

The animation below is the idea. Each amber block is genuinely in the future relative to its indigo training block.

If we write the average score compactly,

\text{Score} = \frac{1}{K} \sum_{k=1}^{K} \frac{1}{|T_k|} \sum_{t\in T_k} L\!\left(y_t,\hat y_{t\mid\mathcal F_{t-h}}\right),

the notation means:

  • K: number of rolling folds;
  • T_k: test timestamps in fold k;
  • |T_k|: number of forecasted observations in that test window;
  • L(\cdot): the error/loss function, such as absolute error;
  • y_t: what actually happened;
  • \hat y_{t\mid\mathcal F_{t-h}}: the forecast for time t made using only information \mathcal F that was available at forecast time t-h;
  • h: forecast horizon.

We use six 28-day folds from July to December 2025 and reproduce the real day-ahead forecasting process inside every fold.

Metrics that fit the decision

Metric Formula Why we use it
MAE \frac{1}{n}\sum \lvert y_t - \hat y_t\rvert In customers per hour; easy to explain
RMSE \sqrt{\frac{1}{n}\sum (y_t - \hat y_t)^2} Punishes the big misses that cause queues
MASE \dfrac{\text{MAE}}{\text{MAE of seasonal naive}} Below 1 means we beat “same hour last week”
Pinball loss \frac{1}{n}\sum \ell_q(y_t, \hat y_{q,t}) Scores the quantile we staff to
80% coverage share of y_t \in [\hat y_{0.1}, \hat y_{0.9}] Should be close to 80% if intervals are honest
👥 Understaffed hours share of hours with fewer cashiers than the actual traffic needed The business metric

We skip MAPE: early-morning hours with very few customers make percentage errors explode.

How to read the table: MAE answers “how many customers off are we, typically?” MASE answers “are we better than simply copying the same hour last week?” Coverage asks whether our uncertainty band is believable. The staffing metrics ask the question the business actually cares about.

A model can therefore have the smallest MAE and still be a poor staffing model if its upper quantiles are badly calibrated.

Show the code: Erlang C staffing rule 
SERVICE_MIN, TARGET_WAIT_MIN, TARGET_SL = 4.0, 3.0, 0.80     # 4-min checkout, 80% of customers wait ≤ 3 min
def erlang_c(c, a):
    """Probability an arriving customer must wait (M/M/c), offered load a = lambda/mu < c."""
    if c <= a: return 1.0
    s = sum(a**k/math.factorial(k) for k in range(c)); top = a**c/math.factorial(c)*c/(c-a)
    return top/(s + top)
def cashiers_needed(lam_per_hour, service=SERVICE_MIN, wait=TARGET_WAIT_MIN, sl=TARGET_SL):
    mu = 60/service; a = lam_per_hour/mu; c = max(1, int(np.ceil(a)))
    while 1 - erlang_c(c, a)*np.exp(-(c*mu - lam_per_hour)*wait/60) < sl: c += 1
    return c
def floor_staff(lam): return max(1, int(np.ceil(lam/30)))    # one floor assistant per 30 customers per hour
CASH = {lam: cashiers_needed(lam) for lam in range(0, 260)}
need = lambda arr: np.array([CASH[int(min(259, max(0, round(v))))] for v in arr])
print("Cashiers needed at 10, 30, 50, 80 customers/hour:", [cashiers_needed(v) for v in (10, 30, 50, 80)])
Cashiers needed at 10, 30, 50, 80 customers/hour: [2, 4, 5, 7]
Show the code: rolling-origin backtest of all models 
origins = pd.date_range("2025-07-01", periods=6, freq="28D")
MODELS = ["Seasonal naive", "ARIMAX", "Prophet-style", "Gradient boosting"]
store = {m: {"mean": [], "p10": [], "p80": [], "p90": []} for m in MODELS}; actual, stamp, fold_id = [], [], []
t0 = time.time()
for k, o in enumerate(origins):
    tr = df[df.ts < o]; te = df[(df.ts >= o) & (df.ts < o + pd.Timedelta(days=28))]
    A, P, G = ARIMAX().fit(tr), ProphetStyle().fit(tr), GBMForecaster().fit(tr)
    a_parts = [A.forecast_day(df[df.day == dd - 1], df[df.day == dd]) for dd in sorted(te.day.unique())]
    outs = {"Seasonal naive": (te.y_lag_week.values, {0.1: .78*te.y_lag_week.values, 0.8: 1.12*te.y_lag_week.values, 0.9: 1.22*te.y_lag_week.values}),
            "ARIMAX": (np.concatenate([p[0] for p in a_parts]), {q: np.concatenate([p[1][q] for p in a_parts]) for q in (0.1, 0.8, 0.9)}),
            "Prophet-style": P.predict(te), "Gradient boosting": G.predict(te)}
    for m, (mean, qd) in outs.items():
        store[m]["mean"].append(mean); store[m]["p10"].append(qd[0.1]); store[m]["p80"].append(qd[0.8]); store[m]["p90"].append(qd[0.9])
    actual.append(te.y.values); stamp.append(te.ts.values); fold_id.append(np.full(len(te), k + 1))
y_bt = np.concatenate(actual); ts_bt = np.concatenate(stamp); fold_bt = np.concatenate(fold_id)
lam_bt = df.set_index("ts").loc[ts_bt, "lam"].values
for m in MODELS:
    for key in store[m]: store[m][key] = np.concatenate(store[m][key])

def pinball(y, f, q): d = y - f; return np.mean(np.maximum(q*d, (q-1)*d))
need_actual = need(y_bt)
mae_naive = np.mean(np.abs(store["Seasonal naive"]["mean"] - y_bt))
rows = []
for m in MODELS:
    s = store[m]; e = s["mean"] - y_bt; plan = need(s["p80"])
    rows.append(dict(Model=m, MAE=np.mean(np.abs(e)), RMSE=np.sqrt(np.mean(e**2)), MASE=np.mean(np.abs(e))/mae_naive,
                     Pinball80=pinball(y_bt, s["p80"], .8), Coverage80=np.mean((y_bt >= s["p10"]) & (y_bt <= s["p90"])),
                     Understaffed=np.mean(plan < need_actual), ExtraStaffHours=np.mean(np.maximum(plan - need_actual, 0))))
res = pd.DataFrame(rows)
oracle_mae = np.mean(np.abs(y_bt - lam_bt))          # best possible: we know the true rate
print(f"Backtest: {len(y_bt):,} forecast hours in {len(origins)} folds · {time.time()-t0:.0f}s · noise floor MAE = {oracle_mae:.2f}")
Backtest: 2,688 forecast hours in 6 folds · 190s · noise floor MAE = 5.36
Show the code 
from IPython.display import Markdown
best = {"MAE": res.MAE.min(), "RMSE": res.RMSE.min(), "MASE": res.MASE.min(), "Pinball80": res.Pinball80.min(),
        "Understaffed": res.Understaffed.min(), "ExtraStaffHours": res.ExtraStaffHours.min()}
cov_best = res.Coverage80.iloc[(res.Coverage80 - .8).abs().argmin()]
def b(v, col, fmt): s = format(v, fmt); return f"**{s}**" if (col == "Coverage80" and v == cov_best) or (col in best and v == best[col]) else s
lines = ["| Model | MAE | RMSE | MASE | Pinball (q=0.8) | 80% coverage | Understaffed hours | Extra staff per hour |", "|---|---|---|---|---|---|---|---|"]
for _, r0 in res.iterrows():
    lines.append(f"| {r0.Model} | {b(r0.MAE,'MAE','.2f')} | {b(r0.RMSE,'RMSE','.2f')} | {b(r0.MASE,'MASE','.3f')} | {b(r0.Pinball80,'Pinball80','.2f')} | "
                 f"{b(r0.Coverage80,'Coverage80','.0%')} | {b(r0.Understaffed,'Understaffed','.1%')} | {b(r0.ExtraStaffHours,'ExtraStaffHours','.2f')} |")
Markdown("\n".join(lines))
Table 1: Rolling-origin results, July–December 2025 (day-ahead forecasts). Best value in each column in bold.
Model MAE RMSE MASE Pinball (q=0.8) 80% coverage Understaffed hours Extra staff per hour
Seasonal naive 11.05 16.12 1.000 5.04 36% 23.2% 0.55
ARIMAX 6.27 8.92 0.568 2.38 83% 7.2% 0.56
Prophet-style 6.34 8.96 0.573 2.40 83% 6.7% 0.59
Gradient boosting 6.87 10.09 0.622 2.68 73% 9.9% 0.49
Show the code 
fig = go.Figure()
for m, c in zip(MODELS, [COLORS["slate"], COLORS["indigo"], COLORS["amber"], COLORS["teal"]]):
    vals = [np.mean(np.abs(store[m]["mean"][fold_bt == k] - y_bt[fold_bt == k])) for k in range(1, 7)]
    fig.add_trace(go.Scatter(x=[f"Fold {k}" for k in range(1, 7)], y=vals, name=m, mode="lines+markers", line=dict(color=c, width=3), marker=dict(size=9), hovertemplate="%{x}<br>MAE: %{y:.2f} customers/hour<extra>" + m + "</extra>"))
fig.add_hline(y=oracle_mae, line=dict(dash="dash", color=COLORS["rose"]), annotation_text="Noise floor (true rate known)", annotation_position="bottom right")
fig.update_layout(height=360, margin=dict(l=60, r=20, t=20, b=40), yaxis_title="MAE (customers per hour)", legend=dict(orientation="h", y=1.12))
fig.show(config=PLOTLY_CONFIG)
Figure 6: MAE per fold. A model you can trust is good in every fold, not just on average. The dashed line is the noise floor: the error even a perfect model would make, because arrivals are random.

The noise floor is the most useful line on this chart. Arrivals are random even when the rate is known: at 40 customers per hour a Poisson count has a standard deviation of about \sqrt{40}\approx 6.3. Our best model’s MAE of 6.27 sits close to the floor of 5.36, so most of the remaining error cannot be modelled away. It has to be handled by staffing to a quantile instead of the mean.

Compare the forecasts for one test week. Use the buttons to switch models.

Show the code 
wk0 = pd.Timestamp("2025-11-17"); sel = (ts_bt >= np.datetime64(wk0)) & (ts_bt < np.datetime64(wk0 + pd.Timedelta(days=7)))
x_lab = pd.to_datetime(ts_bt[sel]).strftime("%a %H:00"); xi = np.arange(sel.sum())
fig = go.Figure(); n_per = 3
for i, (m, c) in enumerate(zip(MODELS[1:], [COLORS["indigo"], COLORS["amber"], COLORS["teal"]])):
    s = store[m]; vis = i == 0
    fig.add_trace(go.Scatter(x=np.r_[xi, xi[::-1]], y=np.r_[s["p90"][sel], s["p10"][sel][::-1]], fill="toself", fillcolor=c, opacity=.18, line=dict(width=0), name="80% interval", visible=vis, hoverinfo="skip"))
    fig.add_trace(go.Scatter(x=xi, y=s["mean"][sel], line=dict(color=c, width=2.5), name=f"{m} forecast", visible=vis, text=x_lab, hovertemplate="%{text}<br>forecast %{y:.0f}<extra></extra>"))
    fig.add_trace(go.Scatter(x=xi, y=y_bt[sel], mode="markers", marker=dict(color=COLORS["rose"], size=5), name="Actual", visible=vis, text=x_lab, hovertemplate="%{text}<br>actual %{y}<extra></extra>"))
btns = [dict(label=m, method="update", args=[{"visible": [j // n_per == i for j in range(3*n_per)]}]) for i, m in enumerate(MODELS[1:])]
fig.update_layout(height=400, margin=dict(l=60, r=20, t=60, b=40), yaxis_title="Customers per hour",
                  xaxis=dict(tickvals=xi[::16], ticktext=[x_lab[i].split()[0] for i in range(0, len(xi), 16)]),
                  updatemenus=[dict(type="buttons", direction="right", x=0, y=1.18, xanchor="left", buttons=btns)], legend=dict(orientation="h", y=-.15))
fig.show(config=PLOTLY_CONFIG)
Figure 7: One test week (Mon–Sun, November 2025): actual arrivals (dots), mean forecast (line) and 80% interval (band) for each model.

🏆 6. Choosing the model

Accuracy is one input. A forecast that feeds a weekly schedule also needs honest uncertainty, stability and something a manager can understand.

Show the code: the selection scorecard 
r_ = res.set_index("Model")
best_acc = r_.drop("Seasonal naive").MAE.idxmin()
best_staff = r_.drop("Seasonal naive").Understaffed.idxmin()
cov_err = (r_.Coverage80 - .8).abs()
champion = ARIMAX().fit(df[df.ts < pd.Timestamp("2025-12-18")])   # refit on all data up to the release date
gap = 100*(r_.loc["ARIMAX", "MAE"]/r_.drop("Seasonal naive").MAE.min() - 1)
print("Most accurate (MAE):     ", best_acc)
print("Fewest understaffed hrs: ", best_staff)
print("Best calibrated (80%):   ", cov_err.drop("Seasonal naive").idxmin())
Most accurate (MAE):      ARIMAX
Fewest understaffed hrs:  Prophet-style
Best calibrated (80%):    Prophet-style
Criterion 📐 ARIMAX 📈 Prophet-style 🌲 Gradient boosting
Accuracy (MAE, this backtest) 6.27 6.34 6.87
Interval calibration (target 80%) 83% 83% 73%
Understaffed hours at P80 7.2% 6.7% 9.9%
Uses the latest surprise (intraday updates) ✅ AR errors ❌ independent errors ⚠️ only through lags we allow
Explains effects to a manager ✅ one coefficient per driver ✅ plots of trend and seasonality ⚠️ needs SHAP or partial dependence
Learns new interactions by itself ❌ hand-engineered ❌ hand-engineered ✅ automatically
Size and serving ✅ about 60 numbers in JSON ✅ small ⚠️ several tree ensembles
Scales to 500 stores ⚠️ one model per store ⚠️ one model per store ✅ one global model
Tip🧑‍🔬 My choice: ARIMAX in production, gradient boosting as the challenger

For a single store, I would ship the ARIMAX (regression with seasonal AR errors). In this backtest its MAE is the best, its staffing record is among the best (7.2% understaffed hours against 6.7% for the best model), its 80% intervals contain 83% of actuals, so the staffing quantiles can be trusted, it can update during the day when the morning is busier than expected, and every coefficient has a business meaning (“rain on an evening costs 23% of traffic”). The whole model is a short JSON file, so the scheduling app does not even need Python.

The Prophet-style model is a fine baseline that analysts can maintain, but it ignores autocorrelation. Gradient boosting becomes my first choice as soon as FreshCart forecasts many stores with one global model, or when new drivers with complex interactions appear. Until then it runs as a challenger in shadow mode, and replaces the champion only if it wins over several months.

The fitted ARIMAX coefficients show why managers like it: on the log scale, e^{\beta}-1 is the percentage effect of each driver.

Show the code 
names = {"promo": "🏷️ Promotion", "promo_we": "🏷️ + on weekends", "promo_eve": "🏷️ + after 17:00", "rain": "🌧️ Rain", "rain_eve": "🌧️ + after 16:00",
         "event_pm": "🏟️ Event day, 16–19h", "event_late": "🏟️ Event day, after 20h", "payday": "💵 Payday", "holiday": "🎄 Public holiday",
         "prehol_pm": "🎄 Day before, after 14h", "comp_wd": "🏪 Competitor (weekdays)", "heat_aft": "🔥 Per +10°C above 20, afternoon", "heat_eve": "🔥 Per +10°C above 20, evening"}
co = pd.DataFrame({"name": [names[k] for k in names], "eff": [100*(np.exp(champion.beta[ARIMAX_FEATURES.index(k)])-1) for k in names]}).sort_values("eff")
fig = go.Figure(go.Bar(x=co.eff, y=co.name, orientation="h", marker_color=[COLORS["rose"] if v < 0 else COLORS["indigo"] for v in co.eff],
                       text=[f"{v:+.0f}%" for v in co.eff], textposition="outside", cliponaxis=False, hovertemplate="%{y}: %{x:+.1f}%<extra></extra>"))
fig.update_layout(height=470, margin=dict(l=220, r=50, t=20, b=40), xaxis_title="Effect on expected customers (%)",
                  xaxis=dict(range=[co.eff.min() - 18, co.eff.max() + 15]))
fig.show(config=PLOTLY_CONFIG)
print(f"phi_1  = {champion.phi[0]:.2f}   memory of the previous hour")
print(f"phi_16 = {champion.phi[1]:.2f}   same hour yesterday")
print(f"sigma  = {champion.sigma:.3f}  innovation s.d. on the log scale")
Figure 8: Effects learned by the ARIMAX model, as % change in expected customers. Hover for exact values. Interaction terms apply on top of the main effect. 
phi_1  = 0.17   memory of the previous hour
phi_16 = 0.03   same hour yesterday
sigma  = 0.284  innovation s.d. on the log scale

👥 7. From forecast to staff: queueing math

A forecast of “42 customers” does not say how many cashiers to schedule. Queueing theory translates demand into capacity.

We approximate the checkout area as an M/M/c queue. You do not need queueing theory beforehand; the notation means:

  • \lambda (“lambda”): arriving customers per hour;
  • \bar s: average checkout time in minutes;
  • \mu=60/\bar s: customers one cashier can serve per hour;
  • c: number of cashiers working in parallel;
  • a=\lambda/\mu: offered workload, measured in cashier-equivalents.

For example, if checkout takes 4 minutes, one cashier can serve \mu=60/4=15 customers/hour. If \lambda=42, the offered load is a=42/15=2.8: fewer than 3 cashiers cannot even keep up on average.

The Erlang C formula then estimates the probability that an arriving customer must wait:

C(c, a) = \frac{\dfrac{a^c}{c!}\dfrac{c}{c-a}}{\displaystyle\sum_{k=0}^{c-1}\frac{a^k}{k!} + \frac{a^c}{c!}\frac{c}{c-a}}, \qquad c > a,

and the service level, the share of customers who wait at most T minutes, is

\operatorname{SL}(c) = 1 - C(c,a)\,e^{-(c\mu - \lambda)\,T/60}.

We choose the smallest c with \operatorname{SL}(c) \ge 80\%.

Which forecast do we plug in? Not the mean. Staffing is a newsvendor problem: if one understaffed hour costs c_u (lost margin and goodwill) and one surplus staff-hour costs c_o (wages), the cost-minimising plan uses the forecast quantile

q^\star = \frac{c_u}{c_u + c_o}.

With c_u = 80 dollars and c_o = 20 dollars, q^\star = 80/100 = 0.8: staff to the P80 forecast. Try the planner below.

Show the code: export tomorrow’s forecast distribution for the planner 
plan_day = df[df.ts.dt.date == pd.Timestamp("2025-12-19").date()]
pm, pq = champion.forecast_day(df[df.day == plan_day.day.iloc[0]-1], plan_day)
zgrid = np.round(np.linspace(-1.6, 1.6, 33), 2)
e_prev = np.log1p(df[df.day == plan_day.day.iloc[0]-1].y.values) - arimax_features(df[df.day == plan_day.day.iloc[0]-1]) @ champion.beta
mu_log = arimax_features(plan_day) @ champion.beta; eta, v, mlog, sd = e_prev[-1], 0, [], []
for j in range(H):
    eta = champion.phi[0]*eta + champion.phi[1]*e_prev[j]; v = champion.phi[0]**2*v + champion.sigma**2; mlog.append(mu_log[j] + eta); sd.append(np.sqrt(v))
ojs_define(plan=[{"hour": int(h), "mlog": float(m), "sd": float(s)} for h, m, s in zip(plan_day.hour, mlog, sd)],
           plan_label=f"{plan_day.ts.iloc[0]:%A %d %B %Y}")

💾 8. Saving the model

A saved model is more than its weights. Six months later someone must be able to answer: what is this, what data trained it, how good was it, and what inputs does it expect? A good artifact therefore contains:

Part Contents Why
🧮 Parameters \boldsymbol\beta, \phi_1, \phi_{16}, \sigma The model itself
🧾 Feature contract ordered feature names, units, the opening hours, the transform \log(1+y) Inputs built the same way in training and serving
🗓️ Lineage training window, data snapshot hash, code version Reproducibility and audits
📏 Validation backtest MAE, MASE, coverage The bar a retrained model must clear
🔐 Integrity SHA-256 checksum of the parameters Detects a corrupted or edited file

For linear and statistical models, store parameters as JSON. It is readable, diff-able, language-neutral and safe to load. For tree ensembles use joblib or skops, and never unpickle files from untrusted sources, because pickle can run arbitrary code.

Show the code: save a versioned, checksummed model artifact 
art_dir = pathlib.Path(tempfile.mkdtemp()) / "models" / "footfall-arimax" / "v2025.12.18"
art_dir.mkdir(parents=True)
params = {"beta": [round(float(b_), 6) for b_ in champion.beta], "phi": [round(float(p_), 6) for p_ in champion.phi], "sigma": round(champion.sigma, 6)}
artifact = {
    "model": "footfall-arimax", "version": "2025.12.18", "store_id": "FC-001",
    "target": {"name": "customers_per_hour", "transform": "log1p"},
    "schedule": {"open_hour": OPEN, "close_hour": CLOSE, "hours_per_day": H},
    "features": ARIMAX_FEATURES,
    "feature_spec": {"temp": "(forecast_high_c - 20) / 10", "trend": "hour_index / (365.25 * 16)", "yr_*": "Fourier of day-of-year, period 365.25"},
    "parameters": params,
    "training": {"from": str(df.ts.min().date()), "to": "2025-12-17", "rows": int((df.ts < pd.Timestamp('2025-12-18')).sum()),
                 "data_sha256": hashlib.sha256(df[df.ts < pd.Timestamp('2025-12-18')].y.values.tobytes()).hexdigest()[:16]},
    "validation": {"scheme": "rolling origin, 6 x 28 days", "mae": round(float(r_.loc["ARIMAX", "MAE"]), 3),
                   "mase": round(float(r_.loc["ARIMAX", "MASE"]), 3), "coverage80": round(float(r_.loc["ARIMAX", "Coverage80"]), 3)},
}
artifact["checksum"] = hashlib.sha256(json.dumps(params, sort_keys=True).encode()).hexdigest()
(art_dir / "model.json").write_text(json.dumps(artifact, indent=2))
(art_dir / "MODEL_CARD.md").write_text(f"# footfall-arimax {artifact['version']}\nIntended use: next-day hourly staffing for store FC-001.\n"
                                       f"Not for: stores without 12 months of history, or forecasting beyond 2 days.\n")
size = (art_dir / "model.json").stat().st_size
print(f"Saved {art_dir.relative_to(art_dir.parents[2])}/model.json  ({size/1024:.1f} KB, {len(params['beta'])} coefficients)")
print(json.dumps({k: artifact[k] for k in ("model", "version", "training", "validation")}, indent=2))
Saved models/footfall-arimax/v2025.12.18/model.json  (2.5 KB, 55 coefficients)
{
  "model": "footfall-arimax",
  "version": "2025.12.18",
  "training": {
    "from": "2024-01-01",
    "to": "2025-12-17",
    "rows": 11472,
    "data_sha256": "8f0fc8c3d405a000"
  },
  "validation": {
    "scheme": "rolling origin, 6 x 28 days",
    "mae": 6.275,
    "mase": 0.568,
    "coverage80": 0.831
  }
}
models/
└── footfall-arimax/
    ├── v2025.12.18/          ← immutable: never edited after release
    │   ├── model.json        ← parameters + feature contract + lineage + checksum
    │   └── MODEL_CARD.md     ← intended use, limits, owner
    └── production -> v2025.12.18   ← a pointer; rollback = move the pointer

Each release gets a new folder; promoting or rolling back a model is just moving the production pointer. The model versioning article covers registries and safe rollouts in more depth.

Production handoff · different role

🚚 9. Shipped: loading the model and forecasting

Now switch roles. You are the developer of FreshCart’s scheduling app. You did not train the model; you received model.json. Your service must load it, verify it, build features exactly as specified, and predict.

Note🧑‍💻 Role change: training team → serving team

Everything above this point was about learning a model. This section is intentionally styled differently because production inference has a different responsibility.

The serving application should not know how ARIMAX was trained. It should know only the artifact contract:

inputs → feature builder → released parameters → forecast → staffing output

This separation prevents accidental retraining, hidden notebook state and training-only features from leaking into production.

Show the code: the inference service, which knows nothing about training 
class FootfallService:
    """Loads a released model.json and produces staffing-ready forecasts. No training code needed."""
    def __init__(self, path):
        self.a = json.loads(pathlib.Path(path).read_text())
        p = self.a["parameters"]
        if hashlib.sha256(json.dumps(p, sort_keys=True).encode()).hexdigest() != self.a["checksum"]:
            raise ValueError("Checksum mismatch: model file is corrupted or was edited")
        if self.a["features"] != ARIMAX_FEATURES:
            raise ValueError("Feature contract changed: refuse to serve")
        self.beta, self.phi, self.sigma = np.array(p["beta"]), np.array(p["phi"]), p["sigma"]
    def forecast(self, yesterday, tomorrow, q=0.8):
        e_prev = np.log1p(yesterday.y.values) - arimax_features(yesterday) @ self.beta
        mu = arimax_features(tomorrow) @ self.beta; eta, v, out = e_prev[-1], 0.0, []
        zq = {0.1: -1.2816, 0.5: 0.0, 0.8: 0.8416, 0.9: 1.2816}[q]
        for j in range(len(tomorrow)):
            eta = self.phi[0]*eta + self.phi[1]*e_prev[j]; v = self.phi[0]**2*v + self.sigma**2; s = np.sqrt(v)
            out.append((mu[j] + eta, s))
        m, s = np.array(out).T
        return pd.DataFrame({"hour": tomorrow.hour.values, "mean": np.expm1(m + s**2/2), "p10": np.expm1(m - 1.2816*s),
                             "p90": np.expm1(m + 1.2816*s), f"p{int(q*100)}": np.expm1(m + zq*s)})

svc = FootfallService(art_dir / "model.json")
print(f"Loaded {svc.a['model']} v{svc.a['version']} · checksum OK · {len(svc.beta)} coefficients")

# forecast the 14 held-out days: 18-31 December 2025 (never seen in training)
hold = []
for dd in sorted(df[df.ts >= pd.Timestamp("2025-12-18")].day.unique()):
    tom, yday = df[df.day == dd], df[df.day == dd - 1]
    f = svc.forecast(yday, tom); f["actual"] = tom.y.values; f["lam"] = tom.lam.values; f["date"] = f"{tom.ts.iloc[0]:%a %d %b}"
    f["need_actual"] = need(f.actual); f["plan"] = need(f.p80)
    hold.append(f)
hold = pd.concat(hold, ignore_index=True)
e = hold["mean"] - hold.actual
print("Held-out 18-31 Dec 2025 (14 days, never seen in training)")
print(f"  MAE               {np.mean(np.abs(e)):.2f} customers/hour")
print(f"  Bias              {e.mean():+.2f} customers/hour")
print(f"  80% coverage      {np.mean((hold.actual >= hold.p10) & (hold.actual <= hold.p90)):.0%}")
print(f"  Understaffed hrs  {np.mean(hold.plan < hold.need_actual):.1%}")
Loaded footfall-arimax v2025.12.18 · checksum OK · 55 coefficients
Held-out 18-31 Dec 2025 (14 days, never seen in training)
  MAE               6.76 customers/hour
  Bias              +1.44 customers/hour
  80% coverage      78%
  Understaffed hrs  6.2%

🎛️ The production console

This console runs the same model.json in your browser: it reads the coefficients, rebuilds the features from the inputs and computes the forecast. Press Play to step through the 14 held-out days, including Christmas. Watch the inputs change, the forecast respond, and the actual arrivals land inside or outside the band. Switch to What-if to set the inputs yourself.

Show the code: hand the artifact and held-out days to the browser 
days_out = []
for dd in sorted(df[df.ts >= pd.Timestamp("2025-12-18")].day.unique()):
    tom, yday = df[df.day == dd], df[df.day == dd - 1]
    e_prev = np.log1p(yday.y.values) - arimax_features(yday) @ svc.beta
    r0 = tom.iloc[0]
    days_out.append({"date": f"{r0.ts:%a %d %b}", "dow": int(r0.dow), "doy": int(r0.doy), "t0": int(r0.t), "temp": round(float(r0.temp), 1),
                     "rain": int(r0.rain), "promo": int(r0.promo), "event": int(r0.event), "payday": int(r0.payday), "holiday": int(r0.holiday),
                     "pre_holiday": int(r0.pre_holiday), "competitor": int(r0.competitor),
                     "e_prev": [round(float(x), 4) for x in e_prev], "actual": [int(x) for x in tom.y.values]})
ojs_define(model_json=json.loads((art_dir / "model.json").read_text()), held_days=days_out)

🔍 Understanding the error of the shipped model

A shipped forecast is judged on fresh data. Three questions matter: how big is the error, is it biased, and where does it come from?

Show the code 
by_h = hold.groupby("hour").apply(lambda g: pd.Series({"mae": np.mean(np.abs(g["mean"] - g.actual)), "floor": np.mean(np.abs(g.actual - g.lam))}), include_groups=False).reset_index()
fig = make_subplots(rows=1, cols=2, subplot_titles=("MAE by hour (held-out)", "Error distribution"), column_widths=[.58, .42])
fig.add_trace(go.Bar(x=by_h.hour, y=by_h.mae, name="Model MAE", marker_color=COLORS["indigo"], hovertemplate="Hour %{x}:00<br>Model MAE: %{y:.2f}<extra></extra>"), 1, 1)
fig.add_trace(go.Scatter(x=by_h.hour, y=by_h.floor, name="Noise floor", mode="lines+markers", line=dict(color=COLORS["rose"], dash="dash"), hovertemplate="Hour %{x}:00<br>Irreducible MAE: %{y:.2f}<extra></extra>"), 1, 1)
fig.add_trace(go.Histogram(x=hold["mean"] - hold.actual, nbinsx=30, marker_color=COLORS["teal"], name="Errors", showlegend=False, hovertemplate="Forecast − actual: %{x:.1f}<br>Hours in bin: %{y}<extra></extra>"), 1, 2)
fig.add_vline(x=0, line=dict(color=COLORS["slate"], dash="dot"), row=1, col=2)
fig.update_xaxes(title_text="Hour", row=1, col=1); fig.update_xaxes(title_text="Forecast − actual", row=1, col=2)
fig.update_yaxes(title_text="Customers", row=1, col=1)
fig.update_layout(height=360, margin=dict(l=60, r=20, t=50, b=50), legend=dict(orientation="h", y=1.18))
fig.show(config=PLOTLY_CONFIG)
xmas = hold[hold.date.str.contains("24 Dec|25 Dec")].groupby("date").apply(lambda g: np.mean(np.abs(g["mean"] - g.actual)), include_groups=False)
print("MAE on 24 and 25 December:", xmas.round(1).to_dict())
Figure 9: Held-out error analysis. Left: MAE by hour versus the noise floor. Right: distribution of forecast errors (forecast − actual); a centred, symmetric shape means no systematic bias. 
MAE on 24 and 25 December: {'Thu 25 Dec': 6.1, 'Wed 24 Dec': 4.4}

The error has three sources, and each calls for a different response:

  1. Irreducible randomness 🎲. Even with the true rate known, counts scatter by roughly \sqrt{\lambda}. Error grows in busy hours for this reason, and it is why we staff to P80 rather than to the mean.
  2. Model error 🧮. The gap between the bars and the dashed floor. This is what better features or models can still win.
  3. Rare events 🎄. Days like Christmas Eve have few examples in two years of history, so their effects are estimated from only one or two occurrences each. For such days, planners should review the forecast and add a manual adjustment.

Errors are also correlated within a day: a cold, quiet morning usually means a quiet afternoon. That is the AR term at work. In production, re-forecasting at noon with the morning’s actual counts tightens the afternoon interval.

📡 10. After launch: monitor and retrain

Shipping is the start, not the end. The monitoring article covers this in depth; for this model the essentials are:

Check Rule Action
📏 Weekly MAE vs. backtest Alert if 4-week MAE > 1.2 × the artifact’s validation.mae Investigate, then retrain
⚖️ Bias Alert if mean error stays beyond ±2 customers/hour for 2 weeks A level shift, such as a new competitor or a road closure
🎯 Coverage Alert if 80% band coverage falls below 70% Intervals too narrow; refit σ
🧾 Input contract Reject requests with missing or out-of-range drivers Fix the upstream feed
🔄 Retraining Monthly on a rolling 2-year window, and only promote if the backtest is not worse New version folder, move the pointer

🧩 Summary

  • Frame the forecast around the decision: horizon, granularity, available inputs and the cost of being wrong come before any model.
  • Ordinary regression with random splits fails on time series: autocorrelated errors, leakage from the future, count data and structural breaks.
  • ARIMAX adds memory to regression, Prophet decomposes trend, seasonality and holidays, and gradient boosting learns interactions from lag features. Each has a clean mathematical core.
  • Validate with rolling origins, report MASE against a seasonal naive baseline, check interval coverage and measure the business metric.
  • Compare against the noise floor to know when to stop modelling and start managing uncertainty.
  • Turn forecasts into staff with Erlang C and the newsvendor quantile q^\star = c_u/(c_u + c_o).
  • Save a contract, not just weights: parameters, feature spec, lineage, validation and checksum, in an immutable version folder.
  • Serve and monitor: the scheduling app loads the artifact, verifies it, forecasts, and alerts when error, bias or coverage drift.
Back to top