time-fcasts-env.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d4467a20",
   "metadata": {},
   "source": [
    "# Forecasting Exercises\n",
    "\n",
    "In this chapter, we're going to do a tour of forecasting exercises: that is, the set of operations, like slicing up time, that you might need to do when performing a forecast. Although there will be some code in this chapter, we're mostly laying the theoretical groundwork.\n",
    "\n",
    "Let's start with some definitions. Features, or regressors, are labelled by $k = 0, \\dots, K-1$ and time by $t = 0, \\dots, T-1$. Note that this means there are $T$ time periods and $K$ features. We'll label different slices of time (or equivalently models trained on different data) by an index $\\mu$ beginning from $\\mu=1$, and we'll give these slices of time different labels: IS for in-sample and OS for out-of-sample. $f_\\mu$ is the model trained on the $\\mu$th set of in-sample data. The target variable (the number we are trying to forecast) is $\\left\\{y_{t+h}\\right\\}_{t=0}^{t=T-1}$, where $h$ is the number of time steps ahead we wish to forecast. Let $\\left\\{x_{tk}\\right\\}_{t=0}^{t=T-1}$ represent feature (or regressor) $k$. In general, a forecast implies $h>0$ (otherwise we're doing a nowcast).\n",
    "\n",
    "Let's first do some setup for the rest of the chapter. We'll be using the **numpy**, **pandas**, **matplotlib**, and **lets-plot** packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6f7024c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from lets_plot import *\n",
    "from lets_plot.mapping import as_discrete\n",
    "\n",
    "LetsPlot.setup_html(no_js=True)\n",
    "width_height_tiles = {\"width\": 0.9, \"height\": 0.5}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51a55374",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "import matplotlib_inline.backend_inline\n",
    "\n",
    "# Plot settings\n",
    "plt.style.use(\n",
    "    \"https://github.com/aeturrell/coding-for-economists/raw/main/plot_style.txt\"\n",
    ")\n",
    "\n",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")\n",
    "# Set pandas max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 15)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7aed965",
   "metadata": {},
   "source": [
    "## Two Period Forecast\n",
    "\n",
    "The most simple forecast you can imagine is if we have two periods (each of which may be made up of more than one value of $t$). The important thing is that the $t$ in the first, \"in-sample\" period do not overlap with the $t$ in the second, \"out-of-sample\" period. That's because we *only want to use information in the in-sample period to forecast the out-of-sample period*. There are some cases where you may want to \"forecast\" the in-sample period, but if you're interested in the real-world performance of your forecast, this is the setup you need.\n",
    "\n",
    "In this world of a single in-sample period and a single out-of-sample period, $\\mu=1$, and we could get drop that notation and just use IS and OS. But we'll leave it in to make the contrast with what's to come clearer.\n",
    "\n",
    "The exercise is to use a model, $f = f_{\\mu=1}$, that is trained to predict $\\{y_{t+h}, t:\\mu_{\\text{IS}}=1\\}$ using $\\{x_{tk}, t:\\mu_{\\text{IS}}=1\\}$ so that it can do an out-of-sample forecast of $\\{y_{t+h}, t:\\mu_\\text{OS}=1\\}$ using $\\{x_{tk}, t:\\mu_{\\text{OS}}=1\\}$. It's actually easier to see it with a picture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a658b0f1",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "s = 1\n",
    "alpha = 4\n",
    "T = 10\n",
    "\n",
    "\n",
    "def in_sample_block_e(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange(0, mu * s + alpha, s)})\n",
    "    xf[\"μ\"] = [mu] * (mu * s + alpha - 0)\n",
    "    xf[\"sample\"] = \"in-sample\"\n",
    "    return xf\n",
    "\n",
    "\n",
    "def os_sample_block_e(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange(mu * s + alpha, T, s)})\n",
    "    xf[\"μ\"] = [mu] * (T - (mu * s + alpha))\n",
    "    xf[\"sample\"] = \"out-of-sample\"\n",
    "    return xf\n",
    "\n",
    "\n",
    "def in_sample_block_r(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange((mu - 1) * s, mu * s + alpha, s)})\n",
    "    xf[\"μ\"] = [mu] * (mu * s + alpha - (mu - 1) * s)\n",
    "    xf[\"sample\"] = \"in-sample\"\n",
    "    return xf\n",
    "\n",
    "\n",
    "def os_sample_block_r(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange(mu * s + alpha, T, s)})\n",
    "    xf[\"μ\"] = [mu] * (T - (mu * s + alpha))\n",
    "    xf[\"sample\"] = \"out-of-sample\"\n",
    "    return xf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce1a7658",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "two_period = pd.concat(\n",
    "    [in_sample_block_e(1, s, alpha, T), os_sample_block_e(1, s, alpha, T)], axis=0\n",
    ")\n",
    "two_period[\"sample\"] = two_period[\"sample\"].astype(\"category\")\n",
    "\n",
    "(\n",
    "    ggplot(two_period, aes(x=\"t\", y=as_discrete(\"sample\"), fill=as_discrete(\"sample\")))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T + 1)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Accent\")\n",
    "    + ggsize(600, 200)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a67b4a8",
   "metadata": {},
   "source": [
    "In this case, the in-sample period runs from $t=0$ to $t=4$, while the out-of-sample period runs from $t=5$ to $t=9$. $T=10$.\n",
    "\n",
    "Let's see how it works with some fake data, and we'll assume $h=3$ in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0fc26f7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import string\n",
    "\n",
    "horizon = 3\n",
    "num = 10 + horizon\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    np.array(\n",
    "        [\n",
    "            range(num),\n",
    "            list(string.ascii_lowercase)[:num],\n",
    "            list(string.ascii_lowercase)[num : 2 * num],\n",
    "        ]\n",
    "    ).T,\n",
    "    columns=[\"Time\", \"y\", \"x\"],\n",
    ")\n",
    "df[\"y at t+h\"] = df[\"y\"].shift(-horizon)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fade8135",
   "metadata": {},
   "source": [
    "What we have here is a dataset with a time index and y and x values. $y_{t+h}$ has been added as a feature by shifting the $y$ variable back by $h=3$ steps. We're not going to be able to predict those NaN values, so the first thing to do is to drop those variables. Second, we'll want to add in which parts are in-sample and which are out-of-sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a542bdf4",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.dropna().copy()\n",
    "df[\"sample\"] = df[\"Time\"].apply(lambda x: \"IS\" if int(x) < 5 else \"OS\")\n",
    "df = df.set_index(\"Time\")\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c364ff4",
   "metadata": {},
   "source": [
    "Note how the sample can either be IS or OS-but not both.\n",
    "\n",
    "Now, if we were doing a forecast, we'd select the IS part to do our initial model building by running something like\n",
    "\n",
    "```python\n",
    "model = ols.fit(X=df.loc[df[\"sample\"]==\"IS\",\"x\"],\n",
    "                y=df.loc[df[\"sample\"]==\"IS\",\"y at t+h\"])\n",
    "```\n",
    "\n",
    "This would then be used to predict the out of sample values:\n",
    "\n",
    "```python\n",
    "y_t_plus_h_os = model.predict(X=df.loc[df[\"sample\"]==\"OS\",\"x\"])\n",
    "```\n",
    "\n",
    "To know how good the estimate was, we'd then compare it with the true value. Typically, a value like the root mean square error is used for this, given by\n",
    "\n",
    "$$\n",
    "{\\displaystyle \\operatorname {RMSE} ={\\sqrt {\\frac {\\sum _{t=0}^{T-1}({\\hat {y}}_{t+h}-y_{t+h})^{2}}{T}}}.}\n",
    "$$\n",
    "\n",
    "where $\\hat{y}_{t+h}$ is the predicted value, equivalent to the object `y_t_plus_h_os` above.\n",
    "\n",
    "### An Example Two Period Forecast\n",
    "\n",
    "Let's see an example of this with some numerical data. Our dataframe needs to have our features, $x$, and target, $y_{t+h}$, in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e389defb",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 10\n",
    "window_size = 5\n",
    "horizon = 3\n",
    "\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"x\": range(0, T + horizon),\n",
    "        \"y\": [0] * horizon\n",
    "        + list(range(0, window_size))\n",
    "        + list(range(T - 1, window_size - 1, -1)),\n",
    "    }\n",
    ")\n",
    "df[\"x\"] = df[\"x\"].astype(float)\n",
    "# Bring y_{t+h} in line with feature x so that models always take data from the same row.\n",
    "# Note that for a variable to be at t+h today, with h the horizon, we need to bring it *back* by h steps\n",
    "df[\"y at t+h\"] = df[\"y\"].shift(-horizon)\n",
    "# drop y so that we don't accidentally put it into model as a feature by mistake\n",
    "df = df.drop(\"y\", axis=1)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25cc68e9",
   "metadata": {},
   "source": [
    "As ever, it's good to look at your data so let's do a quick chart. Note that time and $x$ have the same values here, so this is also the time series of $y_{t+h}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a4244c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot.scatter(x=\"x\", y=\"y at t+h\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19fede0b",
   "metadata": {},
   "source": [
    "Okay, what could possibly go wrong! Let's now do our forecasting exercise. First, we train or fit the model. Note that *pandas* indexing is *inclusive*, so we run to `window_size-1` rather than `window_size` to pick up our 5 in-sample points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d2da857",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "df = df.dropna().copy()\n",
    "model = LinearRegression(fit_intercept=False).fit(\n",
    "    X=df.loc[: window_size - 1, [\"x\"]], y=df.loc[: window_size - 1, [\"y at t+h\"]]\n",
    ")\n",
    "# Make in-sample prediction. Note that the predictions of the model come out in\n",
    "# an array of form [[a], [b], ...] so we flatten this to [a, b, ...]. Also, we pad\n",
    "# the results we didn't predict with [None]\n",
    "df[\"y at t+h IS\"] = list(model.predict(df.loc[: window_size - 1, [\"x\"]]).flatten()) + [\n",
    "    None\n",
    "] * (window_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e37040df",
   "metadata": {},
   "source": [
    "Make out of sample prediction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f4628aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"y at t+h OS\"] = [None] * (window_size) + list(\n",
    "    model.predict(df.loc[window_size:, [\"x\"]]).flatten()\n",
    ")\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e7f2944",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "df.plot.scatter(x=\"x\", y=\"y at t+h\", ax=ax, label=\"y at t+h\", s=80)\n",
    "df.plot.scatter(\n",
    "    x=\"x\",\n",
    "    y=\"y at t+h OS\",\n",
    "    color=\"red\",\n",
    "    ax=ax,\n",
    "    label=\"y at t+h OS\",\n",
    "    s=100,\n",
    ")\n",
    "ax.set_ylabel(r\"$y_{t+h}$ and $\\hat{y}_{t+h}^{\\operatorname{OS}}$\")\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40da0fd8",
   "metadata": {},
   "source": [
    "Oh dear, it doesn't look like our out-of-sample prediction did very well! Let's see what the RMSE error is for both the in-sample and out-of-sample prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8eeac13",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "\n",
    "def get_rmse(df, fcast_name, target_name):\n",
    "    # Get the index of the fcast\n",
    "    index = df[fcast_name].dropna().index\n",
    "    # only compare to values of target for which we have predictions:\n",
    "    return np.sqrt(\n",
    "        mean_squared_error(df.loc[index, fcast_name], df.loc[index, target_name])\n",
    "    )\n",
    "\n",
    "\n",
    "is_rmse = get_rmse(df, \"y at t+h IS\", \"y at t+h\")\n",
    "os_rmse = get_rmse(df, \"y at t+h OS\", \"y at t+h\")\n",
    "\n",
    "print(f\"The rmse for in-sample best predictions is {is_rmse:.2f}\")\n",
    "print(f\"The rmse for out-of-sample best predictions is {os_rmse:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c01368d",
   "metadata": {},
   "source": [
    "In this case - unsurprisingly given the big break in the trend of the time series - even though in-sample forecast is perfect, the *out-of-sample* forecast is terrible! If you ever have a situation where the OS RMSE is *lower* than the IS RMSE, it should ring alarm bells: it suggests something has gone wrong in the code."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "177bc4af",
   "metadata": {},
   "source": [
    "### Adding a Lag\n",
    "\n",
    "Quite frequently, you'll want to include lags in your forecast specification. Let's see how to do that in this simple, two-period forecast example. \n",
    "\n",
    "Just as forecasting a variable at $t+h$ means performing a `.shift(-h)` in order to bring the data together in the same row of a dataframe, so creating a feature (or regressor) that is a lag of another column will involve shifting the *other* way in time. So we're going to use `.shift(lag)`. Let's see it with some data:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9d8fa3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 10\n",
    "window_size = 5\n",
    "horizon = 3\n",
    "\n",
    "\n",
    "df = pd.DataFrame(\n",
    "    {\n",
    "        \"x\": range(-1, T + horizon),\n",
    "        \"y\": list(range(-1 - horizon, window_size))\n",
    "        + list(range(T - 1, window_size - 1, -1)),\n",
    "    }\n",
    ")\n",
    "df[\"x\"] = df[\"x\"].astype(float)\n",
    "# Bring y_{t+h} in line with feature x so that models always take data from the same row.\n",
    "# Note that for a variable to be at t+h today, with h the horizon, we need to bring it *back* by h steps\n",
    "df[\"y at t+h\"] = df[\"y\"].shift(-horizon)\n",
    "\n",
    "# New part: create a single lag of y\n",
    "df[\"lag_y\"] = df[\"y\"].shift(1)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd015aec",
   "metadata": {},
   "source": [
    "Now you can see that there's an extra row we can't use (because we don't have the data for a lag of y from that time). As before, in practice, we'd drop y so that we didn't use it by mistake and also drop the top and tail rows that have nans in. We also reset the index to ensure it starts at 0 (by dropping the first row, we would otherwise lose the first index of zero):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28a2058f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.drop(\"y\", axis=1).dropna().reset_index(drop=True)\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ac0c4c",
   "metadata": {},
   "source": [
    "## Expanding Window\n",
    "\n",
    "In an expanding window forecast exercise, multiple forecasts are made. In each, the in-sample period grows over time. Necessarily (given finite data), the out-of-sample period shrinks over time. The overall out-of-sample forecast is given by the unique union of the out-of-sample forecasts that are trained on the most information. To make this clear, let's first define it mathematically and then with a diagram.\n",
    "\n",
    "We index the different forecasts (synonymous with different models $f_\\mu$) by $\\mu$. The starting size of the window will be $s + \\alpha$ where $s$ is the step size and $\\alpha$ is a parameter that adjusts window size. For arbitrary $z_t$, the $\\mu$th in-sample slice of time series data is: \n",
    "\n",
    "$$\n",
    "I^e_{\\mu}(\\vec{z}) = \\left\\{z_t  \\right\\}_{t=0}^{t=\\mu\\cdot s + \\alpha -1}\n",
    "$$\n",
    "\n",
    "Applied to the features, this is the slice that will be used to train (aka fit) model $f_\\mu$.\n",
    "\n",
    "Likewise, the out-of-sample set is:\n",
    "\n",
    "$$\n",
    "O^e_{\\mu}(\\vec{z}) = \\left\\{z_t \\right\\}^{t=T-1}_{t=\\mu\\cdot s + \\alpha}\n",
    "$$\n",
    "\n",
    "We can visualise this like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccddc831",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "s = 1\n",
    "alpha = 2\n",
    "T = 10\n",
    "max_mu = int((T - 1 + 1 - alpha) / s)\n",
    "\n",
    "multi_period = pd.concat(\n",
    "    (\n",
    "        [in_sample_block_e(mu, s, alpha, T) for mu in range(1, max_mu)]\n",
    "        + [os_sample_block_e(mu, s, alpha, T) for mu in range(1, max_mu)]\n",
    "    ),\n",
    "    axis=0,\n",
    ")\n",
    "multi_period[\"sample\"] = multi_period[\"sample\"].astype(\"category\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e0a7436",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(multi_period, aes(x=\"t\", y=as_discrete(\"μ\"), fill=as_discrete(\"sample\")))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T + 1)), format=\"d\")\n",
    "    + scale_y_continuous(format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Accent\")\n",
    "    + theme(axis_text=element_text(angle=0))\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b309b123",
   "metadata": {},
   "source": [
    "Note that in-sample and out-of-sample never overlap in time for the same value of $\\mu$. This is important because it helps keep our out-of-sample period out-of-sample. We do not want a situation where we *think* we are performing an out-of-sample forecast but actually the model was trained on some of the supposedly out-of-sample data.\n",
    "\n",
    "### The Likely Best Possible Expanding Window Forecast\n",
    "\n",
    "When we are assessing the potential of a model, we want to know i) how it would perform in practice and ii) that we are testing the best version of it possible. To address both of these points, we need to make a choice from the possible out-of-sample predictions because we have several for the same value of $t$:\n",
    "\n",
    "- $|\\mu|$ out-of-sample prediction models for $t=T-1$\n",
    "- $|\\mu|-1$ out-of-sample prediction models for $t=T-2$\n",
    "- and so on...\n",
    "\n",
    "To choose which of our out-of-sample predictions when we have multiple ones to choose from, we create the *union* of step ahead out-of-sample forecasts. Step ahead because that's what will likely be the best forecast (as opposed to many steps ahead) while still being out-of-sample. Because there are no *guarantees* which out-of-sample forecast will be best, it might be more accurate to call this the *preferred* or *one-step ahead* out-of-sample expanding forecast. For each possible value of $t$, this is the out-of-sample prediction with the highest possible value of $\\mu$:\n",
    "\n",
    "$$\n",
    "\\mathcal{O}^{e} = \\bigcup_\\mu  \\bigg\\{ f_\\mu\\left(O_\\mu^e(X)\\right) \\bigg\\}_{t=\\mu s+ \\alpha }^{t=(\\mu+1)s -1+ \\alpha}\n",
    "$$\n",
    "\n",
    "in the diagram above, it's the longest diagonal of out-of-sample predictions. Equivalently, the out-of-sample prediction set is composed of the first step of each test set indexed by $\\mu$. In each case, these out-of-sample predictions are likely to be better than ones for other values of $\\mu$ because they are trained on models that are closer in time to the value that is being predicted.\n",
    "\n",
    "Let's work this out using code too. First, we define a function that, given a data frame with a time index as defined earlier, can return the out-of-sample rows. We know this ranges from $\\mu \\cdot s + \\alpha$ to $T-1$ (inclusive), so will become `range(mu * s + alpha, T, step_size)` inside the function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5d7fb07",
   "metadata": {},
   "outputs": [],
   "source": [
    "step_size = 1\n",
    "alpha = 2\n",
    "window_size = alpha + step_size\n",
    "T = 10\n",
    "max_mu = int((T - 1 + 1 - alpha) / step_size)\n",
    "\n",
    "\n",
    "def os_sample_block_e(mu, step_size, alpha, T):\n",
    "    \"\"\"Returns the times within the $mu$th block of an expanding window (out-of-sample)\"\"\"\n",
    "    xf = pd.DataFrame({\"t\": np.arange(mu * s + alpha, T, step_size)})\n",
    "    xf[\"μ\"] = [mu] * (T - (mu * s + alpha))\n",
    "    xf[\"sample\"] = \"out-of-sample\"\n",
    "    return xf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58c9ea0d",
   "metadata": {},
   "source": [
    "Next, we create the union of best predictions, which will always be the last time step in ascending order of $\\mu$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f569cfc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# iterate over increasing values of mu\n",
    "script_O_e = pd.concat(\n",
    "    [os_sample_block_e(mu, s, alpha, T) for mu in range(1, max_mu)], axis=0\n",
    ")\n",
    "script_O_e = script_O_e.reset_index(drop=True)\n",
    "script_O_e[\"Preferred Prediction\"] = False\n",
    "# Drop any duplicates of timesteps, keeping only the last entry (which will have highest mu)\n",
    "script_O_e.loc[\n",
    "    script_O_e[~script_O_e[\"t\"].duplicated(keep=\"last\")].index, \"Preferred Prediction\"\n",
    "] = True\n",
    "script_O_e[\"μ\"] = script_O_e[\"μ\"].astype(\"int\")\n",
    "script_O_e"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8ff900b",
   "metadata": {},
   "source": [
    "Now we can plot the pattern of predictions we'd like to use as a function of $\\mu$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99b543ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(script_O_e, aes(x=\"t\", y=as_discrete(\"μ\"), fill=\"Preferred Prediction\"))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T)), limits=(-1, T), format=\"d\")\n",
    "    + scale_y_continuous(breaks=list(range(1, max_mu)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Pastel1\")\n",
    "    + ggtitle(\"Union of Expanding Window One-Step \\nAhead Out-of-Sample Predictions\")\n",
    "    + ggsize(700, 300)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7854f9eb",
   "metadata": {},
   "source": [
    "The series in red above is the one that we would use to compare against real data out-of-sample: it is the best *set* of predictions we are likely to make out-of-sample.\n",
    "\n",
    "Of course, it's easy to create a set of in-sample preferred predictions by analogy. These are less interesting but follow the same principle of taking the time steps with the highest numbered value of $\\mu$ available. In this special case of an expanding window, that's just going to be all of the predictions in the final $\\mu$. The expression is\n",
    "\n",
    "$$\n",
    "\\mathcal{I}^{e} = \\bigcup_\\mu \\bigg\\{ f_\\mu\\left(I_\\mu^e(X)\\right)\\bigg\\}_{t=(\\mu -1) s + \\alpha}^{t=\\mu s -1+\\alpha}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9790cd5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def in_sample_block_e(mu, step_size, alpha, T):\n",
    "    \"\"\"Returns the times within the $mu$th block of an expanding window (in-sample)\"\"\"\n",
    "    xf = pd.DataFrame({\"t\": np.arange(0, mu * s + alpha, step_size)})\n",
    "    xf[\"μ\"] = [mu] * (mu * s + alpha - 0)\n",
    "    xf[\"sample\"] = \"in-sample\"\n",
    "    return xf\n",
    "\n",
    "\n",
    "script_I_e = pd.concat(\n",
    "    [in_sample_block_e(mu, s, alpha, T) for mu in range(1, max_mu)], axis=0\n",
    ")\n",
    "script_I_e = script_I_e.reset_index()\n",
    "script_I_e[\"Preferred Prediction\"] = False\n",
    "script_I_e.loc[\n",
    "    script_I_e[~script_I_e[\"t\"].duplicated(keep=\"last\")].index, \"Preferred Prediction\"\n",
    "] = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "738619e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the union of expanding window in-sample predictions:\n",
    "(\n",
    "    ggplot(script_I_e, aes(x=\"t\", y=as_discrete(\"μ\"), fill=\"Preferred Prediction\"))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T)), limits=(-1, T), format=\"d\")\n",
    "    + scale_y_continuous(breaks=list(range(1, max_mu)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Pastel1\")\n",
    "    + ggtitle(\"Union of Expanding Window In-Sample Preferred Predictions\")\n",
    "    + ggsize(700, 300)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2a40111",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# This isn't used\n",
    "# For completeness, this is what the IS and OS time series look like together\n",
    "\n",
    "script_I_e[\"sample\"] = \"In-Sample\"\n",
    "script_O_e[\"sample\"] = \"Out-of-Sample\"\n",
    "script_rolling = pd.concat([script_I_e, script_O_e], axis=0)\n",
    "script_rolling = script_rolling[script_rolling[\"Preferred Prediction\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a48e6d8c",
   "metadata": {},
   "source": [
    "## Rolling Window\n",
    "\n",
    "In a rolling window forecast exercise, multiple forecasts are made. In each, the in-sample period *moves* (rather than grows) over time. Necessarily (given finite data), the out-of-sample period shrinks over time as the exercise progresses through different window configurations.\n",
    "\n",
    "We index the different forecasts (synonymous with different models $f_\\mu$) by $\\mu$. The starting size of the window will be $s + \\alpha$ where $s$ is the step size and $\\alpha$ is a parameter that adjusts window size. For arbitrary $z_t$, the $\\mu$th in-sample slice of time series data is: \n",
    "\n",
    "$$\n",
    "I^r_{\\mu}(\\vec{z}) = \\bigg\\{ z_{t} \\bigg\\}_{t=(\\mu-1)\\cdot s}^{t=\\mu\\cdot s+\\alpha -1}\n",
    "$$\n",
    "\n",
    "Applied to the features, this is the slice that will be used to train (aka fit) model $f_\\mu$.\n",
    "\n",
    "Likewise, the out-of-sample set is:\n",
    "\n",
    "$$\n",
    "O^r_{\\mu}(\\vec{z}) = \\bigg\\{ z_{t} \\bigg\\}_{t=\\mu\\cdot s+\\alpha}^{t=T} \n",
    "$$\n",
    "\n",
    "Below is a visualisation of all the rolling window forecasts in a single forecasting exercise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e63e2e6f",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "s = 1\n",
    "alpha = 2\n",
    "window_size = s + alpha\n",
    "T = 10\n",
    "max_mu = int((T - 1 + 1 - alpha) / s)\n",
    "\n",
    "multi_period_rolling = pd.concat(\n",
    "    (\n",
    "        [in_sample_block_r(mu, s, alpha, T) for mu in range(1, max_mu)]\n",
    "        + [os_sample_block_r(mu, s, alpha, T) for mu in range(1, max_mu)]\n",
    "    ),\n",
    "    axis=0,\n",
    ")\n",
    "multi_period_rolling[\"sample\"] = multi_period_rolling[\"sample\"].astype(\"category\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8080ab90",
   "metadata": {
    "tags": [
     "remove-input"
    ]
   },
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(multi_period_rolling, aes(x=\"t\", y=as_discrete(\"μ\"), fill=\"sample\"))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T)), limits=(-1, T), format=\"d\")\n",
    "    + scale_y_continuous(breaks=list(range(1, max_mu)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Pastel1\")\n",
    "    + ggtitle(\"Rolling Window\")\n",
    "    + ggsize(700, 300)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d87109b6",
   "metadata": {},
   "source": [
    "### The Likely Best Possible Rolling Window Forecast\n",
    "\n",
    "As with the expanding window, we want to know how models will perform in practice and that the best version of it is being assessed. To address both of these points, we need to make a choice from the possible out-of-sample predictions because we have several for the same value of $t$:\n",
    "\n",
    "- $|\\mu|$ out-of-sample prediction models for $t=T-1$\n",
    "- $|\\mu|-1$ out-of-sample prediction models for $t=T-2$\n",
    "- and so on...\n",
    "\n",
    "To do this we create the *union* of step ahead out-of-sample forecasts. Step ahead because that's what will likely be the best forecast (as opposed to many steps ahead) while still being out-of-sample. Because there are no *guarantees* which out-of-sample forecast will be best, it might be more accurate to call this the *preferred* or *one-step ahead* out-of-sample rolling window forecast. For each possible value of $t$, this is the out-of-sample prediction with the highest possible value of $\\mu$:\n",
    "\n",
    "$$\n",
    "\\mathcal{O}^{r} = \\bigcup_\\mu  \\bigg\\{ f_\\mu\\left(O_\\mu^r(X)\\right) \\bigg\\}_{t=\\mu s+ \\alpha }^{t=(\\mu+1)s -1+ \\alpha}\n",
    "$$\n",
    "\n",
    "in the diagram above, it's the longest diagonal of out-of-sample predictions. Equivalently, the out-of-sample prediction set is composed of the first step of each test set indexed by $\\mu$. In each case, these out-of-sample predictions are likely to be better than ones for other values of $\\mu$ because they are trained on models that are closer in time to the value that is being predicted.\n",
    "\n",
    "Let's work this out using code too. First, we define function thats, given a data frame with a time index as defined earlier, can return the in-sample or out-of-sample rows. For out-of-sample, the slices range from $\\mu\\cdot s + \\alpha$ to $T-1$ (inclusive), so will become `range(mu * s + alpha, T, step_size)` inside the function. Likewise, the inclusive/exclusive distinction applies to the in-sample slice, which will have `range((mu - 1) * s, mu * s + alpha, s)` for $(\\mu-1)\\cdot s$ to $\\mu\\cdot s+\\alpha -1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d51c2eaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def in_sample_block_r(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange((mu - 1) * s, mu * s + alpha, s)})\n",
    "    xf[\"μ\"] = [mu] * (mu * s + alpha - (mu - 1) * s)\n",
    "    xf[\"sample\"] = \"in-sample\"\n",
    "    return xf\n",
    "\n",
    "\n",
    "def os_sample_block_r(mu, s, alpha, T):\n",
    "    xf = pd.DataFrame({\"t\": np.arange(mu * s + alpha, T, s)})\n",
    "    xf[\"μ\"] = [mu] * (T - (mu * s + alpha))\n",
    "    xf[\"sample\"] = \"out-of-sample\"\n",
    "    return xf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7395f86b",
   "metadata": {},
   "source": [
    "Having defined these, we can compute the unions of best in- and out-of-sample predictions. For the expanding window case, we ended up repeating some of the same code. This is bad practice. We can do better by putting all of the common elements into a single function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdb4c5c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Config:\n",
    "s = 1\n",
    "alpha = 2\n",
    "window_size = s + alpha\n",
    "T = 10\n",
    "\n",
    "\n",
    "def return_rolling_union_of_slices(s, alpha, T, out_of_sample):\n",
    "    \"\"\"\n",
    "    Returns a dataframe explaining where the best in-sample or out-of-sample forecasts are\n",
    "    for a a given stepsize (s), alpha (=windowsize-s), and total number of time steps of\n",
    "    the data (T).\n",
    "    \"\"\"\n",
    "    if out_of_sample:\n",
    "        slice_function = os_sample_block_r\n",
    "        keep_string = \"last\"\n",
    "    elif not out_of_sample:\n",
    "        slice_function = in_sample_block_r\n",
    "        keep_string = \"first\"\n",
    "    else:\n",
    "        raise ValueError(\"out_of_sample must be True or False (and of type bool).\")\n",
    "\n",
    "    max_mu = int((T - 1 + 1 - alpha) / s)\n",
    "    union = pd.concat(\n",
    "        [slice_function(mu, s, alpha, T) for mu in range(1, max_mu)], axis=0\n",
    "    )\n",
    "    union = union.reset_index(drop=True)\n",
    "    union[\"Preferred Prediction\"] = False\n",
    "    union.loc[\n",
    "        union[~union[\"t\"].duplicated(keep=keep_string)].index, \"Preferred Prediction\"\n",
    "    ] = True\n",
    "    return union\n",
    "\n",
    "\n",
    "script_I_r = return_rolling_union_of_slices(s, alpha, T, out_of_sample=False)\n",
    "script_O_r = return_rolling_union_of_slices(s, alpha, T, out_of_sample=True)\n",
    "script_O_r.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fd8e953",
   "metadata": {},
   "source": [
    "Let's plot the union of best rolling out-of-sample predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74dd9f4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(script_O_r, aes(x=\"t\", y=as_discrete(\"μ\"), fill=\"Preferred Prediction\"))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T)), limits=(-1, T), format=\"d\")\n",
    "    + scale_y_continuous(breaks=list(range(1, max_mu)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Pastel1\")\n",
    "    + ggtitle(\"Union of Rolling Window One Step \\nAhead Out-of-Sample Predictions\")\n",
    "    + ggsize(700, 300)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeea8b8d",
   "metadata": {},
   "source": [
    "And the equivalent for in-sample predictions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77f0c99c",
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(script_I_r, aes(x=\"t\", y=as_discrete(\"μ\"), fill=\"Preferred Prediction\"))\n",
    "    + geom_tile(**width_height_tiles)\n",
    "    + scale_x_continuous(breaks=list(range(T)), limits=(-1, T), format=\"d\")\n",
    "    + scale_y_continuous(breaks=list(range(1, max_mu)), format=\"d\")\n",
    "    + scale_fill_brewer(type=\"qual\", palette=\"Pastel1\")\n",
    "    + ggtitle(\"Union of Rolling Window In-Sample Preferred Predictions\")\n",
    "    + ggsize(700, 300)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6bd3a4a",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Not used\n",
    "# For completeness, this is what the IS and OS time series look like next together\n",
    "\n",
    "script_I_r[\"sample\"] = \"In-Sample\"\n",
    "script_O_r[\"sample\"] = \"Out-of-Sample\"\n",
    "script_rolling = pd.concat([script_I_r, script_O_r], axis=0)\n",
    "script_rolling = script_rolling[script_rolling[\"Preferred Prediction\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d9d64b3",
   "metadata": {},
   "source": [
    "## Transforms of Series in Forecasting Exercises\n",
    "\n",
    "In practice, you often *don't* want to throw a time series straight into a forecast. This is for lots of reasons: stationarity, producing the best performing model, the ranges that certain models will accept, and so on. Instead you want to transform the data first using a general transform $T\\colon X\\to Z$. But when in a forecasting situation, extra care must be taken not to let information leak from the period that a model is *not* supposed to know about. For example, a forecast that is supposed to be out-of-sample, ie it is forecasting a time period $t''$, should not be able to see any information about $X_{t''}$ or $y_{t''}$. Why is this a problem when doing transforms? Because many transforms of, say, a scalar $x_t$ depend on $x_{t'}; \\: t\\neq t'$.\n",
    "\n",
    "First, let's take an example of a common transform with no dependence on any other data: $T(x_t) = \\ln(x_t)$. In this case we do not need to worry.\n",
    "\n",
    "Another example that's quite common is $T(x_t) = \\frac{x_t}{x_{t-1}}$. Again, this transform shouldn't cause us much concern because it only depends on *past* information.\n",
    "\n",
    "Let's now turn to more tricky transforms. Let's imagine we wish to normalise the data, ie take\n",
    "\n",
    "$$\n",
    "T(x_t) = \\frac{1}{\\sigma_{x}} \\left( x_t - \\frac{\\sum_{t'} x_{t'}}{\\sum_{t'} t'} \\right)\n",
    "$$\n",
    "\n",
    "where ${\\displaystyle \\sigma_{x}={\\sqrt {{\\frac {1}{\\sum_{t'} x_t'}} \\sum_{t'}\\left(x_{t'}-{\\frac{\\sum_{t''} x_{t''}}{\\sum_{t''} t''}}\\right)^{2}}}}$. The summations are over *all* values of $t$, but this means that $T(x_t)$ has 'seen' $x_{t'}$, $t'>t$ and $t'$ may well be in the out-of-sample slice. This could be a problem: imagine that the mean or standard deviation changes significantly because of a crisis. Being privy to that shift may make it easier to predict post-crisis values of the target variable.\n",
    "\n",
    "To fix information leakage due to transforms, we redefine transformations to make them a function of the windows being used in the forecasting exercise. For an expanding window transform of a variable $x_t$, this would be:\n",
    "\n",
    "$$\n",
    "T_{\\mu}^e (\\vec{x}) = T_{\\mu}^e (\\vec{x}; I_\\mu^e(\\vec{x})) = \\frac{\\vec{x} - \\langle{I_\\mu^e(\\vec{x})\\rangle}}{\\sigma_{I_\\mu^e(\\vec{x})}} \n",
    "$$\n",
    "\n",
    "where $\\langle x \\rangle$ is the mean of $\\{x_t\\}$. So now the transform is dependent on the *type* of forecasting exercise, the *period* ($\\mu$) of the forecasting exercise, and the *input data*. A similar definition applies to the rolling window configuration.\n",
    "\n",
    "How do we apply this in code? That part's easy: we can re-use the slice functions we made earlier, namely `in_sample_block_r` for in-sample rolling windows and `in_sample_block_e` for out-of-sample expanding windows. A more general way of avoiding data leakage, especially when there are many transformations, is to use **scikit-learn**'s [pipeline](https://scikit-learn.org/stable/modules/compose.html#pipeline) functionality.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4766b5",
   "metadata": {},
   "source": [
    "## Data Vintage and the Real-Time Information Flow\n",
    "\n",
    "There's another complication if you want an accurate assessment of how good your forecast will be in the real world: some data, and especially economic data, are revised over time. But to give an *accurate* picture of how well your forecast approach would have done at the time, you need to use only the data that would have been available at the time the forecast was made: you need to use real-time data vintages. Typically, one only considers different vintages of regressors--but target variables change too and this is worth bearing in mind if you are using lags.\n",
    "\n",
    "A few labelling operations help us use everything defined so far in this chapter in the real-time data vintage case. We introduce an extra label for vintage, $\\tau$, which says when the data were released. Features are now given by $x_{tk\\tau}$. Assuming that vintages are improving over time, models will always be trained using the most up-to-date vintage available up to that time so that a forecast made at $t$ uses the highest available value of $\\tau$ for which $\\tau \\leq t$, eg $\\tau = \\max{ \\{ \\tau' | \\tau' \\leq t \\}}$."
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "formats": "md:myst"
  },
  "kernelspec": {
   "display_name": "codeforecon",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}