{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.11" }, "colab": { "provenance": [] } }, "cells": [ { "cell_type": "markdown", "id": "cell-0000", "metadata": {}, "source": [ "# Forecasting Models\n", "\n", "Split time series without cheating, beat naive baselines with lag features, and meet Holt-Winters and SARIMA — evaluated honestly with MAE and MAPE.\n", "\n", "*Part of the free [Machine Learning](https://ramadnsyh.dev/courses/machine-learning) course by [Muhammad Ramadiansyah](https://ramadnsyh.dev). This notebook is generated from the interactive lesson — [read it online](https://ramadnsyh.dev/courses/machine-learning/time-series-modeling).*" ] }, { "cell_type": "code", "id": "cell-0001", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "%pip install -q statsmodels" ] }, { "cell_type": "markdown", "id": "cell-0002", "metadata": {}, "source": [ "You can now take a series apart; this lesson is about predicting where it\n", "goes next. We'll start with the single most common mistake in forecasting\n", "(shuffled splits), establish baselines that are embarrassingly hard to beat,\n", "turn forecasting into a regression problem scikit-learn can solve, and then\n", "meet the classical specialists: exponential smoothing and ARIMA." ] }, { "cell_type": "markdown", "id": "cell-0003", "metadata": {}, "source": [ "## Never shuffle a time series\n", "\n", "`train_test_split` shuffles by default, and on time series that's not a\n", "small mistake — it's **data leakage**. A shuffled split trains on 1959 and\n", "tests on 1955: the model has literally seen the future, autocorrelation\n", "hands it the answers, and your test score becomes fiction. The honest split\n", "is **temporal**: train on the past, test on the most recent stretch, because\n", "that's exactly the situation the deployed model will face.\n", "\n", "For cross-validation, scikit-learn's `TimeSeriesSplit` respects time: every\n", "fold trains on an expanding window of the past and tests on the block right\n", "after it." ] }, { "cell_type": "code", "id": "cell-0004", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "from sklearn.model_selection import TimeSeriesSplit\n", "\n", "months = pd.date_range(\"2015-01-01\", periods=120, freq=\"MS\")\n", "\n", "tscv = TimeSeriesSplit(n_splits=4, test_size=12)\n", "for fold, (train_idx, test_idx) in enumerate(tscv.split(np.arange(120)), 1):\n", " tr0, tr1 = months[train_idx[0]], months[train_idx[-1]]\n", " te0, te1 = months[test_idx[0]], months[test_idx[-1]]\n", " print(f\"fold {fold}: train {tr0:%Y-%m} to {tr1:%Y-%m} \"\n", " f\"({len(train_idx):3d} mo) | test {te0:%Y-%m} to {te1:%Y-%m}\")" ] }, { "cell_type": "markdown", "id": "cell-0005", "metadata": {}, "source": [ "Each fold answers the question a stakeholder actually asks: \"if I had built\n", "this model a year ago, how would it have done?\" Sliding this scheme forward\n", "one step at a time — refit, predict the next point, repeat — is called\n", "**walk-forward validation**, the gold standard when you can afford the\n", "compute." ] }, { "cell_type": "markdown", "id": "cell-0006", "metadata": {}, "source": [ "## Baselines first: naive and seasonal-naive\n", "\n", "Before any model, establish the score to beat. Two forecasting baselines\n", "are so strong they're humbling:\n", "\n", "- **Naive** — tomorrow equals today. The forecast is a flat line at the\n", " last observed value.\n", "- **Seasonal naive** — this July equals last July. The forecast repeats the\n", " final observed seasonal cycle.\n", "\n", "We'll judge them with two metrics: **MAE** (mean absolute error — average\n", "miss, in the series' own units) and **MAPE** (mean absolute percentage\n", "error — average miss as a percentage, comparable across series)." ] }, { "cell_type": "code", "id": "cell-0007", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from sklearn.metrics import mean_absolute_error, mean_absolute_percentage_error\n", "\n", "rng = np.random.default_rng(42)\n", "months = pd.date_range(\"2015-01-01\", periods=120, freq=\"MS\")\n", "t = np.arange(120)\n", "pattern = np.array([-0.9, -1.2, 0.1, 0.0, 0.2, 1.4,\n", " 2.6, 2.4, 0.7, -0.3, -1.5, -0.4])\n", "y = pd.Series((120 + 1.8 * t) * (1 + 0.13 * pattern[t % 12])\n", " * rng.normal(1, 0.02, 120), index=months)\n", "\n", "train, test = y[:-24], y[-24:]\n", "\n", "naive = pd.Series(train.iloc[-1], index=test.index)\n", "seasonal_naive = pd.Series(np.tile(train[-12:].to_numpy(), 2), index=test.index)\n", "\n", "for name, pred in [(\"naive\", naive), (\"seasonal naive\", seasonal_naive)]:\n", " mae = mean_absolute_error(test, pred)\n", " mape = mean_absolute_percentage_error(test, pred)\n", " print(f\"{name:15s} MAE = {mae:6.1f} MAPE = {mape:.1%}\")\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.plot(y.index[-60:], y[-60:], label=\"actual\")\n", "plt.plot(test.index, naive, \"--\", label=\"naive\")\n", "plt.plot(test.index, seasonal_naive, \"--\", label=\"seasonal naive\")\n", "plt.axvline(test.index[0], color=\"gray\", lw=1)\n", "plt.legend(); plt.title(\"Baselines on the held-out tail\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cell-0008", "metadata": {}, "source": [ "The flat naive line is hopeless on seasonal data, but the seasonal naive\n", "tracks the shape well — its only sin is missing the trend. Any model you\n", "build must beat that dashed line, or it isn't earning its complexity." ] }, { "cell_type": "markdown", "id": "cell-0009", "metadata": {}, "source": [ "## Forecasting as regression: lag features\n", "\n", "Here's the trick that connects this module to everything you've learned:\n", "**turn the series into a supervised table**. Each row's features are its own\n", "past — `lag_1` (last month), `lag_12` (same month last year) — plus calendar\n", "features like the month number. The target is the current value. Then any\n", "regressor you already know can forecast." ] }, { "cell_type": "code", "id": "cell-0010", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "from sklearn.linear_model import LinearRegression\n", "from sklearn.metrics import mean_absolute_error, mean_absolute_percentage_error\n", "\n", "rng = np.random.default_rng(42)\n", "months = pd.date_range(\"2015-01-01\", periods=120, freq=\"MS\")\n", "t = np.arange(120)\n", "pattern = np.array([-0.9, -1.2, 0.1, 0.0, 0.2, 1.4,\n", " 2.6, 2.4, 0.7, -0.3, -1.5, -0.4])\n", "y = pd.Series((120 + 1.8 * t) * (1 + 0.13 * pattern[t % 12])\n", " * rng.normal(1, 0.02, 120), index=months)\n", "\n", "# Build the supervised table: lags + calendar\n", "df = pd.DataFrame({\"y\": y})\n", "for lag in (1, 2, 12):\n", " df[f\"lag_{lag}\"] = df[\"y\"].shift(lag)\n", "df[\"month\"] = df.index.month\n", "df = pd.get_dummies(df, columns=[\"month\"], dtype=float).dropna()\n", "\n", "train, test = df[:-24], df[-24:]\n", "X_train, y_train = train.drop(columns=\"y\"), train[\"y\"]\n", "X_test, y_test = test.drop(columns=\"y\"), test[\"y\"]\n", "\n", "model = LinearRegression().fit(X_train, y_train)\n", "pred = pd.Series(model.predict(X_test), index=y_test.index)\n", "\n", "mae = mean_absolute_error(y_test, pred)\n", "mape = mean_absolute_percentage_error(y_test, pred)\n", "print(f\"lag regression MAE = {mae:6.1f} MAPE = {mape:.1%}\")\n", "\n", "plt.figure(figsize=(9, 4))\n", "plt.plot(y.index[-60:], y[-60:], label=\"actual\")\n", "plt.plot(pred.index, pred, \"--\", lw=2, label=\"lag-feature forecast\")\n", "plt.axvline(pred.index[0], color=\"gray\", lw=1)\n", "plt.legend(); plt.title(\"One-step-ahead forecasts vs actual\")\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "cell-0011", "metadata": {}, "source": [ "The regression crushes both baselines — the lags carry the level and the\n", "season, the month dummies mop up the rest. One honest caveat: these are\n", "**one-step-ahead** predictions, because each test row's lags come from\n", "*actual* observed values. To forecast 24 months into the unknown future\n", "you'd predict one step, feed that prediction back in as `lag_1`, and repeat\n", "— **recursive forecasting** — which lets errors compound as the horizon\n", "grows. This lag-feature recipe is exactly how gradient-boosting models win\n", "most forecasting competitions today; swap `LinearRegression` for\n", "`HistGradientBoostingRegressor` and you have a genuinely modern pipeline." ] }, { "cell_type": "markdown", "id": "cell-0012", "metadata": {}, "source": [ "## Exponential smoothing and Holt-Winters\n", "\n", "The classical specialists take a different route: instead of a feature\n", "table, they maintain running estimates of the components and update them\n", "with each observation. **Simple exponential smoothing** tracks the level as\n", "a weighted average that decays exponentially into the past — good for\n", "series with no trend or season. **Holt's method** adds a second equation\n", "tracking the trend, and **Holt-Winters** adds a third for seasonality —\n", "level, trend, and season, each smoothed with its own parameter. Fit it in\n", "statsmodels (notebook/Colab — not the browser):" ] }, { "cell_type": "code", "id": "cell-0013", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import pandas as pd\n", "from sklearn.metrics import mean_absolute_error\n", "from statsmodels.tsa.holtwinters import ExponentialSmoothing\n", "\n", "url = \"https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv\"\n", "y = pd.read_csv(url, parse_dates=[\"Month\"], index_col=\"Month\")[\"Passengers\"]\n", "train, test = y[:-24], y[-24:]\n", "\n", "hw = ExponentialSmoothing(\n", " train, trend=\"add\", seasonal=\"mul\", seasonal_periods=12\n", ").fit()\n", "\n", "forecast = hw.forecast(24)\n", "print(f\"Holt-Winters MAE on the 24-month holdout: \"\n", " f\"{mean_absolute_error(test, forecast):.1f}\")\n", "\n", "ax = y.plot(figsize=(10, 4), label=\"actual\")\n", "forecast.plot(ax=ax, style=\"--\", label=\"Holt-Winters forecast\")\n", "ax.legend()" ] }, { "cell_type": "markdown", "id": "cell-0014", "metadata": {}, "source": [ "Note `seasonal=\"mul\"` — the airline series is multiplicative, as its\n", "fanning peaks told us last lesson. Unlike the lag-regression above, this is\n", "a true 24-step-ahead forecast made from training data alone, and\n", "Holt-Winters remains a ferociously strong benchmark for seasonal business\n", "data." ] }, { "cell_type": "markdown", "id": "cell-0015", "metadata": {}, "source": [ "## ARIMA and SARIMA, honestly\n", "\n", "**ARIMA(p, d, q)** models a series as a linear function of its own past:\n", "`p` autoregressive terms (past values), `d` rounds of differencing to reach\n", "stationarity, and `q` moving-average terms (past forecast errors). Plain\n", "ARIMA has no notion of seasonality, so in practice you use **SARIMA**, which\n", "adds a seasonal quadruple (P, D, Q, s). The honest summary: ARIMA is\n", "statistically elegant, gives principled confidence intervals, and rewards\n", "expertise — analysts read ACF/PACF plots and compare AIC scores to choose\n", "orders — but it's fiddly to tune, assumes linear dynamics, and on plenty of\n", "real data a seasonal-naive baseline or a lag-feature gradient booster\n", "matches it. Treat it as one candidate to evaluate, not a destination." ] }, { "cell_type": "code", "id": "cell-0016", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import pandas as pd\n", "from sklearn.metrics import mean_absolute_error\n", "from statsmodels.tsa.statespace.sarimax import SARIMAX\n", "\n", "url = \"https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv\"\n", "y = pd.read_csv(url, parse_dates=[\"Month\"], index_col=\"Month\")[\"Passengers\"]\n", "train, test = y[:-24], y[-24:]\n", "\n", "sarima = SARIMAX(train, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12)).fit(disp=False)\n", "forecast = sarima.forecast(24)\n", "print(f\"SARIMA MAE on the 24-month holdout: {mean_absolute_error(test, forecast):.1f}\")\n", "print(f\"AIC: {sarima.aic:.1f}\")" ] }, { "cell_type": "markdown", "id": "cell-0017", "metadata": {}, "source": [ "> **Choosing between all of these**\n", "> \n", "> Evaluate every candidate — baselines, Holt-Winters, SARIMA, lag-feature\n", "> regression — with the same temporal holdout and the same MAE/MAPE, and let\n", "> the numbers decide. On short seasonal business series, Holt-Winters and\n", "> SARIMA are hard to beat; with many related series or rich extra features\n", "> (prices, promotions, weather), the ML route usually wins." ] }, { "cell_type": "markdown", "id": "cell-0018", "metadata": {}, "source": [ "### 🏋️ Exercise — Beat the linear forecaster\n", "\n", "Extend the lag-feature model from the lesson: add `lag_3` and a 3-month\n", "rolling-mean feature (built from `shift(1)` so it never touches the current\n", "value — why does that matter?), and try `Ridge` alongside\n", "`LinearRegression`. Keep the same last-24-months holdout and compare MAE\n", "against the lesson's model. How much did the extra features buy you?" ] }, { "cell_type": "code", "id": "cell-0019", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "# Your solution here" ] }, { "cell_type": "markdown", "id": "cell-0020", "metadata": {}, "source": [ "
✅ Show solution\n", "\n", "```python\n", "import numpy as np\n", "import pandas as pd\n", "from sklearn.linear_model import LinearRegression, Ridge\n", "from sklearn.metrics import mean_absolute_error\n", "\n", "rng = np.random.default_rng(42)\n", "months = pd.date_range(\"2015-01-01\", periods=120, freq=\"MS\")\n", "t = np.arange(120)\n", "pattern = np.array([-0.9, -1.2, 0.1, 0.0, 0.2, 1.4,\n", " 2.6, 2.4, 0.7, -0.3, -1.5, -0.4])\n", "y = pd.Series((120 + 1.8 * t) * (1 + 0.13 * pattern[t % 12])\n", " * rng.normal(1, 0.02, 120), index=months)\n", "\n", "df = pd.DataFrame({\"y\": y})\n", "for lag in (1, 2, 3, 12):\n", " df[f\"lag_{lag}\"] = df[\"y\"].shift(lag)\n", "df[\"roll_3\"] = df[\"y\"].shift(1).rolling(3).mean() # shift first: no leakage\n", "df[\"month\"] = df.index.month\n", "df = pd.get_dummies(df, columns=[\"month\"], dtype=float).dropna()\n", "\n", "train, test = df[:-24], df[-24:]\n", "X_train, y_train = train.drop(columns=\"y\"), train[\"y\"]\n", "X_test, y_test = test.drop(columns=\"y\"), test[\"y\"]\n", "\n", "for name, model in [(\"LinearRegression\", LinearRegression()),\n", " (\"Ridge(alpha=1.0)\", Ridge(alpha=1.0))]:\n", " pred = model.fit(X_train, y_train).predict(X_test)\n", " print(f\"{name:18s} MAE = {mean_absolute_error(y_test, pred):.2f}\")\n", "\n", "# The extra lag and rolling-mean feature shave the MAE a little more;\n", "# Ridge behaves almost identically here since features are few and clean.\n", "```\n", "\n", "
" ] }, { "cell_type": "markdown", "id": "cell-0021", "metadata": {}, "source": [ "Next up: a new module — recommender systems, where the \"past behavior\"\n", "being modeled isn't your own history but everyone else's." ] } ] }