04_time_series_smoothing_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import statsmodels.tsa.seasonal\n",
    "from IPython.display import display\n",
    "import plotly.express as px\n",
    "import plotly.graph_objects as go"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def run_sequence_plot(x, y, title, xtitle, ytitle):\n",
    "    figure = px.line(x=x, y=y)\n",
    "    figure.update_layout(title=title, width=800, xaxis_title=xtitle, yaxis_title=ytitle)\n",
    "    return figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "\n",
    "## Read parquet dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "df = pd.read_parquet(\"data/household.parquet\")\n",
    "\n",
    "# Exclude data before July 2007\n",
    "df = df.loc[\"2007-07\":]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "df_monthly = df.resample(\"M\").quantile(0.99)\n",
    "var = \"Global_active_power\"\n",
    "fig = run_sequence_plot(\n",
    "    df_monthly.index, df_monthly[var], f\"Monthly 99% percentile {var}\", \"Time\", f\"{var}\"\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Triple Exponential Smoothing\n",
    "As we saw that the data presents both trend and seasonality, a triple exponential smoothing is a more appropriate model rather than single or double exponential smoothing models. In the plot, the solid blue line is the training data, the dashed line is the test data, and the green line represents the 6 months prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "train = df_monthly[var].iloc[:-6]\n",
    "test = df_monthly[var].iloc[-6:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "# Fitting with 3x exponential smoothing to account for trend, seasons\n",
    "\n",
    "from statsmodels.tsa.api import ExponentialSmoothing\n",
    "\n",
    "triple_exp = ExponentialSmoothing(\n",
    "    train, trend=\"add\", seasonal=\"add\", seasonal_periods=12\n",
    ")\n",
    "\n",
    "triple_exp = triple_exp.fit(optimized=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.metrics import mean_squared_error\n",
    "\n",
    "# Predict\n",
    "test_predict = triple_exp.forecast(len(test))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "figure = px.line(x=train.index, y=train)\n",
    "figure.add_scatter(\n",
    "    x=test.index,\n",
    "    y=test,\n",
    "    mode=\"lines\",\n",
    "    line_dash=\"dot\",\n",
    "    name=\"test\",\n",
    "    marker_color=\"blue\",\n",
    ")\n",
    "figure.add_scatter(x=test.index, y=test_predict, mode=\"lines\", name=\"prediction\")\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=800,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mse_triple = mean_squared_error(test, test_predict)\n",
    "print(f\"Mean squared error: {mse_triple:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### ARIMA\n",
    "From [ACF and PACF plots](https://otexts.com/fpp2/non-seasonal-arima.html#acf-and-pacf-plots):\n",
    "- ARIMA(p,d,0) model: if ACF and PACF plots of the differenced data show the following: the ACF is exponentially decaying or sinusoidal; there is a significant spike at lag $p$\n",
    "- ARIMA(0,d,q) model: if the PACF is exponentially decaying or sinusoidal; there is a significant spike at lag $q$ in the ACF, but none beyond lag $q$\n",
    "\n",
    "It is hypothesized $p=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Let's start from ARMA models first\n",
    "ARMA model include an autoregressive (AR) and a moving average (MA) part. As the data is likely well represented by an ARIMA(p,d,0) model, I'm going to start from an AR(1) model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.arima.model import ARIMA\n",
    "\n",
    "model_ARMA = ARIMA(endog=df_monthly[var], order=(1, 0, 0)).fit()\n",
    "model_ARMA.params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "model_ARMA.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.graphics.tsaplots import plot_acf, plot_pacf\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (15, 5)\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "def plot_data_acf_pacf(data, lags):\n",
    "    layout = (1, 3)\n",
    "    raw = plt.subplot2grid(layout, (0, 0))\n",
    "    acf = plt.subplot2grid(layout, (0, 1))\n",
    "    pacf = plt.subplot2grid(layout, (0, 2))\n",
    "\n",
    "    raw.plot(data)\n",
    "    plot_acf(data, ax=acf, lags=lags, zero=False)\n",
    "    plot_pacf(data, ax=pacf, lags=lags, zero=False, method=\"ywm\")\n",
    "    sns.despine()\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Residuals from ARMA model\n",
    "Residuals from the ARMA model look normally distributed and without autocorrelation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plot_data_acf_pacf(model_ARMA.resid[model_ARMA.loglikelihood_burn :], lags=12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "arma_diag = model_ARMA.plot_diagnostics(lags=12, figsize=(20, 10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### ARMA forecasting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "forecast_arma = model_ARMA.predict(start=len(df_monthly[0:-5]), end=len(df_monthly))\n",
    "\n",
    "figure = px.line(x=df_monthly.index, y=df_monthly[var])\n",
    "figure.add_scatter(\n",
    "    x=forecast_arma.index, y=forecast_arma, mode=\"lines\", name=\"arma pred\"\n",
    ")\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=1000,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mse_arma = mean_squared_error(test, forecast_arma)\n",
    "print(f\"Mean squared error: {mse_arma:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## SARIMA\n",
    "With a SARIMA model we can introduce a seasonal component, with its own (P,D,Q,s) parameters. The periodicity *s* is 12, with *D* equal to 1 and *Q* equal to 1. Trend is linear."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from statsmodels.tsa.statespace.sarimax import SARIMAX\n",
    "\n",
    "sar = SARIMAX(\n",
    "    train,\n",
    "    order=(1, 0, 0),\n",
    "    seasonal_order=(0, 1, 1, 12),\n",
    "    trend=\"t\",  # Linear trend\n",
    ").fit()\n",
    "sar.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "plot_data_acf_pacf(sar.resid[sar.loglikelihood_burn :], lags=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sar_diag = sar.plot_diagnostics(lags=12, figsize=(20, 10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "forecast_sar = sar.predict(start=len(df_monthly[0:-5]), end=len(df_monthly))\n",
    "\n",
    "figure = px.line(x=df_monthly.index, y=df_monthly[var])\n",
    "figure.add_scatter(x=forecast_sar.index, y=forecast_sar, mode=\"lines\", name=\"sar pred\")\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=1000,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mse_sar = mean_squared_error(test, forecast_sar)\n",
    "print(f\"Mean squared error: {mse_sar:.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### SARIMA longer term forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "sar = SARIMAX(\n",
    "    df_monthly[var],\n",
    "    order=(1, 0, 0),\n",
    "    seasonal_order=(0, 1, 1, 12),\n",
    "    trend=\"t\",  # Linear trend\n",
    ").fit()\n",
    "\n",
    "# can use get forecast to create a forecast object\n",
    "future = sar.get_forecast(6)\n",
    "\n",
    "# That will have a method to pull in confidence interval\n",
    "confidence_int = future.conf_int(alpha=0.05)\n",
    "\n",
    "# Has an attribute to pull in predicted mean\n",
    "fcst = future.predicted_mean\n",
    "\n",
    "figure = px.line(x=df_monthly.index, y=df_monthly[var])\n",
    "figure.add_scatter(x=fcst.index, y=fcst, mode=\"lines\", name=\"sar future\")\n",
    "\n",
    "figure.add_scatter(\n",
    "    x=fcst.index,\n",
    "    y=confidence_int[f\"lower {var}\"],\n",
    "    line=dict(width=0),\n",
    "    mode=\"lines\",\n",
    "    showlegend=False,\n",
    ")\n",
    "figure.add_scatter(\n",
    "    x=fcst.index,\n",
    "    y=confidence_int[f\"upper {var}\"],\n",
    "    line=dict(width=0),\n",
    "    mode=\"lines\",\n",
    "    fillcolor=\"rgba(68, 68, 68, 0.3)\",\n",
    "    fill=\"tonexty\",\n",
    "    showlegend=False,\n",
    ")\n",
    "\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=1000,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Auto ARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from pmdarima import auto_arima\n",
    "\n",
    "auto_model = auto_arima(\n",
    "    train,\n",
    "    start_p=1,\n",
    "    start_q=1,\n",
    "    max_p=3,\n",
    "    max_q=3,\n",
    "    m=12,\n",
    "    start_P=0,\n",
    "    seasonal=True,\n",
    "    D=1,\n",
    "    max_D=3,\n",
    "    trace=True,\n",
    "    error_action=\"ignore\",\n",
    "    suppress_warnings=True,\n",
    "    stepwise=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "auto_model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "auto_diag = auto_model.plot_diagnostics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "print(f\"auto-fit order: {auto_model.order}\")\n",
    "print(f\"auto-fit seasonal_order: {auto_model.seasonal_order}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Auto ARIMA forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "fc, conf_int = auto_model.predict(n_periods=6, return_conf_int=True)\n",
    "\n",
    "# Updates the existing model with a small number of MLE steps\n",
    "auto_model.update(test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "figure = px.line(x=df_monthly.index, y=df_monthly[var])\n",
    "figure.add_scatter(x=test.index, y=fc, mode=\"lines\", name=\"auto arima\")\n",
    "\n",
    "figure.add_scatter(\n",
    "    x=test.index, y=conf_int[:, 1], line=dict(width=0), mode=\"lines\", showlegend=False\n",
    ")\n",
    "figure.add_scatter(\n",
    "    x=test.index,\n",
    "    y=conf_int[:, 0],\n",
    "    line=dict(width=0),\n",
    "    mode=\"lines\",\n",
    "    fillcolor=\"rgba(68, 68, 68, 0.3)\",\n",
    "    fill=\"tonexty\",\n",
    "    showlegend=False,\n",
    ")\n",
    "\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=1000,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "mse_autoar = mean_squared_error(test, fc)\n",
    "print(f\"Mean squared error: {mse_autoar:.5f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "future_auto, conf_int_future = auto_model.predict(n_periods=6, return_conf_int=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "start = df_monthly.index.max() + pd.DateOffset(months=1)\n",
    "new_index = pd.date_range(start, periods=10, freq=\"M\")\n",
    "figure = px.line(x=df_monthly.index, y=df_monthly[var])\n",
    "figure.add_scatter(x=new_index, y=future_auto, mode=\"lines\", name=\"auto arima\")\n",
    "\n",
    "\n",
    "figure.add_scatter(\n",
    "    x=new_index,\n",
    "    y=conf_int_future[:, 1],\n",
    "    line=dict(width=0),\n",
    "    mode=\"lines\",\n",
    "    showlegend=False,\n",
    ")\n",
    "figure.add_scatter(\n",
    "    x=new_index,\n",
    "    y=conf_int_future[:, 0],\n",
    "    line=dict(width=0),\n",
    "    mode=\"lines\",\n",
    "    fillcolor=\"rgba(68, 68, 68, 0.3)\",\n",
    "    fill=\"tonexty\",\n",
    "    showlegend=False,\n",
    ")\n",
    "\n",
    "figure.update_layout(\n",
    "    title=f\"Monthly 99% percentile {var}\",\n",
    "    width=1000,\n",
    "    xaxis_title=\"Time\",\n",
    "    yaxis_title=f\"{var}\",\n",
    ")\n",
    "figure.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}