examples/time_series_cv.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZaGULROK7yS0"
      },
      "source": [
        "\n",
        "# Time Series CV\n",
        "\n",
        "Time series cross validation is deserves special treatment. This is due in part because of autocorrelction factors. `tablespoon` offerts a `TimeSeriesInitialSplit` class with the `split` method. It was designed to fit in naturally with what some users may be used to from `sklearn`. \n",
        "\n",
        "Some may wonder why `sklearn` is not sufficient. It is possible to use the methods in `sklearn` to make a work-around solution. This method is more natural as it allows users to explicitly define their intial time period, increment size, and gap. \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E3YutJBP8aeW"
      },
      "source": [
        "## Rolling Origin Cross Validation\n",
        "\n",
        "![image](https://raw.githubusercontent.com/alexhallam/tablespoon/main/examples/cv1-1.svg)\n",
        "\n",
        "*Source of the above image comes from Rob Hyndman's book fpp3 https://otexts.com/fpp3/tscv.html*\n",
        "\n",
        "Assume the blue dots represents dates to train on and the red dots are for testing.\n",
        "\n",
        "In this picture the initial size of the blue dots are much larger than the size of the testing dots. (6 blue dots vs 1 red). With each iteration train size increases by 1 blue dot.\n",
        "\n",
        "`TimeSeriesInitialSplit` allows this kind of flexibility. A user may set a larger initial training period, which is not bound by the size that the training set is increases by on each iteration. It decouples the initial size from the increment size."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zPjN8IW9-i1b"
      },
      "source": [
        "# Example 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dAyt8T2B-tzh"
      },
      "source": [
        "This example simple shows the functionality `TimeSeriesInitialSplit`. Note the initial time period is set to 7 and increments by 7 as well. There is no gap between the training periods and the test periods. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YirEvGcO8ekh",
        "outputId": "9b575dd0-a0c8-4391-dd5b-2917f568f52b"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "TRAIN: [0 1 2 3 4 5 6] TEST: [ 7  8  9 10 11 12 13]\n",
            "fold number: 1\n",
            "TRAIN: [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13] TEST: [14 15 16 17 18 19 20]\n",
            "fold number: 2\n",
            "TRAIN: [ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20] TEST: [21 22 23 24 25 26 27]\n",
            "fold number: 3\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from tablespoon.model_selection import TimeSeriesInitialSplit\n",
        "\n",
        "X = np.arange(0, 28)\n",
        "tscv = TimeSeriesInitialSplit(initial=7, increment_size=7, gap=0)\n",
        "fold_counter = 0\n",
        "for train_index, test_index in tscv.split(X):\n",
        "    print(\"TRAIN:\", train_index,\"TEST:\", test_index)\n",
        "    fold_counter += 1\n",
        "    print(f\"fold number: {fold_counter}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kikjZjuo_J3v"
      },
      "source": [
        "# Example 2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S_5Ap9Dz_MLO"
      },
      "source": [
        "This is closer to a real world example."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "metadata": {
        "id": "75rJkTaM_TSv"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from tablespoon.model_selection import TimeSeriesInitialSplit\n",
        "from tablespoon.data import WALMART_TX\n",
        "from tablespoon.forecasters import Snaive, Mean, Naive\n",
        "import pandas as pd\n",
        "\n",
        "X = WALMART_TX.assign(ds = lambda df: pd.to_datetime(df.ds))\n",
        "tscv = TimeSeriesInitialSplit(initial= 365, increment_size=7, gap=2)\n",
        "fold_counter = 0\n",
        "vec_mean = np.asarray([])\n",
        "vec_naive = np.asarray([])\n",
        "vec_snaive = np.asarray([])\n",
        "\n",
        "def crps(scalar_y, vec_of_forecast):\n",
        "    x = np.sort(vec_of_forecast)\n",
        "    m = len(vec_of_forecast)\n",
        "    return (2 / m) * np.mean((x - scalar_y) * (m * np.where(scalar_y < x, 1, 0) - np.arange(start=0, stop=m, step=1) + 1 / 2))\n",
        "for train_index, test_index in tscv.split(X):\n",
        "    train, test = X.iloc[train_index], X.iloc[test_index]\n",
        "    h = test.shape[0]\n",
        "    # make some forecasts from different models\n",
        "    df_n = (Naive().predict(df_historical = train, horizon=h, frequency=\"D\", lag = 1, uncertainty_samples = 500).assign(model = 'naive'))    \n",
        "    df_sn = (Snaive().predict(df_historical = train, horizon=h, frequency=\"D\", lag = 7, uncertainty_samples = 500).assign(model = 'snaive'))\n",
        "    df_m = (Mean().predict(df_historical = train, horizon=h, frequency=\"D\", uncertainty_samples = 500).assign(model = 'mean'))\n",
        "    # merge each test set to each model's forecast\n",
        "    df_forecast_actual_n = test.merge(df_n, how='left', on = 'ds')    \n",
        "    df_forecast_actual_sn = test.merge(df_sn, how='left', on = 'ds')\n",
        "    df_forecast_actual_m = test.merge(df_m, how='left', on = 'ds')\n",
        "    # put all the forecast models together\n",
        "    df_forecasts = pd.concat([df_forecast_actual_n, df_forecast_actual_sn, df_forecast_actual_m], axis=0)\n",
        "    # calculate mean crps by model type\n",
        "    df_crps = (df_forecasts\n",
        "            .groupby(by = ['model'], as_index=False)\n",
        "            .apply(lambda df: crps(df[\"y\"].iat[0], df[\"y_sim\"]))\n",
        "            .assign(metric=\"crps\")\n",
        "            .rename(columns={None: \"value\"}))\n",
        "    score_mean = df_crps.value # mean, naive, snaive\n",
        "    vec_mean = np.concatenate((vec_mean, score_mean[0].flatten()), axis = 0)\n",
        "    vec_naive = np.concatenate((vec_naive, score_mean[1].flatten()), axis = 0)\n",
        "    vec_snaive = np.concatenate((vec_snaive, score_mean[2].flatten()), axis = 0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BKINbe8NvTP2"
      },
      "source": [
        "Every forecast method in `tablespoon` is a reference/benchmark method. Find which of the forecast methods perform the best on average. That will be used as the reference.\n",
        "\n",
        "To calculate the skill score the following formulua is used.\n",
        "\n",
        "$$\n",
        "\\frac{\\text{CRPS}_{\\text{Naïve}} - \\text{CRPS}_{\\text{Alternative}}}{\\text{CRPS}_{\\text{Naïve}}}\n",
        "$$\n",
        "\n",
        "With some simplification this is the same as the following\n",
        "\n",
        "$$\n",
        "1 - \\frac{\\text{CRPS}_{\\text{Alternative}}}{\\text{CRPS}_{\\text{Naïve}}}\n",
        "$$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xc916Yjlha22",
        "outputId": "7fee5c59-6519-4d9d-ad54-2ded5f93f8a3"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "avg_mean_score 557.128\n",
            "avg_naive_score 1057.388\n",
            "avg_snaive_score 589.02\n"
          ]
        }
      ],
      "source": [
        "avg_mean_score = np.mean(vec_mean)\n",
        "avg_naive_score = np.mean(vec_naive)\n",
        "avg_snaive_score = np.mean(vec_snaive)\n",
        "print(f'avg_mean_score {np.round(avg_mean_score, 3)}')\n",
        "print(f'avg_naive_score {np.round(avg_naive_score, 3)}')\n",
        "print(f'avg_snaive_score {np.round(avg_snaive_score, 3)}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bRkxJ0MJwDLq"
      },
      "source": [
        "Of all the benchmark method the `Mean` is the best performer. This will be used as a the reference for the skill score calculation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "metadata": {
        "id": "s5OOu5KAr5Yq"
      },
      "outputs": [],
      "source": [
        "def skill_score(best_naive_method, alternative):\n",
        "  return 1 - (alternative / best_naive_method)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Eik7DZ24sdZZ",
        "outputId": "1b75f4f6-0e2a-45ca-aa31-356e6cd998b0"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "mean_skill 0.0\n",
            "naive_skill -0.898\n",
            "snaive_skill -0.057\n"
          ]
        }
      ],
      "source": [
        "mean_skill = skill_score(avg_mean_score, avg_mean_score)\n",
        "naive_skill = skill_score(avg_mean_score, avg_naive_score)\n",
        "snaive_skill = skill_score(avg_mean_score, avg_snaive_score)\n",
        "print(f\"mean_skill {np.round(mean_skill, 3)}\")\n",
        "print(f\"naive_skill {np.round(naive_skill, 3)}\")\n",
        "print(f\"snaive_skill {np.round(snaive_skill, 3)}\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "64CERqbQyzlk"
      },
      "source": [
        "## Notes on Skill Scores\n",
        "\n",
        "1. If the alternative forecast error is equal to the forecast error of the reference then the skill score will be `0`.\n",
        "\n",
        "2. If the alternative forecast error is perfect, score_error = 0, then then the skill score will be `1`.\n",
        "\n",
        "3. If the alternative forecast error is worse than the reference then the skill score will result in a value `< 0`.\n",
        "\n",
        "In summary, skill scores have the range `[-inf, 1]`"
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [],
      "name": "time_series_cv.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3.9.2 64-bit",
      "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.9.2"
    },
    "vscode": {
      "interpreter": {
        "hash": "916dbcbb3f70747c44a77c7bcd40155683ae19c65e1c03b4aa3499c5328201f1"
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}