time-series.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(time-series)=\n",
    "# Time Series\n",
    "\n",
    "In this chapter, we'll look at time series. If you haven't yet looked at the two sections on **pandas**, the [Data Quickstart](data-quickstart) and [Working with Data](working-with-data) chapters, it might be worth taking a quick spin through them first. You may also find it useful to be familiar with a few of the concepts from the previous Chapter, {ref}`time-intro`.\n",
    "\n",
    "While we'll cover the basics here, the full set of time series functionality of **pandas** can be [found here](https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html).\n",
    "\n",
    "This chapter has benefitted from the [Python Data Science Handbook](https://jakevdp.github.io/PythonDataScienceHandbook/) by Jake VanderPlas, Tom Augspurger's [Effective Pandas](https://github.com/TomAugspurger/effective-pandas), [Applied Time Series for Fisheries and Environmental Sciences](https://atsa-es.github.io/), and the documentation of the [**arch**](https://arch.readthedocs.io/) package.\n",
    "\n",
    "Let's imports a few of the packages we'll need first. You may need to install some of these; the Chapter on {ref}`code-preliminaries` covers how to install new packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "from itertools import cycle\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "from rich import inspect\n",
    "from statsmodels.graphics import tsaplots\n",
    "\n",
    "mpl.rcParams.update({\"lines.linewidth\": 1.2})\n",
    "# Set max rows displayed for readability\n",
    "pd.set_option(\"display.max_rows\", 8)\n",
    "\n",
    "# Set seed for random numbers\n",
    "seed_for_prng = 78557\n",
    "prng = np.random.default_rng(\n",
    "    seed_for_prng\n",
    ")  # prng=probabilistic random number generator\n",
    "# Turn off warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "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",
    "matplotlib_inline.backend_inline.set_matplotlib_formats(\"svg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introducing Time Series with **pandas**\n",
    "\n",
    "[**pandas**](https://pandas.pydata.org/) is the workhorse of time series analysis in Python. The basic object is a *timestamp*. The `pd.to_datetime()` function creates timestamps from strings that could reasonably represent datetimes. Let's see an example of using `pd.to_datetime()` to create a timestamp and then inspect all of the methods and attributes of the created timestamp using **rich**'s `inspect()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date = pd.to_datetime(\"16th of February, 2020\")\n",
    "inspect(date)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is of type `Timestamp` and you can see that it has many of the same properties as the built-in Python `datetime.datetime` class from the previous chapter. As with that, the default setting for `tz` (timezone) and `tzinfo` is `None`. There are some extra properties, though, such as `freq` for frequency, which will be very useful when it comes to manipulating time *series* as opposed to just one or two datetimes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating and Using Time Series\n",
    "\n",
    "There are two main scenarios in which you might be creating time series using **pandas**: i) creating one from scratch or ii) reading in data from a file. Let's look at a few ways to do i) first. \n",
    "\n",
    "You can create a time series with **pandas** by taking a date as created above and extending it using **pandas** timedelta function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "date + pd.to_timedelta(np.arange(12), \"D\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This has created a datetime index of type `datetime65[ns]` (remember, an index is a special type of **pandas** column), where \"ns\" stands for nano-second resolution.\n",
    "\n",
    "Another method is to create a range of dates (pass a frequency using the `freq=` keyword argument):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.date_range(start=\"2018/1/1\", end=\"2018/1/8\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to create ranges is to specify the number of periods and the frequency:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.date_range(\"2018-01-01\", periods=3, freq=\"H\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the discussion of the previous chapter on timezones, you can also localise timezones directly in **pandas** dataframes:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dti = pd.date_range(\"2018-01-01\", periods=3, freq=\"H\").tz_localize(\"UTC\")\n",
    "dti.tz_convert(\"US/Pacific\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's see how to turn data that has been read in with a non-datetime type into a vector of datetimes. This happens *all the time* in practice. We'll read in some data on job vacancies for information and communication jobs, ONS code UNEM-JP9P, and then try to wrangle the given \"date\" column into a **pandas** datetime column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "\n",
    "url = \"https://api.beta.ons.gov.uk/data?uri=/employmentandlabourmarket/peopleinwork/employmentandemployeetypes/timeseries/jp9p/lms\"\n",
    "\n",
    "# Get the data from the ONS API:\n",
    "df = pd.DataFrame(pd.json_normalize(requests.get(url).json()[\"months\"]))\n",
    "df[\"value\"] = pd.to_numeric(df[\"value\"])\n",
    "df = df[[\"date\", \"value\"]]\n",
    "df = df.rename(columns={\"value\": \"Vacancies (ICT), thousands\"})\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have the data in. Let's look at the column types that arrived."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the default 'object' type, but we want the date column to have `datetime64[ns]`, which is a datetime type. Again, we use `pd.to_datetime()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"date\"] = pd.to_datetime(df[\"date\"])\n",
    "df[\"date\"].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the conversion from the format of data that was put in of \"2001 MAY\" to datetime worked out-of-the-box. `pd.to_datetime()` will always take an educated guess as to the format, but it won't always work out.\n",
    "\n",
    "What happens if we have a more tricky-to-read-in datetime column? This frequently occurs in practice so it's well worth exploring an example. Let's create some random data with dates in an unusual format with month first, then year, then day, eg \"1, '19, 29\" and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "small_df = pd.DataFrame({\"date\": [\"1, '19, 22\", \"1, '19, 23\"], \"values\": [\"1\", \"2\"]})\n",
    "small_df[\"date\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if we were to run this via `pd.to_datetime()` with no further input, it would misinterpret, for example, the first date as `2022-01-19`. So we must pass a bit more info to `pd.to_datetime()` to help it out. We can pass a `format=` keyword argument with the format that the datetime takes. Here, we'll use `%m` for month in number format, `%y` for year in 2-digit format, and `%d` for 2-digit day. We can also add in the other characters such as `'` and `,`. You can find a list of datetime format identifiers in the previous chapter or over at [https://strftime.org/](https://strftime.org/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.to_datetime(small_df[\"date\"], format=\"%m, '%y, %d\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Datetime Offsets\n",
    "\n",
    "Our data, currently held in `df`, were read in as if they were from the *start* of the month but these data refer to the month that has passed and so should be for the *end* of the month. Fortunately, we can change this using a time offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"date\"] = df[\"date\"] + pd.offsets.MonthEnd()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While we used the `MonthEnd` offset here, there are many different offsets available. You can find a [full table of date offsets here](https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#dateoffset-objects)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The `.dt` accessor\n",
    "\n",
    "When you have a datetime column, you can use the `.dt` accessor to grab lots of useful information from it such as the `minute`, `month`, and so on. Some that are functions, rather than just accessors of underlying properties, are followed by brackets, `()`, because they are functions. Here are a few useful examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rich import print\n",
    "\n",
    "print(\"Using `dt.day_name()`\")\n",
    "print(df[\"date\"].dt.day_name().head())\n",
    "print(\"Using `dt.isocalendar()`\")\n",
    "print(df[\"date\"].dt.isocalendar().head())\n",
    "print(\"Using `dt.month`\")\n",
    "print(df[\"date\"].dt.month.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating a datetime index and setting the frequency\n",
    "\n",
    "For the subsequent parts, we'll set the datetime column to be the index of the dataframe. *This is the standard setup you will likely want to use when dealing with time series.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.set_index(\"date\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if we look at the first few entries of the index of dataframe (a datetime index) using `head()` as above, we'll see that the `freq=` parameter is set as `None`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.index[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This can be set for the whole dataframe using the `asfreq` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.asfreq(\"M\")\n",
    "df.index[:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Although most of the time it doesn't matter about the fact that `freq=None`, some aggregation operations need to know the frequency of the time series in order to work and it's good practice to set it if your data *are* regular. You can use `asfreq` to go from a higher frequency to a lower frequency too: the last entry from the higher frequency that aligns with the lower frequency will be taken, for example in going from months to years, December's value would be used.\n",
    "\n",
    "Note that trying to set the frequency when your datetime index doesn't match up to a particular frequency will cause errors or problems.\n",
    "\n",
    "A few useful frequencies to know about are in the table below; all of these can be used with `pd.to_datetime()` too.\n",
    "\n",
    "| Code  | Represents                                                          |\n",
    "|-------|---------------------------------------------------------------------|\n",
    "| D     | Calendar day                                                        |\n",
    "| W     | Weekly                                                              |\n",
    "| M     | Month end                                                           |\n",
    "| Q     | Quarter end                                                         |\n",
    "| A     | Year end                                                            |\n",
    "| H     | Hours                                                               |\n",
    "| T     | Minutes                                                             |\n",
    "| S     | Seconds                                                             |\n",
    "| B     | Business day                                                        |\n",
    "| BM    | Business month end                                                  |\n",
    "| BQ    | Business quarter end                                                |\n",
    "| BA    | Business year end                                                   |\n",
    "| BH    | Business hours                                                      |\n",
    "| MS    | Month start                                                         |\n",
    "| QS    | Quarter start                                                       |\n",
    "| W-SUN | Weeks beginning with Sunday (similar for other days)                |\n",
    "| 2M    | Every 2 months (works with other combinations of numbers and codes) |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making Quick Time Series Plots\n",
    "\n",
    "Having managed to put your time series into a dataframe, perhaps converting a column of type string into a column of type datetime in the process, you often just want to see the thing! We can achieve this using the `plot()` command, as long as we have a datetime index.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resampling, Rolling, and Shifting\n",
    "\n",
    "Now our data have a *datetime index*, some common time series operations are made very easy for us.\n",
    "\n",
    "### Resampling\n",
    "\n",
    "Quite frequently, there is a situation in which one would like to change the frequency of a given time series. A time index-based dataframe makes this easy via the `resample()` function. `resample()` must be told *how* you'd like to resample the data, for example via the mean or median. Here's an example resampling the monthly data to annual and taking the mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.resample(\"A\").mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As resample is just a special type of aggregation, it can work with all of the usual functions that aggregations do, including in-built functions or user-defined functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.resample(\"5A\").agg([\"mean\", \"std\"]).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resampling can go up in frequency (up-sampling) as well as down, but we no longer need to choose an aggregation function, we must now choose how we'd like to fill in the gaps for the frequencies we didn't have in the original data. In the example below, they are just left as NaNs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.resample(\"D\").asfreq()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Options to fill in missing time series data include using `bfill()` or `ffill()` to fill in the blanks based on the next or last available value, respectively, or `interpolate()` (note how only the first 3 NaNs are replaced using the `limit` keyword argument):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.resample(\"D\").interpolate(method=\"linear\", limit_direction=\"forward\", limit=3)[:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the differences between the filling methods more clearly in this stock market data, following a chart by Jake Vanderplas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get stock market data\n",
    "import pandas_datareader as web\n",
    "\n",
    "xf = web.DataReader(\"AAPL\", \"stooq\", start=\"2017-01-01\", end=\"2019-06-01\")\n",
    "xf = xf.sort_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "data = xf.iloc[:10, 3]\n",
    "\n",
    "\n",
    "colour_list = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "iter_cycle = cycle(colour_list)\n",
    "\n",
    "\n",
    "data.asfreq(\"D\").plot(\n",
    "    ax=ax, marker=\"o\", linestyle=\"None\", zorder=3, color=next(iter_cycle)\n",
    ")\n",
    "data.asfreq(\"D\", method=\"bfill\").plot(ax=ax, style=\"-.o\", lw=1, color=next(iter_cycle))\n",
    "data.asfreq(\"D\", method=\"ffill\").plot(ax=ax, style=\"--o\", lw=1, color=next(iter_cycle))\n",
    "ax.set_ylabel(\"Close ($)\")\n",
    "ax.legend([\"original\", \"back-fill\", \"forward-fill\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rolling Window Functions\n",
    "\n",
    "The `rolling()` and `ewm()` methods are both rolling window functions. The first includes functions of the sequence\n",
    "\n",
    "$$\n",
    "y_t = f(\\{x_{t-i} \\}_{i=0}^{i=R-1})\n",
    "$$\n",
    "\n",
    "where $R$ is the number of periods to use for the rolling window. For example, if the function is the mean, then $f$ takes the form $\\frac{1}{R}\\displaystyle\\sum_{i=0}^{i=R-1} x_{t-i}$.\n",
    "\n",
    "The example below is a 2-period rolling mean:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.rolling(2).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ewm()` includes the class of functions where data point $x_{t-i}$ has a weight $w_i = (1-\\alpha)^i$. As $0 < \\alpha < 1$, points further back in time are given less weight. For example, an exponentially moving average is given by\n",
    "\n",
    "$$\n",
    "y_t = \\frac{x_t + (1 - \\alpha)x_{t-1} + (1 - \\alpha)^2 x_{t-2} + ...\n",
    "+ (1 - \\alpha)^t x_{0}}{1 + (1 - \\alpha) + (1 - \\alpha)^2 + ...\n",
    "+ (1 - \\alpha)^t}\n",
    "$$\n",
    "\n",
    "The example below shows the code for the exponentially weighted moving average:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.ewm(alpha=0.2).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see these methods together on the stock market data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "roll_num = 28\n",
    "alpha = 0.03\n",
    "xf[\"Close\"].plot(label=\"Raw\", alpha=0.5)\n",
    "xf[\"Close\"].expanding().mean().plot(label=\"Expanding Average\", style=\":\")\n",
    "xf[\"Close\"].ewm(alpha=alpha).mean().plot(\n",
    "    label=f\"EWMA ($\\\\alpha=${alpha:.2f})\", style=\"--\"\n",
    ")\n",
    "xf[\"Close\"].rolling(roll_num).mean().plot(label=f\"{roll_num}D MA\", style=\"-.\")\n",
    "ax.legend()\n",
    "ax.set_ylabel(\"Close ($)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more tools to analyse stocks, see the [**Pandas TA**](https://twopirllc.github.io/pandas-ta/) package.\n",
    "\n",
    "We can also use `rolling()` as an intermediate step in creating more than one type of aggregation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "roll = xf[\"Close\"].rolling(50, center=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "m = roll.agg([\"mean\", \"std\"])\n",
    "m[\"mean\"].plot(ax=ax)\n",
    "ax.fill_between(m.index, m[\"mean\"] - m[\"std\"], m[\"mean\"] + m[\"std\"], alpha=0.2)\n",
    "ax.set_ylabel(\"Close ($)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shifting\n",
    "\n",
    "Shifting can move series around in time; it's what we need to create leads and lags of time series. Let's create a lead and a lag in the data. Remember that a lead is going to shift the pattern in the data to the left (ie earlier in time), while the lag is going to shift patterns later in time (ie to the right)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lead = 12\n",
    "lag = 3\n",
    "orig_series_name = df.columns[0]\n",
    "df[f\"lead ({lead} months)\"] = df[orig_series_name].shift(-lead)\n",
    "df[f\"lag ({lag} months)\"] = df[orig_series_name].shift(lag)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.iloc[100:300, :].plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correlations\n",
    "\n",
    "In this section, we'll look at some different ways of comparing time series using simple metrics.\n",
    "\n",
    "### Autocorrelation\n",
    "\n",
    "First, what about comparing a time series to itself? The autocorrelation of a time series tells you how persistent it is. The econometrics package [**statsmodels**]() has some tools for this, most notably `statsmodels.tsa.stattools.acf`. Sometimes what you want is just a visual cue though, in which case the code below produces a nice chart."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = tsaplots.plot_acf(df[\"Vacancies (ICT), thousands\"], lags=24)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rememebr that the autocorrelation can tell you about likely terms to include in a moving average (MA) model.\n",
    "\n",
    "We can also look at the *partial autocorrelation*, which, at a given lag, is the correlation that results after removing the effect of any correlations due to the terms at shorter lags. It can be computer via the `statsmodels.tsa.stattools.pacf` function, but here's a chart:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = tsaplots.plot_pacf(df[\"Vacancies (ICT), thousands\"], lags=24)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember that the partial autocorrelation function can tell you about likely terms to include in an autoregressive (AR) model.\n",
    "\n",
    "By looking at both the ACF and PACF, you can get an idea of what terms would be good to include in *both* an MA and AR model which, combined, make an ARMA model. It's also worth noting that, if the ACF and PACF do not seem to tail off, the series is likely to be non-stationary and differencing will be needed if you are going to model it with an ARMA approach."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Correlations with Other Series\n",
    "\n",
    "To compute correlations with other time series, a popular choice is the cross-correlation function (CCF). Confusingly, there are different definitions of this object in statistics and signal processing. In statistics, it's given by\n",
    "\n",
    "$$\n",
    "r_k = \\frac{g_k}{\\sqrt{\\sigma_x\\sigma_y}};\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "g_k = \\frac{1}{n}\\sum_{t=1}^{n-k} \\left(y_t-\\bar{y}\\right) \\left(x_{t+k}-\\bar{x}\\right),\n",
    "$$\n",
    "\n",
    "and $k$ is set by a rule of thumb, usually something like $10 \\log_{10}\\left(\\min(|x|, |y|)/2\\right)$. $g_k$ is the 'sliding' dot product of $y$, the response, with $x$, the predictor, up to the $k$ th lag.\n",
    "\n",
    "The default cross-correlation function in **statsmodels** is different. It only looks forward, and it computes *all* possible lags. You can call it via `sm.ccf` if you have already run `import statsmodels.api as sm`.\n",
    "\n",
    "We'll use the standard stats treatment here, but we'll need to define a function to do the computation we need. Then we'll apply the CCF to a classic example: the response of the number of Lynx trapped (in Canada) to sunspot activity.\n",
    "\n",
    "First, defining the function. We'll make use of the signal processing `ss.correlate()` function but modify it so that adopts the definition of CCF used in statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal as ss\n",
    "\n",
    "\n",
    "def ccf_stats(x, y, lag_max=None):\n",
    "    if not lag_max:\n",
    "        lag_max = int(10 * np.log10(min(len(x), len(y)) / 2))\n",
    "    result = ss.correlate(x - np.mean(x), y - np.mean(y), method=\"direct\") / (\n",
    "        np.std(y) * np.std(x) * len(y)\n",
    "    )\n",
    "    length = (len(result) - 1) // 2\n",
    "    lo = length - lag_max\n",
    "    hi = length + (lag_max + 1)\n",
    "    return result[lo:hi]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's grab the two different datasets and merge them, taking a quick look while we're at it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lynx_df = pd.read_csv(\n",
    "    \"https://vincentarelbundock.github.io/Rdatasets/csv/datasets/lynx.csv\", index_col=0\n",
    ")\n",
    "spot_df = pd.read_csv(\n",
    "    \"https://vincentarelbundock.github.io/Rdatasets/csv/datasets/sunspot.year.csv\",\n",
    "    index_col=0,\n",
    ")\n",
    "df = pd.merge(lynx_df, spot_df, on=\"time\", how=\"inner\", suffixes=(\"_lynx\", \"_sunspots\"))\n",
    "df = df.set_index(\"time\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.plot(subplots=True, ylim=(0, None));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use our defined function on the data. We'll add in a rough confidence interval of $\\pm \\frac{2}{\\sqrt{n}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ccf_res_stats = ccf_stats(df[\"value_sunspots\"], np.log(df[\"value_lynx\"]))\n",
    "conf_int = 2.0 / np.sqrt(len(df))\n",
    "ax.stem(\n",
    "    range(-int(len(ccf_res_stats) / 2), int(len(ccf_res_stats) / 2) + 1),\n",
    "    ccf_res_stats,\n",
    "    basefmt=\" \",\n",
    "    markerfmt=\" \",\n",
    ")\n",
    "ax.axhline(-conf_int, ls=\"--\", lw=0.9, color=\"k\")\n",
    "ax.axhline(conf_int, ls=\"--\", lw=0.9, color=\"k\")\n",
    "ax.axhline(0, color=\"k\", lw=1)\n",
    "ax.set_ylabel(\"Cross-correlation\")\n",
    "ax.set_xlabel(\"Lag\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The significant correlations at lags of -3 to -5 seem to suggest that lynx numbers fall in the 3 to 5 years after high sunspot activity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transformations\n",
    "\n",
    "There are many reasons why you might want to perform transformations on series, including to normalise data and to adjust for changes (such as inflation). In this section, we'll look at some of the most common transforms for time series.\n",
    "\n",
    "### Mathematical Transformations\n",
    "\n",
    "The usual transformations that we know and love from other series of numerical data can equally be applied to time series, for example the log and power transforms. Let's see some of these in action using data on US GDP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime\n",
    "\n",
    "import pandas_datareader.data as web\n",
    "\n",
    "ts_start_date = pd.to_datetime(\"2009-01-01\")\n",
    "\n",
    "df = web.DataReader(\"NA000334Q\", \"fred\", start=ts_start_date, end=datetime.now())\n",
    "df.columns = [\"gdp\"]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's create new columns for each of log, square, and square root transforms. All of these mathematical transformations can be found in **numpy**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"log\"] = np.log(df[\"gdp\"])\n",
    "df[\"square\"] = np.power(df[\"gdp\"], 2)\n",
    "df[\"sqrt\"] = np.sqrt(df[\"gdp\"])\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A very common transformation is to change to growth space. In time series maths, this is $\\frac{y_t}{y_{t-1}}$ for period-on-period growth or, say, $\\frac{y_t}{y_{t-4}}$ for growth relative to 4 periods ago.\n",
    "\n",
    "Let's add in quarter-on-quarter-a-year-ago and quarter-on-quarter growth rates to our growing list of transforms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"q_on_q\"] = df[\"gdp\"].pct_change(1) * 100\n",
    "df[\"q_on_q_yr_ago\"] = df[\"gdp\"].pct_change(4) * 100\n",
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another classic transform you might need is to de-mean your data. We are transforming such that\n",
    "\n",
    "$$\n",
    "y_t' = y_t - \\langle y_t \\rangle\n",
    "$$\n",
    "\n",
    "where the angular brackets denote a mean. Here's an example of that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"demean\"] = df[\"gdp\"] - df[\"gdp\"].mean()\n",
    "df[[\"gdp\", \"demean\"]].tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be wary though! The mean is a function that, by default, uses all time periods, so you need to be careful if you're working on a real-time data problem such as forecasting. In this case you only want to use the mean *up to that point in time* to demean a series. Mathematically, this is\n",
    "\n",
    "$$\n",
    "y_t' = y_t - \\frac{1}{t+1}\\displaystyle{\\sum_{\\tau=0}^{\\tau=t} y_\\tau}\n",
    "$$\n",
    "\n",
    "To do this, we can make use of the **pandas** *expanding* functionality we saw earlier. Here's an example of real-time demeaning:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"real_time_demean\"] = df[\"gdp\"] - df[\"gdp\"].expanding().mean()\n",
    "df[[\"gdp\", \"real_time_demean\"]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note how the second real-time demeaned entry is the second entry for GDP minus the average of the second and first entries of GDP. A more general way of avoiding data leakage is to use **scikit-learn**'s [pipeline](https://scikit-learn.org/stable/modules/compose.html#pipeline) functionality.\n",
    "\n",
    "Let's move on to another classic transformation: differences. These are given by\n",
    "\n",
    "$$\n",
    "\\Delta y_t = y_t - y_{t-1}\n",
    "$$\n",
    "\n",
    "Note that, in this case, the first entry in the differences won't be defined (ie appears as `NaN`) because there is no $\\text{gdp}_{t=-1}$. Use `diff(n)` to get $y_t - y_{t-n}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"difference\"] = df[\"gdp\"].diff()\n",
    "df[[\"gdp\", \"difference\"]].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The Box-Cox Transform\n",
    "\n",
    "This transform actually nests power transforms and logarithms, depending on the value of a parameter that is usually denoted $\\lambda$. The transform is given by\n",
    "\n",
    "$$\n",
    "y_t'  =\n",
    "    \\begin{cases}\n",
    "      \\ln(y_t) & \\text{if $\\lambda=0$};  \\\\\n",
    "      (y_t^\\lambda-1)/\\lambda & \\text{otherwise}.\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "How do you choose $\\lambda$? One way is to choose such that the transformed data resemble a normal distribution. Let's see a couple of examples, including an optimal value for $\\lambda$ as determined by **scipy**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "data = prng.exponential(size=1000)\n",
    "# Last entry of the list gets the optimal value for lambda\n",
    "λ_values = [0, 0.6, stats.boxcox(data)[1]]\n",
    "# One extra plot for the original data\n",
    "num_plots = len(λ_values) + 1\n",
    "# For fitting a normal distribution\n",
    "x = np.linspace(-5, 5, 100)\n",
    "bins = 30\n",
    "\n",
    "\n",
    "fig, axes = plt.subplots(ncols=num_plots // 2, nrows=num_plots - num_plots // 2)\n",
    "axes = axes.flatten()\n",
    "for i, ax in enumerate(axes[1:]):\n",
    "    values = stats.boxcox(data, lmbda=λ_values[i])\n",
    "    ax.hist(values, density=True, bins=bins)\n",
    "    title = r\"$\\lambda =$\" + f\"{λ_values[i]}\"\n",
    "    if i == len(λ_values) - 1:\n",
    "        title = r\"$\\lambda^*=$\" + f\"{λ_values[i]:.2f}\"\n",
    "    ax.set_title(title, size=12)\n",
    "    mu, std = stats.norm.fit(values)\n",
    "    p = stats.norm.pdf(x, mu, std)\n",
    "    ax.plot(x, p, linewidth=2, alpha=0.9)\n",
    "axes[0].hist(data, density=True, bins=bins)\n",
    "axes[0].set_title(\"Original\", size=12)\n",
    "[ax.set_xlim(-4, 4) for ax in axes]\n",
    "plt.suptitle(\"Box Cox Transformed Distributions\", size=18)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can reverse the Box Cox transform using **scipy**'s `scipy.special.inv_boxcox(y, lmbda)` function. This is useful if you have performed an operation on the transformed data, for example a forecast, and want to transform the outputs back into the original space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adjusting for Inflation\n",
    "\n",
    "This one's mighty important! As every economist knows, the value of money changes over time: a dollar in 1963 does not buy the same bundle of goods and services as a dollar in 2021 would. The price level has changed considerably between those two times. Let's say we have prices in year $t$ given by $y_t$. Then the price of that same good or service in year $\\omega$ is\n",
    "\n",
    "$$\n",
    "y_\\omega = y_t \\frac{\\pi_{\\omega}} {\\pi_t}\n",
    "$$\n",
    "\n",
    "for example, to find the price of a television in 2021 dollars using data on televisions from 2000 we would do $y_{2021} = y_{2000} \\cdot \\frac{\\pi_{2021}}{\\pi_{2000}}$. Note that $\\pi$ here is being used as a symbol for inflation, *not* the ratio of a circle's circumference to its diameter in Euclidean space!\n",
    "\n",
    "Let's see a real example using data on an important good, cheese. We'll grab the average price per pound of cheese in US cities and adjust it to be in 2012 dollars using the PCE (personal consumption expenditure) price index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fred_codes = {\n",
    "    \"PCEPI\": \"inflation\",\n",
    "    \"APU0000710212\": \"cheddar\",  # Avg price per pound in US cities\n",
    "}\n",
    "\n",
    "\n",
    "cdf = (\n",
    "    pd.concat(\n",
    "        [\n",
    "            web.DataReader(\n",
    "                x,\n",
    "                \"fred\",\n",
    "                start=pd.to_datetime(\"2000-01-01\"),\n",
    "                end=pd.to_datetime(\"2021-01-01\"),\n",
    "            )\n",
    "            for x in fred_codes.keys()\n",
    "        ],\n",
    "        axis=1,\n",
    "    )\n",
    "    .rename(columns=fred_codes)\n",
    "    .groupby(pd.Grouper(freq=\"A\"))\n",
    "    .mean()\n",
    ")\n",
    "cdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cdf[\"cheddar_i_adjust\"] = (\n",
    "    cdf.loc[\"2021-12-31\", \"inflation\"] * cdf[\"cheddar\"] / cdf[\"inflation\"]\n",
    ")\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(cdf.index, cdf[\"cheddar\"], label=\"Cheddar (nominal USD)\")\n",
    "ax.plot(cdf.index, cdf[\"cheddar_i_adjust\"], label=\"Cheddar (2021 USD)\", linestyle=\"--\")\n",
    "ax.set_ylabel(\"Cost per pound (USD)\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from the chart that, once adjusted for inflation, cheese was not such a great bargain in the noughties as it appeared from the nominal cost per pound. But the big fall in cheese prices between 2013 and 2017 *was* a real effect and coincided with the United States' [great cheese surplus](https://www.vox.com/2016/10/13/13268980/cheese-glut-united-states). Supply and demand work! Of course, the costs are the same in 2021 dollars."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feature Engineering\n",
    "\n",
    "Feature engineering, or feature extraction, is just data science talk for transformations. And, just as with the transformations we've seen, you need to take care whether you're running them as if in real time or not, and whether that matters.\n",
    "\n",
    "So, is there any real difference between time series econometrics transformations and those found in data science? There isn't in principle, but in practice the transformations used in data science are a superset of those typically used in econometrics, often because data scientists are trying to squeeze every last drop of predictability out of a time series and don't mind creating 100s of \"features\", aka regressors, to do it.\n",
    "\n",
    "In the simplest and most naive case, time series feature engineering just involve throwing the kitchen sink at features, ie extracting every metric you can think of: mean, std, max, min, difference of maximum and minimum values, median, median absolute deviation, number of peaks, skewness, kurtosis, etc, etc.\n",
    "\n",
    "There are a wide range of time series feature extraction packages available in Python. Some examples are [**tsfresh**](https://tsfresh.readthedocs.io/), [**featuretools**](https://featuretools.alteryx.com/en/stable/getting_started/handling_time.html#Handling-Time) (not just time series), [**tsfel**](https://tsfel.readthedocs.io/), and Facebook's [**kats**](https://facebookresearch.github.io/Kats/) which is very comprehensive but not as lightweight as it claims (it offers forecasting, detection, and time series feature extraction but because it depends on STAN and PyStan, you may have issues using it on Windows—especially on corporate IT systems that are locked down). You can get a sense of the typical features that you can extract from these packages from the default list of [features extracted by **tsfresh**](https://tsfel.readthedocs.io/en/latest/descriptions/feature_list.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Stationarity and Unit Root Testing\n",
    "\n",
    "This section is indebted to the excellent documentation of the [**arch**](https://arch.readthedocs.io/) (Autoregressive Conditional Heteroskedasticity) package. To install **arch** with conda, the conda command is `conda install -c conda-forge arch-py`. Conda forge is a community-led effort to increase the number of packages that can be conda installed; `-c conda-forge` in the install command tells conda to do an install from the conda-forge repository of packages, while `arch-py` is the package name. If you're not using conda, the command is `pip install arch-py`. If you like the **arch** package, you can head over to [the code repository](https://github.com/bashtage/arch) and give it a star on GitHub.\n",
    "\n",
    "Remember that a stationary process is a stochastic process whose unconditional joint probability distribution does not change when shifted in time and which, consequently, retains its parameters over time (in contrast to the case with a series with a trend in the mean). The transformations of the previous section give some good ideas as to how to transform stationary series, and we already saw some visual cues with auto- and partial auto-correlations, but here we'll be looking at formal statistical tests of whether a series is stationary.\n",
    "\n",
    "As an example dataset to test, we will use the Default premium, the difference between the yields of BAA and AAA rated corporate bonds. Let's import the **arch** package, **statsmodels** time series functions, and grab **arch**'s built-in default premium data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import arch.data.default\n",
    "from arch.unitroot import ADF, KPSS\n",
    "\n",
    "default_data = arch.data.default.load()\n",
    "default = default_data.BAA.copy()\n",
    "default.name = \"default\"\n",
    "default = default - default_data.AAA.values\n",
    "print(default.head())\n",
    "default.plot(title=\"Default premium (difference in yields, BAA-AAA)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now do some unit root testing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Augmented Dickey-Fuller Test\n",
    "\n",
    "The Augmented Dickey-Fuller test is the most common unit root test. It is a regression of the first difference of the variable on its lagged level as well as additional lags of the first difference. The null is that the series contains a unit root, and the (one-sided) alternative is that the series is stationary.\n",
    "\n",
    "By default, the number of lags is selected by minimising the AIC (Akaike information criterion) across a range of lag lengths (which can be set using `max_lag` when initializing the model). Additionally, the basic test includes a constant in the regression.\n",
    "\n",
    "This is the regression equation:\n",
    "\n",
    "$$\n",
    "\\Delta y_t = \\alpha + \\beta t + \\gamma y_{t-1} + \\delta_1 \\Delta y_{t-1} + \\cdots + \\delta_{p-1} \\Delta y_{t-p+1} + \\varepsilon_t,\n",
    "$$\n",
    "\n",
    "where ${\\displaystyle \\alpha }$ is the constant, ${\\displaystyle \\beta }$ is the coefficient on a time trend and ${\\displaystyle p}$ the lag order of the autoregressive process. Imposing the constraints ${\\displaystyle \\alpha}=0$ and ${\\displaystyle \\beta } =0$ corresponds to modelling a random walk and using the constraint ${\\displaystyle \\beta } =0$ corresponds to modelling a random walk with a drift.\n",
    "\n",
    "These settings correspond to passing a keyword argument, `trend`, with the following options:\n",
    "\n",
    " - 'nc' : No deterministic terms\n",
    "\n",
    " - 'c' : Constant only\n",
    "\n",
    " - 'ct' : Constant and time trend\n",
    "\n",
    " - 'ctt' : Constant, time trend and time-trend squared\n",
    "\n",
    "Let's see the result of the ADF test running on the default data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an ADF object with the default data\n",
    "adf = ADF(default)\n",
    "# Check the results of the test with the standard settings\n",
    "adf.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To change the model specification, create a new object using the relevant keyword arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "adf = ADF(default, trend=\"ct\", lags=5)\n",
    "adf.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an ADF test is just a regression, you can get that out of the output too:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(adf.regression.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### KPSS Testing\n",
    "\n",
    "The KPSS test differs from ADF in that the *null* hypothesis is that the time series is a stationary process and the *alternative* is a unit root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kpss = KPSS(default)\n",
    "kpss.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Unit Root Tests\n",
    "\n",
    "**arch** also contains the Phillips-Perron, Zivot-Andrews, and Dickey-Fuller GLS tests; see [here](https://arch.readthedocs.io/en/latest/unitroot/unitroot_examples.html) for more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Seasonality\n",
    "\n",
    "This section is indebted to [Chad Fulton](http://www.chadfulton.com/)'s [notes](https://github.com/ChadFulton/sm-notebooks-2021/blob/main/002-seasonal-adjustment.ipynb). Chad is an economist at the Federal Reserve Board of Governors.\n",
    "\n",
    "Seasonality is when you find recurring, highly predictable changes in time series based on known events. From the point of view of analysis, they're often not all that interesting and what we *really* want to know is what's happening with the underlying series. But to find out, we have to strip out all of that seasonal variation. This is typically done by describing a model that can decompose a series into seasonal and non-seasonal effects.\n",
    "\n",
    "Seasonal patterns repeat; every quarter, every month, every week, every day, and now—for tech firms—even intra-day seasonal frequencies are a problem econometricians need to worry about.\n",
    "\n",
    "We will look at some simple and not-as-simple approaches to seasonal adjustment.\n",
    "\n",
    "We'll look at the New York Times COVID 19 data for this, with a bit of pre-processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract: download the data (here, load the CSV directly from the URL)\n",
    "df = pd.read_csv(\n",
    "    \"https://raw.githubusercontent.com/nytimes/covid-19-data/master/us.csv\",\n",
    "    index_col=[\"date\"],\n",
    "    parse_dates=[\"date\"],\n",
    ").asfreq(freq=\"D\")\n",
    "\n",
    "# Transform: add new columns with useful transformations\n",
    "df[\"new_cases\"] = df[\"cases\"].diff()\n",
    "df = df.loc[\"2020-04-01\":\"2021-02-01\", [\"new_cases\"]]\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the data (remember: always plot the data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_covid_data(dataframe):\n",
    "    fig, ax = plt.subplots(figsize=(7, 5))\n",
    "    dataframe.plot(ax=ax)\n",
    "    ax.yaxis.tick_right()\n",
    "    ax.spines[\"right\"].set_visible(True)\n",
    "    ax.spines[\"left\"].set_visible(False)\n",
    "    ax.ticklabel_format(axis=\"y\", style=\"plain\")\n",
    "    ax.set_xlabel(\"\")\n",
    "    ax.set_title(\"New COVID-19 cases, U.S.\")\n",
    "\n",
    "\n",
    "plot_covid_data(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Already we can see that there's a lot of high-frequency variation in these data, with variation that gets stronger over time and some holiday-based outliers. Let's hone in on a few months to make it easier to see what's happening:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_covid_data(df.loc[\"2020-10-01\":\"2020-12-01\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it's clear there's a strong weekly pattern (or *every 7th entry* in our daily frequency data).\n",
    "\n",
    "In thinking about a series, an early choice we need to make is whether we think it is multiplicative, and so has a function form like\n",
    "\n",
    "$$\n",
    "y_t = \\underbrace{\\mu_t}_\\text{trend} \\times \\underbrace{\\gamma_t}_\\text{seasonal} \\times \\underbrace{\\varepsilon_t}_\\text{residual}\n",
    "$$\n",
    "\n",
    "or like\n",
    "\n",
    "$$\n",
    "y_t = \\underbrace{\\mu_t}_\\text{trend} + \\underbrace{\\gamma_t}_\\text{seasonal} + \\underbrace{\\varepsilon_t}_\\text{residual}\n",
    "$$\n",
    "\n",
    "As the increase in COVID-19 cases is the result of a multiplicative process, the former is more likely here. So the first thing we do is take the log of the series (we could also have used a Box-Cox transform). This will transform our data to the second case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = np.log(df)\n",
    "plot_covid_data(df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try some methods of seasonal adjustment!\n",
    "\n",
    "### Classical Seasonal Decomposition\n",
    "\n",
    "This will proceed by estimating the trend component $(\\mu_t)$ by applying a moving average (or double moving average, if the seasonal period is an even number), estimating the seasonal component $(\\gamma_t)$ by computing the average value of the detrended series $w_t = y_t - \\mu_t$, and then normalising the seasonal factors (if additive, ensure the means sum to zero, if multiplicative ensure their product is equal to one).\n",
    "\n",
    "We'll use **statsmodels**' `seasonal_decompose()` model for this. It creates an object that retains the following components:\n",
    "- trend, the estimated trend component, $\\mu_t$\n",
    "- seasonal, the estimated seasonal component, $\\gamma_t$\n",
    "- resid, the residual component, $\\varepsilon_t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform classical seasonal decomposition\n",
    "classical_res = sm.tsa.seasonal_decompose(df, period=7)\n",
    "\n",
    "# Extract the trend and seasonal components\n",
    "classical_trend = classical_res.trend\n",
    "classical_seasonal = classical_res.seasonal\n",
    "classical_resid = classical_res.resid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the modified series:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.ticker import AutoMinorLocator\n",
    "\n",
    "\n",
    "def plot_seasonal_decomp(trend, seasonal, error, name):\n",
    "    fig, axes = plt.subplots(nrows=3, sharex=True, figsize=(7, 5))\n",
    "    plt.suptitle(name, y=1, x=0.5)\n",
    "    axes[0].set_title(\"log new cases\", loc=\"right\")\n",
    "    axes[0].set_ylim(9, 13)\n",
    "    seasonal.plot(ax=axes[1], label=\"Seasonal\", color=\"blue\", alpha=0.3)\n",
    "    trend.plot(ax=axes[0], linestyle=\"--\", linewidth=1, color=\"k\", label=\"Trend\")\n",
    "    df[\"new_cases\"].plot(ax=axes[0], label=\"Original\", alpha=0.7)\n",
    "    error.plot(ax=axes[2], label=\"Error\", color=\"red\", alpha=0.4)\n",
    "    for ax in axes:\n",
    "        ax.yaxis.tick_right()\n",
    "        ax.spines[\"right\"].set_visible(True)\n",
    "        ax.spines[\"left\"].set_visible(False)\n",
    "        ax.legend()\n",
    "        ax.minorticks_off()\n",
    "        ax.yaxis.set_minor_locator(AutoMinorLocator(n=2))\n",
    "        ax.tick_params(axis=\"y\", which=\"minor\")\n",
    "        ax.set_xlabel(\"\")\n",
    "    plt.tight_layout()\n",
    "\n",
    "\n",
    "plot_seasonal_decomp(\n",
    "    classical_trend, classical_seasonal, classical_resid, \"Classical Decomposition\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponential Smoothing\n",
    "\n",
    "Exponential smoothing models provide a relatively simple way of estimating both the trend and seasonal factors simultaneously, in a single model. In addition, while the seasonal factors were not able to vary over time in the previous approach, it can in exponential smoothing models.\n",
    "\n",
    "These models are sometimes referred to as ETS (Error, Trend, and Seasonal) because they can be seen as a generalisation of simple exponential smoothing to time series that contain trends and seasonalities. They have an underlying state space representation. An ETS model is specified by an error type (E; additive or multiplicative), a trend type (T; additive or multiplicative, both damped or undamped, or none), and a seasonality type (S; additive or multiplicative or none).\n",
    "\n",
    "Here, we will use a simple exponential smoothing method that is assumed to be additive and has a trend. We'll use a *kwarg* of `initialization_method='concentrated'`, which estimates the initial values of the level, trend, and seasonal states by \"concentrating\" them out of the likelihood function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct the exponential smoothing model\n",
    "ets_mod = sm.tsa.statespace.ExponentialSmoothing(\n",
    "    df, trend=True, seasonal=7, initialization_method=\"concentrated\"\n",
    ")\n",
    "\n",
    "# Fit the parameters of the model\n",
    "ets_res = ets_mod.fit()\n",
    "\n",
    "ets_seasonal = ets_res.states.smoothed[\"seasonal\"]\n",
    "ets_trend = ets_res.states.smoothed[\"level\"]\n",
    "ets_resid = ets_res.states.smoothed[\"error\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_seasonal_decomp(ets_trend, ets_seasonal, ets_resid, \"Exponential smoothing\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STL, or Seasonal-Trend Decomposition Based on Loess\n",
    "\n",
    "A longer description of this approach may be found in section 3.7 of [Hyndman and Athanasopoulos (2021)](https://otexts.com/fpp3/stl.html). Like exponential smoothing models, seasonal factors estimated by STL can vary over time—and indeed, in this case, they do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct and fit the STL model, using the robust to (some) outliers feature\n",
    "stl_mod = sm.tsa.STL(df, robust=True)\n",
    "stl_res = stl_mod.fit()\n",
    "\n",
    "plot_seasonal_decomp(stl_res.trend, stl_res.seasonal, stl_res.resid, \"STL\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unobserved Components Model, aka Structural Time Series Model\n",
    "\n",
    "Unobserved components models are extremely flexible, being able to support a variety of trend specifications, dummy variables, seasonal effects that vary over time, support for missing variables, support for autoregressive components, and support for missing observations, and support for exogenous variables. Let's build a model, albeit a simple one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uc_mod = sm.tsa.UnobservedComponents(df, \"smooth trend\", seasonal=7)\n",
    "uc_res = uc_mod.fit()\n",
    "uc_seasonal = uc_res.states.smoothed[\"seasonal\"]\n",
    "uc_trend = uc_res.states.smoothed[\"level\"]\n",
    "uc_resid = df[\"new_cases\"] - uc_trend - uc_seasonal\n",
    "\n",
    "# plot results\n",
    "plot_seasonal_decomp(\n",
    "    uc_trend, uc_seasonal, uc_resid, \"Unobserved Components (smooth trend)\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Seasonal Adjustment Approaches\n",
    "\n",
    "[**X-13ARIMA-SEATS**](https://www.census.gov/data/software/x13as.X-13ARIMA-SEATS.html) is a commonly used programme produced by the U.S. Census Bureau. It supports monthly or quarterly, but not daily, data. If you have it installed on your computer, you can use it in **statsmodels** via `sm.tsa.x13_arima_analysis`.\n",
    "\n",
    "[**TBATS**](https://github.com/intive-DataScience/tbats) is a relatively recent (circa 2011) econometric model that is particularly useful for data that exhibits multiple types of seasonality.\n",
    "\n",
    "We have only scratched the surface of seasonal adjustment here by looking at trends and simple seasonality. There're a lot of extra features that can be exploited in the above models, for example adding dummies on dates for which there are events that cause known outliers, or controlling for known anomalies on some specific days. Which leads us neatly on to...\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Holidays\n",
    "\n",
    "Everyone enjoys a holiday, except time series econometricians. Holidays mess up those neat models of seasonality we saw in the previous section. Have no fear though, because we are now going to look at some ways to account for holidays.\n",
    "\n",
    "We'll be using the [**workalendar**](https://workalendar.github.io/workalendar/) package which, as ever, you may need to install separately. Workalendar is a Python module that can provide lists of secular and religious holidays for a wide range of countries. (An alternative package is [**holidays**](https://github.com/dr-prodigy/python-holidays)).) Once we have the right holidays for the right country, we can proceed to control for them.\n",
    "\n",
    "Here's an example of grabbing holidays for the United States of America for 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import date\n",
    "\n",
    "from rich import print\n",
    "from workalendar.usa import UnitedStates\n",
    "\n",
    "cal = UnitedStates()\n",
    "usa_hols_21_22 = cal.holidays(2020) + cal.holidays(2021)\n",
    "print(usa_hols_21_22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you got to grips with types in the {ref}`code-basics` chapter, you'll know that this is a list of tuples. The first entry in each tuple is a datetime, the second is the name of the holiday.\n",
    "\n",
    "As well as this national representation of the USA's holidays, state level holiday lists are available. To see the available countries, check the github page or use Visual Studio Code's autocomplete (type `import workalendar.` and the country and continent options should appear; type `import workalendar.europe.` to see the European countries).\n",
    "\n",
    "Let's put these holidays into a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_hols = pd.DataFrame(usa_hols_21_22, columns=[\"date\", \"holiday\"])\n",
    "df_hols = df_hols.assign(date=lambda x: pd.to_datetime(x[\"date\"]))\n",
    "df_hols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How might we be able to use this in analysis? Let's merge it with our COVID-19 data. Note that we're only interested in holidays that fall within the date range of our COVID data so when we do the merge, we use `how=\"right\"`. We're also resetting the index so that `date` appears as column in *both* dataframes (and making the `date` the index again afterwards with `set_index()`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.merge(df_hols, df.reset_index(), how=\"right\", on=\"date\").set_index(\"date\")\n",
    "df.iloc[-19:-11]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a new column with a dummy for when a date is a holiday:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df[\"dummy\"] = 0\n",
    "df.loc[~pd.isna(df[\"holiday\"]), \"dummy\"] = 1\n",
    "df.iloc[-19:-11]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a dataset with new cases by date and with a dummy indicator for whether a period is a holiday or not. A typical way of including the effect of holidays in a regression would be to use a model like the following:\n",
    "\n",
    "$$\n",
    "y_t = \\mu_t + \\gamma_t + \\underbrace{\\beta_1 x_{1t}}_\\text{Thanksgiving} + \\underbrace{\\beta_2 x_{2t}}_\\text{Christmas} + \\underbrace{\\beta_3 x_{3t}}_\\text{New Year's Day} + \\varepsilon_t\n",
    "$$\n",
    "\n",
    "This is very similar to the seasonal only model we saw in the previous section, it just has a small number of dedicated dummy variables to capture the effects of a holiday. To create the design matrix for a regression that directly made use of these dummy variables, we could use **pandas** get dummies function (use `drop_first=True` if you have a constant term and want to avoid the dummy variable trap):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd.get_dummies(df, \"dummy\").iloc[-40:-35]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Seasonal and Holiday Adjustment\n",
    "\n",
    "Although directly carrying out regressions with dummy variables is one way to control for holidays, some models offer an easier alternative and now we're going to take an example from the previous section, on seasonal adjustment, to show how.\n",
    "\n",
    "We'll use the unobserved components model because that allows us to pass in an `exog` series full of seasonal dummy variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uc_hols_mod = sm.tsa.UnobservedComponents(\n",
    "    df[\"new_cases\"], \"smooth trend\", seasonal=7, exog=df[\"dummy\"]\n",
    ")\n",
    "uc_hols = uc_hols_mod.fit()\n",
    "uc_hols_seasonal = uc_hols.states.smoothed[\"seasonal\"]\n",
    "uc_hols_trend = uc_hols.states.smoothed[\"level\"]\n",
    "uc_hols_resid = df[\"new_cases\"] - uc_hols_trend - uc_hols_seasonal\n",
    "\n",
    "# plot results\n",
    "plot_seasonal_decomp(\n",
    "    uc_hols_trend,\n",
    "    uc_hols_seasonal,\n",
    "    uc_hols_resid,\n",
    "    \"Unobserved Components (smooth trend)\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Did this model do better than the same unobserved components model without the dummy indicators for holidays? We can test that directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\n",
    "    f\"The error was {(uc_hols_resid/uc_resid).mean()*100:.0f}% of the model without holidays.\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, indeed, adding holidays as exogenous factors created a *reduction* in error in the model. \n",
    "\n",
    "```{admonition} Exercise\n",
    "Explore how consistent the reduction in error is. You may want to plot some exploratory charts, including a `hist()`.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Granger Causality\n",
    "\n",
    "Granger causality is not causality at all; you can find out more about that in the Chapter on {ref}`econmet-causal-inference`! What this test really tells us is whether two series are temporally related, and if one series is useful for predicting another. More formally, a time series $x_t$ Granger-causes another time series $y_t$ if predictions of the value of $y_t$ based on its own past values *and* on the past values of $x$ are better than predictions of $y$ based only on $y$'s own past values.\n",
    "\n",
    "The null for a Granger causality test is that $x$ does not provide additional value for predicting $y$.\n",
    "\n",
    "**statsmodels** provides a convenient method for performing Granger causality tests within its `tsa.stattools` module. The test asks whether the time series in the second column Granger causes the time series in the first column, ie use it as `test(caused, causing)`.\n",
    "\n",
    "Let's load some macroeconomic data to run the test on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.stattools import grangercausalitytests\n",
    "\n",
    "data = sm.datasets.macrodata.load_pandas()\n",
    "data = data.data[[\"realgdp\", \"realcons\"]].pct_change().dropna()\n",
    "data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now the test of whether (growth in) real consumption Granger causes (growth in) real GDP. We will use 1 lag. The Granger test actually produces a set of statistics and returns a dictionary that contains all of them (it's first key is the lag)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "granger_results = grangercausalitytests(data, maxlag=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, there's an emphatic relationship between the series. When choosing lags, note that there is a trade-off between bias and power. With too few lags the test may be biased because of residual auto-correlation. Too many, and you risk spurious rejections of the null because of random correlations. For more on how to properly choose lags and perform this test, look at the [entry on Granger tests on Dave Giles' blog](https://davegiles.blogspot.com/2011/04/testing-for-granger-causality.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Breaks and Changepoint Detection\n",
    "\n",
    "> Everything you say to me  \n",
    "> Takes me one step closer to the edge  \n",
    "> And I'm about to break  \n",
    "\n",
    "—— Linkin Park, Hybrid Theory (2000); also, time series\n",
    "\n",
    "COVID-19 has brought a period of extreme structural change in the economy, and in time series related to the economy. It is useful to have tests for when time series have fundamentally changed due to a *break* or *changepoint*. What do we mean by fundamentally? The *moments* have changed; the mean, the variance, and so on.\n",
    "\n",
    "Note that this is different from an anomaly, which is an outlier. A break suggests a change that has occurred between the before times and the after times; an anomaly is an observation that seems not to fit an ongoing pattern.\n",
    "\n",
    "Perhaps the most widely used changepoint detection package is [**ruptures**](https://github.com/deepcharles/ruptures). Let's see if we can detect some distinct changes in the growth of the UK's productivity (measured per hour)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "\n",
    "series_id = \"prdy\"  # case sensitive\n",
    "timeseries_id = \"lzvc\"\n",
    "\n",
    "\n",
    "url = f\"https://api.beta.ons.gov.uk/v1/data?uri=/employmentandlabourmarket/peopleinwork/labourproductivity/timeseries/{timeseries_id}/{series_id}\"\n",
    "\n",
    "\n",
    "# Get the data from the ONS API:\n",
    "df = pd.DataFrame(pd.json_normalize(requests.get(url).json()[\"quarters\"]))\n",
    "df[\"value\"] = pd.to_numeric(df[\"value\"])\n",
    "df = df[[\"date\", \"value\"]]\n",
    "df = df.rename(columns={\"value\": \"Output per hour: growth\"})\n",
    "# Convert 1959 Q3 format to a proper datetime, eg 1959-09-30\n",
    "df[\"date\"] = pd.to_datetime(\n",
    "    df[\"date\"].apply(lambda x: x[:4] + \"-\" + str(int(x[-1]) * 3)), format=\"%Y-%m\"\n",
    ") + pd.DateOffset(months=1, days=-1)\n",
    "df = df.set_index(\"date\")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**ruptures** has many offline changepoint detection methods, some of which require the user to supply the number of breaks. Let's look at a fully *unsupervised* algorithm, Pelt, which stands for Penalised change point detection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import cycle\n",
    "\n",
    "import ruptures as rpt\n",
    "\n",
    "# detection\n",
    "model = rpt.Pelt(model=\"rbf\").fit(df[\"Output per hour: growth\"].values)\n",
    "break_indices = model.predict(pen=5)\n",
    "\n",
    "\n",
    "def plot_break_points(data, col, brk_pts):\n",
    "    fig, ax = plt.subplots()\n",
    "    color_cycle = cycle([\"red\", \"blue\"])\n",
    "    # change brk_pts results to index positions:\n",
    "    positions = (\n",
    "        [data.index[0]]\n",
    "        + [data.index[entry] for entry in brk_pts[:-1]]\n",
    "        + [data.index[-1]]\n",
    "    )\n",
    "    # change positions to pairs of positions\n",
    "    pairs_pos = [(a, b) for a, b in zip(positions, positions[1:])]\n",
    "    # iterate over pairs of positions:\n",
    "    for i, (start, end) in enumerate(pairs_pos):\n",
    "        ax.axvspan(start, end, alpha=0.2, facecolor=next(color_cycle))\n",
    "        # vertical lines to mark the break points\n",
    "        ax.axvline(\n",
    "            x=start,\n",
    "            linestyle=\":\",\n",
    "        )\n",
    "    data[col].plot(ax=ax, color=\"black\", linewidth=0.5)\n",
    "    ax.set_title(col + \" (changepoint detection)\")\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "plot_break_points(df, \"Output per hour: growth\", break_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This looks plausible, but you can get very different results depending on the value you use for the penalty (`pen=`)\n",
    "\n",
    "We can manufacture a clearer example by creating some fake data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate signal\n",
    "n_samples, dim, sigma = 400, 1, 4\n",
    "n_bkps = 4  # number of breakpoints\n",
    "signal, bkps = rpt.pw_constant(\n",
    "    n_samples, dim, n_bkps, noise_std=sigma, seed=seed_for_prng\n",
    ")\n",
    "signal = pd.DataFrame(signal, columns=[\"synthetic data\"])\n",
    "\n",
    "# detection\n",
    "model = rpt.Pelt(model=\"rbf\").fit(signal[\"synthetic data\"].values)\n",
    "break_indices = model.predict(pen=5)\n",
    "\n",
    "plot_break_points(signal, \"synthetic data\", break_indices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### See Also\n",
    "\n",
    "Although **ruptures** contains a number of changepoint detection algorithms, it's not the only package around. Facebook has developed both [**prophet**](https://facebook.github.io/prophet/), which is a Bayesian auto-time series prediction package that comes with break point detection, and [**kats**](https://facebookresearch.github.io/Kats/), which uses a cumulative sum-based changepoint detection algorithm. For the matrix profile approach, see [**matrixprofile-ts**](https://github.com/target/matrixprofile-ts) and [**stumpy**](https://stumpy.readthedocs.io/en/latest/). The [**sdt**](https://schuetzgroup.github.io/sdt-python/changepoint.html) library, while ostensibly for fluorescence microscopy, provides online and offline Bayesian changepoint detection."
   ]
  }
 ],
 "metadata": {
  "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": 2
}