Forecasting Sales
What happens next?.
How many coffee beans should Harbor Coffee order for next month? Order too little and shelves go empty. Order too much, and it goes stale. Good forecasts directly save money, and that's what makes this one of the most practical skills in data science.
- How many beans for next month?
- Too little or too much
- Forecasts save money
Splitting a series into parts
A time series is usually made of a trend, a repeating pattern, and noise. Seasonal decompose splits them apart. We use a seven day period. Run it. There's the steady upward trend. There's the weekly rhythm, with its strong weekends. And there's what's left. Seeing the parts makes it clear what a forecast has to capture.
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
sales = pd.read_csv("sales.csv", parse_dates=["date"])
daily = sales.groupby("date")["revenue"].sum().asfreq("D")
parts = seasonal_decompose(daily, period=7, model="multiplicative")
fig, ax = plt.subplots(3, 1, figsize=(9, 5), sharex=True)
for a, s, n in zip(ax, [parts.trend, parts.seasonal, parts.resid], ["Trend", "Weekly", "Noise"]):
a.plot(s); a.set_ylabel(n)
plt.tight_layout(); plt.savefig("parts.png", dpi=110)
print("Strongest weekday effect:", round(parts.seasonal.max(), 2))Strongest weekday effect: 1.17
Baselines first
Before any model, set some simple baselines. Tomorrow equals today. The recent average. And the same day last week. Each one is scored by average percentage error. We hold back the last eight weeks. Run it. The same day last week does far better, as it captures the weekly pattern. This is the score that any real model must beat.
import numpy as np
import pandas as pd
sales = pd.read_csv("sales.csv", parse_dates=["date"])
y = sales.groupby("date")["revenue"].sum().asfreq("D")
train, test = y[:-56], y[-56:]
mape = lambda p: (np.abs(test.values - p) / test.values).mean()
print(f"Last value: {mape(np.full(56, train.iloc[-1])):.1%}")
print(f"Mean of last 28: {mape(np.full(56, train.iloc[-28:].mean())):.1%}")
seasonal = np.tile(train.iloc[-7:].values, 8)
print(f"Same day last week: {mape(seasonal):.1%}")Last value: 11.2% Mean of last 28: 10.1% Same day last week: 7.1%
Exponential smoothing
Holt Winters is a classic forecasting method, which tracks a level, a trend, and a season. It gives more weight to recent days. Run it. It beats our baselines, and it needs only one line. Classic methods like this are fast, stable, and easy to explain, so they're often all that you need.
import numpy as np
import pandas as pd
from statsmodels.tsa.holtwinters import ExponentialSmoothing
sales = pd.read_csv("sales.csv", parse_dates=["date"])
y = sales.groupby("date")["revenue"].sum().asfreq("D")
train, test = y[:-56], y[-56:]
model = ExponentialSmoothing(train, trend="add", seasonal="mul", seasonal_periods=7).fit()
pred = model.forecast(56)
print(f"Error: {(np.abs(test - pred) / test).mean():.1%}")
print("Next week total:", round(pred.iloc[:7].sum()))Error: 2.8% Next week total: 74804
Boosting with lags
Machine learning can do this too, by turning time into a table of lag features. We use last week, four weeks ago, a recent average, and calendar columns. Run it. It's a good score, but the simpler model wins here, because trees can't project a trend beyond what they've seen. Boosting really shines with many series, and with extra data, such as promotions or weather.
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
sales = pd.read_csv("sales.csv", parse_dates=["date"])
y = sales.groupby("date")["revenue"].sum().asfreq("D").rename("y")
df = y.to_frame()
df["lag_7"], df["lag_28"] = df["y"].shift(7), df["y"].shift(28)
df["mean_28"] = df["y"].shift(7).rolling(28).mean()
df["dow"], df["doy"] = df.index.dayofweek, df.index.dayofyear
df = df.dropna()
train, test = df[:-56], df[-56:]
cols = ["lag_7", "lag_28", "mean_28", "dow", "doy"]
m = HistGradientBoostingRegressor(random_state=0).fit(train[cols], train["y"])
print(f"Boosting error: {(np.abs(test['y'] - m.predict(test[cols])) / test['y']).mean():.1%}")Boosting error: 4.0%
Reporting a range
A forecast without a range is a guess that's dressed up. We use real past errors, to build a range for next week. Run it. The business can now order for the top of the range, or plan for the bottom. That's far more useful than one number.
import numpy as np
import pandas as pd
sales = pd.read_csv("sales.csv", parse_dates=["date"])
y = sales.groupby("date")["revenue"].sum().asfreq("D")
week = y.rolling(7).sum()
errors = (week / week.shift(7) - 1).dropna()
low, high = errors.iloc[-180:].quantile([0.05, 0.95])
forecast = week.iloc[-1]
print(f"Next week forecast: {forecast:,.0f}")
print(f"90% range: {forecast * (1 + low):,.0f} to {forecast * (1 + high):,.0f}")Next week forecast: 82,345 90% range: 81,963 to 85,520
Summary
Start with baselines, try classic methods, then boosting, and always give a range. Now that we've got answers, how do we get people to believe them? Next, the skill that separates great data scientists from good ones. Communication.
Key takeaways- Trend, season, noise
- Beat the baselines
- Forecast with a range