notebooks/scalable-forecasting-with-functime.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2828f80c",
   "metadata": {
    "papermill": {
     "duration": 0.010867,
     "end_time": "2023-10-18T00:04:29.849976",
     "exception": false,
     "start_time": "2023-10-18T00:04:29.839109",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "![image.png](https://github.com/descendant-ai/functime/raw/main/docs/img/banner_dark_bg.png)\n",
    "\n",
    "[functime](functime.ai) is a lightweight time series library that enables scaling forecasts of large panel datasets effortless, thanks to Polars and global forecasting.\n",
    "\n",
    "In other words: with functime, you can **train forecasting models on thousands of time series on your laptop** (no PySpark required!), and feature engineering is blazingly fast thanks to Polars' parralelisation capabilities.\n",
    "\n",
    "# What you'll see\n",
    "\n",
    "1. How functime supports you in displaying panel datasets.\n",
    "2. How functions natively operate on panel datasets, abstracting away the aggregations.\n",
    "3. Provides naive models to benchmark your forecasts.\n",
    "4. Has a sklearn-compatible API.\n",
    "5. Makes it easy to perform time-series cross-validation.\n",
    "6. Provides utilities to rank your forecasts and residuals to improve model diagnostics.\n",
    "7. Provides plotting utilities to spot poorly-behaved time series in the panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f704b48",
   "metadata": {
    "_cell_guid": "b1076dfc-b9ad-4769-8c92-a6c4dae69d19",
    "_uuid": "8f2839f25d086af736a60e9eeb907d3b93b6e0e5",
    "papermill": {
     "duration": 39.411306,
     "end_time": "2023-10-18T00:05:09.271461",
     "exception": false,
     "start_time": "2023-10-18T00:04:29.860155",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import functime\n",
    "except ModuleNotFoundError:\n",
    "    # NOTE: this might fail on Kaggle: it requires an account and phone verification.\n",
    "    # If it happens, go to the \"sidebar > Notebook options\" and verify your phone number.\n",
    "    !pip install --quiet -- \"functime[plot,lgb]\"\n",
    "\n",
    "import polars as pl\n",
    "import functime\n",
    "\n",
    "print(\n",
    "    f\"polars: {pl.__version__}\",\n",
    "    f\"functime: {functime.__version__}\",\n",
    "    sep=\"\\n\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd3c435d",
   "metadata": {
    "papermill": {
     "duration": 0.01008,
     "end_time": "2023-10-18T00:05:09.293099",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.283019",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 1. Data Description and Problem Statement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e70262a",
   "metadata": {
    "papermill": {
     "duration": 0.010924,
     "end_time": "2023-10-18T00:05:09.314194",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.30327",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The dataset contains monthly observations of some 70 commodities, between 1960 and 2023."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7221b7ce",
   "metadata": {
    "papermill": {
     "duration": 0.236079,
     "end_time": "2023-10-18T00:05:09.560478",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.324399",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "url = \"https://github.com/neocortexdb/functime/raw/main/data/commodities.parquet\"\n",
    "y = pl.scan_parquet(url).with_columns(pl.col(\"time\").cast(pl.Date))\n",
    "\n",
    "y.group_by(\"commodity_type\").tail().collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "505b50f4",
   "metadata": {
    "papermill": {
     "duration": 0.012478,
     "end_time": "2023-10-18T00:05:09.583463",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.570985",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "A panel dataset always contains a column that identifies the time series, the datetime indexes, and the target values. There can also be additional covariates, such as holiday effects or macroeconomic variables.\n",
    "\n",
    "functime assumes that the first column is always the time series identifier (in the documentation, it is referred to as _entity column_), and is followed by the time column and the target column.\n",
    "\n",
    "> RFC: Do you like these defaults? What do you think of a `PanelDataFrame` class that inherits from a `polars.DataFrame`, whose constructor allows you to specify the names of the entity, time and target cols?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfaaa2b1",
   "metadata": {
    "papermill": {
     "duration": 0.019087,
     "end_time": "2023-10-18T00:05:09.613348",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.594261",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "entity_col, time_col, target_col = y.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74cd2dbb",
   "metadata": {
    "papermill": {
     "duration": 0.010058,
     "end_time": "2023-10-18T00:05:09.633751",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.623693",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The data has a monthly frequency. Our goal is to forecast the price of the next quarter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2d55a6b",
   "metadata": {
    "papermill": {
     "duration": 0.018149,
     "end_time": "2023-10-18T00:05:09.662129",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.64398",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "freq = \"1mo\"\n",
    "forecasting_horizon = 4 # a quarter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "650f3a5c",
   "metadata": {
    "papermill": {
     "duration": 0.010431,
     "end_time": "2023-10-18T00:05:09.682845",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.672414",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 2. Exploratory Data Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f530c1ff",
   "metadata": {
    "papermill": {
     "duration": 0.010013,
     "end_time": "2023-10-18T00:05:09.766248",
     "exception": false,
     "start_time": "2023-10-18T00:05:09.756235",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The number of observations per entity is not constant: as we can see from the plot below, the number ranges from approximately 250 to 750."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5177b569",
   "metadata": {},
   "outputs": [],
   "source": [
    "from functime import plotting\n",
    "\n",
    "plotting.plot_entities(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57fd2045",
   "metadata": {
    "papermill": {
     "duration": 0.012366,
     "end_time": "2023-10-18T00:05:13.590863",
     "exception": false,
     "start_time": "2023-10-18T00:05:13.578497",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "How do these series look like? Since there are 70 series (or thousands, potentially), we cannot inspect each of them. We could sample some:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48f4463d",
   "metadata": {
    "papermill": {
     "duration": 0.379802,
     "end_time": "2023-10-18T00:05:13.981761",
     "exception": false,
     "start_time": "2023-10-18T00:05:13.601959",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Here is a little bit of magic to display a sample of time series.\n",
    "# You won't have to do this in functime 0.10.0 coming out December 2023!\n",
    "\n",
    "k = 6\n",
    "seed = 42\n",
    "\n",
    "commodities = (\n",
    "    y\n",
    "    .select(pl.col(entity_col).unique())\n",
    ")\n",
    "\n",
    "# sample K entity IDs to plot\n",
    "sample = (\n",
    "    commodities\n",
    "    .collect()\n",
    "    .to_series()\n",
    "    .sample(k, seed=seed)\n",
    ")\n",
    "\n",
    "(\n",
    "    y\n",
    "    .filter(pl.col(entity_col).is_in(sample))\n",
    "    .pipe(\n",
    "        plotting.plot_panel,\n",
    "        title=f\"Sample of {k} time series\",\n",
    "        height=125*k,\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "635d63e7",
   "metadata": {
    "papermill": {
     "duration": 0.010931,
     "end_time": "2023-10-18T00:05:14.00399",
     "exception": false,
     "start_time": "2023-10-18T00:05:13.993059",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "> In the next release, there will be a `ids: int | Container[int] | str | Container[str]` parameter to plot specific entities; there will also be a `shuffle` parameter to sample `ids` if it is an int.\n",
    "\n",
    "For this reason, we should look at other time series summary statistics. One of such is the coefficient of variation, measured as the variance of a series divided by its mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc06fb85",
   "metadata": {
    "papermill": {
     "duration": 0.289798,
     "end_time": "2023-10-18T00:05:14.306018",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.01622",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# With the next release, you will also be able to show the top/bottom K \n",
    "# series in the panel ranked according a scoring function\n",
    "\n",
    "def coefficient_of_variation(col):\n",
    "    return pl.col(target_col).std() / pl.col(target_col).mean() \n",
    "\n",
    "top_k_cv = (\n",
    "    y\n",
    "    .group_by(entity_col)\n",
    "    .agg(coefficient_of_variation(time_col).alias(\"cv\"))\n",
    "    .top_k(k=k, by=\"cv\")\n",
    "    .collect()\n",
    "    .get_column(entity_col)\n",
    ")\n",
    "\n",
    "(\n",
    "    y\n",
    "    .filter(pl.col(entity_col).is_in(top_k_cv))\n",
    "    .pipe(\n",
    "        plotting.plot_panel,\n",
    "        title=f\"Top {k} series by coefficient of variation\",\n",
    "        height=125*k,\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e2243c",
   "metadata": {
    "papermill": {
     "duration": 0.012269,
     "end_time": "2023-10-18T00:05:14.636525",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.624256",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 3. Train-test split"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe102f5d",
   "metadata": {
    "papermill": {
     "duration": 0.011943,
     "end_time": "2023-10-18T00:05:14.660729",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.648786",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "functime features the familiar scikit-learn API:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915e9077",
   "metadata": {
    "papermill": {
     "duration": 0.021784,
     "end_time": "2023-10-18T00:05:14.694638",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.672854",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.cross_validation import train_test_split\n",
    "\n",
    "splitter = train_test_split(forecasting_horizon)\n",
    "y_train, y_test = splitter(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "759d4aa6",
   "metadata": {
    "papermill": {
     "duration": 0.011647,
     "end_time": "2023-10-18T00:05:14.718689",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.707042",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "functime functions, however, are designed to work with panel data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "747e85cf",
   "metadata": {
    "papermill": {
     "duration": 0.128961,
     "end_time": "2023-10-18T00:05:14.859611",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.73065",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "y_test.head(12).collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4301e12",
   "metadata": {
    "papermill": {
     "duration": 0.011962,
     "end_time": "2023-10-18T00:05:14.883961",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.871999",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 4. Benchmark forecast"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2e392fb",
   "metadata": {
    "papermill": {
     "duration": 0.011839,
     "end_time": "2023-10-18T00:05:14.907995",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.896156",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Some forecasting methods are extremely simple and surprisingly effective. Naive and seasonal naive forecasters are surprisingly hard to beat! You should always consider using the naive and seasonal naive forecasts as benchmarks. The naive forecaster repeats the last value and, for a multi-step forecast, returns a flat line. The seasonal naive returns the same pattern of length `sp` (seasonal periodicity). [reference](https://otexts.com/fpp3/simple-methods.html#simple-methods)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69278c55",
   "metadata": {
    "papermill": {
     "duration": 7.273909,
     "end_time": "2023-10-18T00:05:22.19466",
     "exception": false,
     "start_time": "2023-10-18T00:05:14.920751",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.forecasting import snaive\n",
    "\n",
    "forecaster_naive = snaive(\n",
    "    freq=freq,\n",
    "    sp=12,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "255ce487",
   "metadata": {
    "papermill": {
     "duration": 0.01201,
     "end_time": "2023-10-18T00:05:22.219224",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.207214",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Now we train the seasonal naive model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92b539f7",
   "metadata": {
    "papermill": {
     "duration": 0.147947,
     "end_time": "2023-10-18T00:05:22.379522",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.231575",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "_ = forecaster_naive.fit(y=y_train)\n",
    "y_bench_pred = forecaster_naive.predict(fh=forecasting_horizon)\n",
    "\n",
    "y_bench_pred.head(12)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bddd96d",
   "metadata": {},
   "source": [
    "## 4.1 Diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dca2677f",
   "metadata": {
    "papermill": {
     "duration": 0.201988,
     "end_time": "2023-10-18T00:05:22.593718",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.39173",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_forecasts(\n",
    "    y_test.filter(pl.col(entity_col).is_in(sample)).collect(),\n",
    "    y_bench_pred.filter(pl.col(entity_col).is_in(sample)),\n",
    "    title=f\"Sample of {k} naive forecasts\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ebf1d11",
   "metadata": {
    "papermill": {
     "duration": 0.013231,
     "end_time": "2023-10-18T00:05:22.619745",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.606514",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "Once again, how can we evaluate the goodness of fit across dozens of series? functime offers the `rank_*` functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da547da9",
   "metadata": {
    "papermill": {
     "duration": 0.134968,
     "end_time": "2023-10-18T00:05:22.768144",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.633176",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.evaluation import rank_point_forecasts\n",
    "\n",
    "rank_bench_point_forecast = rank_point_forecasts(\n",
    "    y_true=y_test,\n",
    "    y_pred=y_bench_pred,\n",
    "    sort_by=\"mape\",\n",
    "    descending=True,\n",
    ")\n",
    "\n",
    "rank_bench_point_forecast.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89aca476",
   "metadata": {
    "papermill": {
     "duration": 0.207695,
     "end_time": "2023-10-18T00:05:22.988393",
     "exception": false,
     "start_time": "2023-10-18T00:05:22.780698",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "worst_bench_point_forecast = rank_bench_point_forecast.head(k).get_column(entity_col)\n",
    "\n",
    "plotting.plot_forecasts(\n",
    "    y_test.filter(pl.col(entity_col).is_in(worst_bench_point_forecast)).collect(),\n",
    "    y_bench_pred.filter(pl.col(entity_col).is_in(worst_bench_point_forecast)),\n",
    "    title=f\"Worst {k} naive forecasts by mean absolute percentage error (MAPE)\",\n",
    "    height=125*k,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d53c4713",
   "metadata": {
    "papermill": {
     "duration": 0.013408,
     "end_time": "2023-10-18T00:05:23.238334",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.224926",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "There are, however, better approaches to glance the goodness of fit - which we shall see after we fit another model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ae38984",
   "metadata": {
    "papermill": {
     "duration": 0.013342,
     "end_time": "2023-10-18T00:05:23.265401",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.252059",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 5. Forecasting with global models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3aca3fe2",
   "metadata": {
    "papermill": {
     "duration": 0.013832,
     "end_time": "2023-10-18T00:05:23.292953",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.279121",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "What makes functime so fast is also the fact that only one model is fit on the whole dataset - hence why the name _global forecasting_. This strategy first proved successful with the [M5 competition](https://www.sciencedirect.com/science/article/pii/S0169207021001874?via%3Dihub), while there is a growing amount of [literature](https://www.sciencedirect.com/science/article/abs/pii/S0169207021000558?via%3Dihub) (mostly from Hyndman and his coauthors). This enables a [\"new paradigm of forecasting\"](https://blogs.sas.com/content/forecasting/2016/10/25/changing-the-paradigm-for-business-forecasting-part-10/):\n",
    "\n",
    "> [...] [T]he amount of time, effort, and resources spent on forecasting is not commensurate with the benefit achieved – the improvement in accuracy.\n",
    "\n",
    "> We spend far too many resources generating, reviewing, adjusting, and approving our forecasts, while almost invariably failing to achieve the level of accuracy desired. The evidence now shows that a large proportion of typical business forecasting efforts fail to improve the forecast, or even make it worse. So the conversation needs to change. The focus needs to change.\n",
    "\n",
    "> We need to shift our attention from esoteric model building to the forecasting process itself – its efficiency and its effectiveness."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd93ffdf",
   "metadata": {
    "papermill": {
     "duration": 0.013278,
     "end_time": "2023-10-18T00:05:23.319577",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.306299",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## 5.1 Linear models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4053d88a",
   "metadata": {
    "papermill": {
     "duration": 0.013106,
     "end_time": "2023-10-18T00:05:23.346019",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.332913",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "A simple linear regression can be made a global forecasting model simply by fitting in on the whole dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed88df51",
   "metadata": {
    "papermill": {
     "duration": 0.021961,
     "end_time": "2023-10-18T00:05:23.381376",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.359415",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.forecasting import linear_model\n",
    "from functime.preprocessing import scale, add_fourier_terms, roll\n",
    "from functime.seasonality import add_calendar_effects\n",
    "\n",
    "target_transforms = scale()\n",
    "feature_transforms = roll(window_sizes=(6, 12, 24), stats=(\"mean\", \"std\"), freq=freq)\n",
    "\n",
    "forecaster_linear = linear_model(\n",
    "    freq=freq,\n",
    "    lags=12,\n",
    "    target_transform=target_transforms,\n",
    "    feature_transform=feature_transforms,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2978a8e4",
   "metadata": {
    "papermill": {
     "duration": 0.014961,
     "end_time": "2023-10-18T00:05:23.41023",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.395269",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "In this case, we scale the target variable and perform a couple of other common steps in feature engineering for time series:\n",
    "\n",
    "1. Add 12 lagged values of the target variable.\n",
    "2. Compute the rolling mean and standard deviation with window size of 6 and 12 and 24 (6 total).\n",
    "\n",
    "functime also supports de-trending and de-seasonalisation, boxcox transform but also fractional differential transform as explained in the influential Advances in Financial Machine Learning by Marcos Lopez de Prado. See them all [here](https://docs.functime.ai/ref/preprocessing).\n",
    "\n",
    "To improve our model evaluation, we also perform time-series cross validation with 5 splits of length 180 (i.e., 15 years worth of data)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "621c7bd2",
   "metadata": {
    "papermill": {
     "duration": 0.02262,
     "end_time": "2023-10-18T00:05:23.449892",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.427272",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "backtesting_opts = dict(\n",
    "    y=y_train,\n",
    "    window_size=180, # 15 years of training data in each fold\n",
    "    test_size=forecasting_horizon,\n",
    "    step_size=1,\n",
    "    n_splits=5\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "988bf3e0",
   "metadata": {
    "papermill": {
     "duration": 4.911324,
     "end_time": "2023-10-18T00:05:28.375353",
     "exception": false,
     "start_time": "2023-10-18T00:05:23.464029",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "y_linear_preds, y_linear_resids = forecaster_linear.backtest(**backtesting_opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "475cb3ee",
   "metadata": {},
   "source": [
    "### 5.1.1 Diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d93ce9f3",
   "metadata": {
    "papermill": {
     "duration": 0.215697,
     "end_time": "2023-10-18T00:05:28.630112",
     "exception": false,
     "start_time": "2023-10-18T00:05:28.414415",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_backtests(\n",
    "    y_train.filter(pl.col(entity_col).is_in(sample)).collect(),\n",
    "    y_linear_preds.filter(pl.col(entity_col).is_in(sample)),\n",
    "    title=f\"Sample of {k} cross-validated linear forecasts\",\n",
    "    last_n=36,\n",
    "    height=125*k\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9938870c",
   "metadata": {
    "papermill": {
     "duration": 0.014708,
     "end_time": "2023-10-18T00:05:29.315183",
     "exception": false,
     "start_time": "2023-10-18T00:05:29.300475",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "We can even compare residuals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56cd9b38",
   "metadata": {
    "papermill": {
     "duration": 0.340875,
     "end_time": "2023-10-18T00:05:29.671086",
     "exception": false,
     "start_time": "2023-10-18T00:05:29.330211",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.evaluation import rank_residuals\n",
    "\n",
    "rank_linear_residuals = rank_residuals(\n",
    "    y_resids=y_linear_resids,\n",
    "    sort_by=\"abs_bias\",\n",
    "    descending=True\n",
    ")\n",
    "\n",
    "best_linear_residuals = rank_linear_residuals.tail(k).get_column(entity_col)\n",
    "\n",
    "plotting.plot_residuals(\n",
    "    y_resids=y_linear_resids.filter(pl.col(entity_col).is_in(best_linear_residuals)),\n",
    "    n_bins=200,\n",
    "    height=800,\n",
    "    width=1000,\n",
    "    title=f\"Top {k} linear forecast residuals\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cfd0bb8",
   "metadata": {
    "papermill": {
     "duration": 0.035161,
     "end_time": "2023-10-18T00:05:29.741605",
     "exception": false,
     "start_time": "2023-10-18T00:05:29.706444",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "## 5.2 Gradient Boosted Trees"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52d00526",
   "metadata": {
    "papermill": {
     "duration": 0.041639,
     "end_time": "2023-10-18T00:05:29.816015",
     "exception": false,
     "start_time": "2023-10-18T00:05:29.774376",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.forecasting import lightgbm\n",
    "\n",
    "forecaster_gbm = lightgbm(\n",
    "    freq=freq,\n",
    "    lags=12,\n",
    "    target_transform=target_transforms,\n",
    "    feature_transform=feature_transforms\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e96289c",
   "metadata": {
    "papermill": {
     "duration": 10.757739,
     "end_time": "2023-10-18T00:05:40.606568",
     "exception": false,
     "start_time": "2023-10-18T00:05:29.848829",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "y_lgb_preds, y_lgb_resids = forecaster_gbm.backtest(**backtesting_opts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab49aa28",
   "metadata": {
    "papermill": {
     "duration": 0.034065,
     "end_time": "2023-10-18T00:05:40.674214",
     "exception": false,
     "start_time": "2023-10-18T00:05:40.640149",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "### 5.1.2 Diagnostics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "116e378b",
   "metadata": {
    "papermill": {
     "duration": 0.21501,
     "end_time": "2023-10-18T00:05:40.926778",
     "exception": false,
     "start_time": "2023-10-18T00:05:40.711768",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_backtests(\n",
    "    y_train.filter(pl.col(entity_col).is_in(sample)).collect(),\n",
    "    y_lgb_preds.filter(pl.col(entity_col).is_in(sample)),\n",
    "    title=f\"Sample of {k} cross-validated LightGBM forecasts\",\n",
    "    last_n=36,\n",
    "    height=125*k\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51981125",
   "metadata": {
    "papermill": {
     "duration": 0.221433,
     "end_time": "2023-10-18T00:05:41.462778",
     "exception": false,
     "start_time": "2023-10-18T00:05:41.241345",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "rank_lgb_residuals = rank_residuals(\n",
    "    y_resids=y_lgb_resids,\n",
    "    sort_by=\"abs_bias\",\n",
    "    descending=True\n",
    ")\n",
    "\n",
    "best_lgb_residuals = rank_lgb_residuals.tail(k).get_column(entity_col)\n",
    "\n",
    "plotting.plot_residuals(\n",
    "    y_resids=y_lgb_resids.filter(pl.col(entity_col).is_in(best_lgb_residuals)),\n",
    "    n_bins=200,\n",
    "    height=800,\n",
    "    width=1000,\n",
    "    title=f\"Top {k} LightGBM forecast residuals\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa2d3c45",
   "metadata": {
    "papermill": {
     "duration": 0.057841,
     "end_time": "2023-10-18T00:05:41.577215",
     "exception": false,
     "start_time": "2023-10-18T00:05:41.519374",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "# 6. Comparing different models\n",
    "\n",
    "Currently naive forecasters raise an error when `.backtest` is called on them. The error is caused by computing the residuals; it will be fixed before the next minor release."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c0adec0",
   "metadata": {
    "_kg_hide-input": true,
    "papermill": {
     "duration": 1.316606,
     "end_time": "2023-10-18T00:05:42.946977",
     "exception": false,
     "start_time": "2023-10-18T00:05:41.630371",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from functime.backtesting import backtest\n",
    "from functime.cross_validation import sliding_window_split, train_test_split\n",
    "\n",
    "splitter = train_test_split(forecasting_horizon)\n",
    "y_train, y_test = splitter(y)\n",
    "\n",
    "# make a new forecast each month\n",
    "step_size = 1\n",
    "\n",
    "# backtest: refit one year\n",
    "n_splits = 12\n",
    "\n",
    "cv_sliding = sliding_window_split(\n",
    "    test_size=backtesting_opts[\"test_size\"],\n",
    "    step_size=backtesting_opts[\"step_size\"],\n",
    "    n_splits=backtesting_opts[\"n_splits\"],\n",
    "    window_size=backtesting_opts[\"window_size\"],\n",
    ")\n",
    "\n",
    "y_bench_preds = backtest(\n",
    "    forecaster=forecaster_naive,\n",
    "    fh=backtesting_opts[\"test_size\"],\n",
    "    y=backtesting_opts[\"y\"],\n",
    "    cv=cv_sliding,\n",
    "    residualize=False # currently unsupported, will be fixed in time for 0.10.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b36188",
   "metadata": {
    "papermill": {
     "duration": 0.05408,
     "end_time": "2023-10-18T00:05:43.05265",
     "exception": false,
     "start_time": "2023-10-18T00:05:42.99857",
     "status": "completed"
    },
    "tags": []
   },
   "source": [
    "The comet plot uses sMAPE (currently, will be fixed):\n",
    "\n",
    "$\\text{sMAPE} = \\text{mean}\\left(200|y_{t} - \\hat{y}_{t}|/(y_{t}+\\hat{y}_{t})\\right)$\n",
    "\n",
    "This metric has its weaknesses, despite being used (e.g. in M3):\n",
    "\n",
    "> However, if $y_t$ is close to zero, $\\hat{y}_{t}$ is also likely to be close to zero. Thus, the measure still involves division by a number close to zero, making the calculation unstable. Also, the value of sMAPE can be negative, so it is not really a measure of “absolute percentage errors” at all. [ref](https://otexts.com/fpp3/accuracy.html#percentage-errors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5688d28",
   "metadata": {},
   "source": [
    "## 6.1 Forecast Value Added Plot (FVA)\n",
    "\n",
    "Compares the performance of a foreaster's prediction against a naive foreaster, for all time series in the panel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a517fe9",
   "metadata": {
    "papermill": {
     "duration": 0.398366,
     "end_time": "2023-10-18T00:05:43.505431",
     "exception": false,
     "start_time": "2023-10-18T00:05:43.107065",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_fva(\n",
    "    y_train,\n",
    "    y_lgb_preds,\n",
    "    y_bench_preds,\n",
    "    height=900,\n",
    "    width=900,\n",
    "    title=\"SMAPE: Seasonal Naive (x axis) and LightGBM (y axis)\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e797427",
   "metadata": {
    "papermill": {
     "duration": 0.368752,
     "end_time": "2023-10-18T00:05:43.924812",
     "exception": false,
     "start_time": "2023-10-18T00:05:43.55606",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_fva(\n",
    "    y_train,\n",
    "    y_linear_preds,\n",
    "    y_bench_preds,\n",
    "    height=900,\n",
    "    width=900,\n",
    "    title=\"SMAPE: Seasonal Naive (x axis) and Linear (y axis)\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a6dca8c1",
   "metadata": {},
   "source": [
    "## 6.2 Comet Plot\n",
    "\n",
    "Plot the forecast errors against the coefficient of variation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fea9419",
   "metadata": {
    "papermill": {
     "duration": 0.354278,
     "end_time": "2023-10-18T00:05:44.331066",
     "exception": false,
     "start_time": "2023-10-18T00:05:43.976788",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_comet(\n",
    "    y_train=y_train.collect(),\n",
    "    y_test=y_train.collect(),\n",
    "    y_pred=y_lgb_preds,\n",
    "    height=900,\n",
    "    width=900,\n",
    "    title=\"LightGBM: Coefficient of Variation (x) vs SMAPE (y)\",\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bbe177c",
   "metadata": {
    "papermill": {
     "duration": 0.355876,
     "end_time": "2023-10-18T00:05:44.738241",
     "exception": false,
     "start_time": "2023-10-18T00:05:44.382365",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plotting.plot_comet(\n",
    "    y_train=y_train.collect(),\n",
    "    y_test=y_train.collect(),\n",
    "    y_pred=y_linear_preds,\n",
    "    height=900,\n",
    "    width=900,\n",
    "    title=\"Linear Model: Coefficient of Variation (x) vs SMAPE (y)\",\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kaggle": {
   "accelerator": "none",
   "dataSources": [],
   "dockerImageVersionId": 30587,
   "isGpuEnabled": false,
   "isInternetEnabled": true,
   "language": "python",
   "sourceType": "notebook"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 81.374164,
   "end_time": "2023-10-18T00:05:47.008724",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2023-10-18T00:04:25.63456",
   "version": "2.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}