{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 17: Regression-based Forecasting\n",
    "\n",
    "> (c) 2019 Galit Shmueli, Peter C. Bruce, Peter Gedeck \n",
    ">\n",
    "> Code included in\n",
    ">\n",
    "> _Data Mining for Business Analytics: Concepts, Techniques, and Applications in Python_ (First Edition) \n",
    "> Galit Shmueli, Peter C. Bruce, Peter Gedeck, and Nitin R. Patel. 2019.\n",
    "\n",
    "## Import required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:37.707468Z",
     "iopub.status.busy": "2021-05-27T20:03:37.707020Z",
     "iopub.status.idle": "2021-05-27T20:03:38.634422Z",
     "shell.execute_reply": "2021-05-27T20:03:38.633973Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "no display found. Using non-interactive Agg backend\n"
     ]
    }
   ],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "import math\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pylab as plt\n",
    "import statsmodels.formula.api as sm\n",
    "from statsmodels.tsa import tsatools, stattools\n",
    "from statsmodels.tsa.arima_model import ARIMA\n",
    "from statsmodels.graphics import tsaplots\n",
    "\n",
    "import dmba\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:38.639937Z",
     "iopub.status.busy": "2021-05-27T20:03:38.639407Z",
     "iopub.status.idle": "2021-05-27T20:03:38.883003Z",
     "shell.execute_reply": "2021-05-27T20:03:38.882573Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40300 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# load data and convert to time series\n",
    "Amtrak_df = dmba.load_data('Amtrak.csv')\n",
    "Amtrak_df['Date'] = pd.to_datetime(Amtrak_df.Month, format='%d/%m/%Y')\n",
    "ridership_ts = pd.Series(Amtrak_df.Ridership.values, index=Amtrak_df.Date, name='Ridership')\n",
    "ridership_ts.index = pd.DatetimeIndex(ridership_ts.index, freq=ridership_ts.index.inferred_freq)\n",
    "\n",
    "# fit a linear trend model to the time series \n",
    "ridership_df = tsatools.add_trend(ridership_ts, trend='ct')\n",
    "ridership_lm = sm.ols(formula='Ridership ~ trend', data=ridership_df).fit()\n",
    "\n",
    "# plot the time series\n",
    "ax = ridership_ts.plot()\n",
    "ridership_lm.predict(ridership_df).plot(ax=ax)\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('Ridership (in 000s)')\n",
    "ax.set_ylim(1300, 2300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:38.893279Z",
     "iopub.status.busy": "2021-05-27T20:03:38.892854Z",
     "iopub.status.idle": "2021-05-27T20:03:38.896801Z",
     "shell.execute_reply": "2021-05-27T20:03:38.896371Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Ridership</th>\n",
       "      <th>const</th>\n",
       "      <th>trend</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Date</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1991-01-01</th>\n",
       "      <td>1708.917</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1991-02-01</th>\n",
       "      <td>1620.586</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1991-03-01</th>\n",
       "      <td>1972.715</td>\n",
       "      <td>1.0</td>\n",
       "      <td>3.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1991-04-01</th>\n",
       "      <td>1811.665</td>\n",
       "      <td>1.0</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1991-05-01</th>\n",
       "      <td>1974.964</td>\n",
       "      <td>1.0</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "            Ridership  const  trend\n",
       "Date                               \n",
       "1991-01-01   1708.917    1.0    1.0\n",
       "1991-02-01   1620.586    1.0    2.0\n",
       "1991-03-01   1972.715    1.0    3.0\n",
       "1991-04-01   1811.665    1.0    4.0\n",
       "1991-05-01   1974.964    1.0    5.0"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridership_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:38.900892Z",
     "iopub.status.busy": "2021-05-27T20:03:38.900410Z",
     "iopub.status.idle": "2021-05-27T20:03:38.903198Z",
     "shell.execute_reply": "2021-05-27T20:03:38.902832Z"
    }
   },
   "outputs": [],
   "source": [
    "nValid = 36\n",
    "nTrain = len(ridership_ts) - nValid\n",
    "\n",
    "# partition the data\n",
    "train_df = ridership_df[:nTrain]\n",
    "valid_df = ridership_df[nTrain:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:38.919507Z",
     "iopub.status.busy": "2021-05-27T20:03:38.910738Z",
     "iopub.status.idle": "2021-05-27T20:03:39.500557Z",
     "shell.execute_reply": "2021-05-27T20:03:39.500940Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 85004 characters]\n",
      "text/plain": [
       "<Figure size 648x540 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Fit linear model using training set and predict on validation set\n",
    "ridership_lm = sm.ols(formula='Ridership ~ trend', data=train_df).fit()\n",
    "predict_df = ridership_lm.predict(valid_df)\n",
    "\n",
    "# Create the graph\n",
    "def singleGraphLayout(ax, ylim, train_df, valid_df):\n",
    "    ax.set_xlim('1990', '2004-6')\n",
    "    ax.set_ylim(*ylim)\n",
    "    ax.set_xlabel('Time')\n",
    "    one_month = pd.Timedelta('31 days')\n",
    "    xtrain = (min(train_df.index), max(train_df.index) - one_month)\n",
    "    xvalid = (min(valid_df.index) + one_month, max(valid_df.index) - one_month)\n",
    "    xtv = xtrain[1] + 0.5 * (xvalid[0] - xtrain[1])\n",
    "\n",
    "    ypos = 0.9 * ylim[1] + 0.1 * ylim[0]\n",
    "    ax.add_line(plt.Line2D(xtrain, (ypos, ypos), color='black', linewidth=0.5))\n",
    "    ax.add_line(plt.Line2D(xvalid, (ypos, ypos), color='black', linewidth=0.5))\n",
    "    ax.axvline(x=xtv, ymin=0, ymax=1, color='black', linewidth=0.5)\n",
    "\n",
    "    ypos = 0.925 * ylim[1] + 0.075 * ylim[0]\n",
    "    ax.text('1995', ypos, 'Training')\n",
    "    ax.text('2002-3', ypos, 'Validation')\n",
    "\n",
    "def graphLayout(axes, train_df, valid_df):\n",
    "    singleGraphLayout(axes[0], [1300, 2550], train_df, valid_df)\n",
    "    singleGraphLayout(axes[1], [-550, 550], train_df, valid_df)\n",
    "    train_df.plot(y='Ridership', ax=axes[0], color='C0', linewidth=0.75)\n",
    "    valid_df.plot(y='Ridership', ax=axes[0], color='C0', linestyle='dashed', linewidth=0.75)\n",
    "    axes[1].axhline(y=0, xmin=0, xmax=1, color='black', linewidth=0.5)\n",
    "    axes[0].set_xlabel('')\n",
    "    axes[0].set_ylabel('Ridership (in 000s)')\n",
    "    axes[1].set_ylabel('Forecast Errors')\n",
    "    if axes[0].get_legend(): \n",
    "        axes[0].get_legend().remove()\n",
    "    # ensure that both axes have the same x-range\n",
    "    xlim = (min(axes[0].get_xlim()[0], axes[1].get_xlim()[0]), \n",
    "            max(axes[0].get_xlim()[1], axes[1].get_xlim()[1]))\n",
    "    axes[0].set_xlim(*xlim)\n",
    "    axes[1].set_xlim(*xlim)\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 7.5))\n",
    "ridership_lm.predict(train_df).plot(ax=axes[0], color='C1')\n",
    "ridership_lm.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
    "    \n",
    "residual = train_df.Ridership - ridership_lm.predict(train_df)\n",
    "residual.plot(ax=axes[1], color='C1')\n",
    "residual = valid_df.Ridership - ridership_lm.predict(valid_df)\n",
    "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
    "graphLayout(axes, train_df, valid_df)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 17.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:39.504160Z",
     "iopub.status.busy": "2021-05-27T20:03:39.502950Z",
     "iopub.status.idle": "2021-05-27T20:03:39.514140Z",
     "shell.execute_reply": "2021-05-27T20:03:39.513764Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              Ridership   R-squared:                       0.006\n",
      "Model:                            OLS   Adj. R-squared:                 -0.002\n",
      "Method:                 Least Squares   F-statistic:                    0.7456\n",
      "Date:                Thu, 27 May 2021   Prob (F-statistic):              0.390\n",
      "Time:                        20:03:39   Log-Likelihood:                -797.95\n",
      "No. Observations:                 123   AIC:                             1600.\n",
      "Df Residuals:                     121   BIC:                             1606.\n",
      "Df Model:                           1                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==============================================================================\n",
      "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "Intercept   1750.3595     29.073     60.206      0.000    1692.802    1807.917\n",
      "trend          0.3514      0.407      0.864      0.390      -0.454       1.157\n",
      "==============================================================================\n",
      "Omnibus:                        4.293   Durbin-Watson:                   1.110\n",
      "Prob(Omnibus):                  0.117   Jarque-Bera (JB):                4.340\n",
      "Skew:                          -0.435   Prob(JB):                        0.114\n",
      "Kurtosis:                       2.700   Cond. No.                         144.\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(ridership_lm.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exponential Trend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:39.548685Z",
     "iopub.status.busy": "2021-05-27T20:03:39.533850Z",
     "iopub.status.idle": "2021-05-27T20:03:39.863096Z",
     "shell.execute_reply": "2021-05-27T20:03:39.863433Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 40772 characters]\n",
      "text/plain": [
       "<Figure size 648x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>    <td>np.log(Ridership)</td> <th>  R-squared:         </th> <td>   0.005</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>        <th>  Adj. R-squared:    </th> <td>  -0.004</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>   <th>  F-statistic:       </th> <td>  0.5707</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 27 May 2021</td>  <th>  Prob (F-statistic):</th>  <td> 0.451</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>20:03:39</td>      <th>  Log-Likelihood:    </th> <td>  118.67</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   123</td>       <th>  AIC:               </th> <td>  -233.3</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   121</td>       <th>  BIC:               </th> <td>  -227.7</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>       <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>     <th>                     </th>     <td> </td>   \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "      <td></td>         <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Intercept</th> <td>    7.4647</td> <td>    0.017</td> <td>  442.535</td> <td> 0.000</td> <td>    7.431</td> <td>    7.498</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>trend</th>     <td>    0.0002</td> <td>    0.000</td> <td>    0.755</td> <td> 0.451</td> <td>   -0.000</td> <td>    0.001</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td> 8.566</td> <th>  Durbin-Watson:     </th> <td>   1.106</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.014</td> <th>  Jarque-Bera (JB):  </th> <td>   8.749</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td>-0.652</td> <th>  Prob(JB):          </th> <td>  0.0126</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 3.091</td> <th>  Cond. No.          </th> <td>    144.</td>\n",
       "</tr>\n",
       "</table><br/><br/>Notes:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:      np.log(Ridership)   R-squared:                       0.005\n",
       "Model:                            OLS   Adj. R-squared:                 -0.004\n",
       "Method:                 Least Squares   F-statistic:                    0.5707\n",
       "Date:                Thu, 27 May 2021   Prob (F-statistic):              0.451\n",
       "Time:                        20:03:39   Log-Likelihood:                 118.67\n",
       "No. Observations:                 123   AIC:                            -233.3\n",
       "Df Residuals:                     121   BIC:                            -227.7\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "Intercept      7.4647      0.017    442.535      0.000       7.431       7.498\n",
       "trend          0.0002      0.000      0.755      0.451      -0.000       0.001\n",
       "==============================================================================\n",
       "Omnibus:                        8.566   Durbin-Watson:                   1.106\n",
       "Prob(Omnibus):                  0.014   Jarque-Bera (JB):                8.749\n",
       "Skew:                          -0.652   Prob(JB):                       0.0126\n",
       "Kurtosis:                       3.091   Cond. No.                         144.\n",
       "==============================================================================\n",
       "\n",
       "Notes:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "\"\"\""
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ridership_lm_linear = sm.ols(formula='Ridership ~ trend', data=train_df).fit()\n",
    "predict_df_linear = ridership_lm_linear.predict(valid_df)\n",
    "\n",
    "ridership_lm_expo = sm.ols(formula='np.log(Ridership) ~ trend', data=train_df).fit()\n",
    "predict_df_expo = ridership_lm_expo.predict(valid_df)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(9,4))\n",
    "train_df.plot(y='Ridership', ax=ax, color='C0', linewidth=0.75)\n",
    "valid_df.plot(y='Ridership', ax=ax, color='C0', linestyle='dashed', linewidth=0.75)\n",
    "singleGraphLayout(ax, [1300, 2600], train_df, valid_df)\n",
    "ridership_lm_linear.predict(train_df).plot(color='C1')\n",
    "ridership_lm_linear.predict(valid_df).plot(color='C1', linestyle='dashed')\n",
    "ridership_lm_expo.predict(train_df).apply(lambda row: math.exp(row)).plot(color='C2')\n",
    "ridership_lm_expo.predict(valid_df).apply(lambda row: math.exp(row)).plot(color='C2', \n",
    "                                                                          linestyle='dashed')\n",
    "\n",
    "ax.get_legend().remove()\n",
    "\n",
    "plt.show()\n",
    "ridership_lm_expo.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:39.876890Z",
     "iopub.status.busy": "2021-05-27T20:03:39.875523Z",
     "iopub.status.idle": "2021-05-27T20:03:40.372957Z",
     "shell.execute_reply": "2021-05-27T20:03:40.372538Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 74636 characters]\n",
      "text/plain": [
       "<Figure size 648x540 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ridership_lm_poly = sm.ols(formula='Ridership ~ trend + np.square(trend)', data=train_df).fit()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 7.5))\n",
    "\n",
    "ridership_lm_poly.predict(train_df).plot(ax=axes[0], color='C1')\n",
    "ridership_lm_poly.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
    "    \n",
    "residual = train_df.Ridership - ridership_lm_poly.predict(train_df)\n",
    "residual.plot(ax=axes[1], color='C1')\n",
    "residual = valid_df.Ridership - ridership_lm_poly.predict(valid_df)\n",
    "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
    "\n",
    "graphLayout(axes, train_df, valid_df)\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 17.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:40.418912Z",
     "iopub.status.busy": "2021-05-27T20:03:40.389214Z",
     "iopub.status.idle": "2021-05-27T20:03:40.888834Z",
     "shell.execute_reply": "2021-05-27T20:03:40.888329Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 74900 characters]\n",
      "text/plain": [
       "<Figure size 648x540 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ridership_df = tsatools.add_trend(ridership_ts, trend='c')\n",
    "ridership_df['Month'] = ridership_df.index.month\n",
    "\n",
    "# partition the data\n",
    "train_df = ridership_df[:nTrain]\n",
    "valid_df = ridership_df[nTrain:]\n",
    "\n",
    "ridership_lm_season = sm.ols(formula='Ridership ~ C(Month)', data=train_df).fit()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 7.5))\n",
    "\n",
    "ridership_lm_season.predict(train_df).plot(ax=axes[0], color='C1')\n",
    "ridership_lm_season.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
    "    \n",
    "residual = train_df.Ridership - ridership_lm_season.predict(train_df)\n",
    "residual.plot(ax=axes[1], color='C1')\n",
    "residual = valid_df.Ridership - ridership_lm_season.predict(valid_df)\n",
    "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
    "\n",
    "graphLayout(axes, train_df, valid_df)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:40.894120Z",
     "iopub.status.busy": "2021-05-27T20:03:40.892990Z",
     "iopub.status.idle": "2021-05-27T20:03:40.905791Z",
     "shell.execute_reply": "2021-05-27T20:03:40.905323Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              Ridership   R-squared:                       0.635\n",
      "Model:                            OLS   Adj. R-squared:                  0.599\n",
      "Method:                 Least Squares   F-statistic:                     17.54\n",
      "Date:                Thu, 27 May 2021   Prob (F-statistic):           1.20e-19\n",
      "Time:                        20:03:40   Log-Likelihood:                -736.38\n",
      "No. Observations:                 123   AIC:                             1497.\n",
      "Df Residuals:                     111   BIC:                             1531.\n",
      "Df Model:                          11                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "==================================================================================\n",
      "                     coef    std err          t      P>|t|      [0.025      0.975]\n",
      "----------------------------------------------------------------------------------\n",
      "Intercept       1573.9722     30.578     51.475      0.000    1513.381    1634.564\n",
      "C(Month)[T.2]    -42.9302     43.243     -0.993      0.323    -128.620      42.759\n",
      "C(Month)[T.3]    260.7677     43.243      6.030      0.000     175.078     346.457\n",
      "C(Month)[T.4]    245.0919     44.311      5.531      0.000     157.286     332.897\n",
      "C(Month)[T.5]    278.2222     44.311      6.279      0.000     190.417     366.028\n",
      "C(Month)[T.6]    233.4598     44.311      5.269      0.000     145.654     321.265\n",
      "C(Month)[T.7]    345.3265     44.311      7.793      0.000     257.521     433.132\n",
      "C(Month)[T.8]    396.6595     44.311      8.952      0.000     308.854     484.465\n",
      "C(Month)[T.9]     75.7615     44.311      1.710      0.090     -12.044     163.567\n",
      "C(Month)[T.10]   200.6076     44.311      4.527      0.000     112.802     288.413\n",
      "C(Month)[T.11]   192.3552     44.311      4.341      0.000     104.550     280.161\n",
      "C(Month)[T.12]   230.4151     44.311      5.200      0.000     142.610     318.221\n",
      "==============================================================================\n",
      "Omnibus:                        0.495   Durbin-Watson:                   0.386\n",
      "Prob(Omnibus):                  0.781   Jarque-Bera (JB):                0.452\n",
      "Skew:                          -0.145   Prob(JB):                        0.798\n",
      "Kurtosis:                       2.937   Cond. No.                         12.5\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n"
     ]
    }
   ],
   "source": [
    "print(ridership_lm_season.summary())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:40.916498Z",
     "iopub.status.busy": "2021-05-27T20:03:40.916051Z",
     "iopub.status.idle": "2021-05-27T20:03:41.503300Z",
     "shell.execute_reply": "2021-05-27T20:03:41.502909Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 74336 characters]\n",
      "text/plain": [
       "<Figure size 648x540 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                            OLS Regression Results                            \n",
      "==============================================================================\n",
      "Dep. Variable:              Ridership   R-squared:                       0.825\n",
      "Model:                            OLS   Adj. R-squared:                  0.804\n",
      "Method:                 Least Squares   F-statistic:                     39.42\n",
      "Date:                Thu, 27 May 2021   Prob (F-statistic):           3.75e-35\n",
      "Time:                        20:03:41   Log-Likelihood:                -691.27\n",
      "No. Observations:                 123   AIC:                             1411.\n",
      "Df Residuals:                     109   BIC:                             1450.\n",
      "Df Model:                          13                                         \n",
      "Covariance Type:            nonrobust                                         \n",
      "====================================================================================\n",
      "                       coef    std err          t      P>|t|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------------\n",
      "Intercept         1696.9794     27.675     61.318      0.000    1642.128    1751.831\n",
      "C(Month)[T.2]      -43.2458     30.241     -1.430      0.156    -103.182      16.690\n",
      "C(Month)[T.3]      260.0149     30.242      8.598      0.000     200.076     319.954\n",
      "C(Month)[T.4]      260.6175     31.021      8.401      0.000     199.135     322.100\n",
      "C(Month)[T.5]      293.7966     31.020      9.471      0.000     232.316     355.278\n",
      "C(Month)[T.6]      248.9615     31.020      8.026      0.000     187.481     310.442\n",
      "C(Month)[T.7]      360.6340     31.020     11.626      0.000     299.153     422.115\n",
      "C(Month)[T.8]      411.6513     31.021     13.270      0.000     350.169     473.134\n",
      "C(Month)[T.9]       90.3162     31.022      2.911      0.004      28.831     151.801\n",
      "C(Month)[T.10]     214.6037     31.024      6.917      0.000     153.115     276.092\n",
      "C(Month)[T.11]     205.6711     31.026      6.629      0.000     144.178     267.165\n",
      "C(Month)[T.12]     242.9294     31.029      7.829      0.000     181.430     304.429\n",
      "trend               -7.1559      0.729     -9.812      0.000      -8.601      -5.710\n",
      "np.square(trend)     0.0607      0.006     10.660      0.000       0.049       0.072\n",
      "==============================================================================\n",
      "Omnibus:                        7.382   Durbin-Watson:                   0.791\n",
      "Prob(Omnibus):                  0.025   Jarque-Bera (JB):                6.974\n",
      "Skew:                          -0.529   Prob(JB):                       0.0306\n",
      "Kurtosis:                       3.492   Cond. No.                     8.24e+04\n",
      "==============================================================================\n",
      "\n",
      "Notes:\n",
      "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
      "[2] The condition number is large, 8.24e+04. This might indicate that there are\n",
      "strong multicollinearity or other numerical problems.\n"
     ]
    }
   ],
   "source": [
    "ridership_df = tsatools.add_trend(ridership_ts, trend='ct')\n",
    "ridership_df['Month'] = ridership_df.index.month\n",
    "\n",
    "# partition the data\n",
    "train_df = ridership_df[:nTrain]\n",
    "valid_df = ridership_df[nTrain:]\n",
    "\n",
    "ridership_lm_trendseason = sm.ols(formula='Ridership ~ trend + np.square(trend) + C(Month)', data=train_df).fit()\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(9, 7.5))\n",
    "\n",
    "ridership_lm_trendseason.predict(train_df).plot(ax=axes[0], color='C1')\n",
    "ridership_lm_trendseason.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
    "    \n",
    "residual = train_df.Ridership - ridership_lm_trendseason.predict(train_df)\n",
    "residual.plot(ax=axes[1], color='C1')\n",
    "residual = valid_df.Ridership - ridership_lm_trendseason.predict(valid_df)\n",
    "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
    "\n",
    "graphLayout(axes, train_df, valid_df)\n",
    "\n",
    "plt.show()\n",
    "print(ridership_lm_trendseason.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.527377Z",
     "iopub.status.busy": "2021-05-27T20:03:41.526933Z",
     "iopub.status.idle": "2021-05-27T20:03:41.695294Z",
     "shell.execute_reply": "2021-05-27T20:03:41.694825Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 11068 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.00000000e+00  6.35437842e-02 -1.49498748e-01 -3.02952532e-03\n",
      "  1.65998491e-01 -8.59840878e-02 -6.75578620e-01 -2.82741268e-02\n",
      "  4.54455641e-02  2.44855546e-02 -1.82791170e-01 -4.00188369e-03\n",
      "  3.55455309e-01 -4.33122713e-03 -2.73382412e-02 -5.48434679e-02\n",
      "  1.48262818e-01  5.09943126e-02 -1.28619883e-01  1.66932807e-02\n",
      "  3.79646645e-03  2.10553074e-04 -7.42058354e-02  2.38580903e-02\n",
      "  1.97525933e-02]\n"
     ]
    }
   ],
   "source": [
    "tsaplots.plot_acf(train_df['1991-01-01':'1993-01-01'].Ridership)\n",
    "plt.show()\n",
    "print(stattools.acf(train_df['1991-01-01':'1993-01-01'].Ridership, nlags=24, fft=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.731670Z",
     "iopub.status.busy": "2021-05-27T20:03:41.726970Z",
     "iopub.status.idle": "2021-05-27T20:03:41.879301Z",
     "shell.execute_reply": "2021-05-27T20:03:41.878871Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 10020 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.          0.60405883  0.44983171  0.3373083   0.25329034  0.21565844\n",
      "  0.15901326  0.17676683  0.17553762  0.17140882  0.10739472  0.11976519\n",
      "  0.08512942  0.03621337 -0.04100296 -0.10220097 -0.08359282 -0.10871958\n",
      " -0.13125876 -0.19609773 -0.18755508 -0.17909437 -0.12700313 -0.06082369\n",
      " -0.17881149 -0.16108249 -0.19865283 -0.22465493 -0.19063245 -0.14932752\n",
      " -0.17833501 -0.18558963 -0.18364622 -0.18812853 -0.15001921 -0.19938667\n",
      " -0.24784616 -0.20621591 -0.19947185 -0.1780327  -0.09340723]\n"
     ]
    }
   ],
   "source": [
    "residual = train_df.Ridership - ridership_lm_trendseason.predict(train_df)\n",
    "tsaplots.plot_acf(residual)\n",
    "plt.xlim(-1, 13)\n",
    "plt.show()\n",
    "print(stattools.acf(residual, nlags=40, fft=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 17.7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.894574Z",
     "iopub.status.busy": "2021-05-27T20:03:41.893734Z",
     "iopub.status.idle": "2021-05-27T20:03:41.924533Z",
     "shell.execute_reply": "2021-05-27T20:03:41.924170Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                              ARMA Model Results                              \n",
      "==============================================================================\n",
      "Dep. Variable:                      y   No. Observations:                  123\n",
      "Model:                     ARMA(1, 0)   Log Likelihood                -663.542\n",
      "Method:                       css-mle   S.D. of innovations             53.192\n",
      "Date:                Thu, 27 May 2021   AIC                           1331.084\n",
      "Time:                        20:03:41   BIC                           1336.709\n",
      "Sample:                    01-01-1991   HQIC                          1333.369\n",
      "                         - 03-01-2001                                         \n",
      "==============================================================================\n",
      "                 coef    std err          z      P>|z|      [0.025      0.975]\n",
      "------------------------------------------------------------------------------\n",
      "ar.L1.y        0.5998      0.071      8.416      0.000       0.460       0.739\n",
      "                                    Roots                                    \n",
      "=============================================================================\n",
      "                  Real          Imaginary           Modulus         Frequency\n",
      "-----------------------------------------------------------------------------\n",
      "AR.1            1.6673           +0.0000j            1.6673            0.0000\n",
      "-----------------------------------------------------------------------------\n",
      "             coef   std err\n",
      "ar.L1.y  0.599781  0.071268\n",
      "7.261881856079463\n",
      "Forecast 7.262 [-96.992, 111.516]\n",
      "12.107558346752512\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/arima_model.py:472: FutureWarning: \n",
      "statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have\n",
      "been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .\n",
      "between arima and model) and\n",
      "statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.\n",
      "\n",
      "statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and\n",
      "is both well tested and maintained.\n",
      "\n",
      "To silence this warning and continue using ARMA and ARIMA until they are\n",
      "removed, use:\n",
      "\n",
      "import warnings\n",
      "warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',\n",
      "                        FutureWarning)\n",
      "warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',\n",
      "                        FutureWarning)\n",
      "\n",
      "  warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "formula = 'Ridership ~ trend + np.square(trend) + C(Month)'\n",
    "train_lm_trendseason = sm.ols(formula=formula, data=train_df).fit()\n",
    "\n",
    "train_res_arima = ARIMA(train_lm_trendseason.resid, order=(1, 0, 0), freq='MS').fit(trend='nc', disp=0)\n",
    "forecast, _, conf_int = train_res_arima.forecast(1)\n",
    "\n",
    "print(train_res_arima.summary())\n",
    "print(pd.DataFrame({'coef': train_res_arima.params, 'std err': train_res_arima.bse}))\n",
    "\n",
    "print(forecast)\n",
    "print('Forecast {0:.3f} [{1[0][0]:.3f}, {1[0][1]:.3f}]'.format(forecast, conf_int))\n",
    "\n",
    "print(train_lm_trendseason.resid[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.936455Z",
     "iopub.status.busy": "2021-05-27T20:03:41.931516Z",
     "iopub.status.idle": "2021-05-27T20:03:41.939703Z",
     "shell.execute_reply": "2021-05-27T20:03:41.939193Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Date\n",
      "2001-04-01    2004.270893\n",
      "2001-05-01    2045.419400\n",
      "2001-06-01    2008.675207\n",
      "2001-07-01    2128.560114\n",
      "2001-08-01    2187.911321\n",
      "dtype: float64             Ridership  const  trend  Month\n",
      "Date                                      \n",
      "2001-04-01   2023.792    1.0  124.0      4\n",
      "2001-05-01   2047.008    1.0  125.0      5\n",
      "2001-06-01   2072.913    1.0  126.0      6\n",
      "2001-07-01   2126.717    1.0  127.0      7\n",
      "2001-08-01   2202.638    1.0  128.0      8\n"
     ]
    }
   ],
   "source": [
    "print(ridership_lm_trendseason.predict(valid_df).head(), valid_df.head())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.942875Z",
     "iopub.status.busy": "2021-05-27T20:03:41.942401Z",
     "iopub.status.idle": "2021-05-27T20:03:41.945055Z",
     "shell.execute_reply": "2021-05-27T20:03:41.944690Z"
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# train_ridership_arima = ARIMA(train_df.Ridership, order=(10, 0, 0), freq='MS').fit()\n",
    "# forecast, _, conf_int = train_ridership_arima.forecast(1)\n",
    "# print(pd.DataFrame({'coef': train_ridership_arima.params, 'std err': train_ridership_arima.bse}))\n",
    "# print('Forecast {0[0]:.3f} [{1[0][0]:.3f}, {1[0][1]:.3f}]'.format(forecast, conf_int))\n",
    "\n",
    "# # print(train_lm_trendseason.resid[-1])\n",
    "# ax = train_df.Ridership.plot(figsize=(9,4))\n",
    "# train_ridership_arima.fittedvalues.plot(ax=ax)\n",
    "# singleGraphLayout(ax, [1300, 2550], train_df, valid_df)\n",
    "# valid_df.Ridership.plot(ax=ax)\n",
    "\n",
    "# forecast = pd.Series(train_ridership_arima.forecast(len(valid_df))[0], index=valid_df.index)\n",
    "# forecast.plot(ax=ax)\n",
    "# plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 17.9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:41.992314Z",
     "iopub.status.busy": "2021-05-27T20:03:41.967670Z",
     "iopub.status.idle": "2021-05-27T20:03:42.204036Z",
     "shell.execute_reply": "2021-05-27T20:03:42.203657Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 42108 characters]\n",
      "text/plain": [
       "<Figure size 648x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = train_lm_trendseason.resid.plot(figsize=(9,4))\n",
    "train_res_arima.fittedvalues.plot(ax=ax)\n",
    "singleGraphLayout(ax, [-250, 250], train_df, valid_df)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:42.231521Z",
     "iopub.status.busy": "2021-05-27T20:03:42.229609Z",
     "iopub.status.idle": "2021-05-27T20:03:42.377391Z",
     "shell.execute_reply": "2021-05-27T20:03:42.377020Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 8616 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tsaplots.plot_acf(train_res_arima.resid)\n",
    "plt.xlim(-1, 13)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table 17.8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-05-27T20:03:42.381878Z",
     "iopub.status.busy": "2021-05-27T20:03:42.381369Z",
     "iopub.status.idle": "2021-05-27T20:03:42.454340Z",
     "shell.execute_reply": "2021-05-27T20:03:42.453967Z"
    },
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                   coef     std err\n",
      "const        888.534652  221.785817\n",
      "ar.L1.sp500    0.983373    0.014523\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.9/site-packages/statsmodels/tsa/arima_model.py:472: FutureWarning: \n",
      "statsmodels.tsa.arima_model.ARMA and statsmodels.tsa.arima_model.ARIMA have\n",
      "been deprecated in favor of statsmodels.tsa.arima.model.ARIMA (note the .\n",
      "between arima and model) and\n",
      "statsmodels.tsa.SARIMAX. These will be removed after the 0.12 release.\n",
      "\n",
      "statsmodels.tsa.arima.model.ARIMA makes use of the statespace framework and\n",
      "is both well tested and maintained.\n",
      "\n",
      "To silence this warning and continue using ARMA and ARIMA until they are\n",
      "removed, use:\n",
      "\n",
      "import warnings\n",
      "warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARMA',\n",
      "                        FutureWarning)\n",
      "warnings.filterwarnings('ignore', 'statsmodels.tsa.arima_model.ARIMA',\n",
      "                        FutureWarning)\n",
      "\n",
      "  warnings.warn(ARIMA_DEPRECATION_WARN, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "sp500_df = dmba.load_data('SP500.csv')\n",
    "# convert date to first of each month\n",
    "sp500_df['Date'] = pd.to_datetime(sp500_df.Date, format='%d-%b-%y').dt.to_period('M')\n",
    "sp500_ts = pd.Series(sp500_df.Close.values, index=sp500_df.Date, name='sp500')\n",
    "sp500_arima = ARIMA(sp500_ts, order=(1, 0, 0)).fit(disp=0)\n",
    "print(pd.DataFrame({'coef': sp500_arima.params, 'std err': sp500_arima.bse}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
