!pip install -q chronos-forecasting statsforecast timesfm tirex-ts
1 Template: your own forecasting experiment
School of Digital Business 5.0 — Time Series in the AI Era.
The same toolkit as the worked example, for your series. Change only the settings cell, run everything (Runtime → Run all), then write your conclusions in the last section.
Public or anonymised data only. Everything runs inside your Colab session.
1.1 0. Setup
import os, time, warnings
os.environ["NIXTLA_ID_AS_COL"] = "true"
os.environ["HF_HUB_DISABLE_PROGRESS_BARS"] = "1"
warnings.filterwarnings("ignore")
import transformers
transformers.utils.logging.disable_progress_bar()
transformers.utils.logging.set_verbosity_error()
import numpy as np, pandas as pd, torch, matplotlib.pyplot as plt
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS, AutoTheta, MSTL, Naive, SeasonalNaive
plt.rcParams.update({"figure.figsize": (12, 3.6), "axes.grid": True, "grid.alpha": 0.3})
PINK, NAVY, TEAL, ORANGE, SLATE = "#e64173", "#181485", "#20B2AA", "#FFA500", "#314f4f"
CLASSICAL = {"Naive": lambda s: Naive(alias="Naive"),
"Seasonal naive": lambda s: SeasonalNaive(s, alias="Seasonal naive"),
"ETS": lambda s: AutoETS(season_length=s, alias="ETS"),
"ARIMA": lambda s: AutoARIMA(season_length=s, alias="ARIMA"),
"Theta": lambda s: AutoTheta(season_length=s, alias="Theta"),
"MSTL": lambda s: MSTL(season_length=[s, 365], alias="MSTL")}
_loaded = {}
def foundation(name):
# Load each foundation model once; the first call downloads its weights.
if name not in _loaded:
if name == "Chronos-2":
from chronos import BaseChronosPipeline
_loaded[name] = BaseChronosPipeline.from_pretrained("amazon/chronos-2", device_map="cpu")
elif name == "TimesFM":
import timesfm
m = timesfm.TimesFM_2p5_200M_torch.from_pretrained("google/timesfm-2.5-200m-pytorch")
m.compile(timesfm.ForecastConfig(max_context=1024, max_horizon=128, normalize_inputs=True,
use_continuous_quantile_head=True))
_loaded[name] = m
elif name == "TiRex":
import tirex
_loaded[name] = tirex.load_model("NX-AI/TiRex", device="cpu")
return _loaded[name]
def forecast_fm(name, history, h):
# Median and 80% interval (P10, P90) from a foundation model, for one series.
ctx = torch.tensor(np.asarray(history, dtype=np.float32))
model = foundation(name)
if name == "Chronos-2":
q, _ = model.predict_quantiles([ctx], prediction_length=h, quantile_levels=[0.1, 0.5, 0.9])
q = np.asarray(q[0]).reshape(h, 3)
return q[:, 1], q[:, 0], q[:, 2]
if name == "TimesFM":
point, q = model.forecast(horizon=h, inputs=[ctx.numpy()])
q = np.asarray(q)[0]
return np.asarray(point)[0], q[:, 1], q[:, 9]
q, mean = model.forecast(context=[ctx], prediction_length=h)
q = np.asarray(q)[0]
return np.asarray(mean)[0], q[:, 0], q[:, 8]
def backtest(y, h, windows, season, classical, foundations, freq="D"):
# Rolling-origin backtest: `windows` cuts of `h` steps at the end of y. One row per date and cut.
df = pd.DataFrame({"unique_id": "series", "ds": y.index, "y": y.values})
timing = {}
start = time.time()
models = [CLASSICAL[name](season) for name in classical]
res = StatsForecast(models=models, freq=freq).cross_validation(df=df, h=h, n_windows=windows,
step_size=h, level=[80])
for name in classical:
timing[name] = (time.time() - start) / len(classical)
for name in foundations:
start, parts = time.time(), []
for cut in sorted(res["cutoff"].unique()):
mid, lo, hi = forecast_fm(name, y[y.index <= cut].values, h)
parts.append(pd.DataFrame({"ds": res.loc[res["cutoff"] == cut, "ds"].values, "cutoff": cut,
name: mid, f"{name}-lo-80": lo, f"{name}-hi-80": hi}))
res = res.merge(pd.concat(parts), on=["ds", "cutoff"])
timing[name] = time.time() - start
res["h"] = res.groupby("cutoff").cumcount() + 1
return res, timing
def score(res, y, season, names, timing=None):
# MASE (scaled by the seasonal-naive error in the year before the test), MAPE, bias, 80% coverage.
past = y[y.index <= res["cutoff"].min()].values[-365:]
scale = np.mean(np.abs(past[season:] - past[:-season]))
rows = {}
for name in names:
e = res["y"] - res[name]
row = {"MASE": np.mean(np.abs(e)) / scale, "MAPE, %": 100 * np.mean(np.abs(e) / res["y"]),
"bias, %": 100 * np.mean(-e / res["y"])}
if f"{name}-lo-80" in res:
row["inside P10-P90"] = np.mean((res["y"] >= res[f"{name}-lo-80"]) & (res["y"] <= res[f"{name}-hi-80"]))
if timing and name in timing:
row["seconds"] = timing[name]
rows[name] = row
return pd.DataFrame(rows).T, scale1.2 1. Settings
PRODUCT— one of the lecture’s ten products, orNoneto upload your own CSV (two columns: date, value; one row per day or per week).SEASON— 7 for daily data with a weekly rhythm, 52 for weekly data.H— how far ahead;WINDOWS— how many backtest cuts.CLASSICAL,FOUNDATION— which forecasters to compare. ARIMA and MSTL are slow on long daily series.
PRODUCT = "Duolingo" # ChatGPT, Duolingo, Spotify, Netflix, Tinder, Strava, TikTok, Notion, Threads, Black Friday, or None
SEASON = 7 # 7 for daily data, 52 for weekly data
H = 28 # forecast horizon, in steps
WINDOWS = 6 # backtest windows
CLASSICAL_MODELS = ["Seasonal naive", "ETS", "Theta"] # also: "Naive", "ARIMA", "MSTL" (slow)
FOUNDATION = ["Chronos-2", "TiRex"] # also: "TimesFM"1.3 2. Look at your data
if PRODUCT:
data = pd.read_csv("https://raw.githubusercontent.com/Aranaur/aranaur.rbind.io/main/lectures/genesis/sdb/data/pageviews.csv", parse_dates=["ds"])
y = data[data["unique_id"] == PRODUCT].set_index("ds")["y"].dropna()
TITLE = PRODUCT
else:
from google.colab import files
uploaded = files.upload()
raw = pd.read_csv(next(iter(uploaded)))
raw.columns = ["ds", "y"] + list(raw.columns[2:])
y = raw.assign(ds=pd.to_datetime(raw["ds"])).set_index("ds")["y"].sort_index().dropna()
TITLE = "Uploaded"
FREQ = pd.infer_freq(y.index[-30:]) or "D"
y = y.asfreq(FREQ).interpolate()
print(len(y), y.index.min().date(), y.index.max().date(), "frequency:", FREQ)
fig, ax = plt.subplots()
ax.plot(y.index, y.values, color="lightgray", lw=0.8, label="Daily views")
ax.plot(y.index, y.rolling(7).mean(), color=SLATE, lw=1.5, label="7-day average")
ax.legend(loc="upper left"); ax.set_title(TITLE); plt.show()
recent = y[y.index >= y.index.max() - pd.Timedelta(days=3 * 365)]
rel = recent / recent.groupby(recent.index.year).transform("mean") - 1
fig, axes = plt.subplots(1, 2, figsize=(12, 3.4))
for year, g in rel.groupby(rel.index.year):
axes[0].plot(range(7), g.groupby(g.index.dayofweek).mean().reindex(range(7)).values, "o-", label=str(year))
axes[1].plot(range(1, 13), g.groupby(g.index.month).mean().reindex(range(1, 13)).values, "o-", label=str(year))
axes[0].set_xticks(range(7), ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']); axes[0].set_title("Average day of week, by year (% of the year's mean)", fontsize=10)
axes[1].set_xticks(range(1, 13)); axes[1].set_title("Average by month, by year (% of the year's mean)", fontsize=10)
for ax in axes:
ax.yaxis.set_major_formatter(plt.matplotlib.ticker.PercentFormatter(1.0)); ax.legend(fontsize=8)
plt.show()1.4 3. Backtest
res, timing = backtest(y, H, WINDOWS, SEASON, CLASSICAL_MODELS, FOUNDATION, freq=FREQ)
NAMES = CLASSICAL_MODELS + FOUNDATION
print("cut-offs:", [f"{c:%Y-%m-%d}" for c in sorted(res["cutoff"].unique())])1.5 4. Scoreboard
board, scale = score(res, y, SEASON, NAMES, timing)
board.round(3).style.background_gradient(subset=["MASE"], cmap="RdYlGn_r").format(precision=3)1.6 5. By horizon and by window
SHOW = NAMES
by_h = res.groupby((res["h"] - 1) // 7).apply(
lambda g: pd.Series({m: np.mean(np.abs(g["y"] - g[m])) / scale for m in SHOW}), include_groups=False)
by_w = res.groupby("cutoff").apply(
lambda g: pd.Series({m: np.mean(np.abs(g["y"] - g[m])) / scale for m in NAMES}), include_groups=False)
fig, axes = plt.subplots(1, 2, figsize=(13, 3.8), width_ratios=[1, 1.4])
for m in SHOW:
axes[0].plot(by_h.index + 1, by_h[m], "o-", label=m)
axes[0].set_xticks(by_h.index + 1); axes[0].set_xlabel("week"); axes[0].set_title("MASE by week of the horizon", fontsize=10)
axes[0].legend(fontsize=8)
im = axes[1].imshow(by_w.T.values, cmap="RdYlGn_r", vmin=0.2, vmax=1.6, aspect="auto")
axes[1].set_yticks(range(len(NAMES)), NAMES); axes[1].set_xticks(range(len(by_w)), [f"{c:%d.%m}" for c in by_w.index])
for i in range(len(NAMES)):
for j in range(len(by_w)):
axes[1].text(j, i, f"{by_w.iloc[j, i]:.2f}", ha="center", va="center", fontsize=7)
axes[1].set_title("MASE by backtest window", fontsize=10)
plt.show()
covered = board["inside P10-P90"].dropna()
fig, ax = plt.subplots(figsize=(12, 2.8))
ax.bar(covered.index, covered.values, color=[PINK if m in FOUNDATION else SLATE for m in covered.index])
ax.axhline(0.8, color="black", lw=1, ls="--")
ax.yaxis.set_major_formatter(plt.matplotlib.ticker.PercentFormatter(1.0)); ax.set_ylim(0, 1.05)
ax.set_title("Share of days inside the 80% interval (target: 80%)", fontsize=10); plt.show()1.7 6. Forecast with the backtest winner
best = board["MASE"].idxmin()
print("backtest winner:", best)
future = pd.date_range(y.index[-1], periods=H + 1, freq=FREQ)[1:]
if best in FOUNDATION:
mid, lo, hi = forecast_fm(best, y.values, H)
else:
out = StatsForecast(models=[CLASSICAL[best](SEASON)], freq=FREQ).forecast(
df=pd.DataFrame({"unique_id": "s", "ds": y.index, "y": y.values}), h=H, level=[80])
mid, lo, hi = out[best].values, out[f"{best}-lo-80"].values, out[f"{best}-hi-80"].values
fig, ax = plt.subplots()
ax.plot(y.index[-90:], y.values[-90:], color=SLATE, lw=1, label="actual")
ax.fill_between(future, lo, hi, color=PINK, alpha=0.2, label="P10–P90")
ax.plot(future, mid, color=PINK, lw=2, label=best)
ax.legend(loc="upper left"); ax.set_title("Next 28 days"); plt.show()
pd.DataFrame({"P10": lo, "P50": mid, "P90": hi}, index=future).round(0).head(14)1.8 7. Your conclusions
Write two or three sentences for each — the numbers above are your evidence:
- Which forecaster won, and by how much over “same period last week”? In how many windows?
- Is the range calibrated? Which model’s P10–P90 is closest to 80%?
- How does the error grow with the horizon? How often would you re-forecast?
- What would you plan with this forecast — and what does it not know about the coming weeks?