# Colab: install the one library the rubric needs that isn't preinstalled.
!pip -q install pmdarima[?25l [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m0.0/689.1 kB[0m [31m?[0m eta [36m-:--:--[0m [2K [91m━━━━━━━━━━━━━━━━[0m[90m╺[0m[90m━━━━━━━━━━━━━━━━━━━━━━━[0m [32m276.5/689.1 kB[0m [31m8.1 MB/s[0m eta [36m0:00:01[0m [2K [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m689.1/689.1 kB[0m [31m11.3 MB/s[0m eta [36m0:00:00[0m [?25h
!pip -q install keras-tuner
import warnings; warnings.filterwarnings("ignore")
import numpy as np, pandas as pd, matplotlib.pyplot as plt
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.tsa.stattools import adfuller, kpss
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import pmdarima as pm
from xgboost import XGBRegressor
from sklearn.model_selection import GridSearchCV, TimeSeriesSplit
from sklearn.preprocessing import MinMaxScaler
import tensorflow as tf
from tensorflow import keras
from tensorflow.keras import layers
import keras_tuner as kt[?25l [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m0.0/129.4 kB[0m [31m?[0m eta [36m-:--:--[0m [2K [90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━[0m [32m129.4/129.4 kB[0m [31m4.1 MB/s[0m eta [36m0:00:00[0m [?25h
# Both workbooks have four tabs; read every tab and concatenate.
weekly = pd.concat(pd.read_excel("UK Weekly Trended Timeline from 200101_202429.xlsx",
sheet_name=None, dtype={"ISBN": str}).values(), ignore_index=True)
isbns = pd.concat(pd.read_excel("ISBN List.xlsx",
sheet_name=None, dtype={"ISBN": str}).values(), ignore_index=True)
weekly.columns = weekly.columns.str.replace(" ", "_")
print(weekly.shape, "rows;", weekly["ISBN"].nunique(), "ISBNs")
weekly.head(3)(227224, 13) rows; 500 ISBNs
| ISBN | Title | Author | Interval | End_Date | Volume | Value | ASP | RRP | Binding | Imprint | Publisher_Group | Product_Class | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 9780002261821 | One For My Baby | Parsons, Tony | 200513 | 2005-04-02 | 1 | 15.99 | 15.99 | 15.99 | Hardback | HarperCollins Publishers | HarperCollins Grp | F1.1 General & Literary Fiction |
| 1 | 9780002261821 | One For My Baby | Parsons, Tony | 200503 | 2005-01-22 | 1 | 15.99 | 15.99 | 15.99 | Hardback | HarperCollins Publishers | HarperCollins Grp | F1.1 General & Literary Fiction |
| 2 | 9780002261821 | One For My Baby | Parsons, Tony | 200422 | 2004-05-29 | 1 | 11.19 | 11.19 | 15.99 | Hardback | HarperCollins Publishers | HarperCollins Grp | F1.1 General & Literary Fiction |
FREQ = "W-SAT" # weekly grid: week ending Saturday
weekly["ISBN"] = weekly["ISBN"].astype(str)
weekly["End_Date"] = pd.to_datetime(weekly["End_Date"])
weekly["Volume"] = pd.to_numeric(weekly["Volume"], errors="coerce").fillna(0)
def series(isbn, start=None):
# weekly units for one ISBN on a fixed Saturday grid, silent weeks = 0
s = weekly.loc[weekly["ISBN"] == isbn].set_index("End_Date")["Volume"].sort_index()
s = s.groupby(level=0).sum().resample(FREQ).sum().fillna(0.0)
return s.loc[s.index > start] if start else slast_sale = weekly.groupby("ISBN")["End_Date"].max()
active = sorted(last_sale[last_sale >= "2024-07-01"].index) # still-selling backlist (used in 2.3)
# Show ALL qualifying ISBNs with title/author (captured for the report), not just the count.
active_tbl = (isbns[isbns["ISBN"].isin(active)][["ISBN", "Title", "Author"]]
.drop_duplicates("ISBN").sort_values("Title").reset_index(drop=True))
active_tbl["Title"] = active_tbl["Title"].str.slice(0, 45) # trim one 150-char outlier for a tidy listing
print(f"{len(active)} ISBNs still selling after 2024-07-01 (all {len(active_tbl)} listed):")
print(active_tbl.to_string(index=False))61 ISBNs still selling after 2024-07-01 (all 61 listed):
ISBN Title Author
9780747268161 After You'd Gone O'Farrell, Maggie
9780722532935 Alchemist,The Coelho, Paulo
9780752844299 Arthur: The Seeing Stone:Book 1:Arthur Crossley-Holland, Kevin
9780752846576 Asterix: Asterix and The Actress:Album 31:Ast Uderzo, Albert
9781841150437 Bad Blood:A Memoir Sage, Lorna
9780552998444 Best a Man Can Get,The O'Farrell, John
9780552145954 Between Extremes Keenan, Brian & McCarthy, John
9780552998000 Blackberry Wine:from Joanne Harris, the bests Harris, Joanne
9780006550433 Bonesetter’s Daughter,The Tan, Amy
9780749397548 Captain Corelli's Mandolin:AS SEEN ON BBC BET de Bernières, Louis
9780552998482 Chocolat:(Chocolat 1) Harris, Joanne
9780140276619 Consolations of Philosophy,The de Botton, Alain
9780552997034 Down Under Bryson, Bill
9780099286578 Elizabeth Starkey, Dr David
9780552997348 Emotionally Weird Atkinson, Kate
9780140285215 English Passengers Kneale, Matthew
9780099285823 Experience Amis, Martin
9780140294231 Extra Virgin Hawes, Annie
9780349112763 Fortune's Rocks Shreve, Anita
9780349114033 Four Blondes Bushnell, Candace
9780099286387 Geisha Dalby, Liza
9780006514091 Glass Palace,The Ghosh, Amitav
9780099428558 Howard Marks' Book Of Dope Stories Marks, Howard
9780006531203 In the Heart of the Sea:The Epic True Story T Philbrick, Nathaniel
9780552145053 Irresistible Forces Steel, Danielle
9780224060875 It's Not About the Bike:My Journey Back to Li Armstrong, Lance
9780440864141 Jacqueline Wilson Double Decker Wilson, Jacqueline
9780440864554 Jacqueline Wilson's Superstars Wilson, Jacqueline
9780552145060 Journey Steel, Danielle
9781841461502 KS2 English Study Book - Ages 7-11:CGP KS2 En CGP Books
9781841462509 KS2 Science Study Book:CGP KS2 Science CGP Books
9781841460406 KS3 Maths Revision Guide – Foundation (includ CGP Books
9781841460307 KS3 Maths Revision Guide – Higher (includes O CGP Books
9781841462400 KS3 Science Revision Guide – Foundation (incl CGP Books
9781841462301 KS3 Science Revision Guide – Higher (includes CGP Books
9780349113609 Last Time They Met,The Shreve, Anita
9780099422587 London:The Biography Ackroyd, Peter
9780261103252 Lord of the Rings,The Tolkien, J. R. R.
9780006512134 Man and Boy Parsons, Tony
9780552998727 Marrying The Mistress:an irresistible and gri Trollope, Joanna
9780340766057 McCarthy's Bar:A Journey of Discovery in Ire McCarthy, Pete
9780099771517 Memoirs of a Geisha:The Literary Sensation an Golden, Arthur
9780006514213 Miss Garnet’s Angel Vickers, Salley
9780006647553 Mog’s Bad Thing Kerr, Judith
9780749395698 Mr Nice Marks, Howard
9780140259506 My East End:Memories of Life in Cockney Londo O'Neill, Gilda
9780340696767 Nathaniel's Nutmeg:How One Man's Courage Chan Milton, Giles
9780719559792 Now We Are Sixty:20th Anniversary Edition Matthew, Christopher
9780340786055 One Hundred Ways for a Cat to Train Its Human Haddon, Celia
9780140295962 One-hit Wonder Jewell, Lisa
9780140275421 Other Side of the Dale,The Phinn, Gervase
9780140281293 Over Hill and Dale Phinn, Gervase
9780091867775 Round Ireland with a Fridge Hawks, Tony
9780007101887 Sky is Falling,The Sheldon, Sidney
9780330355667 Strange Places, Questionable People Simpson, John
9780099244721 Timeline Crichton, Michael
9780593048153 Universe In A Nutshell,The:the beautifully il Hawking, Stephen
9780241003008 Very Hungry Caterpillar,The:The Very Hungry C Carle, Eric
9780744523232 We're Going on a Bear Hunt Rosen, Michael
9780140276336 White Teeth Smith, Zadie
9780091816971 Who Moved My Cheese Johnson, Dr Spencer
plt.rcParams["figure.figsize"] = (12, 4)
ncols = 4; nrows = int(np.ceil(len(active)/ncols))
fig, axes = plt.subplots(nrows, ncols, figsize=(16, 2.2*nrows))
for ax, isbn in zip(axes.ravel(), active):
s = series(isbn)
ax.plot(s.index, s.values, linewidth=0.5); ax.set_title(isbn, fontsize=7)
ax.tick_params(labelsize=6)
for ax in axes.ravel()[len(active):]: ax.axis("off")
fig.tight_layout(); plt.show()
import plotly.express as px
# Per-book market view over the active backlist: speed (copies/week), price (ASP)
# and total revenue (bubble size). Silent weeks are already zero from the resample.
mk = weekly[weekly["ISBN"].isin(active)].copy()
mk["Value"] = pd.to_numeric(mk["Value"].astype(str).str.replace("£", "", regex=False)
.str.replace(",", "", regex=False), errors="coerce").fillna(0)
stats = (mk.groupby("ISBN")
.agg(Volume=("Volume", "sum"), Revenue=("Value", "sum"), Weeks=("End_Date", "nunique"))
.reset_index())
stats["Copies_per_week"] = stats["Volume"] / stats["Weeks"]
stats["ASP"] = stats["Revenue"] / stats["Volume"].replace(0, np.nan)
stats = stats.merge(active_tbl[["ISBN", "Title", "Author"]], on="ISBN", how="left")
fig = px.scatter(stats, x="Copies_per_week", y="ASP", size="Revenue", hover_name="Title",
hover_data={"Author": True, "Revenue": ":,.0f",
"Copies_per_week": ":.0f", "ASP": ":.2f"},
size_max=55, height=560, template="plotly_white",
color_discrete_sequence=["#2a9d8f"],
labels={"Copies_per_week": "average copies sold per week",
"ASP": "average selling price, £", "Revenue": "total revenue, £"},
title="The active backlist as a market: speed, price and total size")
fig.update_traces(marker=dict(opacity=0.7, line=dict(color="white", width=0.7)))
fig.show()
# Top five active titles on each measure.
for label, col in [("copies per week", "Copies_per_week"), ("average price (ASP)", "ASP"),
("total revenue", "Revenue")]:
print(f"\nTop 5 by {label}:")
print(stats.nlargest(5, col)[["Title", "Author", col]].round(2).to_string(index=False))
ALCH, CAT = "9780722532935", "9780241003008"
raw = {"Alchemist": series(ALCH, "2012-01-01"), "Caterpillar": series(CAT, "2012-01-01")}
for n, s in raw.items(): print(f"{n}: {len(s)} weeks, {s.index.min().date()} to {s.index.max().date()}")Alchemist: 655 weeks, 2012-01-07 to 2024-07-20 Caterpillar: 655 weeks, 2012-01-07 to 2024-07-20
M = 52 # weekly data: the annual season is 52 weeks
def stat(s):
out = {}
for label, x in [("raw", s), ("d=1", s.diff().dropna()), ("D=1 (lag52)", s.diff(M).dropna())]:
try: adf_p = round(adfuller(x)[1], 4)
except Exception: adf_p = np.nan
try: kpss_p = round(kpss(x, regression="c", nlags="auto")[1], 4)
except Exception: kpss_p = np.nan
out[label] = {"ADF p": adf_p, "KPSS p": kpss_p}
return pd.DataFrame(out).T
eras = {"pre_covid": ("2012-01-01", "2020-03-01"),
"during_covid": ("2020-03-01", "2021-06-01"),
"post_covid": ("2021-06-01", "2025-01-01")}
for name, s in raw.items():
print(f"\n=== {name}: whole series ===" ); print(stat(s))
for era, (a, b) in eras.items():
seg = s.loc[(s.index >= a) & (s.index < b)]
if len(seg) > M + 5:
print(f"-- {name} {era} (n={len(seg)}) --"); print(stat(seg))
=== Alchemist: whole series ===
ADF p KPSS p
raw 0.0 0.022
d=1 0.0 0.100
D=1 (lag52) 0.0 0.100
-- Alchemist pre_covid (n=426) --
ADF p KPSS p
raw 0.0000 0.0100
d=1 0.0000 0.1000
D=1 (lag52) 0.0096 0.0483
-- Alchemist during_covid (n=65) --
ADF p KPSS p
raw 0.0026 0.1000
d=1 0.0000 0.1000
D=1 (lag52) 0.0000 0.0412
-- Alchemist post_covid (n=164) --
ADF p KPSS p
raw 0.0 0.10
d=1 0.0 0.10
D=1 (lag52) 0.0 0.01
=== Caterpillar: whole series ===
/tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned.
ADF p KPSS p
raw 0.0296 0.01
d=1 0.0000 0.10
D=1 (lag52) 0.0000 0.10
-- Caterpillar pre_covid (n=426) --
ADF p KPSS p
raw 0.2293 0.01
d=1 0.0000 0.10
D=1 (lag52) 0.0145 0.10
-- Caterpillar during_covid (n=65) --
ADF p KPSS p
raw 0.0123 0.1000
d=1 0.0000 0.1000
D=1 (lag52) 0.0105 0.0655
-- Caterpillar post_covid (n=164) --
ADF p KPSS p
raw 0.0 0.0116
d=1 0.0 0.1000
D=1 (lag52) 0.0 0.1000
/tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is smaller than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned. /tmp/ipykernel_2115/2987534856.py:7: InterpolationWarning: The test statistic is outside of the range of p-values available in the look-up table. The actual p-value is greater than the p-value returned.
COVID = [("2020-03-21","2020-06-27"), ("2020-10-03","2021-04-03")]
def in_covid(idx):
m = np.zeros(len(idx), bool)
for a, b in COVID: m |= (idx >= a) & (idx <= b)
return m
def impute_covid(s):
s = s.copy(); woy = s.index.isocalendar().week.values
zero_covid = (s.values == 0) & in_covid(s.index)
for i in np.where(zero_covid)[0]:
past = s.values[(woy == woy[i]) & (s.values > 0) & ~in_covid(s.index)]
if len(past): s.iloc[i] = float(np.mean(past))
return s
books = {n: impute_covid(s) for n, s in raw.items()}
for n in books:
print(f"{n}: imputed {int(((raw[n]==0) & in_covid(raw[n].index)).sum())} COVID-zero weeks; "
f"{int((books[n]<=0).sum())} non-positive weeks remain")Alchemist: imputed 27 COVID-zero weeks; 0 non-positive weeks remain Caterpillar: imputed 27 COVID-zero weeks; 0 non-positive weeks remain
def amp_vs_level_corr(s):
g = s.groupby(s.index.year)
full = g.count() >= 40
amp, lvl = (g.max()-g.min())[full], g.mean()[full]
return float(np.corrcoef(lvl, amp)[0, 1])
for name, s in books.items():
print(f"{name}: amplitude-vs-level correlation = {amp_vs_level_corr(s):.2f} -> multiplicative")
d = seasonal_decompose(s, model="multiplicative", period=M, extrapolate_trend="freq")
d.plot(); plt.gcf().set_size_inches(11, 6); plt.suptitle(name); plt.tight_layout(); plt.show()Alchemist: amplitude-vs-level correlation = 0.85 -> multiplicative
Caterpillar: amplitude-vs-level correlation = 0.67 -> multiplicative
fig, axes = plt.subplots(2, 2, figsize=(14, 6))
for j, (name, s) in enumerate(books.items()):
plot_acf(s, lags=120, ax=axes[0, j], title=f"{name}: ACF")
plot_pacf(s, lags=120, method="ywm", ax=axes[1, j], title=f"{name}: PACF")
plt.tight_layout(); plt.show()
H = 32 # forecast horizon: the final 32 weeks (rubric)
arima, arima_fc = {}, {}
for name, s in books.items():
train, test = s.iloc[:-H], s.iloc[-H:]
m = pm.auto_arima(train, seasonal=True, m=M, max_p=2, max_q=2, max_d=1,
max_P=1, max_Q=1, max_D=1, information_criterion="aicc",
stepwise=True, suppress_warnings=True, error_action="ignore")
fc = pd.Series(m.predict(H), index=test.index)
arima[name], arima_fc[name] = m, fc
print(f"{name}: {m.order} x {m.seasonal_order} AIC={m.aic():.0f}")Alchemist: (0, 1, 2) x (1, 0, 1, 52) AIC=7424 Caterpillar: (2, 1, 1) x (1, 0, 1, 52) AIC=8683
fig, axes = plt.subplots(1, 2, figsize=(14, 3))
for ax, (name, m) in zip(axes, arima.items()):
pd.Series(m.resid()).plot(ax=ax, linewidth=0.6, title=f"{name}: SARIMA residuals")
plt.tight_layout(); plt.show()
fig, axes = plt.subplots(1, 2, figsize=(15, 4))
for ax, (name, s) in zip(axes, books.items()):
train, test = s.iloc[:-H], s.iloc[-H:]
mean, ci = arima[name].predict(H, return_conf_int=True)
ax.plot(train.index[-104:], train.values[-104:], linewidth=0.7, label="train")
ax.plot(test.index, test.values, label="actual", linewidth=1.2)
ax.plot(test.index, mean, label="forecast", linewidth=1.2)
ax.fill_between(test.index, ci[:, 0], ci[:, 1], alpha=0.2)
ax.set_title(name); ax.legend(loc="upper left", fontsize=8)
plt.tight_layout(); plt.show()
def scores(y, f):
y, f = np.asarray(y, float), np.asarray(f, float)
nz = y != 0
return {"MAE": float(np.mean(np.abs(y-f))),
"MAPE": float(np.mean(np.abs((y[nz]-f[nz])/y[nz]))*100),
"RMSE": float(np.sqrt(np.mean((y-f)**2)))}
results = [] # one row per model: book, model, MAE, MAPE, RMSE
for name, s in books.items():
sc = scores(s.iloc[-H:], arima_fc[name])
results.append({"book": name, "model": "SARIMA", **sc})
print(f"{name} SARIMA MAE={sc['MAE']:.1f} MAPE={sc['MAPE']:.1f}%")Alchemist SARIMA MAE=127.8 MAPE=21.0% Caterpillar SARIMA MAE=345.4 MAPE=17.3%
SEED = 42; np.random.seed(SEED) # reproducible from the first ML model on
def lags(s, w):
f = pd.concat([s.shift(i).rename(f"lag{i}") for i in range(1, w+1)], axis=1)
f["y"] = s.values
return f.dropna()
grid = {"n_estimators": [200, 400], "max_depth": [3, 5], "learning_rate": [0.05, 0.1]}
xgb_fc = {}
for name, s in books.items():
train, test = s.iloc[:-H], s.iloc[-H:]
best = None
for w in [12, 26, 52]: # window_length is part of the search
fr = lags(train, w)
gs = GridSearchCV(XGBRegressor(random_state=SEED, n_jobs=-1), grid,
cv=TimeSeriesSplit(3), scoring="neg_mean_absolute_error")
gs.fit(fr.drop(columns="y"), fr["y"])
if best is None or gs.best_score_ > best[0]:
best = (gs.best_score_, w, gs.best_estimator_)
_, w, model = best
hist = list(s.iloc[:-H].values) # recursive multi-step forecast
for _ in range(H): # features are lag1..lagw = most-recent first
hist.append(float(model.predict(np.array(hist[-w:][::-1]).reshape(1, -1))[0]))
fc = pd.Series(np.clip(hist[-H:], 0, None), index=test.index)
xgb_fc[name] = fc
sc = scores(test, fc); results.append({"book": name, "model": "XGBoost", **sc})
print(f"{name} XGBoost (window={w}) MAE={sc['MAE']:.1f} MAPE={sc['MAPE']:.1f}%")
fig, axes = plt.subplots(1, 2, figsize=(15, 4))
for ax, (name, s) in zip(axes, books.items()):
ax.plot(s.index[-104:], s.values[-104:], linewidth=0.7, label="actual")
ax.plot(xgb_fc[name].index, xgb_fc[name].values, label="XGBoost"); ax.set_title(name); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()Alchemist XGBoost (window=52) MAE=134.1 MAPE=27.6% Caterpillar XGBoost (window=52) MAE=368.1 MAPE=18.5%
tf.random.set_seed(SEED) # seed TensorFlow before building the LSTM
WIN = 26
def seqs(v, w): return np.stack([v[i-w:i] for i in range(w, len(v))]), v[w:]
lstm_fc = {}
for name, s in books.items():
train, test = s.iloc[:-H], s.iloc[-H:]
sc_ = MinMaxScaler(); z = sc_.fit_transform(train.values.reshape(-1, 1)).ravel()
X, y = seqs(z, WIN); X = X[..., None]
Xf, yf, Xv, yv = X[:-H], y[:-H], X[-H:], y[-H:]
def build(hp):
mdl = keras.Sequential([layers.Input((WIN, 1)),
layers.LSTM(hp.Int("units", 16, 96, step=16), dropout=hp.Float("dropout", 0.0, 0.3, step=0.1)),
layers.Dense(1)])
mdl.compile(optimizer=keras.optimizers.Adam(hp.Choice("lr", [1e-3, 5e-4])), loss="mse")
return mdl
tuner = kt.RandomSearch(build, objective="val_loss", max_trials=5, overwrite=True,
seed=SEED, directory="kt", project_name=f"lstm_{name}")
stop = keras.callbacks.EarlyStopping(monitor="val_loss", patience=5, restore_best_weights=True)
tuner.search(Xf, yf, validation_data=(Xv, yv), epochs=30, batch_size=32, callbacks=[stop], verbose=0)
model = tuner.get_best_models(1)[0]
model.fit(X, y, validation_split=0.1, epochs=30, batch_size=32, callbacks=[stop], verbose=0)
hist = list(z)
for _ in range(H):
hist.append(float(model.predict(np.array(hist[-WIN:]).reshape(1, WIN, 1), verbose=0)[0, 0]))
fc = pd.Series(np.clip(sc_.inverse_transform(np.array(hist[-H:]).reshape(-1, 1)).ravel(), 0, None), index=test.index)
lstm_fc[name] = fc
sc = scores(test, fc); results.append({"book": name, "model": "LSTM", **sc})
print(f"{name} LSTM MAE={sc['MAE']:.1f} MAPE={sc['MAPE']:.1f}%")
fig, axes = plt.subplots(1, 2, figsize=(15, 4))
for ax, (name, s) in zip(axes, books.items()):
ax.plot(s.index[-104:], s.values[-104:], linewidth=0.7, label="actual")
ax.plot(lstm_fc[name].index, lstm_fc[name].values, label="LSTM"); ax.set_title(name); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()Alchemist LSTM MAE=207.8 MAPE=26.8% Caterpillar LSTM MAE=901.5 MAPE=38.6%
hyb_seq = {}
for name, s in books.items():
test = s.iloc[-H:]
resid = pd.Series(arima[name].resid(), index=s.iloc[:-H].index).dropna()
sc_ = MinMaxScaler(); z = sc_.fit_transform(resid.values.reshape(-1, 1)).ravel()
X, y = seqs(z, WIN); X = X[..., None]
def build(hp):
mdl = keras.Sequential([layers.Input((WIN, 1)),
layers.LSTM(hp.Int("units", 16, 64, step=16), dropout=hp.Float("dropout", 0.0, 0.3, step=0.1)),
layers.Dense(1)]); mdl.compile(optimizer="adam", loss="mse"); return mdl
tuner = kt.RandomSearch(build, objective="val_loss", max_trials=5, overwrite=True,
seed=SEED, directory="kt", project_name=f"hyb_{name}")
stop = keras.callbacks.EarlyStopping(monitor="val_loss", patience=5, restore_best_weights=True)
tuner.search(X[:-H], y[:-H], validation_data=(X[-H:], y[-H:]), epochs=30, batch_size=32, callbacks=[stop], verbose=0)
model = tuner.get_best_models(1)[0]; model.fit(X, y, epochs=30, batch_size=32, callbacks=[stop], verbose=0)
hist = list(z)
for _ in range(H):
hist.append(float(model.predict(np.array(hist[-WIN:]).reshape(1, WIN, 1), verbose=0)[0, 0]))
resid_fc = sc_.inverse_transform(np.array(hist[-H:]).reshape(-1, 1)).ravel()
fc = pd.Series(np.clip(arima_fc[name].values + resid_fc, 0, None), index=test.index)
hyb_seq[name] = fc
sc = scores(test, fc); results.append({"book": name, "model": "Hybrid-sequential", **sc})
print(f"{name} Hybrid-seq MAE={sc['MAE']:.1f} MAPE={sc['MAPE']:.1f}%")Alchemist Hybrid-seq MAE=127.7 MAPE=21.4% Caterpillar Hybrid-seq MAE=339.5 MAPE=17.7%
hyb_par = {}
fig, axes = plt.subplots(1, 2, figsize=(14, 3.5))
for ax, (name, s) in zip(axes, books.items()):
test = s.iloc[-H:]
sweep = [(w, scores(test, w*arima_fc[name].values + (1-w)*lstm_fc[name].values)["MAE"]) for w in np.linspace(0, 1, 21)]
ax.plot(*zip(*sweep), marker="o"); ax.set_title(f"{name}: weight sweep"); ax.set_xlabel("SARIMA weight"); ax.set_ylabel("MAE")
fc = pd.Series((0.5*arima_fc[name].values + 0.5*lstm_fc[name].values).clip(0), index=test.index) # submitted = fixed 50/50
hyb_par[name] = fc
sc = scores(test, fc); results.append({"book": name, "model": "Hybrid-parallel(0.5)", **sc})
print(f"{name} Hybrid-par 50/50 MAE={sc['MAE']:.1f} MAPE={sc['MAPE']:.1f}% (best sweep weight={min(sweep, key=lambda t: t[1])[0]:.2f})")
plt.tight_layout(); plt.show()Alchemist Hybrid-par 50/50 MAE=164.7 MAPE=23.4% (best sweep weight=1.00) Caterpillar Hybrid-par 50/50 MAE=548.8 MAPE=24.0% (best sweep weight=1.00)
HM = 8
fig, axes = plt.subplots(1, 2, figsize=(15, 4))
for ax, (name, s) in zip(axes, books.items()):
m = s.resample("MS").sum(); tr, te = m.iloc[:-HM], m.iloc[-HM:]
sar = pd.Series(pm.auto_arima(tr, seasonal=True, m=12, max_p=2, max_q=2, max_d=1, max_P=1, max_Q=1, max_D=1,
suppress_warnings=True, error_action="ignore").predict(HM), index=te.index)
fr = lags(tr, 12); xg = XGBRegressor(n_estimators=300, max_depth=3, learning_rate=0.05, random_state=SEED).fit(fr.drop(columns="y"), fr["y"])
hist = list(tr.values)
for _ in range(HM): hist.append(float(xg.predict(np.array(hist[-12:][::-1]).reshape(1, -1))[0]))
xgm = pd.Series(np.clip(hist[-HM:], 0, None), index=te.index)
for label, fc in [("SARIMA", sar), ("XGBoost", xgm)]:
sc = scores(te, fc); print(f"{name} monthly {label} MAE={sc['MAE']:.0f} MAPE={sc['MAPE']:.1f}%")
ax.plot(te.index, te.values, label="actual"); ax.plot(te.index, sar, label="SARIMA"); ax.plot(te.index, xgm, label="XGBoost")
ax.set_title(f"{name}: monthly (8-mo)"); ax.legend(fontsize=8)
plt.tight_layout(); plt.show()Alchemist monthly SARIMA MAE=777 MAPE=37.6% Alchemist monthly XGBoost MAE=693 MAPE=36.5% Caterpillar monthly SARIMA MAE=1948 MAPE=20.7% Caterpillar monthly XGBoost MAE=2903 MAPE=29.1%
summary = pd.DataFrame(results).sort_values(["book", "RMSE"]).round(1)
summary| book | model | MAE | MAPE | RMSE | |
|---|---|---|---|---|---|
| 2 | Alchemist | XGBoost | 134.1 | 27.6 | 188.4 |
| 6 | Alchemist | Hybrid-sequential | 127.7 | 21.4 | 192.1 |
| 0 | Alchemist | SARIMA | 127.8 | 21.0 | 193.4 |
| 8 | Alchemist | Hybrid-parallel(0.5) | 164.7 | 23.4 | 265.1 |
| 4 | Alchemist | LSTM | 207.8 | 26.8 | 352.1 |
| 7 | Caterpillar | Hybrid-sequential | 339.5 | 17.7 | 417.6 |
| 1 | Caterpillar | SARIMA | 345.4 | 17.3 | 424.3 |
| 3 | Caterpillar | XGBoost | 368.1 | 18.5 | 480.1 |
| 9 | Caterpillar | Hybrid-parallel(0.5) | 548.8 | 24.0 | 662.7 |
| 5 | Caterpillar | LSTM | 901.5 | 38.6 | 1022.1 |
val = {ALCH: "Alchemist", CAT: "Caterpillar"}
weekly["Value"] = pd.to_numeric(weekly["Value"].astype(str).str.replace("£","",regex=False).str.replace(",","",regex=False), errors="coerce")
rows = []
for isbn, name in val.items():
s = books[name]; test = s.iloc[-H:]
v = weekly.loc[weekly["ISBN"]==isbn].set_index("End_Date")["Value"].resample(FREQ).sum().reindex(test.index).fillna(0).values
price = np.divide(v, test.values, out=np.zeros_like(v), where=test.values>0)
err = test.values - arima_fc[name].values # +ve = under-forecast (units)
festive = test.index.month.isin([12]) | ((test.index.month==1) & (test.index.day<=15))
rows.append({"book": name, "peak_week_£_missed": round(float(err[np.argmax(test.values)]*price[np.argmax(test.values)])),
"festive_£_under_ordered": round(float(np.clip(err[festive],0,None) @ price[festive])),
"holdout_£_under_ordered": round(float(np.clip(err,0,None) @ price))})
pd.DataFrame(rows)| book | peak_week_£_missed | festive_£_under_ordered | holdout_£_under_ordered | |
|---|---|---|---|---|
| 0 | Alchemist | 7472 | 12622 | 22903 |
| 1 | Caterpillar | 2620 | 2620 | 30137 |