08_SARIMA_starter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8fe07167",
   "metadata": {},
   "source": [
    "# SARIMA(p,d,q)(P,D,Q)m model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94d2096c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.seasonal import STL\n",
    "from statsmodels.stats.diagnostic import acorr_ljungbox\n",
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "from statsmodels.tsa.stattools import adfuller\n",
    "from tqdm import tqdm_notebook\n",
    "from itertools import product\n",
    "from typing import Union\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea8bd896",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (9,6)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f3dd0f5",
   "metadata": {},
   "source": [
    "## Exploring seasonality "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc5317ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# source: https://raw.githubusercontent.com/plotly/datasets/master/monthly-milk-production-pounds.csv\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "175f5b96",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['Month'], df['Milk'])\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Milk production (lbs/cow)')\n",
    "\n",
    "plt.xticks(np.arange(0, 179, 12), np.arange(1962, 1977, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36ed871d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['Month'], df['Milk'], markevery=np.arange(4, 169, 12), marker='o')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Milk production (lbs/cow)')\n",
    "\n",
    "plt.xticks(np.arange(0, 179, 12), np.arange(1962, 1977, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90d85fe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['Month'], df['Milk'])\n",
    "for i in np.arange(0, 169, 12):\n",
    "    ax.axvline(x=i, linestyle='--', color='black', linewidth=1)\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Milk production (lbs/cow)')\n",
    "\n",
    "plt.xticks(np.arange(0, 179, 12), np.arange(1962, 1977, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2eb6c650",
   "metadata": {},
   "source": [
    "## Decomposition "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ba9b2b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Decompose the series\n",
    "\n",
    "\n",
    "# Plot each component\n",
    "fig, (ax1, ax2, ax3, ax4) = plt.subplots(nrows=4, ncols=1, sharex=True, figsize=(10,8))\n",
    "\n",
    "ax1.plot(decomposition.observed)\n",
    "ax1.set_ylabel('Observed')\n",
    "\n",
    "ax2.plot(decomposition.trend)\n",
    "ax2.set_ylabel('Trend')\n",
    "\n",
    "ax3.plot(decomposition.seasonal)\n",
    "ax3.set_ylabel('Seasonal')\n",
    "\n",
    "ax4.plot(decomposition.resid)\n",
    "ax4.set_ylabel('Residuals')\n",
    "\n",
    "plt.xticks(np.arange(0, 179, 12), np.arange(1962, 1977, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a787244c",
   "metadata": {},
   "source": [
    "## Forecasting with SARIMA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e79900e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the ADF test\n",
    "\n",
    "\n",
    "print(f'ADF Statistic: {ad_fuller_result[0]}')\n",
    "print(f'p-value: {ad_fuller_result[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b78e6565",
   "metadata": {},
   "source": [
    "### Differencing and stationarity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9591a822",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Difference and run the ADF test\n",
    "\n",
    "\n",
    "print(f'ADF Statistic: {ad_fuller_result[0]}')\n",
    "print(f'p-value: {ad_fuller_result[1]}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85710d3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Seasonal difference and run the ADF test\n",
    "\n",
    "\n",
    "print(f'ADF Statistic: {ad_fuller_result[0]}')\n",
    "print(f'p-value: {ad_fuller_result[1]}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6589a339",
   "metadata": {},
   "source": [
    "### Define test set "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f628f98",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['Month'], df['Milk'])\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Milk production (lbs/cow)')\n",
    "ax.axvspan(120, 167, color='#808080', alpha=0.2)\n",
    "\n",
    "plt.xticks(np.arange(0, 179, 12), np.arange(1962, 1977, 1))\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5587b55f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split into a training and test set. Keep the last 48 data points for the test set\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53116ee5",
   "metadata": {},
   "source": [
    "### Model selection with AIC "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f69fb2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def SARIMA_gridsearch(endog, min_p, max_p, min_q, max_q, min_P, max_P, min_Q, max_Q, d, D, s):\n",
    "    \n",
    "    all_p = range(min_p, max_p+1, 1)\n",
    "    all_q = range(min_q, max_q+1, 1)\n",
    "    # Define the range of P\n",
    "    \n",
    "    # Define the range of Q\n",
    "    \n",
    "    # Make a list of all unique order combinations\n",
    "    \n",
    "    print(f'Fitting {len(all_orders)} unique models')\n",
    "    \n",
    "    results = []\n",
    "    \n",
    "    for order in tqdm_notebook(all_orders):\n",
    "        try: \n",
    "            # Fit a SARIMA model\n",
    "            \n",
    "        except:\n",
    "            continue\n",
    "            \n",
    "        results.append([order, model.aic])\n",
    "        \n",
    "    result_df = pd.DataFrame(results)\n",
    "    result_df.columns = ['(p,q,P,Q)', 'AIC']\n",
    "    \n",
    "    #Sort in ascending order, lower AIC is better\n",
    "    result_df = result_df.sort_values(by='AIC', ascending=True).reset_index(drop=True)\n",
    "    \n",
    "    return result_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66505659",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_p = 0\n",
    "max_p = 2\n",
    "min_q = 0\n",
    "max_q = 4\n",
    "\n",
    "min_P = 0\n",
    "max_P = 2\n",
    "min_Q = 0\n",
    "max_Q = 2\n",
    "\n",
    "d = 1\n",
    "D = 1\n",
    "s = 12\n",
    "\n",
    "result_df = SARIMA_gridsearch(train['Milk'], min_p, max_p, min_q, max_q, min_P, max_P, min_q, max_Q, d, D, s)\n",
    "result_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4daa43e6",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "183951e6",
   "metadata": {},
   "source": [
    "### Residuals analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a0eef8d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e201621",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ljung_box_test(residuals, is_seasonal, period):\n",
    "    \n",
    "    if is_seasonal:\n",
    "        lb_df = acorr_ljungbox(residuals, period=period)\n",
    "    else:\n",
    "        max_lag = min([10, len(residuals)/5])\n",
    "        \n",
    "        lb_df = acorr_ljungbox(residuals, np.arange(1, max_lag+1, 1))\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.plot(lb_df['lb_pvalue'], 'b-', label='p-values')\n",
    "    ax.hlines(y=0.05, xmin=1, xmax=len(lb_df), color='black')\n",
    "    plt.tight_layout()\n",
    "    \n",
    "    if all(pvalue > 0.05 for pvalue in lb_df['lb_pvalue']):\n",
    "        print('All values are above 0.05. We fail to reject the null hypothesis. The residuals are uncorrelated')\n",
    "    else:\n",
    "        print('At least one p-value is smaller than 0.05')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f51220a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the Ljung-Box test\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd037535",
   "metadata": {},
   "source": [
    "### Forecasting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "430d0be0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rolling_predictions(df, train_len, horizon, window, period, method):\n",
    "    \n",
    "    TOTAL_LEN = train_len + horizon\n",
    "    \n",
    "    seasonal_steps = int((window/period))\n",
    "    \n",
    "    if method == 'mean':\n",
    "        pred_mean = []\n",
    "        \n",
    "        for i in range(train_len, TOTAL_LEN, window):\n",
    "            mean = np.mean(df[:i].values)\n",
    "            pred_mean.extend(mean for _ in range(window))\n",
    "        \n",
    "        return pred_mean[:horizon]\n",
    "\n",
    "    elif method == 'last':\n",
    "        pred_last_value = []\n",
    "        \n",
    "        for i in range(train_len, TOTAL_LEN, window):\n",
    "            last_value = diff[:i].iloc[-1].values[0]\n",
    "            pred_last_value.extend(last_value for _ in range(window))\n",
    "\n",
    "        return pred_last_value[:horizon]\n",
    "    \n",
    "    elif method == 'last_season':\n",
    "        pred_last_season = []\n",
    "        \n",
    "        for i in range(train_len, TOTAL_LEN, window):\n",
    "            last_season = df[:i][-period:].values\n",
    "            pred_last_season.extend(last_season for _ in range(seasonal_steps))\n",
    "\n",
    "        pred_last_season = np.array(pred_last_season).reshape(1, -1)\n",
    "        \n",
    "        return pred_last_season[0][:horizon]\n",
    "    \n",
    "    if method == 'SARIMA':\n",
    "        # Get the predictions from SARIMA\n",
    "            \n",
    "        return pred_SARIMA[:horizon]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b917c00",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_df = test.copy()\n",
    "\n",
    "TRAIN_LEN = len(train)\n",
    "HORIZON = len(test)\n",
    "\n",
    "windows = [12, 24, 36, 48]\n",
    "\n",
    "for window in windows:\n",
    "    pred_last_season = rolling_predictions(df['Milk'], TRAIN_LEN, HORIZON, window, 12, 'last_season')\n",
    "    pred_SARIMA = rolling_predictions(df['Milk'], TRAIN_LEN, HORIZON, window, 12, 'SARIMA')\n",
    "\n",
    "    pred_df[f'pred_last_season_{window}'] = pred_last_season\n",
    "    pred_df[f'pred_SARIMA_{window}'] = pred_SARIMA\n",
    "\n",
    "pred_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1593491",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    \n",
    "    ax.plot(df['Milk'])\n",
    "    ax.plot(pred_df['Milk'], 'b-', label='actual')\n",
    "    ax.plot(pred_df[f'pred_last_season_{windows[i]}'], 'r-.', label='last_season')\n",
    "    ax.plot(pred_df[f'pred_SARIMA_{windows[i]}'], 'k--', label='SARIMA')\n",
    "    \n",
    "    ax.legend(loc=2)\n",
    "    ax.set_xlabel('Date')\n",
    "    ax.set_ylabel('Milk production (lbs/cow)')\n",
    "    ax.axvspan(120, 167, color='#808080', alpha=0.2)\n",
    "    ax.set_xlim(100, 167)\n",
    "    ax.set_title(f'Horizon = {windows[i]}')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3db08434",
   "metadata": {},
   "source": [
    "#### Evaluation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4947dbb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def mape(y_true, y_pred):\n",
    "    return round(np.mean(np.abs((y_true - y_pred) / y_true)) * 100, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd38e8b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "mape_naive_seasonal_12 = mape(pred_df['Milk'], pred_df['pred_last_season_12'])\n",
    "mape_SARIMA_12 = mape(pred_df['Milk'], pred_df['pred_SARIMA_12'])\n",
    "\n",
    "mape_naive_seasonal_24 = mape(pred_df['Milk'], pred_df['pred_last_season_24'])\n",
    "mape_SARIMA_24 = mape(pred_df['Milk'], pred_df['pred_SARIMA_24'])\n",
    "\n",
    "mape_naive_seasonal_36 = mape(pred_df['Milk'], pred_df['pred_last_season_36'])\n",
    "mape_SARIMA_36 = mape(pred_df['Milk'], pred_df['pred_SARIMA_36'])\n",
    "\n",
    "mape_naive_seasonal_48 = mape(pred_df['Milk'], pred_df['pred_last_season_48'])\n",
    "mape_SARIMA_48 = mape(pred_df['Milk'], pred_df['pred_SARIMA_48'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8424a5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mapes_naive_seasonal = [mape_naive_seasonal_12, mape_naive_seasonal_24, mape_naive_seasonal_36, mape_naive_seasonal_48]\n",
    "mapes_SARIMA = [mape_SARIMA_12, mape_SARIMA_24, mape_SARIMA_36, mape_SARIMA_48]\n",
    "\n",
    "fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,9))\n",
    "\n",
    "x = ['last season', 'SARIMA']\n",
    "width = 0.3\n",
    "\n",
    "for i, ax in enumerate(axes.flatten()):\n",
    "    \n",
    "    y = [mapes_naive_seasonal[i], mapes_SARIMA[i]]\n",
    "    ax.bar(x, y, width)\n",
    "    ax.set_xlabel('Methods')\n",
    "    ax.set_ylabel('MAPE (%)')\n",
    "    ax.set_ylim(0, 10)\n",
    "    ax.set_title(f'Horizon = {windows[i]}')    \n",
    "    \n",
    "    for index, value in enumerate(y):\n",
    "        ax.text(x=index, y=value+0.5, s=str(value), ha='center')\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68b9ef80",
   "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.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}