15_theta_starter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a24421d1",
   "metadata": {},
   "source": [
    "# Theta model "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65c5226b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import statsmodels.api as sm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "419acc7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams[\"figure.figsize\"] = (9,6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4a81c6e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc50336f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['co2'])\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('CO2 concentration (ppmw)')\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33c7bd00",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "080dd290",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2b2da04",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['co2'])\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('CO2 concentration (ppmw)')\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c501b13",
   "metadata": {},
   "source": [
    "## Modeling "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53b54b32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train/test split (keep the last 2 years of data for the test set)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19864972",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.forecasting.theta import ThetaModel\n",
    "\n",
    "def rolling_forecast(df: pd.DataFrame, train_len: int, horizon: int, window: int, period: int, method: str) -> list:\n",
    "    \n",
    "    seasonal_steps = int((window/period))\n",
    "    \n",
    "    TOTAL_LEN = train_len + horizon\n",
    "\n",
    "    if 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",
    "    elif method == 'theta':\n",
    "        # Get predictions from Theta\n",
    "\n",
    "        return pred_theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57f6667b",
   "metadata": {},
   "outputs": [],
   "source": [
    "TRAIN_LEN = len(train)\n",
    "HORIZON = len(test)\n",
    "PERIOD = 52\n",
    "WINDOW = 52\n",
    "\n",
    "pred_last_season = rolling_forecast(df, TRAIN_LEN, HORIZON, WINDOW, PERIOD, 'last_season')\n",
    "pred_theta = rolling_forecast(df, TRAIN_LEN, HORIZON, WINDOW, PERIOD, 'theta')\n",
    "\n",
    "test = test.copy()\n",
    "\n",
    "test.loc[:, 'pred_last_season'] = pred_last_season\n",
    "test.loc[:, 'pred_theta'] = pred_theta\n",
    "\n",
    "test.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0139d13e",
   "metadata": {},
   "source": [
    "## Testing triple exponential smoothing "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8c9332b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.holtwinters import ExponentialSmoothing\n",
    "\n",
    "def rolling_forecast(df: pd.DataFrame, train_len: int, horizon: int, window: int, period: int, method: str) -> list:\n",
    "    \n",
    "    seasonal_steps = int((window/period))\n",
    "    \n",
    "    TOTAL_LEN = train_len + horizon\n",
    "\n",
    "    if 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",
    "    elif method == 'theta':\n",
    "        pred_theta = []\n",
    "        \n",
    "        for i in range(train_len, TOTAL_LEN, window):\n",
    "            tm = ThetaModel(endog=df[:i], period=52)\n",
    "            res = tm.fit()\n",
    "            predictions = res.forecast(window)\n",
    "            pred_theta.extend(predictions)\n",
    "            \n",
    "\n",
    "        return pred_theta\n",
    "            \n",
    "    elif method == 'tes':\n",
    "        # Get predictions from TES\n",
    "\n",
    "        return pred_tes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "121bfae9",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_tes = rolling_forecast(df, TRAIN_LEN, HORIZON, WINDOW, PERIOD, 'tes')\n",
    "\n",
    "test.loc[:, 'pred_tes'] = pred_tes\n",
    "\n",
    "test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a1af69f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['co2'])\n",
    "ax.plot(test['co2'], 'b-', label='actual')\n",
    "ax.plot(test['pred_last_season'], 'r:', label='baseline')\n",
    "ax.plot(test['pred_theta'], 'g-.', label='Theta')\n",
    "ax.plot(test['pred_tes'], 'k--', label='TES')\n",
    "\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('CO2 concentration (ppmv)')\n",
    "ax.axvspan('2000-01-08', '2001-12-29', color='#808080', alpha=0.2)\n",
    "\n",
    "ax.legend(loc='best')\n",
    "\n",
    "ax.set_xlim('1998-03-07', '2001-12-29')\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b03a6c5a",
   "metadata": {},
   "source": [
    "## Evaluation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22aac8a2",
   "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": "ef22a629",
   "metadata": {},
   "outputs": [],
   "source": [
    "mape_baseline = mape(test['co2'], test['pred_last_season'])\n",
    "mape_theta = mape(test['co2'], test['pred_theta'])\n",
    "mape_tes = mape(test['co2'], test['pred_tes'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d9831bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "x = ['Baseline', 'Theta', 'TES']\n",
    "y = [mape_baseline, mape_theta, mape_tes]\n",
    "\n",
    "ax.bar(x, y, width=0.4)\n",
    "ax.set_xlabel('Exponential smoothing models')\n",
    "ax.set_ylabel('MAPE (%)')\n",
    "ax.set_ylim(0, 1)\n",
    "\n",
    "for index, value in enumerate(y):\n",
    "    plt.text(x=index, y=value + 0.05, s=str(value), ha='center')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5250a371",
   "metadata": {},
   "source": [
    "## Testing SARIMA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d1f718c",
   "metadata": {},
   "outputs": [],
   "source": [
    "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"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cad2d44",
   "metadata": {},
   "source": [
    "### Test for stationarity "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eeb4617",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run the ADF test\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": "19a2414f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Difference the series 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": "56c4648f",
   "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",
    "    all_P = range(min_P, max_P+1, 1)\n",
    "    all_Q = range(min_Q, max_Q+1, 1)\n",
    "    \n",
    "    all_orders = list(product(all_p, all_q, all_P, all_Q))\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",
    "            model = SARIMAX(\n",
    "                endog, \n",
    "                order=(order[0], d, order[1]),\n",
    "                seasonal_order=(order[2], D, order[3], s)).fit(disp=False)\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": "43298bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_p = 1\n",
    "max_p = 3\n",
    "min_q = 1\n",
    "max_q = 3\n",
    "\n",
    "min_P = 1\n",
    "max_P = 1\n",
    "min_Q = 1\n",
    "max_Q = 1\n",
    "\n",
    "d = 1\n",
    "D = 0\n",
    "s = 52\n",
    "\n",
    "result_df = SARIMA_gridsearch(train, 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": "7ee5c8e0",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9b559ebf",
   "metadata": {},
   "source": [
    "### Residuals analysis "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c193cef8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17eccb2e",
   "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": "a8973ee5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "d1391226",
   "metadata": {},
   "source": [
    "### Forecasting "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d95c7d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from statsmodels.tsa.holtwinters import ExponentialSmoothing\n",
    "\n",
    "def rolling_forecast(df: pd.DataFrame, train_len: int, horizon: int, window: int, period: int, method: str) -> list:\n",
    "    \n",
    "    seasonal_steps = int((window/period))\n",
    "    \n",
    "    TOTAL_LEN = train_len + horizon\n",
    "\n",
    "    if 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",
    "    elif method == 'theta':\n",
    "        pred_theta = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            tm = ThetaModel(endog=df[:i], period=52)\n",
    "            res = tm.fit()\n",
    "            predictions = res.forecast(window)\n",
    "            pred_theta.extend(predictions)\n",
    "            \n",
    "\n",
    "        return pred_theta\n",
    "            \n",
    "    elif method == 'tes':\n",
    "        pred_tes = []\n",
    "        \n",
    "        for i in range(train_len, total_len, window):\n",
    "            tes = ExponentialSmoothing(\n",
    "                df[:i],\n",
    "                trend='add',\n",
    "                seasonal='add',\n",
    "                seasonal_periods=52,\n",
    "                initialization_method='estimated'\n",
    "            ).fit()\n",
    "            \n",
    "            predictions = tes.forecast(window)\n",
    "            pred_tes.extend(predictions)\n",
    "\n",
    "        return pred_tes\n",
    "\n",
    "    elif method == 'SARIMA':\n",
    "        # Get predictions from SARIMA\n",
    "            \n",
    "        return pred_SARIMA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "047fec03",
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_SARIMA = rolling_forecast(df, TRAIN_LEN, HORIZON, WINDOW, PERIOD, 'SARIMA')\n",
    "\n",
    "test.loc[:, 'pred_SARIMA'] = pred_SARIMA\n",
    "\n",
    "test.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "132d3daa",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(df['co2'])\n",
    "ax.plot(test['co2'], 'b-', label='actual')\n",
    "ax.plot(test['pred_last_season'], 'r:', label='baseline')\n",
    "ax.plot(test['pred_theta'], 'g-.', label='Theta')\n",
    "ax.plot(test['pred_tes'], 'k--', label='TES')\n",
    "ax.plot(test['pred_SARIMA'], label='SARIMA')\n",
    "\n",
    "ax.set_xlabel('Time')\n",
    "ax.set_ylabel('CO2 concentration (ppmv)')\n",
    "ax.axvspan('2000-01-08', '2001-12-29', color='#808080', alpha=0.2)\n",
    "\n",
    "ax.legend(loc='best')\n",
    "\n",
    "ax.set_xlim('1998-03-07', '2001-12-29')\n",
    "\n",
    "fig.autofmt_xdate()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fbe1d85",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "x = ['Baseline', 'Theta', 'TES', 'SARIMA']\n",
    "y = [mape_baseline, mape_theta, mape_tes, mape_SARIMA]\n",
    "\n",
    "ax.bar(x, y, width=0.4)\n",
    "ax.set_xlabel('Exponential smoothing models')\n",
    "ax.set_ylabel('MAPE (%)')\n",
    "ax.set_ylim(0, 1)\n",
    "\n",
    "for index, value in enumerate(y):\n",
    "    plt.text(x=index, y=value + 0.05, s=str(value), ha='center')\n",
    "    \n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15bfa4ad",
   "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
}