PROJECT 002 · PYTHON SOURCE
Next-day bike demand forecasting
Study-specific code. Shared modules, dependency versions, and reproduction instructions are included in all project files.
Download Python file ↓"""Next-day hourly demand forecast with calendar-time lags and rolling daily origins."""
from __future__ import annotations
import json
from pathlib import Path
import joblib
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.ensemble import HistGradientBoostingRegressor
from sklearn.impute import SimpleImputer
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.pipeline import make_pipeline
PROJECT = Path(__file__).resolve().parents[1]
OUTPUT = PROJECT / "outputs"
FIGURES = PROJECT / "reports" / "figures"
TEST_START = pd.Timestamp("2012-09-01")
FEATURES = ["hr", "weekday", "workingday", "holiday", "month_sin", "month_cos", "hour_sin", "hour_cos", "trend_days", "lag_24h", "lag_48h", "lag_168h", "lag_336h", "past_7_same_hour_mean", "previous_day_mean"]
FOLDS = [("2012-03-01", "2012-05-01"), ("2012-05-01", "2012-07-01"), ("2012-07-01", "2012-09-01")]
CANDIDATES = [
{"max_leaf_nodes": 15, "l2_regularization": 1.0},
{"max_leaf_nodes": 31, "l2_regularization": 10.0},
]
def load_data() -> pd.DataFrame:
path = PROJECT / "data" / "raw" / "hour.csv"
if not path.exists():
raise FileNotFoundError("Run python scripts/download_data.py bikes from the repository root.")
return pd.read_csv(path, parse_dates=["dteday"])
def feature_table(raw: pd.DataFrame) -> tuple[pd.DataFrame, dict]:
source = raw.copy()
source["timestamp"] = pd.to_datetime(source["dteday"]) + pd.to_timedelta(source["hr"], unit="h")
if source["timestamp"].duplicated().any():
raise ValueError("Duplicate hourly timestamps")
if source.isna().any().any():
raise ValueError("Missing source values; review data policy before modeling")
if not (source["cnt"] == source["casual"] + source["registered"]).all():
raise ValueError("Counts fail source reconciliation")
if (source["cnt"] < 0).any():
raise ValueError("Negative demand")
source = source.set_index("timestamp").sort_index()
timeline = pd.date_range(source.index.min(), source.index.max(), freq="h")
counts = source["cnt"].reindex(timeline)
# Shift on a complete time grid. A row shift on the source would jump across missing hours.
frame = pd.DataFrame(index=timeline)
frame["cnt"] = counts
frame["hr"] = timeline.hour
frame["weekday"] = timeline.dayofweek
frame["workingday"] = source["workingday"].reindex(timeline)
frame["holiday"] = source["holiday"].reindex(timeline)
frame["month_sin"] = np.sin(2 * np.pi * timeline.month / 12)
frame["month_cos"] = np.cos(2 * np.pi * timeline.month / 12)
frame["hour_sin"] = np.sin(2 * np.pi * timeline.hour / 24)
frame["hour_cos"] = np.cos(2 * np.pi * timeline.hour / 24)
frame["trend_days"] = (timeline - pd.Timestamp("2011-01-01")).total_seconds() / 86400
for lag in [24, 48, 168, 336]:
frame[f"lag_{lag}h"] = counts.shift(lag)
frame["past_7_same_hour_mean"] = pd.concat([counts.shift(24 * d) for d in range(1, 8)], axis=1).mean(axis=1)
previous_day = counts.resample("D").mean().shift(1)
frame["previous_day_mean"] = frame.index.normalize().map(previous_day)
audit = {
"source_hourly_rows": len(source), "calendar_hours": len(timeline),
"missing_calendar_hours": int(counts.isna().sum()),
"source_start": str(source.index.min()), "source_end": str(source.index.max()),
"total_rentals": int(source["cnt"].sum()),
"rows_with_missing_lag_24h": int(frame.loc[frame["cnt"].notna(), "lag_24h"].isna().sum()),
}
# No invented targets: evaluate observed hours only; allow missing historical features.
frame = frame.loc[frame["cnt"].notna() & (frame.index >= timeline.min() + pd.Timedelta(days=14))].copy()
audit["modeling_rows_after_warmup"] = len(frame)
if {"cnt", "casual", "registered", "temp", "atemp", "hum", "windspeed", "weathersit"} & set(FEATURES):
raise AssertionError("Target leakage or future observed weather in features")
return frame, audit
def model_for(params: dict):
# Imputation learns only the training distribution, including in every validation fold.
return make_pipeline(
SimpleImputer(strategy="median", add_indicator=True),
HistGradientBoostingRegressor(loss="poisson", learning_rate=.07, max_iter=180,
early_stopping=False, random_state=42, **params),
)
def metrics(y, prediction) -> dict:
y, prediction = np.asarray(y), np.asarray(prediction)
return {
"mae": float(mean_absolute_error(y, prediction)),
"rmse": float(np.sqrt(mean_squared_error(y, prediction))),
"wape": float(np.abs(y - prediction).sum() / y.sum()),
"bias": float((prediction - y).mean()),
}
def baseline_predictions(train: pd.DataFrame, target: pd.DataFrame) -> dict:
calendar = train.groupby(["weekday", "hr"])["cnt"].mean()
fallback = np.array([calendar.get((row.weekday, row.hr), train["cnt"].mean()) for row in target.itertuples()])
return {
"Calendar average": fallback,
"Previous day": target["lag_24h"].fillna(pd.Series(fallback, index=target.index)).to_numpy(),
"Previous week": target["lag_168h"].fillna(pd.Series(fallback, index=target.index)).to_numpy(),
}
def evaluate(frame: pd.DataFrame, audit: dict) -> dict:
records = []
baseline_records = []
for start, end in FOLDS:
train = frame[frame.index < start]
validation = frame[(frame.index >= start) & (frame.index < end)]
assert train.index.max() < validation.index.min() < TEST_START
for name, pred in baseline_predictions(train, validation).items():
baseline_records.append({"fold_start": start, "model": name, "rows": len(validation), **metrics(validation["cnt"], pred)})
for candidate, params in enumerate(CANDIDATES):
model = model_for(params)
model.fit(train[FEATURES], train["cnt"])
pred = model.predict(validation[FEATURES])
records.append({"fold_start": start, "candidate": candidate, "rows": len(validation), **params, **metrics(validation["cnt"], pred)})
validation_scores = pd.DataFrame(records)
scores = validation_scores.assign(weighted_error=lambda d: d["mae"] * d["rows"]).groupby("candidate").agg(error=("weighted_error", "sum"), rows=("rows", "sum"))
selected = int((scores["error"] / scores["rows"]).idxmin())
baseline_validation = pd.DataFrame(baseline_records)
baseline_scores = baseline_validation.assign(weighted_error=lambda d: d["mae"] * d["rows"]).groupby("model").agg(error=("weighted_error", "sum"), rows=("rows", "sum"))
selected_baseline = str((baseline_scores["error"] / baseline_scores["rows"]).idxmin())
train = frame[frame.index < TEST_START]
test = frame[frame.index >= TEST_START]
model = model_for(CANDIDATES[selected])
model.fit(train[FEATURES], train["cnt"])
predictions = baseline_predictions(train, test)
predictions["Gradient boosting"] = model.predict(test[FEATURES])
test_scores = pd.DataFrame([{"model": name, "rows": len(test), **metrics(test["cnt"], pred)} for name, pred in predictions.items()])
pred = pd.DataFrame({"actual": test["cnt"], **predictions}, index=test.index)
pred.index.name = "timestamp"
selected_mae = test_scores.set_index("model").loc["Gradient boosting", "mae"]
baseline_mae = test_scores.set_index("model").loc[selected_baseline, "mae"]
summary = {
**audit, "train_rows": len(train), "test_rows": len(test),
"train_end": str(train.index.max()), "test_start": str(test.index.min()),
"selected_candidate": selected, "selected_parameters": CANDIDATES[selected],
"baseline_selected_on_validation": selected_baseline,
"test_mae_improvement_vs_selected_baseline": float(1 - selected_mae / baseline_mae),
"features": FEATURES,
"forecast_origin": "00:00 each target day; assumes prior-day counts finalized by midnight",
"evaluation": "Rolling next-day prediction with fixed model; realized past test counts become available on subsequent days",
}
# Exploratory uncertainty: resample whole days so the 24 within-day errors stay together.
errors = pd.DataFrame({"model": np.abs(pred["Gradient boosting"] - pred["actual"]),
"baseline": np.abs(pred[selected_baseline] - pred["actual"])}).resample("D").agg(["sum", "count"])
rng = np.random.default_rng(42)
deltas = []
for _ in range(1000):
sample = errors.iloc[rng.integers(0, len(errors), len(errors))]
deltas.append(float((sample["baseline"]["sum"].sum() - sample["model"]["sum"].sum()) / sample["model"]["count"].sum()))
summary["daily_bootstrap_mae_reduction_95_interval"] = np.quantile(deltas, [.025, .975]).tolist()
summary["bootstrap_limitation"] = "Exploratory interval; day resampling does not preserve dependence across successive days."
diagnostic = []
for month, group in pred.groupby(pred.index.strftime("%Y-%m")):
diagnostic.append({"slice": month, "rows": len(group), **metrics(group["actual"], group["Gradient boosting"])})
for name, mask in {"Commute hours (all days)": pred.index.hour.isin([7, 8, 9, 16, 17, 18]), "Other hours": ~pred.index.hour.isin([7, 8, 9, 16, 17, 18])}.items():
group = pred[mask]
diagnostic.append({"slice": name, "rows": len(group), **metrics(group["actual"], group["Gradient boosting"])})
return {"summary": summary, "validation": validation_scores, "baseline_validation": baseline_validation,
"test_scores": test_scores, "predictions": pred, "diagnostics": pd.DataFrame(diagnostic), "model": model}
def make_figures(raw: pd.DataFrame, result: dict) -> None:
FIGURES.mkdir(parents=True, exist_ok=True)
plt.rcParams.update({"font.size": 11, "axes.spines.top": False, "axes.spines.right": False, "axes.titleweight": "bold"})
fig, ax = plt.subplots(figsize=(10, 4.5), layout="constrained")
for working, label, color in [(1, "Working days", "#164e99"), (0, "Non-working days", "#b94e18")]:
# Historical exploration is restricted to pre-test data.
profile = raw[(raw["workingday"] == working) & (raw["dteday"] < TEST_START)].groupby("hr")["cnt"].mean()
ax.plot(profile.index, profile, marker="o", color=color, label=label)
ax.set(title="Demand follows different daily rhythms", xlabel="Hour of day", ylabel="Mean observed hourly rentals", xticks=range(0, 24, 3))
ax.legend(frameon=False)
fig.savefig(FIGURES / "hourly_demand.png", dpi=160)
plt.close(fig)
scores = result["test_scores"].sort_values("mae")
fig, ax = plt.subplots(figsize=(10, 4.5), layout="constrained")
bars = ax.barh(scores["model"], scores["mae"], color=["#164e99" if x == "Gradient boosting" else "#9caabd" for x in scores["model"]])
ax.bar_label(bars, fmt="%.1f", padding=5)
ax.set_xlim(0, scores["mae"].max() * 1.15)
ax.set(title="Next-day forecast error on unseen dates", xlabel="Mean absolute error · rentals per observed hour · Sep–Dec 2012")
ax.invert_yaxis()
fig.savefig(FIGURES / "model_comparison.png", dpi=160)
plt.close(fig)
week = result["predictions"].loc["2012-09-03":"2012-09-09"]
fig, ax = plt.subplots(figsize=(12, 4.5), layout="constrained")
ax.plot(week.index, week["actual"], color="#1a2638", label="Observed", linewidth=1.6)
ax.plot(week.index, week["Gradient boosting"], color="#2574cf", label="Predicted", linewidth=1.5)
ax.set(title="First full Monday–Sunday week in the test period", ylabel="Hourly rentals", xlabel="Predictions issued separately at midnight each day")
ax.legend(frameon=False)
fig.savefig(FIGURES / "forecast_week.png", dpi=160)
plt.close(fig)
def write_outputs(frame: pd.DataFrame, result: dict) -> None:
OUTPUT.mkdir(parents=True, exist_ok=True)
for key in ["validation", "baseline_validation", "test_scores", "diagnostics"]:
result[key].to_csv(OUTPUT / f"{key}.csv", index=False)
result["predictions"].to_csv(OUTPUT / "test_predictions.csv")
(OUTPUT / "metrics.json").write_text(json.dumps(result["summary"], indent=2) + "\n", encoding="utf-8")
joblib.dump({"model": result["model"], "features": FEATURES, "training_end": result["summary"]["train_end"]}, OUTPUT / "demand_model.joblib")
def main() -> dict:
print("Bikes: creating historical features and evaluating chronological splits", flush=True)
raw = load_data()
frame, audit = feature_table(raw)
# Independently reconcile the provided daily and hourly source files.
day = pd.read_csv(PROJECT / "data" / "raw" / "day.csv", parse_dates=["dteday"])
daily = raw.groupby("dteday")["cnt"].sum().sort_index()
np.testing.assert_array_equal(daily.to_numpy(), day.sort_values("dteday")["cnt"].to_numpy())
result = evaluate(frame, audit)
write_outputs(frame, result)
make_figures(raw, result)
print(result["test_scores"].to_string(index=False), flush=True)
print(json.dumps(result["summary"], indent=2), flush=True)
return result
if __name__ == "__main__":
main()