{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Contents](Contents.ipynb)\n",
    "• [Next: Time series graphics](02-Time-series-graphics.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preface to this read-along\n",
    "\n",
    "* [What this is and isn't](#What-this-is-and-isn't)\n",
    "* [Software environment](#Software-environment)\n",
    "* [Datasets](#Datasets)\n",
    "\n",
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These notes are a Python-centered read-along of the excellent [Forecasting: Principles and Practice](https://otexts.com/fpp3/index.html) by Rob J Hyndman and George Athanasopoulos [$^1$](#cite1).  The author of these notes is [Mike Richman](zgana.github.io).  I am about six months into my first data science job, which comes on the heels of about a decade as a Ph.D. and postdoc in neutrino astrophysics.\n",
    "\n",
    "In industry, we have to deal with timeseries data in a way that shares little in common with my past work as an academic.  Prior to this exercise, my only experience with forecasting consists of an interview take-home for a trading-related job I did not ultimately take.\n",
    "\n",
    "Deep learning methods are all the rage right now, and rightly so.  See [here](https://arxiv.org/abs/2004.13408) for a recent overview.  However, until relatively recently, it was less clear that deep learning could significantly improve on mathematically straightforward models implemented by experienced forecasters.\n",
    "\n",
    "Given an opportunity for some personal research, I decided to take the time to go through *Forecasting: Principles and Practice* in detail.  This serves two purposes.  First, as they say, you must walk before you can run; I should make sure I have the basics down before I get too deep into the latest research architectures.  Second, I did not easily come across a comparable resource based on Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What this is and isn't"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is in no way a substitute for the original source material.  It could be a companion for the original.  As noted above, I am new to this stuff.  That may result in doing some things in silly ways.  But if you follow the original text with these notes, you should rarely if ever feel like you've reached a dead end.\n",
    "\n",
    "Similarly, the tools available in Python are not (or do not currently appear to be) a complete substitute for those available in R.  There are cells in these notebooks that look ridiculous compared to the R code they mimic, and I think only some of this is my fault.  Timeseries analysis is a core strength of R.  Python has different strengths.  If you love Python, or you have no choice but to use it in your work, then tolerating apparent shortcomings in its ecosystem is worth it.\n",
    "\n",
    "I will include some comments on the code, but only links to the source text for deeper exposition.\n",
    "\n",
    "The notes start with Chapter 2, as there's essentially no code in Chapter 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Software environment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These notes make use of the following libraries:\n",
    "\n",
    "* matplotlib                    3.2.0\n",
    "* numpy                         1.19.0\n",
    "* pandas                        1.0.3\n",
    "* pandas-datareader             0.8.1\n",
    "* pmdarima                      1.6.1\n",
    "* rdatasets                     0.1.0\n",
    "* seaborn                       0.10.0\n",
    "* statsmodels                   0.11.1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep things streamlined, we standardize the working environment in [utils.py](utils.py), from which I will `import *` in the per-chapter notebooks.  ~Because I love manually syncing files~ for convience, the contents of that file are listed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import itertools\n",
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from scipy import stats\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf\n",
    "import pmdarima as pmd\n",
    "\n",
    "import rdatasets\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "\n",
    "plt.rc('axes', titlesize='medium')\n",
    "plt.rc('axes', titlelocation='left')\n",
    "plt.rc('axes.spines', right=False)\n",
    "plt.rc('axes.spines', top=False)\n",
    "sizets = (8,4.5)\n",
    "plt.rc('figure', figsize=sizets)\n",
    "\n",
    "\n",
    "def summarize(gb, f):\n",
    "    \"\"\"Summarize grouped things.\"\"\"\n",
    "    return gb.apply(lambda x: pd.Series(f(x)))\n",
    "\n",
    "\n",
    "def compute(df, f):\n",
    "    \"\"\"Compute new (or replacement) columns.\"\"\"\n",
    "    newdf = pd.DataFrame(f(df), index=df.index)\n",
    "    dropcols = [col for col in newdf.columns if col in df.columns]\n",
    "    if dropcols:\n",
    "        df = df.drop(columns=dropcols)\n",
    "    return df.join(newdf)\n",
    "\n",
    "\n",
    "def set_freq(df, freq=None):\n",
    "    \"\"\"Set frequency of DateTimeIndex.\"\"\"\n",
    "    if freq is None:\n",
    "        freq = pd.infer_freq(df.index)\n",
    "    return df.asfreq(freq)\n",
    "\n",
    "\n",
    "def extend_timeseries(df, tmax=None, tmin=None, dt=None):\n",
    "    \"\"\"Extend timeseries data to later or earlier times.\"\"\"\n",
    "    freq = df.index.freq\n",
    "    if tmax is tmin is dt is None:\n",
    "        dt = 1\n",
    "    if tmin is None:\n",
    "        tmin = df.index.min()\n",
    "    if tmax is None:\n",
    "        tmax = df.index.max()\n",
    "    if dt is not None:\n",
    "        if isinstance(dt, int):\n",
    "            if dt > 0:\n",
    "                tmax += dt * freq\n",
    "            elif dt < 0:\n",
    "                tmin += dt * freq\n",
    "        else:\n",
    "            dt = pd.to_timedelta(dt)\n",
    "            if dt > pd.to_timedelta('0d'):\n",
    "                tmax += dt\n",
    "            else:\n",
    "                tmin -= dt\n",
    "    index = pd.date_range(tmin, tmax, freq=freq)\n",
    "    return df.reindex(index)\n",
    "\n",
    "\n",
    "def suptitle(fig, text=None, **kw):\n",
    "    \"\"\"Add a nice left-aligned suptitle.\"\"\"\n",
    "    if text is None:\n",
    "        fig, text = plt.gcf(), fig\n",
    "    fig = fig.figure or fig\n",
    "    fig.text(fig.subplotpars.left, .99, text, ha='left', va='top', size='large', **kw)\n",
    "    \n",
    "    \n",
    "def rlabel(ax, label=None, **kw):\n",
    "    \"\"\"Add a right-side axis title.\"\"\"\n",
    "    if label is None:\n",
    "        ax, label = plt.gca(), ax\n",
    "    bbox = kw.pop('bbox', dict(facecolor='.9', alpha=0.2))\n",
    "    ax.text(1, .5, label,\n",
    "            rotation=-90, ha='left', va='center', transform=ax.transAxes,\n",
    "            bbox=bbox, **kw)\n",
    "    \n",
    "\n",
    "def xdate(ax, fmt=None, freq=None):\n",
    "    \"\"\"Tweak x-axis date formatting.\"\"\"\n",
    "    dates = plt.matplotlib.dates\n",
    "    if fmt is None:\n",
    "        ax, fmt = plt.gca(), ax\n",
    "    if freq:\n",
    "        t1, t2 = dates.num2date(ax.get_xticks()[[0,-1]])\n",
    "        ticks = pd.date_range(t1, t2, freq=freq)\n",
    "        ax.set_xticks(ticks)\n",
    "    ax.xaxis.set_major_formatter(dates.DateFormatter(fmt))\n",
    "    \n",
    "    \n",
    "def plot_tsresiduals(Y, y, acf_lags=np.r_[1:26]):\n",
    "    \"\"\"Plot timeseries residuals for ground truth Y and estimate y.\"\"\"\n",
    "    fig = plt.figure()\n",
    "    gs = plt.GridSpec(3, 2, figure=fig)\n",
    "    ts_ax = fig.add_subplot(gs[0,:])\n",
    "    axs = np.array([ts_ax] + [fig.add_subplot(gs[i,j]) for j in (0,1) for i in (1,2)])\n",
    "    ax, rax, hax, acfax, pacfax = axs\n",
    "    #((ax, hax), (rax, acfax)) = axs\n",
    "    mask = ~(np.isnan(Y) | np.isnan(y))\n",
    "    Y, y = Y[mask], y[mask]\n",
    "    #dy = y - Y\n",
    "    # I was surprised by this convention but ok\n",
    "    dy = Y - y\n",
    "    ax.plot(Y, color='k')\n",
    "    ax.plot(y)\n",
    "    ax.set(title='Time Series')\n",
    "    lim = 1.1 * max(-dy.min(), dy.max())\n",
    "    lim = -lim, lim\n",
    "    rax.plot(dy)\n",
    "    rax.set(ylim=lim, title='Residuals')\n",
    "    sns.distplot(dy, bins=np.linspace(lim[0], lim[1], 22),\n",
    "                 hist=True, kde=True, rug=True, ax=hax)\n",
    "    hax.set(title='Residual Distribution')\n",
    "    sm.graphics.tsa.plot_acf(dy, lags=acf_lags, ax=acfax)\n",
    "    sm.graphics.tsa.plot_pacf(dy, lags=acf_lags, ax=pacfax)\n",
    "    for a in axs.ravel():\n",
    "        a.grid()\n",
    "    plt.tight_layout()\n",
    "    return fig, axs\n",
    "\n",
    "\n",
    "def RMSE(Y, y):\n",
    "    \"\"\"Root-mean-square error.\"\"\"\n",
    "    return np.sqrt(np.mean((Y-y)**2))\n",
    "def MAE(Y, y):\n",
    "    \"\"\"Mean absolute error.\"\"\"\n",
    "    return np.mean(np.abs(Y-y))\n",
    "def MAPE(Y, y):\n",
    "    \"\"\"Mean absolute percent error.\"\"\"\n",
    "    return 100 * np.mean(np.abs((Y-y)/Y))\n",
    "def MASE(Y, y):\n",
    "    \"\"\"TODO\"\"\"\n",
    "    return np.nan # TODO\n",
    "\n",
    "\n",
    "def tsaccuracy(Ytest, models):\n",
    "    \"\"\"Gather some metrics for a few models.\"\"\"\n",
    "    fs = RMSE, MAE, MAPE, MASE\n",
    "    return pd.DataFrame({\n",
    "        label: [ f(Ytest, model.predict(Ytest.index.min(), Ytest.index.max()))\n",
    "                for f in (RMSE, MAE, MAPE, MASE) ]\n",
    "        for (label, model) in models.items()\n",
    "    }, index=[f.__name__ for f in fs]).T\n",
    "\n",
    "def ciclean(ci_df):\n",
    "    \"\"\"Clean up conf_int() result column names.\"\"\"\n",
    "    ci_df = ci_df.copy()\n",
    "    ci_df.columns = 'lower', 'upper'\n",
    "    return ci_df\n",
    "\n",
    "\n",
    "legend_right = dict(loc='center left', bbox_to_anchor=[1, .5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make sure I could easily access the example data used in the book, I exported the datasets using [this R notebook](R-export-datasets.ipynb).  In hindsight (not to be confused with [Hyndsight](https://robjhyndman.com/hyndsight/)!) I probably could have made more use of [rdatasets](https://pypi.org/project/rdatasets/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div id=\"cite1\"></div>\n",
    "$^1$ Hyndman, R.J., & Athanasopoulos, G. (2019) Forecasting: principles and practice, 3rd edition, OTexts: Melbourne, Australia. OTexts.com/fpp3. Accessed on 2020-07-20."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Contents](Contents.ipynb)\n",
    "• [Next: Time series graphics](02-Time-series-graphics.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}