Pharma Sales Time Series Forecasting_stage 14 - without outputs.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "faEKk8GqjpBW"
      },
      "source": [
        "Data Source: https://www.kaggle.com/datasets/milanzdravkovic/pharma-sales-data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZWsehj6Z7AZm"
      },
      "source": [
        "## Libraries"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "BVlEsru-6QH-"
      },
      "outputs": [],
      "source": [
        "import pandas as pd\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import re\n",
        "import math\n",
        "import calendar\n",
        "import calendar\n",
        "from datetime import datetime"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tZzc3EWn8bXb"
      },
      "outputs": [],
      "source": [
        "from sklearn import linear_model\n",
        "from sklearn.cluster import KMeans\n",
        "from sklearn.decomposition import PCA\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.preprocessing import StandardScaler\n",
        "from sklearn.preprocessing import PolynomialFeatures\n",
        "from sklearn.svm import SVR\n",
        "from sklearn.ensemble import VotingRegressor"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mMrgw5YAIttm",
        "outputId": "b0b922b7-4807-4e52-cec7-ba2e23ff009a"
      },
      "outputs": [],
      "source": [
        "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",
        "from statsmodels.tsa.holtwinters import ExponentialSmoothing"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Jnssy5jk8bUf"
      },
      "outputs": [],
      "source": [
        "from pylab import rcParams\n",
        "rcParams['figure.figsize'] = 12, 8\n",
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X8BsDslQ8bRm"
      },
      "outputs": [],
      "source": [
        "np.set_printoptions(precision=2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PdcdK_Au4m0w",
        "outputId": "287db865-3b1f-410b-a8d8-bf529a2e553f"
      },
      "outputs": [],
      "source": [
        "!pip install mlxtend"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oX45WDXG4oNd"
      },
      "outputs": [],
      "source": [
        "from pathlib import Path\n",
        "import heapq\n",
        "from collections import defaultdict\n",
        "import pandas as pd\n",
        "import matplotlib.pylab as plt\n",
        "from mlxtend.frequent_patterns import apriori\n",
        "from mlxtend.frequent_patterns import association_rules"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mnm8okxwSOwY",
        "outputId": "c32e1198-0cf2-41d8-dc69-621c88b17210"
      },
      "outputs": [],
      "source": [
        "!pip install dmba"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "L31e29P-SNx7",
        "outputId": "a5495cd7-1fe1-427f-e14d-7f897379033d"
      },
      "outputs": [],
      "source": [
        "from dmba import regressionSummary"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8joMEHanejE8",
        "outputId": "eb898ae3-6ff9-4c86-e5c0-b6cbf450a720"
      },
      "outputs": [],
      "source": [
        "!pip install --upgrade plotly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_VXhRC9xep_E",
        "outputId": "4ea9e39f-52fa-4ba5-8700-de30d6105533"
      },
      "outputs": [],
      "source": [
        "!pip install prophet"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xvIdEkxreW3C"
      },
      "outputs": [],
      "source": [
        "%load_ext rpy2.ipython\n",
        "%matplotlib inline\n",
        "\n",
        "from prophet import Prophet\n",
        "import pandas as pd\n",
        "import logging\n",
        "import warnings\n",
        "\n",
        "logging.getLogger('prophet').setLevel(logging.ERROR)\n",
        "logging.getLogger('numexpr').setLevel(logging.ERROR)\n",
        "warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8uFvWh0V8o9c",
        "outputId": "f5eda15a-3dd3-4f00-f6fd-809d9635a0ce"
      },
      "outputs": [],
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YWkrKY2R_1D9"
      },
      "source": [
        "## Load data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "sVMlmT9F8rnI"
      },
      "outputs": [],
      "source": [
        "df_sales_daily = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesdaily.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ppQnh_AY_2Ya"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesmonthly.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KcYUg2m__2TL"
      },
      "outputs": [],
      "source": [
        "df_sales_hourly = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/saleshourly.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DB-EBlIgluBF"
      },
      "outputs": [],
      "source": [
        "df_sales_weekly = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesweekly.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "cKUqbaLO-ePk",
        "outputId": "8a2e0a35-979a-4cee-b8da-5941017905af"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r1BRKMVjdNAN"
      },
      "source": [
        "## Data Analysis"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gammoab1ZqFp"
      },
      "source": [
        "### M01AB drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Fb0JLuc3Zpu_"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB = df_sales_daily[['M01AB', 'Weekday Name']]\n",
        "result_M01AB = df_sales_daily_M01AB.groupby(['Weekday Name'], as_index=False).sum().sort_values('M01AB', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "pCdkUMtXFluo",
        "outputId": "2753f24c-8d8f-4dbf-9c0b-f5a28ad02c93"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9EiUnxOpFoh1"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2wuPxmaGZr0O"
      },
      "outputs": [],
      "source": [
        "result_day_M01AB = result_M01AB.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-4RvHCAzaDK1"
      },
      "outputs": [],
      "source": [
        "result_value_M01AB = round(result_M01AB.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3h4Bi6DFZrw_",
        "outputId": "4626f55b-ecf6-4032-9dfd-593b6a1fa074"
      },
      "outputs": [],
      "source": [
        "print('The 1st drug, M01AB, was most frequently sold on ' + str(result_day_M01AB) \n",
        "      + ' with the volume of ' + str(result_value_M01AB))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a-nCqoYL_pGA"
      },
      "source": [
        "### M01AE drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3BwQppR--hIB"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AE = df_sales_daily[['M01AE', 'Weekday Name']]\n",
        "result_M01AE = df_sales_daily_M01AE.groupby(['Weekday Name'], as_index=False).sum().sort_values('M01AE', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p4zO4USQ-9At"
      },
      "outputs": [],
      "source": [
        "result_day_M01AE = result_M01AE.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zLxy5DHv_LKV"
      },
      "outputs": [],
      "source": [
        "result_value_M01AE = round(result_M01AE.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uVe0SZHQ_RlJ",
        "outputId": "362cad76-fe19-45ce-e075-256bb81a4b7f"
      },
      "outputs": [],
      "source": [
        "print('The 2nd drug, M01AE, was most frequently sold on ' + str(result_day_M01AE) \n",
        "      + ' with the volume of ' + str(result_value_M01AE))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nXIky_QeaQJ8"
      },
      "source": [
        "### N02BA drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c0406ljYaP6F"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_N02BA = df_sales_daily[['N02BA', 'Weekday Name']]\n",
        "result_N02BA = df_sales_daily_N02BA.groupby(['Weekday Name'], as_index=False).sum().sort_values('N02BA', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-KncC88WaYzI"
      },
      "outputs": [],
      "source": [
        "result_day_N02BA = result_N02BA.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6ewcA0vCaP3O"
      },
      "outputs": [],
      "source": [
        "result_value_N02BA = round(result_N02BA.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "B6-ijUrjaP0V",
        "outputId": "3cc3a258-b2f4-4394-bf94-c854c2b6e307"
      },
      "outputs": [],
      "source": [
        "print('The 3rd drug, N02BA, was most frequently sold on ' + str(result_day_N02BA) \n",
        "      + ' with the volume of ' + str(result_value_N02BA))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8u4_MvtBanQB"
      },
      "source": [
        "### N02BE drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OnCm9l-naPu0"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_N02BE = df_sales_daily[['N02BE', 'Weekday Name']]\n",
        "result_N02BE = df_sales_daily_N02BE.groupby(['Weekday Name'], as_index=False).sum().sort_values('N02BE', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "427vZVuSaPsJ"
      },
      "outputs": [],
      "source": [
        "result_day_N02BE = result_N02BE.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qAZleuphaPpZ"
      },
      "outputs": [],
      "source": [
        "result_value_N02BE = round(result_N02BE.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kvV1YWEXaPmg",
        "outputId": "776f8fd0-7321-4471-c0df-55e73351d960"
      },
      "outputs": [],
      "source": [
        "print('The 4th drug, N02BE, was most frequently sold on ' + str(result_day_N02BE) \n",
        "      + ' with the volume of ' + str(result_value_N02BE))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gm4kYdefassQ"
      },
      "source": [
        "### N05B drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KclwI96raPkA"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_N05B = df_sales_daily[['N05B', 'Weekday Name']]\n",
        "result_N05B = df_sales_daily_N05B.groupby(['Weekday Name'], as_index=False).sum().sort_values('N05B', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9ZBjp2q5aPhQ"
      },
      "outputs": [],
      "source": [
        "result_day_N05B = result_N05B.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5DWuKV1daPer"
      },
      "outputs": [],
      "source": [
        "result_value_N05B = round(result_N05B.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ESfK9i_daPbg",
        "outputId": "3459b80f-e541-4fa0-c16c-3fab299473dc"
      },
      "outputs": [],
      "source": [
        "print('The 5th drug, N05B, was most frequently sold on ' + str(result_day_N05B) \n",
        "      + ' with the volume of ' + str(result_value_N05B))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qTkf2tKQa9KK"
      },
      "source": [
        "### N05C drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tNw4xAiCaPYw"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_N05C = df_sales_daily[['N05C', 'Weekday Name']]\n",
        "result_N05C = df_sales_daily_N05C.groupby(['Weekday Name'], as_index=False).sum().sort_values('N05C', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LtISgRo_aPWG"
      },
      "outputs": [],
      "source": [
        "result_day_N05C = result_N05C.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lqqkFIZsbL8R"
      },
      "outputs": [],
      "source": [
        "result_value_N05C = round(result_N05C.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jEcBArcJbMZx",
        "outputId": "72e84fea-9617-4c56-86da-20e86728a71e"
      },
      "outputs": [],
      "source": [
        "print('The 6th drug, N05C, was most frequently sold on ' + str(result_day_N05C) \n",
        "      + ' with the volume of ' + str(result_value_N05C))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K1MXmubHbWL6"
      },
      "source": [
        "### R03 drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dgQu-6SVbWfY"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_R03 = df_sales_daily[['R03', 'Weekday Name']]\n",
        "result_R03 = df_sales_daily_R03.groupby(['Weekday Name'], as_index=False).sum().sort_values('R03', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-PNMrCBhbXD3"
      },
      "outputs": [],
      "source": [
        "result_day_R03 = result_R03.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SAVCmGHjbW-F"
      },
      "outputs": [],
      "source": [
        "result_value_R03 = round(result_R03.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dpo5WwbxbW2g",
        "outputId": "eb515c7f-f647-4869-c41c-7c6ddeed1b0a"
      },
      "outputs": [],
      "source": [
        "print('The 7th drug, R03, was most frequently sold on ' + str(result_day_R03) \n",
        "      + ' with the volume of ' + str(result_value_R03))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PINqqJ9CbXec"
      },
      "source": [
        "### R06 drug by Weekday"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ctWjt566bXyO"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_R06 = df_sales_daily[['R06', 'Weekday Name']]\n",
        "result_R06 = df_sales_daily_R06.groupby(['Weekday Name'], as_index=False).sum().sort_values('R06', ascending=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-iR-IQ5PbYdl"
      },
      "outputs": [],
      "source": [
        "result_day_R06 = result_R06.iloc[0,0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gFllkZtpbYQJ"
      },
      "outputs": [],
      "source": [
        "result_value_R06 = round(result_R06.iloc[0,1], 2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dM7lmp5zbYG2",
        "outputId": "fff00a3d-e155-4d28-8fec-281685c9d4a7"
      },
      "outputs": [],
      "source": [
        "print('The 8th drug, R06, was most frequently sold on ' + str(result_day_R06) \n",
        "      + ' with the volume of ' + str(result_value_R06))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "70ifrnu5_9w5"
      },
      "source": [
        "## Top 3 drugs gaving the largest sales in January 2015, July 2016, September 2017"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "yR0507Kl_mBU",
        "outputId": "b2c0bb2a-f560-494d-bd3c-06ec4d2cfe1b"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1PueuyiIAcD6"
      },
      "outputs": [],
      "source": [
        "def top_3_drugs_by_month(month, year):\n",
        "    month = str(month) if (month > 9) else '0' + str(month)\n",
        "    year = str(year)\n",
        "\n",
        "    # Filter by Date\n",
        "    sales = df_sales_monthly.loc[df_sales_monthly['datum'].str.contains('^' + year + '\\-' + month + '' , flags=re.I, regex=True)]\n",
        "    \n",
        "    # Reset index\n",
        "    sales = sales.reset_index()\n",
        "    \n",
        "    # Filter chosen drugs\n",
        "    top_Sales_by_Product = sales[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']]\n",
        "    \n",
        "    # Sort values horizontally in descending order\n",
        "    top_Sales_by_Product = top_Sales_by_Product.sort_values(by=0, ascending=False, axis=1)\n",
        "    \n",
        "    # Print results\n",
        "    print('Top 3 Drugs by Sale in ' + calendar.month_name[int(month)] + ' ' + year)\n",
        "    for drug in top_Sales_by_Product.columns.values[0:3]:\n",
        "        print(' - Product: ' + str(drug) + ', Volume sold: ' + str(round(top_Sales_by_Product[drug].iloc[0], 2)))\n",
        "    print(\"\\n\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "VME9cFUeEGiv",
        "outputId": "e98124de-7f2c-4434-f09b-cd13a9f17627"
      },
      "outputs": [],
      "source": [
        "# Top 3 Drugs by Sale in January 2014\n",
        "top_3_drugs_by_month(1, 2014)\n",
        "\n",
        "# Top 3 Drugs by Sale in January 2015\n",
        "top_3_drugs_by_month(1, 2015)\n",
        "\n",
        "# Top 3 Drugs by Sale in July 2016\n",
        "top_3_drugs_by_month(7, 2016)\n",
        "\n",
        "# Top 3 Drugs by Sale in September 2017\n",
        "top_3_drugs_by_month(9, 2017)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oxMhSdmBFS-e"
      },
      "source": [
        "## Top 1 Drug sold on Mondays in years"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dsEa8UjaEIOc"
      },
      "outputs": [],
      "source": [
        "df_sales_daily = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesdaily.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "UD3UeKb1Hk0U",
        "outputId": "9065c947-8f3d-44f8-938b-d48e875a6cb5"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4XqkHqYDdbDj"
      },
      "source": [
        "### 2014"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7EiZV9KsdeY7"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2014 = df_sales_daily.loc[df_sales_daily['datum'].str.contains('2014', flags=re.I, regex=True) & (df_sales_daily['Weekday Name'] == 'Monday')]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v4dd02G3deJk"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2014 = df_sales_daily_2014.groupby(['Weekday Name'], as_index=False).sum()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "xkfkfc25deGj",
        "outputId": "5401ceae-889a-4320-e5ca-d75495b667b2"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2014 = df_sales_daily_2014[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']]\n",
        "result_2014 = df_sales_daily_2014.sort_values(by=0, ascending=False, axis=1)\n",
        "for field in result_2014.columns.values[0:1]:\n",
        "    print('The most frequently sold Drug on Mondays in 2014 is ' + str(field) + ' with the Volume of ' + str(round(result_2014[field].iloc[0], 2)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "afxWYMXxd2vn"
      },
      "source": [
        "### 2015"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HSBkAdjcd3Ap"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2015 = df_sales_daily.loc[df_sales_daily['datum'].str.contains('2015', flags=re.I, regex=True) & (df_sales_daily['Weekday Name'] == 'Monday')]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "oAJbELw1d3cc"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2015 = df_sales_daily_2015.groupby(['Weekday Name'], as_index=False).sum()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "J_DbWGOdd3VD",
        "outputId": "3f795ae6-10ae-483c-d610-52a94ca09d35"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2015 = df_sales_daily_2015[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']]\n",
        "result_2015 = df_sales_daily_2015.sort_values(by=0, ascending=False, axis=1)\n",
        "for field in result_2015.columns.values[0:1]:\n",
        "    print('The most frequently sold Drug on Mondays in 2015 is ' + str(field) + ' with the Volume of ' + str(round(result_2015[field].iloc[0], 2)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xqI1jHCyeHMV"
      },
      "source": [
        "### 2016"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5Ep94IrpeHfg"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2016 = df_sales_daily.loc[df_sales_daily['datum'].str.contains('2016', flags=re.I, regex=True) & (df_sales_daily['Weekday Name'] == 'Monday')]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-n2Qd2rdeH54"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2016 = df_sales_daily_2016.groupby(['Weekday Name'], as_index=False).sum()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Qxz6UEvveHxZ",
        "outputId": "abf06da3-3b54-4305-f49b-a1f369fedad4"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2016 = df_sales_daily_2016[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']]\n",
        "result_2016 = df_sales_daily_2016.sort_values(by=0, ascending=False, axis=1)\n",
        "for field in result_2016.columns.values[0:1]:\n",
        "    print('The most frequently sold Drug on Mondays in 2016 is ' + str(field) + ' with the Volume of ' + str(round(result_2016[field].iloc[0], 2)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LkwiTjXDc023"
      },
      "source": [
        "### 2017"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lm_mP0G_G7tH"
      },
      "source": [
        "#### Filter out all except Mondays in 2017"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bWJc00DEGykE"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2017 = df_sales_daily.loc[df_sales_daily['datum'].str.contains('2017', flags=re.I, regex=True) & (df_sales_daily['Weekday Name'] == 'Monday')]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nEX8-3TDHrwA"
      },
      "source": [
        "#### Group by Weekday and Sum up"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SZ4CcPV-Hd3n"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2017 = df_sales_daily_2017.groupby(['Weekday Name'], as_index=False).sum()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bev0ZfLoH1vu"
      },
      "source": [
        "#### Filter by the chosen Drugs and Sort by descending order"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dLbx4eJbHxHe"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_2017 = df_sales_daily_2017[['M01AB', 'M01AE', 'N02BA', 'N02BE', 'N05B', 'N05C', 'R03', 'R06']]\n",
        "result_2017 = df_sales_daily_2017.sort_values(by=0, ascending=False, axis=1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UvQV3pnGH0-e",
        "outputId": "b9d465bc-9248-4237-ad94-0f16f93243cb"
      },
      "outputs": [],
      "source": [
        "for field in result_2017.columns.values[0:1]:\n",
        "    print('The most frequently sold Drug on Mondays in 2017 is ' + str(field) + ' with the Volume of ' + str(round(result_2017[field].iloc[0], 2)))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OdGaFs68vACg"
      },
      "source": [
        "## Pharmaceutical Portfolio Management Analysis"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-be-IzABvIte"
      },
      "source": [
        "### Volatility and Risk\n",
        "\n",
        "The volatility is measured by the average squared deviation from the mean, which is the standard deviation."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 250
        },
        "id": "IRgJjcCSvKXX",
        "outputId": "e93d8091-cf4e-4aa5-e436-5c389084d535"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "IXTwxw4svKUU",
        "outputId": "0bc3be4a-8c55-460f-d3a6-8453c50da69c"
      },
      "outputs": [],
      "source": [
        "returns = df_sales_daily[[\"M01AB\", \"M01AE\", \"N02BA\", \"N02BE\", \"N05B\", \"N05C\", \"R03\", \"R06\"]].pct_change()\n",
        "returns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "ypyukTEwvKQx",
        "outputId": "43d36547-9dcb-4a87-f07f-95115fd3e6d7"
      },
      "outputs": [],
      "source": [
        "returns = returns.dropna()\n",
        "returns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 533
        },
        "id": "UkcfmBeuvqTu",
        "outputId": "8158fa07-d23b-4028-a65d-ee7445e28a21"
      },
      "outputs": [],
      "source": [
        "returns.replace([np.inf, -np.inf], 0, inplace=True)\n",
        "returns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "uDFB5d5T3kff",
        "outputId": "026be13d-3f0f-4b7c-954f-e77e5e512dac"
      },
      "outputs": [],
      "source": [
        "returns.info()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 364
        },
        "id": "BMMm3qVv3mJT",
        "outputId": "a2d56dc9-5256-4970-fa9b-6c13fe30c259"
      },
      "outputs": [],
      "source": [
        "returns.describe()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "luxiT13lv7TA",
        "outputId": "56d470a4-9d7f-4839-c0f8-f8dab13dc0f0"
      },
      "outputs": [],
      "source": [
        "deviations = returns - returns.mean()\n",
        "squared_deviations = deviations**2\n",
        "mean_squared_deviations = squared_deviations.mean()\n",
        "\n",
        "volatility = np.sqrt(mean_squared_deviations)\n",
        "volatility"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WwzkUfcEv7Py",
        "outputId": "318d5eef-14d6-451f-bcbe-308291d22c0b"
      },
      "outputs": [],
      "source": [
        "returns.std()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "wCyho24hv7Me",
        "outputId": "53eb641c-9d28-49c3-a48e-c20595feb1dd"
      },
      "outputs": [],
      "source": [
        "returns.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "D-a51niuv7JH",
        "outputId": "c4974971-6232-464c-a066-07c962167de1"
      },
      "outputs": [],
      "source": [
        "number_of_obs = returns.shape[0]\n",
        "mean_squared_deviations = squared_deviations.sum()/(number_of_obs-1)\n",
        "volatility = np.sqrt(mean_squared_deviations)\n",
        "volatility"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HPwSvHvcwFU0",
        "outputId": "4519d06d-e91c-45e4-ac24-ba303fecc8ed"
      },
      "outputs": [],
      "source": [
        "returns.std()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "C37YyW4Z3PvO"
      },
      "source": [
        "### Annualizing Volatility\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PsL4R0mO3UhV"
      },
      "source": [
        "I annualize volatility by scaling (multiplying) it by the square root of the number of periods per observation\n",
        "\n",
        "Therefore, to annualize the volatility of a monthly series, I muiltiply it by the square root of 12. Instead of using the `np.sqrt()` I can raise it to the power of $0.5$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IHStZCQlwHJS",
        "outputId": "dac1154d-d318-452c-d654-ae367dd795b6"
      },
      "outputs": [],
      "source": [
        "annualized_vol = returns.std()*(12**0.5)\n",
        "annualized_vol"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "naI-KkyK3FX3"
      },
      "source": [
        "### Risk Adjusted Returns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "2oBWWiXHx41l"
      },
      "outputs": [],
      "source": [
        "returns_plot = returns/100"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "zpUX70yVyDyT",
        "outputId": "0414828a-be9b-4a0f-f013-607989a21e70"
      },
      "outputs": [],
      "source": [
        "returns_plot"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 364
        },
        "id": "pFbG-cr_yFHL",
        "outputId": "6feb5600-02d3-4200-94f7-a960ff0264d7"
      },
      "outputs": [],
      "source": [
        "returns_plot.describe()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 358
        },
        "id": "mhLz2aILwHF2",
        "outputId": "f074f4f8-7c1f-4452-cf9a-16a6e3e4d7ea"
      },
      "outputs": [],
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (45,5)\n",
        "returns_plot.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pKXr45gMwHCq",
        "outputId": "a711b553-1dc4-4499-e8ae-d8b317d44cb0"
      },
      "outputs": [],
      "source": [
        "annualized_vol = returns.std()*np.sqrt(12)\n",
        "annualized_vol"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nZpqgIQwz4Rs",
        "outputId": "79327e5b-3ad9-4e93-e31e-f94e07e05e96"
      },
      "outputs": [],
      "source": [
        "returns.shape[0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gqz8yzp60Suf",
        "outputId": "4bc6c84d-82b5-41b9-e6ea-4f87b6abf79a"
      },
      "outputs": [],
      "source": [
        "annualized_vol = returns.std()*(12**0.5)\n",
        "annualized_vol"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "YjStRxDrwG_O",
        "outputId": "5da879ae-4082-43a1-b7a8-53d3ba127662"
      },
      "outputs": [],
      "source": [
        "n_months = 3\n",
        "return_per_month = (returns+1).prod()**(1/n_months) - 1\n",
        "return_per_month"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "foBZVkGNzy3c"
      },
      "outputs": [],
      "source": [
        "annualized_return = (return_per_month + 1)**12-1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "IP8gC7zzzy0s",
        "outputId": "c3bdbc71-abd4-46d0-85cc-a9cb55aa855b"
      },
      "outputs": [],
      "source": [
        "annualized_return = (returns+1).prod()**(12/n_months) - 1\n",
        "annualized_return"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0FUe5r4pzyxd",
        "outputId": "afc3dc2b-26d2-4e3b-c08c-a29f0635d35c"
      },
      "outputs": [],
      "source": [
        "annualized_return/annualized_vol"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "mFp8Me2J02Qj",
        "outputId": "4ed8c4ef-a8c3-41bc-d6d6-00bcde40ad96"
      },
      "outputs": [],
      "source": [
        "riskfree_rate = 0.03\n",
        "excess_return = annualized_return - riskfree_rate\n",
        "sharpe_ratio = excess_return/annualized_vol\n",
        "sharpe_ratio"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 367
        },
        "id": "6Knjri2Q2IMb",
        "outputId": "5d7cf33a-3a44-4f30-bd08-4c9a299b682d"
      },
      "outputs": [],
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (10,5)\n",
        "sharpe_ratio.plot.bar()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "diGN1Nah3wMj"
      },
      "source": [
        "### Computing Drawdowns\n",
        "\n",
        "1. Convert the time series of returns to a time series that represents a wealth index\n",
        "2. Compute a time series of the previous peaks\n",
        "3. Compute the Drawdown as the difference between the previous peak and the current value\n",
        "\n",
        "Let's do this for M01AB stocks."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 468
        },
        "id": "Q20sZ-N74b4Q",
        "outputId": "38a55602-6389-4949-9ed4-d25651af25ad"
      },
      "outputs": [],
      "source": [
        "df_sales_daily"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "nuo1qQ1d4Dgq",
        "outputId": "1d0280fe-d909-4be4-c86d-80aa6b616f37"
      },
      "outputs": [],
      "source": [
        "rets = df_sales_daily[[\"datum\", \"M01AB\", \"M01AE\", \"N02BA\", \"N02BE\", \"N05B\", \"N05C\", \"R03\", \"R06\"]]\n",
        "rets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 490
        },
        "id": "JREbtyO141Tq",
        "outputId": "d86100c1-97fc-4d56-d220-b5c95aa60893"
      },
      "outputs": [],
      "source": [
        "rets['Date'] = pd.to_datetime(rets.datum, format='%m/%d/%Y')\n",
        "rets.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 551
        },
        "id": "teD6dLF-66Au",
        "outputId": "05eec619-fa64-48a9-e7a6-5b4027054787"
      },
      "outputs": [],
      "source": [
        "rets['Date 1'] = rets['Date'].dt.strftime('%Y%m%d')\n",
        "rets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "X5UcmBrm658D"
      },
      "outputs": [],
      "source": [
        "date_temp = rets[[\"Date 1\"]]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "UgCbb5iM4bEa",
        "outputId": "2fdabcc8-7bd6-4f55-a81f-c876ceb1a64c"
      },
      "outputs": [],
      "source": [
        "rets = rets[[\"M01AB\", \"M01AE\", \"N02BA\", \"N02BE\", \"N05B\", \"N05C\", \"R03\", \"R06\"]]\n",
        "rets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "b0P6lzc_7p0O"
      },
      "outputs": [],
      "source": [
        "rets = rets.pct_change()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "2etG2WgO7vL2",
        "outputId": "5d1e9abf-32fa-4056-b827-6c70739af2b7"
      },
      "outputs": [],
      "source": [
        "rets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RLatOnEt9CBp"
      },
      "outputs": [],
      "source": [
        "rets = rets.iloc[1: , :]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 533
        },
        "id": "gyqACfnE8AkQ",
        "outputId": "7cb12be9-048a-4635-cacb-945181217b28"
      },
      "outputs": [],
      "source": [
        "rets.replace([np.inf, -np.inf, np.nan], 0, inplace=True)\n",
        "rets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "cveSzhTD8RiK",
        "outputId": "76887bae-010a-4e3b-9984-78d195a0ff1f"
      },
      "outputs": [],
      "source": [
        "date_temp"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "YcnVLq7S8Cqd",
        "outputId": "6811a83e-251a-45ef-e28a-bda2d300b201"
      },
      "outputs": [],
      "source": [
        "rets2 = pd.concat([rets, date_temp], axis=1, join=\"inner\")\n",
        "rets2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4fxcXlUw9pPW",
        "outputId": "837b5217-bb18-4716-ede8-275ab99890b7"
      },
      "outputs": [],
      "source": [
        "rets2.info()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 641
        },
        "id": "GKDc1SVb-Coo",
        "outputId": "715077d5-43d9-48cb-c379-b44136a04672"
      },
      "outputs": [],
      "source": [
        "rets2['Date']=pd.to_datetime(rets2['Date 1'],format='%Y%m%d')\n",
        "rets2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 693
        },
        "id": "wVkhcM9E_amh",
        "outputId": "4a76349f-c849-4173-ba7b-7fd37f1b50ae"
      },
      "outputs": [],
      "source": [
        "rets2 = rets2.groupby([rets2.Date.dt.year, rets2.Date.dt.month]).last()\n",
        "rets2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 499
        },
        "id": "lyaBzmJM-Chw",
        "outputId": "618dd26a-c642-472b-afef-c4c7c7d2b39d"
      },
      "outputs": [],
      "source": [
        "rets2 = rets2.set_index('Date')\n",
        "rets2"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eyHOSYt9_91U"
      },
      "source": [
        "#### M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "94JZYZRL3v1M",
        "outputId": "7fc795f4-a027-4c73-ab12-f2d9b82e235a"
      },
      "outputs": [],
      "source": [
        "wealth_index = 1000*(1+rets2[\"M01AB\"]).cumprod()\n",
        "wealth_index.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "qff2yxer3vx7",
        "outputId": "de28b3f9-3c16-4810-cf89-76a728b0ef29"
      },
      "outputs": [],
      "source": [
        "previous_peaks = wealth_index.cummax()\n",
        "previous_peaks.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "MfY3Jd413vvH",
        "outputId": "6353c256-8ba5-45a4-d6f2-a38b92a952b0"
      },
      "outputs": [],
      "source": [
        "drawdown = (wealth_index - previous_peaks)/previous_peaks\n",
        "drawdown.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ZVaULXQp3vr3",
        "outputId": "e69a83f2-c291-4fa9-d439-c2be0471575c"
      },
      "outputs": [],
      "source": [
        "drawdown.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "mEd3qR1V3von",
        "outputId": "c36e2ee6-46ab-4ff0-c084-c34e48427559"
      },
      "outputs": [],
      "source": [
        "drawdown[\"2014\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "acKpF6KgAhYV",
        "outputId": "b5b1ed42-ad06-455c-a509-d5e98ea30591"
      },
      "outputs": [],
      "source": [
        "drawdown_M01AE[\"2014\":].min()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dxvlPmUSABYI"
      },
      "source": [
        "#### M01AE"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 347
        },
        "id": "K4m_2yrB_80j",
        "outputId": "b684a507-70a2-4650-c577-ceca74f66282"
      },
      "outputs": [],
      "source": [
        "wealth_index_M01AE = 1000*(1+rets2[\"M01AE\"]).cumprod()\n",
        "wealth_index_M01AE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 347
        },
        "id": "rzvZVzNn_8xe",
        "outputId": "a0146c66-0453-4e41-a0f1-7535798220ea"
      },
      "outputs": [],
      "source": [
        "previous_peaks_M01AE = wealth_index_M01AE.cummax()\n",
        "previous_peaks_M01AE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "K8vGPEzs_8uc",
        "outputId": "e479e949-68f5-40e1-ea32-b8906acb7140"
      },
      "outputs": [],
      "source": [
        "drawdown_M01AE = (wealth_index_M01AE - previous_peaks_M01AE)/previous_peaks_M01AE\n",
        "drawdown_M01AE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ptRn0pF5AP7s",
        "outputId": "92d94755-67df-42db-a45e-b9e1b2ba6e24"
      },
      "outputs": [],
      "source": [
        "drawdown_M01AE.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "K_4coutXAPy3",
        "outputId": "41180316-229e-4015-985a-5f1a9a027f38"
      },
      "outputs": [],
      "source": [
        "drawdown_M01AE[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "7CSrpgWrAXvV",
        "outputId": "13f53751-5a26-407e-d015-dc12b883d426"
      },
      "outputs": [],
      "source": [
        "drawdown_M01AE[\"2016\":].min()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tSWTWczMAoFK"
      },
      "source": [
        "#### N02BA"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "zL61I6BLAwC4",
        "outputId": "0c9a11a0-8afd-446e-8bf5-a3cebc083eb7"
      },
      "outputs": [],
      "source": [
        "wealth_index_N02BA = 1000*(1+rets2[\"N02BA\"]).cumprod()\n",
        "wealth_index_N02BA.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "JT_Be6wNA1tR",
        "outputId": "403490d6-8ca0-40db-ae3d-b5f7285efb2c"
      },
      "outputs": [],
      "source": [
        "previous_peaks_N02BA = wealth_index_N02BA.cummax()\n",
        "previous_peaks_N02BA.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "UFTGz7g8A7yj",
        "outputId": "7cb3a91f-3a76-4c1f-b66c-039251295a05"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BA = (wealth_index_N02BA - previous_peaks_N02BA)/previous_peaks_N02BA\n",
        "drawdown_N02BA.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dpIXJdj7BAMt",
        "outputId": "26a47f90-9745-4412-d985-fddc2b2eb83c"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BA.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "t8JTfNjSBHtK",
        "outputId": "aa453474-3b75-4776-9e6e-2f2345514f9d"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BA[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TZTOKSQ4BJ5n",
        "outputId": "5956ef68-c0d8-4464-967f-ceea2e33aeb3"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BA[\"2017\":].max()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AiJnBEn_Bjw1"
      },
      "source": [
        "#### N02BE"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "gkvXL_CxAxze",
        "outputId": "d1bc093b-adb6-42d6-c3d0-d789202d4b75"
      },
      "outputs": [],
      "source": [
        "wealth_index_N02BE = 1000*(1+rets2[\"N02BE\"]).cumprod()\n",
        "wealth_index_N02BE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "HC2IxLDHA2E7",
        "outputId": "8315a9c6-4a78-4955-d862-f5b184ef208b"
      },
      "outputs": [],
      "source": [
        "previous_peaks_N02BE = wealth_index_N02BE.cummax()\n",
        "previous_peaks_N02BE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "j4nFwLSTA8Na",
        "outputId": "dc416aca-298b-4148-c705-a7274507fadb"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BE = (wealth_index_N02BE - previous_peaks_N02BE)/previous_peaks_N02BE\n",
        "drawdown_N02BE.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pXuGhz9yBBur",
        "outputId": "f0f49afd-4560-4566-df23-49489fe55241"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BE.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "RcvFou8DBGJe",
        "outputId": "42361883-1585-4513-beb1-ef5ac30b8344"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BE[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "6HRpENODBKrv",
        "outputId": "9bc00202-0930-4e61-e7a8-1a60bb0fcd91"
      },
      "outputs": [],
      "source": [
        "drawdown_N02BE[\"2018\":].min()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KshXT7ClCALd"
      },
      "source": [
        "#### N05B"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 347
        },
        "id": "Zjgz_jRqAxur",
        "outputId": "0adb8d81-2f24-4cda-df44-9aa18cdaf30c"
      },
      "outputs": [],
      "source": [
        "wealth_index_N05B = 1000*(1+rets2[\"N05B\"]).cumprod()\n",
        "wealth_index_N05B.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 347
        },
        "id": "SOLMVulXA2ce",
        "outputId": "ac2c7453-d0dc-461d-9c27-4e1c4b6df87c"
      },
      "outputs": [],
      "source": [
        "previous_peaks_N05B = wealth_index_N05B.cummax()\n",
        "previous_peaks_N05B.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "Q4IlsMNpA8oP",
        "outputId": "6547c018-eb5a-43a9-beb8-c81111ccb53d"
      },
      "outputs": [],
      "source": [
        "drawdown_N05B = (wealth_index_N05B - previous_peaks_N05B)/previous_peaks_N05B\n",
        "drawdown_N05B.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "a--scTV8BCRm",
        "outputId": "a9586b61-9b3f-44cd-fd3f-2974fa059d5f"
      },
      "outputs": [],
      "source": [
        "drawdown_N05B.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "kSq13lMgBFXX",
        "outputId": "83042881-0aac-4ced-c4f8-d723bead69d9"
      },
      "outputs": [],
      "source": [
        "drawdown_N05B[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "KBnZCaGIBLRx",
        "outputId": "6cb975d6-b531-465e-e6b5-c73b468abb5e"
      },
      "outputs": [],
      "source": [
        "drawdown_N05B[\"2017\":].max()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "A8uZC3piCYy0"
      },
      "source": [
        "#### R03"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "UxO5fhDRAxq5",
        "outputId": "3824c648-5ec1-4402-f9e0-8e610adae16f"
      },
      "outputs": [],
      "source": [
        "wealth_index_R03 = 1000*(1+rets2[\"R03\"]).cumprod()\n",
        "wealth_index_R03.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "54pEmVddA3Li",
        "outputId": "4cb6f838-c309-4a6d-e53b-1a02f6c7a4c1"
      },
      "outputs": [],
      "source": [
        "previous_peaks_R03 = wealth_index_R03.cummax()\n",
        "previous_peaks_R03.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "DmrBzdL9A9BW",
        "outputId": "f0f808bc-4d72-462d-c029-bf28956491b8"
      },
      "outputs": [],
      "source": [
        "drawdown_R03 = (wealth_index_R03 - previous_peaks_R03)/previous_peaks_R03\n",
        "drawdown_R03.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HNxqFiX1BC4J",
        "outputId": "dab8ac26-b8be-472f-fa6b-c08b52f0417f"
      },
      "outputs": [],
      "source": [
        "drawdown_R03.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "NdtTw5YrBE5u",
        "outputId": "16d151d6-8b49-4d0a-84a7-ff9f93ff2966"
      },
      "outputs": [],
      "source": [
        "drawdown_R03[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "FZS6pVF0BL-F",
        "outputId": "3afc53fc-a9cc-4c97-965a-6bee8cac852e"
      },
      "outputs": [],
      "source": [
        "drawdown_R03[\"2015\":].min()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PYsHVHV6Cnd4"
      },
      "source": [
        "#### R06"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "ZPeYRNiQAxnu",
        "outputId": "3d72b659-cb70-4671-ab24-ed45e58f1d1a"
      },
      "outputs": [],
      "source": [
        "wealth_index_R06 = 1000*(1+rets2[\"R06\"]).cumprod()\n",
        "wealth_index_R06.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "4WSXQcZ6Axkd",
        "outputId": "29581d19-691a-40b0-8afa-0b02761e61a5"
      },
      "outputs": [],
      "source": [
        "previous_peaks_R06 = wealth_index_R06.cummax()\n",
        "previous_peaks_R06.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "Db6RPfkoAxhM",
        "outputId": "b87f09d1-b79c-4034-a8ea-98cfcdb316af"
      },
      "outputs": [],
      "source": [
        "drawdown_R06 = (wealth_index_R06 - previous_peaks_R06)/previous_peaks_R06\n",
        "drawdown_R06.plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "lWmU67e0Axdz",
        "outputId": "5a3245ef-abf9-43e7-c1db-d8806fa702a9"
      },
      "outputs": [],
      "source": [
        "drawdown_R06.min()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 336
        },
        "id": "n-8OFiH0AxYk",
        "outputId": "a90bd710-2ce3-4a9d-922c-197d686797d2"
      },
      "outputs": [],
      "source": [
        "drawdown_R06[\"2015\":].plot()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "kzOQLsr8BMkE",
        "outputId": "d591c8e4-d46d-4ed7-f88d-f82c8fa842ec"
      },
      "outputs": [],
      "source": [
        "drawdown_R06[\"2015\":].min()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WwWYY5pYfsJZ"
      },
      "source": [
        "## Time Series Analysis (part 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ohFefAtmIFJ_"
      },
      "outputs": [],
      "source": [
        "df_sales_daily = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesdaily.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "QxGZDDcegmgh",
        "outputId": "2f2935a6-e384-43b5-a589-64568c4b7aa3"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "coHbg383wYNH"
      },
      "source": [
        "### Seasonality Analysis\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "CSh5Y-NBfyoG",
        "outputId": "17305262-1e23-4d54-b0f3-8fb45efc87e2"
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(8, 1, figsize=(15, 50), sharex=True)\n",
        "for name, ax in zip(['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06'], axes):\n",
        "    sns.boxplot(data=df_sales_daily, x='Month', y=name, ax=ax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "GAP3MgVQgAIg",
        "outputId": "0ef35e6b-df86-418c-86c1-007fa79701b1"
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(8, 1, figsize=(15, 50), sharex=True)\n",
        "for name, ax in zip(['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06'], axes):\n",
        "    sns.boxplot(data=df_sales_daily, x='Weekday Name', y=name, ax=ax)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "HIvjVCOXjq7y",
        "outputId": "677aaa8f-e4ee-41ef-c5e3-529ddea13875"
      },
      "outputs": [],
      "source": [
        "cols_plot = ['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06']\n",
        "dfatc_365d = df_sales_daily[cols_plot].rolling(window=365, center=True).mean()\n",
        "dfatc_30d = df_sales_daily[cols_plot].rolling(30, center=True).mean()\n",
        "dfatc_std = df_sales_daily[cols_plot].rolling(30, center=True).std()\n",
        "subplotindex = 0\n",
        "numrows = 4\n",
        "numcols = 2\n",
        "fig, ax = plt.subplots(numrows, numcols, figsize=(30, 50))\n",
        "plt.subplots_adjust(wspace=0.1, hspace=0.3)\n",
        "\n",
        "for x in cols_plot:\n",
        "    rowindex=math.floor(subplotindex/numcols)\n",
        "    colindex=subplotindex-(rowindex*numcols)\n",
        "    ax[rowindex,colindex].plot(df_sales_daily.loc[:,x], linewidth=0.5, label='Daily sales')\n",
        "    ax[rowindex,colindex].plot(dfatc_30d.loc[:,x], label='30-d Rolling Mean')\n",
        "    ax[rowindex,colindex].plot(dfatc_365d.loc[:,x], color='0.2', linewidth=3, label='365-d Rolling Mean')\n",
        "    ax[rowindex,colindex].plot(dfatc_std.loc[:,x], color='0.5', linewidth=3, label='30-d Rolling Std')\n",
        "    ax[rowindex,colindex].set_ylabel('Sales')\n",
        "    ax[rowindex,colindex].legend()\n",
        "    ax[rowindex,colindex].set_title('Trends in ' + x + ' drugs sales');   \n",
        "    subplotindex = subplotindex + 1\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 340
        },
        "id": "VyF-_qQ-j9Gj",
        "outputId": "e1605c7b-1812-49b0-a9a1-f603254713fe"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(30, 10))\n",
        "for nm in cols_plot:\n",
        "    ax.plot(dfatc_365d[nm], label=nm, marker='.', linestyle='-', linewidth=0.5)\n",
        "    ax.legend()\n",
        "    ax.set_ylabel('Drug sales')\n",
        "    ax.set_title('Trends in Drug Sales for different groups (365-d Rolling Means)');"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "rrqa2LKMlTL7",
        "outputId": "b5f6bbf4-efb8-4115-f2c3-42749f7f7623"
      },
      "outputs": [],
      "source": [
        "from statsmodels.tsa.seasonal import seasonal_decompose\n",
        "result = seasonal_decompose(df_sales_daily['M01AB'].rolling(30, center=True).mean().dropna(), freq=365, filt=None)\n",
        "plt.rcParams[\"figure.figsize\"] = (20,20)\n",
        "result.plot()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2jEfP45OlU7Y",
        "outputId": "80eab7f4-9665-43ac-8dbb-9baf5363b02e"
      },
      "outputs": [],
      "source": [
        "for x in ['M01AB','M01AE','N02BA','N02BE','N05B','N05C','R03','R06']:\n",
        "    result = seasonal_decompose(df_sales_weekly[x], freq=52, model = 'additive')\n",
        "    dfs = pd.concat([result.trend, result.seasonal, result.resid, result.observed], axis=1)\n",
        "    dfs.columns = ['trend', 'seasonal', 'residuals', 'observed']\n",
        "    dfs = dfs.dropna()\n",
        "    res = dfs['residuals'].values\n",
        "    obs = dfs['observed'].values\n",
        "    resmean = np.mean(np.abs(res))\n",
        "    obsmean = np.mean(np.abs(obs))\n",
        "    perc = resmean*100/obsmean\n",
        "    print(x + ' RESMEAN:' + str(resmean) + ', OBSMEAN:' + str(obsmean) + ', PERC:' + str(perc) + '%')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rPbTr_-kl-g9"
      },
      "source": [
        "### Stationarity Analysis"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jVouR1t9wrtM",
        "outputId": "15b1ad00-98fd-4e89-90dd-71291c08c21a"
      },
      "outputs": [],
      "source": [
        "from statsmodels.tsa.stattools import adfuller"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Ld9W9Woqwrdt",
        "outputId": "03443d69-52e2-41ff-adcc-062f7587b543"
      },
      "outputs": [],
      "source": [
        "for x in ['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06']:\n",
        "    dftest = adfuller(df_sales_weekly[x], regression='ct', autolag='AIC')\n",
        "    print(\"ADF test for \"+x)\n",
        "    print(\"-----------------------------\")\n",
        "    print(\"Test statistic = {:.3f}\".format(dftest[0]))\n",
        "    print(\"P-value = {:.3f}\".format(dftest[1]))\n",
        "    print(\"Critical values :\")\n",
        "    for k, v in dftest[4].items():\n",
        "        print(\"\\t{}: {} - The data is {} stationary with {}% confidence\".format(k, v, \"not\" if v<dftest[0] else \"\", 100-int(k[:-1])))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "qtZWrMrQwraY"
      },
      "outputs": [],
      "source": [
        "from statsmodels.tsa.stattools import kpss\n",
        "#warnings.filterwarnings(\"ignore\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TkNZTRKfyHIB",
        "outputId": "be167143-ff87-48f5-8255-e81a6226c1bf"
      },
      "outputs": [],
      "source": [
        "for x in ['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06']:\n",
        "    print(\" > Is \" + x + \" data stationary ?\")\n",
        "    dftest = kpss(np.log(df_sales_weekly[x]), 'ct')\n",
        "    print(\"Test statistic = {:.3f}\".format(dftest[1]))\n",
        "    print(\"Critical values: \")\n",
        "    for k, v in dftest[3].items():\n",
        "        print(\"\\t{}: {}\".format(k, v))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "diH1vW3O0B_K"
      },
      "source": [
        "### Regularity Analysis"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LDBQI2zd0Kj_"
      },
      "outputs": [],
      "source": [
        "def ApEn(U, m, r):\n",
        "    def _maxdist(x_i, x_j):\n",
        "        return max([abs(ua - va) for ua, va in zip(x_i, x_j)])\n",
        "    def _phi(m):\n",
        "        x = [[U[j] for j in range(i, i + m - 1 + 1)] for i in range(N - m + 1)]\n",
        "        C = [len([1 for x_j in x if _maxdist(x_i, x_j) <= r]) / (N - m + 1.0) for x_i in x]\n",
        "        return (N - m + 1.0)**(-1) * sum(np.log(C))\n",
        "    N = len(U)\n",
        "    return abs(_phi(m+1) - _phi(m))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Mr0Rndcr1H0i",
        "outputId": "28d5bcfd-096a-4a71-99e2-c1c5116a2333"
      },
      "outputs": [],
      "source": [
        "for x in ['M01AB','M01AE','N02BA','N02BE', 'N05B','N05C','R03','R06']:\n",
        "    print(x + ': ' + str(ApEn(df_sales_weekly[x].values, m=2, r=0.2*np.std(df_sales_weekly[x].values))))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QYg6GIPo21EN"
      },
      "source": [
        "### Autocorrelation analysis"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "T3gwCfCX22JF"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "whzOY1Un22Fi"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Bei_9Pcvl9l1"
      },
      "source": [
        "## Sales Prediction with Machine Learning"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OTlQnImtcJwj"
      },
      "source": [
        "### Predict the pharma sales in January 2020\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "cSO1jFPcl3Dc"
      },
      "outputs": [],
      "source": [
        "def scatterData(X_train, y_train, X_test, y_test, title):\n",
        "    plt.title('Prediction using ' + title)\n",
        "    plt.xlabel('Month sequence', fontsize=20)\n",
        "    plt.ylabel('Sales', fontsize=20)\n",
        "\n",
        "    # Use Matplotlib Scatter Plot\n",
        "    plt.scatter(X_train, y_train, color='blue', label='Training observation points')\n",
        "    plt.scatter(X_test, y_test, color='cyan', label='Testing observation points')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cSUDWCFEc5_y"
      },
      "source": [
        "### Linear Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PzhP8D2gc21z"
      },
      "outputs": [],
      "source": [
        "def predictLinearRegression(X_train, y_train, X_test, y_test):\n",
        "\n",
        "    y_train = y_train.reshape(-1, 1)\n",
        "    y_test = y_test.reshape(-1, 1)\n",
        "\n",
        "    scatterData(X_train, y_train, X_test, y_test, 'Linear Regression')\n",
        "\n",
        "    reg = linear_model.LinearRegression()\n",
        "    reg.fit(X_train, y_train)\n",
        "    plt.plot(X_train, reg.predict(X_train), color='red', label='Linear regressor')\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "\n",
        "    # LINEAR REGRESSION - Predict/Test model\n",
        "    y_predict_linear = reg.predict(X_test)\n",
        "\n",
        "    # LINEAR REGRESSION - Predict for January 2020\n",
        "    linear_predict = reg.predict([[predictFor]])\n",
        "    # linear_predict = reg.predict([[predictFor]])[0]\n",
        "\n",
        "    # LINEAR REGRESSION - Accuracy\n",
        "    accuracy = reg.score(X_train, y_train)\n",
        "\n",
        "    # LINEAR REGRESSION - Error\n",
        "    # error = round(np.mean((y_predict_linear-y_test)**2), 2)\n",
        "    \n",
        "    # Results\n",
        "    print('Linear Regression: ' + str(linear_predict) + ' (Accuracy: ' + str(round(accuracy*100)) + '%)')\n",
        "\n",
        "    return {'regressor':reg, 'values':linear_predict}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-f7mHsjXc_ns"
      },
      "source": [
        "### Polynomial Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "k3QvemHdc2y3"
      },
      "outputs": [],
      "source": [
        "def predictPolynomialRegression(X_train, y_train, X_test, y_test):\n",
        "\n",
        "    y_train = y_train.reshape(-1, 1)\n",
        "    y_test = y_test.reshape(-1, 1)\n",
        "\n",
        "    scatterData(X_train, y_train, X_test, y_test, 'Polynomial Regression')\n",
        "    \n",
        "    poly_reg = PolynomialFeatures(degree = 2)\n",
        "    X_poly = poly_reg.fit_transform(X_train)\n",
        "    poly_reg_model = linear_model.LinearRegression()\n",
        "    poly_reg_model.fit(X_poly, y_train)\n",
        "    plt.plot(X_train, poly_reg_model.predict(poly_reg.fit_transform(X_train)), color='green', label='Polynomial regressor')\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "\n",
        "    # Polynomial Regression - Predict/Test model\n",
        "    y_predict_polynomial = poly_reg_model.predict(X_poly)\n",
        "\n",
        "    # Polynomial Regression - Predict for January 2020\n",
        "    polynomial_predict = poly_reg_model.predict(poly_reg.fit_transform([[predictFor]]))\n",
        "\n",
        "    # Polynomial Regression - Accuracy\n",
        "    # X_poly_test = poly_reg.fit_transform(X_test)\n",
        "    accuracy = poly_reg_model.score(X_poly, y_train)\n",
        "\n",
        "    # Polynomial Regression - Error\n",
        "    # error = round(np.mean((y_predict_polynomial-y_train)**2), 2)\n",
        "\n",
        "    # Result\n",
        "    print('Polynomial Regression: ' + str(polynomial_predict) + ' (Accuracy: ' + str(round(accuracy*100)) + '%)')\n",
        "    return {'regressor':poly_reg_model, 'values':polynomial_predict}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Pu0rbzC5dGRM"
      },
      "source": [
        "### Simple Vector Regression (SVR)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "usMcjs0mc2wN"
      },
      "outputs": [],
      "source": [
        "def predictSVR(X_train, y_train, X_test, y_test):\n",
        "\n",
        "    y_train = y_train.reshape(-1, 1)\n",
        "    y_test = y_test.reshape(-1, 1)\n",
        "\n",
        "    scatterData(X_train, y_train, X_test, y_test, 'Simple Vector Regression (SVR)')\n",
        "\n",
        "    svr_regressor = SVR(kernel='rbf', gamma='auto')\n",
        "    svr_regressor.fit(X_train, y_train.ravel())\n",
        "\n",
        "    # plt.scatter(X_train, y_train, color='red', label='Actual observation points')\n",
        "    plt.plot(X_train, svr_regressor.predict(X_train), label='SVR regressor')\n",
        "    plt.legend()\n",
        "    plt.show()\n",
        "\n",
        "    # Simple Vector Regression (SVR) - Predict/Test model\n",
        "    y_predict_svr = svr_regressor.predict(X_test)\n",
        "\n",
        "    # Simple Vector Regression (SVR) - Predict for January 2020\n",
        "    svr_predict = svr_regressor.predict([[predictFor]])\n",
        "\n",
        "    # Simple Vector Regression (SVR) - Accuracy\n",
        "    accuracy = svr_regressor.score(X_train, y_train)\n",
        "\n",
        "    # Simple Vector Regression (SVR) - Error\n",
        "    # error = round(np.mean((y_predict_svr-y_train)**2), 2)\n",
        "    \n",
        "    # Result\n",
        "    print('Simple Vector Regression (SVR): ' + str(svr_predict) + ' (Accuracy: ' + str(round(accuracy*100)) + '%)')\n",
        "    return {'regressor':svr_regressor, 'values':svr_predict}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bIahGGl0dMbz"
      },
      "source": [
        "### Define a pharma product"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "W_dk6GZcdFxY"
      },
      "outputs": [],
      "source": [
        "product_N02BA = 'N02BA'"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "z2Aw0GKRdOLe"
      },
      "outputs": [],
      "source": [
        "regResults_N02BA = pd.DataFrame(columns=('Linear', 'Polynomial', 'SVR', 'Voting Regressor'), index=[product_N02BA])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CtchTzbTdT8K"
      },
      "outputs": [],
      "source": [
        "rcParams['figure.figsize'] = 12, 8"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1MOvbYb_dapc"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly = pd.read_csv('/content/drive/My Drive/Colab Notebooks/Projects/Pharma Sales/salesmonthly.csv')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "DgvDdbXQdl1n"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly = df_sales_monthly.loc[df_sales_monthly['datum'].str.contains(\"2014\") | df_sales_monthly['datum'].str.contains(\"2015\") | df_sales_monthly['datum'].str.contains(\"2016\") | df_sales_monthly['datum'].str.contains(\"2017\") | df_sales_monthly['datum'].str.contains(\"2018\") | df_sales_monthly['datum'].str.contains(\"2019\")]\n",
        "df_sales_monthly = df_sales_monthly.reset_index()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "ZUc4aLt2do7B",
        "outputId": "78a87215-3c86-4be6-d964-708bfb1769e8"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zPqblJ5zdqZD"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly['datumNumber'] = 1\n",
        "for index, row in df_sales_monthly.iterrows():\n",
        "    df_sales_monthly.loc[index, 'datumNumber'] = index+1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p3pCWQJOduMn"
      },
      "outputs": [],
      "source": [
        "# the first and the last available month is quite low which may indicate that it might be incomplete\n",
        "# and skewing results so we're dropping it\n",
        "df_sales_monthly.drop(df_sales_monthly.head(1).index,inplace=True)\n",
        "df_sales_monthly.drop(df_sales_monthly.tail(1).index,inplace=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bJE-csabdx00"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly = df_sales_monthly[df_sales_monthly[product_N02BA] != 0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "szk0cJKDd2M1",
        "outputId": "9b8307e9-1b92-4310-e9bb-34fd0e4baa42"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "v9HKkhxzd4J-"
      },
      "outputs": [],
      "source": [
        "predictFor = len(df_sales_monthly) + 5"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TOAssK2LeILk",
        "outputId": "729adefc-089c-4faf-efce-36932d1b1150"
      },
      "outputs": [],
      "source": [
        "print('Predictions for the product ' + str(product_N02BA) + ' sales in January 2020')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "0Gova2sYeOVh"
      },
      "outputs": [],
      "source": [
        "regValues_N02BA = {}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6x45yM64eP1p"
      },
      "outputs": [],
      "source": [
        "dfSplit_N02BA = df_sales_monthly[['datumNumber', product_N02BA]]\n",
        "\n",
        "# We are going to keep 30% of the dataset in test dataset\n",
        "train, test = train_test_split(dfSplit_N02BA, test_size=3/10, random_state=0)\n",
        "\n",
        "trainSorted = train.sort_values('datumNumber', ascending=True)\n",
        "testSorted = test.sort_values('datumNumber', ascending=True)\n",
        "\n",
        "X_train = trainSorted[['datumNumber']].values\n",
        "y_train = trainSorted[product_N02BA].values\n",
        "X_test = testSorted[['datumNumber']].values\n",
        "y_test = testSorted[product_N02BA].values"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WplLU24YedoM"
      },
      "source": [
        "### Feature Scaling\n",
        "\n",
        "To improve the performance of the model"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NWouUzaDeY-3"
      },
      "outputs": [],
      "source": [
        "# scale_X = StandardScaler()\n",
        "# scale_y = StandardScaler()\n",
        "\n",
        "# X_train = scale_X.fit_transform(X_train)\n",
        "# y_train = scale_y.fit_transform(y_train.reshape(-1, 1))\n",
        "\n",
        "# X_test = scale_X.fit_transform(X_test)\n",
        "# y_test = scale_y.fit_transform(y_test.reshape(-1, 1))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4MJ1X1aPelE0"
      },
      "source": [
        "### Performance"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0_TS_6VPenJs"
      },
      "source": [
        "### Linear Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 540
        },
        "id": "EMyBTYYFebxi",
        "outputId": "870b1a85-f7d1-4f27-c199-50a07cd9c749"
      },
      "outputs": [],
      "source": [
        "linearResult_N02BA = predictLinearRegression(X_train, y_train, X_test, y_test)\n",
        "reg_N02BA = linearResult_N02BA['regressor']\n",
        "regValues_N02BA['Linear'] = round(linearResult_N02BA['values'][0][0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sOeY7ykUeyk_"
      },
      "source": [
        "#### Polynomial Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 540
        },
        "id": "XqfEQqboebu-",
        "outputId": "53db9cb1-90ff-4954-a02a-b6eab996006c"
      },
      "outputs": [],
      "source": [
        "polynomialResult_N02BA = predictPolynomialRegression(X_train, y_train, X_test, y_test)\n",
        "polynomial_regressor_N02BA = polynomialResult_N02BA['regressor']\n",
        "regValues_N02BA['Polynomial'] = round(polynomialResult_N02BA['values'][0][0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0mDJJhXBfBO8"
      },
      "source": [
        "#### Simple Vector Regression (SVR)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 540
        },
        "id": "XKhnx2Yvebsh",
        "outputId": "f345224e-7e35-4cf5-8227-69789043afe2"
      },
      "outputs": [],
      "source": [
        "svrResult_N02BA = predictSVR(X_train, y_train, X_test, y_test)\n",
        "svr_regressor_N02BA = svrResult_N02BA['regressor']\n",
        "regValues_N02BA['SVR'] = round(svrResult_N02BA['values'][0])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7AKvlxtFfSM3"
      },
      "source": [
        "#### Voting Regressor"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "XCn0PDIVfQKe",
        "outputId": "91d79025-86d5-48bb-8d46-507fdc315f19"
      },
      "outputs": [],
      "source": [
        "vRegressor_N02BA = VotingRegressor(estimators=[('reg', reg_N02BA), ('polynomial_regressor', polynomial_regressor_N02BA), ('svr_regressor', svr_regressor_N02BA)])\n",
        "\n",
        "vRegressorRes = vRegressor_N02BA.fit(X_train, y_train.ravel())\n",
        "\n",
        "# VotingRegressor - Predict for January 2020\n",
        "vRegressor_predict_N02BA = vRegressor_N02BA.predict([[predictFor]])\n",
        "regValues_N02BA['Voting Regressor'] = round(vRegressor_predict_N02BA[0])\n",
        "print('Voting Regressor January 2020 predicted value: ' + str(round(vRegressor_predict_N02BA[0])))\n",
        "regResults_N02BA.loc[product_N02BA] = regValues_N02BA"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 81
        },
        "id": "4suD8fcQfg_t",
        "outputId": "6555f071-1396-4021-f8c6-13ffd2fd3c4f"
      },
      "outputs": [],
      "source": [
        "regResults_N02BA"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vIbjPzfhyygW"
      },
      "source": [
        "## Association Rules (Market Basket Analysis)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "JnyAk2HQfw0p",
        "outputId": "c04466c1-f988-4aa1-db23-d495ac8ec1c1"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products = df_sales_daily.copy()\n",
        "df_sales_daily_products.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "-CaiEcj-0Kbl",
        "outputId": "0b102937-6924-4c11-8fff-0ef36744fea3"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products = df_sales_daily_products[['datum',\t'M01AB',\t'M01AE',\t'N02BA', 'N02BE', 'N05B',\t'N05C',\t'R03',\t'R06']]\n",
        "df_sales_daily_products.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "Ln6leWML3UhV",
        "outputId": "e0161fca-0f4b-486c-a26c-8ee0339ea333"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products.set_index('datum', inplace=True)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "yVZ0yjhi5LsV",
        "outputId": "4bf5fd46-0453-4328-dd37-7acdef97fd44"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['M01AB'] = (df_sales_daily_products['M01AB'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "xZULR4QJ8Omn",
        "outputId": "b9eedba5-92fc-438f-85f4-ad776e414694"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['M01AE'] = (df_sales_daily_products['M01AE'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "YOmwkkYQ6iQm",
        "outputId": "0cf7b45f-37cc-4e8e-d717-6eb734ce0746"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['N02BA'] = (df_sales_daily_products['N02BA'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "F-mUO_Cu6iGy",
        "outputId": "8614f319-6d67-4ae6-86a4-2570afc92f2a"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['N02BE'] = (df_sales_daily_products['N02BE'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "RUWM7eud6h_3",
        "outputId": "134aaab9-c0a8-4150-83e9-aaf424f7d505"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['N05B'] = (df_sales_daily_products['N05B'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "e9YJRrZl6h4f",
        "outputId": "348a11ae-44ee-4f4e-807b-ff04a06efd6e"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['N05C'] = (df_sales_daily_products['N05C'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "MgOMGZiV6hxt",
        "outputId": "868db8ed-33b9-4b28-8e01-9946f2a76f37"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['R03'] = (df_sales_daily_products['R03'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "AHllk1V86hsW",
        "outputId": "f9f36a1e-33d4-4117-e05b-d4a5df4d2fbf"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_products['R06'] = (df_sales_daily_products['R06'] > 0).astype(int)\n",
        "df_sales_daily_products"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "SPbQwDk13eCI",
        "outputId": "546f80b0-17a2-4fed-9099-66165f432890"
      },
      "outputs": [],
      "source": [
        "itemsets = apriori(df_sales_daily_products, min_support=0.2, use_colnames=True)\n",
        "itemsets"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "o2k1bTnY4dnv",
        "outputId": "d464b79e-2f94-4d71-a18a-13b1814d43a8"
      },
      "outputs": [],
      "source": [
        "rules = association_rules(itemsets, metric='confidence', min_threshold=0.5)\n",
        "rules.sort_values(by=['lift'], ascending=False).head(6)\n",
        "print(rules.sort_values(by=['lift'], ascending=False).drop(columns=['antecedent support', 'consequent support', 'conviction']).head(6))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 316
        },
        "id": "ceB_cDvo4dlF",
        "outputId": "a4889a0b-ffe3-4fb9-e4f5-582ce6d6a8ff"
      },
      "outputs": [],
      "source": [
        "rules[[len(c) == 1 for c in rules.consequents]].sort_values(by=['lift'], ascending=False).head(6)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OcHb2CXvBYHZ"
      },
      "source": [
        "If `M01AB, N02BA, R06` is purchased, then with confidence **78.66%** `R03` will also be purchased. This rule has a lift ratio of\n",
        "**1.0214**.\n",
        "\n",
        "If `M01AB, N02BA, R06, N02BE` is purchased, then with confidence **78.66%** `R03` will also be purchased. This rule has a lift ratio of\n",
        "**1.0214**.\n",
        "\n",
        "If `N05B, M01AB, N02BA, R06` is purchased, then with confidence **78.65%** `R03` will also be purchased. This rule has a lift ratio of\n",
        "**1.0212**.\n",
        "\n",
        "If `R06, M01AB, N02BA, N05B, N02BE` is purchased, then with confidence **78.65%** `R03` will also be purchased. This rule has a lift ratio of\n",
        "**1.0212**.\n",
        "\n",
        "If `M01AE, R06, M01AB, N02BA, N05B, N02BE` is purchased, then with confidence **78.62%** `R03` will also be purchased. This rule has a lift ratio of **1.0208**.\n",
        "\n",
        "If `M01AE, R06, M01AB, N02BA, N05B` is purchased, then with confidence **78.62%** `R03` will also be purchased. This rule has a lift ratio of\n",
        "**1.0208**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qwthK2M3A_5W"
      },
      "source": [
        "The support for the rule indicates its impact in terms of overall size: How many transactions are affected? If\n",
        "only a small number of transactions are affected, the rule may be of little use (unless the consequent is very\n",
        "valuable and/or the rule is very efficient in finding it).\n",
        "\n",
        "• The lift ratio indicates how efficient the rule is in finding consequents, compared to random selection. A very\n",
        "efficient rule is preferred to an inefficient rule, but we must still consider support: A very efficient rule that\n",
        "has very low support may not be as desirable as a less efficient rule with much greater support.\n",
        "\n",
        "• The confidence tells us at what rate consequents will be found and is useful in determining the business or\n",
        "operational usefulness of a rule: A rule with low confidence may find consequents at too low a rate to be\n",
        "worth the cost of (say) promoting the consequent in all the transactions that involve the antecedent."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u6sHrg6hEjMK"
      },
      "source": [
        "## Time Series Forecasting (part 2)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V7neAYogSho-"
      },
      "source": [
        "#### M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "g9QyGypUD8g5",
        "outputId": "c1a243a1-ad63-4bc0-a756-2f94763e86e3"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum = df_sales_daily[['datum', 'M01AB']]\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 490
        },
        "id": "dZFyU0gYGPl2",
        "outputId": "6c9b7464-54c3-45fc-c813-0536971ff20e"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum['Date'] = pd.to_datetime(df_sales_daily_M01AB_datum.datum, format='%m/%d/%Y')\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "px73mXmcGPZq"
      },
      "outputs": [],
      "source": [
        "M01AB_ts = pd.Series(df_sales_daily_M01AB_datum.M01AB.values, index=df_sales_daily_M01AB_datum.Date, name='M01AB_Sales')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 300
        },
        "id": "hQeWiIngHBcB",
        "outputId": "33a6b7fc-3890-4c99-8b33-329b47a7aafa"
      },
      "outputs": [],
      "source": [
        "# Pandas Version\n",
        "M01AB_ts.plot(ylim=[0, 20], legend=False)\n",
        "plt.xlabel('Year'); plt.ylabel('M01AB')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CBpGA4bRHBZ3"
      },
      "outputs": [],
      "source": [
        "# Create short time series from 2018 to 2018 using a slice \n",
        "M01AB_ts_1yr = M01AB_ts['2018':'2018']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "yFfJfjEvHBXP",
        "outputId": "34756467-a700-4779-cd80-1a350334f408"
      },
      "outputs": [],
      "source": [
        "# create a data frame with additional predictors from time series\n",
        "# the following command adds a constant term, a trend term and a quadratic trend term\n",
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ctt')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 394
        },
        "id": "oL_unXZaIqrh",
        "outputId": "adadd26f-9600-449d-b397-3f7304b4ce22"
      },
      "outputs": [],
      "source": [
        "M01AB_df.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "A10biB3MIqk6"
      },
      "outputs": [],
      "source": [
        "# fit a linear trend model to the time series\n",
        "M01AB_lm = sm.ols(formula='M01AB_Sales ~ trend + trend_squared', data=M01AB_df).fit()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3h2Ru5s4KEiF"
      },
      "source": [
        "##### Plots that enhance the different components of the time series\n",
        "\n",
        "(a) Zoom-in to 1 year of data \n",
        "\n",
        "(b) Original series with overlaid quadratic trendline"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 287
        },
        "id": "LsQ3Rm9KIqiA",
        "outputId": "de20d7ec-0fe9-4315-c591-125b935132c4"
      },
      "outputs": [],
      "source": [
        "# shorter and longer time series\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(35, 15))\n",
        "M01AB_ts_1yr.plot(ax=axes[0])\n",
        "M01AB_ts.plot(ax=axes[1])\n",
        "for ax in axes:\n",
        "    ax.set_xlabel('Time')\n",
        "    ax.set_ylabel('M01AB')\n",
        "    ax.set_ylim(0, 20)\n",
        "M01AB_lm.predict(M01AB_df).plot(ax=axes[1])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 168
        },
        "id": "00nofGsoKLet",
        "outputId": "ef734465-3623-4972-9776-02014058028e"
      },
      "outputs": [],
      "source": [
        "# plot the time series\n",
        "ax = M01AB_ts.plot(figsize= (35, 8))\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB (in 000s)')\n",
        "ax.set_ylim(0, 20)\n",
        "M01AB_lm.predict(M01AB_df).plot(ax=ax)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8BlPtJPZKLb8"
      },
      "outputs": [],
      "source": [
        "nValid = 36\n",
        "nTrain = len(M01AB_ts) - nValid\n",
        "# partition the data\n",
        "train_ts_M01AB = M01AB_ts[:nTrain]\n",
        "valid_ts_M01AB = M01AB_ts[nTrain:]\n",
        "# generate the naive and seasonal naive forecast\n",
        "naive_pred_M01AB = pd.Series(train_ts_M01AB[-1], index=valid_ts_M01AB.index)\n",
        "last_season_M01AB = train_ts_M01AB[-12:]\n",
        "seasonal_pred_M01AB = pd.Series(pd.concat([last_season_M01AB]*5)[:len(valid_ts_M01AB)].values, index=valid_ts_M01AB.index)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0IQihTHBKcQ7",
        "outputId": "fc366ac4-5b8a-4cd6-a625-15c41ab2153f"
      },
      "outputs": [],
      "source": [
        "train_ts_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5PV5mCCvKLZH",
        "outputId": "59bba95d-4e51-47c9-ae33-ce7de91dce5d"
      },
      "outputs": [],
      "source": [
        "valid_ts_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WbnrBky7KzdE",
        "outputId": "b03fa90b-19b9-4128-ca5f-db9a0dfeb195"
      },
      "outputs": [],
      "source": [
        "naive_pred_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "bOaqa1VNKzZQ",
        "outputId": "17b0b18d-c1e6-4d38-cc31-09db356417d8"
      },
      "outputs": [],
      "source": [
        "last_season_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nhijxym9KzVl",
        "outputId": "2ba77c5d-0f82-4fed-b52a-2788f5b4eec9"
      },
      "outputs": [],
      "source": [
        "seasonal_pred_M01AB"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EvW_R0VzNxrx"
      },
      "source": [
        "##### Naive and Seasonal Naive Forecasts in a 3-year validation set for M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 305
        },
        "id": "K1qG6n_2RKQ2",
        "outputId": "d2cc9f04-3b87-49be-d446-45f947bccf73"
      },
      "outputs": [],
      "source": [
        "# plot forecasts and actual in the training and validation sets\n",
        "ax = train_ts_M01AB.plot(color='C0', linewidth=0.75, figsize= (35,15))\n",
        "valid_ts_M01AB.plot(ax=ax, color='C0', linestyle='dashed', linewidth=0.75)\n",
        "#ax.set_xlim('2015', '2016')\n",
        "ax.set_ylim(0, 20)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('Ridership (in 000s)')\n",
        "naive_pred_M01AB.plot(ax=ax, color='green')\n",
        "seasonal_pred_M01AB.plot(ax=ax, color='orange')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 162
        },
        "id": "6snZI9vCL8zz",
        "outputId": "b901f0f7-6f44-4017-ebd8-34ff125d276a"
      },
      "outputs": [],
      "source": [
        "# plot forecasts and actual in the training and validation sets\n",
        "ax = train_ts_M01AB.plot(color='C0', linewidth=0.75, figsize= (35,7))\n",
        "valid_ts_M01AB.plot(ax=ax, color='C0', linestyle='dashed', linewidth=0.75)\n",
        "#ax.set_xlim('2015', '2016')\n",
        "ax.set_ylim(0, 20)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "naive_pred_M01AB.plot(ax=ax, color='green')\n",
        "seasonal_pred_M01AB.plot(ax=ax, color='orange')\n",
        "\n",
        "# determine coordinates for drawing the arrows and lines\n",
        "one_month = pd.Timedelta('31 days')\n",
        "xtrain = (min(train_ts_M01AB.index), max(train_ts_M01AB.index) - one_month)\n",
        "xvalid = (min(valid_ts_M01AB.index) + one_month, max(valid_ts_M01AB.index) - one_month)\n",
        "xfuture = (max(valid_ts_M01AB.index) + one_month, '2020')\n",
        "xtv = xtrain[1] + 0.5 * (xvalid[0] - xtrain[1])\n",
        "xvf = xvalid[1] + 0.5 * (xfuture[0] - xvalid[1])\n",
        "ax.add_line(plt.Line2D(xtrain, (2450, 2450), color='black', linewidth=0.5))\n",
        "ax.add_line(plt.Line2D(xvalid, (2450, 2450), color='black', linewidth=0.5))\n",
        "ax.add_line(plt.Line2D(xfuture, (2450, 2450), color='black', linewidth=0.5))\n",
        "ax.text('2014', 22, 'Training')\n",
        "ax.text('2019-09-03', 22, 'Validation')\n",
        "ax.text('2019-10-09', 22, 'Future')\n",
        "ax.axvline(x=xtv, ymin=0, ymax=1, color='black', linewidth=0.5)\n",
        "ax.axvline(x=xvf, ymin=0, ymax=1, color='black', linewidth=0.5)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nv_YhygfRvrg"
      },
      "source": [
        "##### Predictive accuracy of naive and seasonal naive forecasts in the validation and training set for M01AB"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zLX6x0arR1DK"
      },
      "source": [
        "Table below compares the accuracies of these two naive forecasts.\n",
        "Because M01AB has monthly seasonality, the seasonal\n",
        "naive forecast is the clear winner on both training and validation and\n",
        "on all popular measures. In choosing between the two models, the\n",
        "accuracy on the validation set is more relevant than the accuracy on\n",
        "the training set. Performance on the validation set is more indicative\n",
        "of how the models will perform in the future."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d-iXsDoRR7P9"
      },
      "source": [
        "###### Validation Set"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "HechN5jtL8ry",
        "outputId": "e76654f3-c012-4dbc-b76e-cbdfa2c0f263"
      },
      "outputs": [],
      "source": [
        "regressionSummary(valid_ts_M01AB, naive_pred_M01AB)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "f8UHdkV4K2LO",
        "outputId": "a4167606-a74c-4b3b-f159-45a4eb87d0df"
      },
      "outputs": [],
      "source": [
        " regressionSummary(valid_ts_M01AB, seasonal_pred_M01AB)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "u_fmI-KNSTVe"
      },
      "source": [
        "###### Training set"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5ncelHUMK2IN",
        "outputId": "f82a5b6a-7e43-4ac5-8ba7-f909d2f0cd5b"
      },
      "outputs": [],
      "source": [
        "# Calculate naive metrics for training set (shifted by 1 month)\n",
        "regressionSummary(train_ts_M01AB[1:], train_ts_M01AB[:-1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "JEElY10CSAAz",
        "outputId": "5477e347-1b47-4093-888d-f65185f3cc56"
      },
      "outputs": [],
      "source": [
        " # Calculate seasonal naive metrics for training set (shifted by 12 months)\n",
        "regressionSummary(train_ts_M01AB[12:], train_ts_M01AB[:-12])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lL63Gy7OEmLK"
      },
      "source": [
        "### Regression Based Forecasting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 510
        },
        "id": "rzApn618En8-",
        "outputId": "f5dccac6-4df5-4e13-978e-5d14c2aeb438"
      },
      "outputs": [],
      "source": [
        "# fit a linear trend model to the time series\n",
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ct')\n",
        "M01AB_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cfvWV6VrWR0A"
      },
      "source": [
        "#### Linear Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tNSLmQ1DX8wl"
      },
      "outputs": [],
      "source": [
        "nValid = 36\n",
        "nTrain = len(M01AB_df) - nValid\n",
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "BJkEwd8_X_0A",
        "outputId": "a5086e12-b6d3-4cdd-c523-273f235cc4f0"
      },
      "outputs": [],
      "source": [
        "train_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "8Tit_8DlYBAu",
        "outputId": "acd7b3ea-1180-4300-b314-0350e547b90c"
      },
      "outputs": [],
      "source": [
        "valid_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Dj8nMWWjVSOw"
      },
      "outputs": [],
      "source": [
        "M01AB_lm = sm.ols(formula='M01AB_Sales ~ trend', data=M01AB_df).fit()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eMz8vuz7WQaW"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_linear = sm.ols(formula='M01AB_Sales ~ trend', data=train_df).fit()\n",
        "predict_df_linear = M01AB_lm_linear.predict(valid_df)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 287
        },
        "id": "NmWGMnwaVSK5",
        "outputId": "742e5e1c-943c-4070-9b67-200a467f7747"
      },
      "outputs": [],
      "source": [
        "# plot the time series\n",
        "ax = M01AB_ts.plot(figsize= (35,15))\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.set_ylim(0, 20)\n",
        "M01AB_lm.predict(M01AB_df).plot(ax=ax)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zkOXbH-ZVbpA"
      },
      "source": [
        "##### A linear trend fit to M01AB in the training period and forecasted in the validation period\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "c6JfzVfrVU28"
      },
      "outputs": [],
      "source": [
        "# fit linear model using training set and predict on validation set\n",
        "M01AB_lm = sm.ols(formula='M01AB_Sales ~ trend', data=train_df).fit()\n",
        "predict_df = M01AB_lm.predict(valid_df)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9S5S_YZZeUiN"
      },
      "source": [
        "### Single Graph Layout"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "1Yx5wcoMVenx"
      },
      "outputs": [],
      "source": [
        "def singleGraphLayout(ax, ylim, train_df, valid_df):\n",
        "    ax.set_xlim('2014', '2019-09')\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",
        "    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",
        "    ypos = 0.925 * ylim[1] + 0.075 * ylim[0]\n",
        "    ax.text('2014', ypos, 'Training')\n",
        "    ax.text('2019-09-3', ypos, 'Validation')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IDJU3GfCeXQC"
      },
      "source": [
        "### Graph Layout"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Llg9SJ3fVek6"
      },
      "outputs": [],
      "source": [
        "def graphLayout(axes, train_df, valid_df):\n",
        "    singleGraphLayout(axes[0], [0, 20], train_df, valid_df)\n",
        "    singleGraphLayout(axes[1], [-7, 13], train_df, valid_df)\n",
        "    train_df.plot(y='M01AB_Sales', ax=axes[0], color='C0', linewidth=0.75)\n",
        "    valid_df.plot(y='M01AB_Sales', ax=axes[0], color='C0', linestyle='dashed',\n",
        "    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('M01AB_Sales')\n",
        "    axes[1].set_ylabel('Forecast Errors')\n",
        "    if axes[0].get_legend():\n",
        "        axes[0].get_legend().remove()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "tt9Qv4AkZuVT",
        "outputId": "6ce6f992-4b18-4f99-ec9b-55c07d4a924c"
      },
      "outputs": [],
      "source": [
        "train_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T8i4ohsyVnia"
      },
      "source": [
        "##### A linear trend fit to M01AB in the training period (a) and forecasted in the validation period (b)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 375
        },
        "id": "X5BRx8VfVehs",
        "outputId": "c15463b2-138d-4e84-9d1b-0b298814088f"
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(35, 20))\n",
        "M01AB_lm.predict(train_df).plot(ax=axes[0], color='C1')\n",
        "M01AB_lm.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
        "residual = train_df.M01AB_Sales - M01AB_lm.predict(train_df)\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_df.M01AB_Sales - M01AB_lm.predict(valid_df)\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_df, valid_df)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eBWCyqHAejij"
      },
      "source": [
        "##### Summary"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 478
        },
        "id": "gpLWpO25Veet",
        "outputId": "693155ac-af48-4e17-f60b-dde5afba6bdd"
      },
      "outputs": [],
      "source": [
        "M01AB_lm.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-dVOxTzWWXKI"
      },
      "source": [
        "#### Exponential Regression"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AAZZrFfhVebn",
        "outputId": "7ad92e80-a2ff-42d5-a8c2-78dc00dfbd02"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_expo = sm.ols(formula='np.log(M01AB_Sales) ~ trend', data=train_df).fit()\n",
        "predict_df_expo = M01AB_lm_expo.predict(valid_df)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_9DGFIHHWbYK"
      },
      "source": [
        "##### Exponential (green) and linear (orange) trend used to forecast M01AB\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 192
        },
        "id": "T3VGRMOrV6ah",
        "outputId": "07beb44e-65b3-44c1-e3bb-117e5d28af8a"
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(55,15))\n",
        "train_df.plot(y='M01AB_Sales', ax=ax, color='C0', linewidth=0.75)\n",
        "valid_df.plot(y='M01AB_Sales', ax=ax, color='C0', linestyle='dashed', linewidth=0.75)\n",
        "singleGraphLayout(ax, [0, 20], train_df, valid_df)\n",
        "M01AB_lm_linear.predict(train_df).plot(color='C1')\n",
        "M01AB_lm_linear.predict(valid_df).plot(color='C1', linestyle='dashed')\n",
        "M01AB_lm_expo.predict(train_df).apply(lambda row: math.exp(row)).plot(color='C2')\n",
        "M01AB_lm_expo.predict(valid_df).apply(lambda row: math.exp(row)).plot(color='C2',  linestyle='dashed')\n",
        "ax.get_legend().remove()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 478
        },
        "id": "dfEf2iY4egZF",
        "outputId": "cff7a050-3bc0-4e42-a04b-82230162671a"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_expo.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "evmjsu1YbN5v",
        "outputId": "22c31282-2064-4526-d3a9-de1fad2dbf3b"
      },
      "outputs": [],
      "source": [
        "M01AB_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 510
        },
        "id": "op3a6lmLW3an",
        "outputId": "724968cc-873b-44d6-9a73-6bc6fed48526"
      },
      "outputs": [],
      "source": [
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ct')\n",
        "M01AB_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "aVCsWYcsW6M8",
        "outputId": "4091341d-c921-41c1-c694-260e1d8b938f"
      },
      "outputs": [],
      "source": [
        "M01AB_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "wdUQWeMwW7La"
      },
      "outputs": [],
      "source": [
        "nValid = 36\n",
        "nTrain = len(M01AB_df) - nValid\n",
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "ccRcWKDVW7Il",
        "outputId": "e02d594a-eb1a-467c-d5b3-7f2558001dd6"
      },
      "outputs": [],
      "source": [
        "train_df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mxcQ7gEddhmX"
      },
      "source": [
        "## Polynomial Regression"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-vglDJ6IXDDH"
      },
      "source": [
        "##### Quadratic trend model used to forecast M01AB. Plots of fitted, forecasted, and actual values (a) and forecast errors (b)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 242
        },
        "id": "cxnwgQVWW7Fw",
        "outputId": "cb05ab66-cd47-4b67-851f-6e8ecd5ddb63"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_poly = sm.ols(formula='M01AB_Sales ~ trend + np.square(trend)', data=train_df).fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(55, 20))\n",
        "M01AB_lm_poly.predict(train_df).plot(ax=axes[0], color='C1')\n",
        "M01AB_lm_poly.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
        "residual = train_df.M01AB_Sales - M01AB_lm_poly.predict(train_df)\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_df.M01AB_Sales - M01AB_lm_poly.predict(valid_df)\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_df, valid_df)\n",
        "plt.show()\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jXHLrXJjeB34"
      },
      "source": [
        "##### Summary of Polynomial Regression with Quadratic Trend"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 499
        },
        "id": "NRqzFxQRd-QT",
        "outputId": "0b6475a7-5fdc-46e8-ff5e-00fc5e0f7c8e"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_poly.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v6I5jJaDXHik"
      },
      "source": [
        "##### Summary of output from fitting additive seasonality to the M01AB data in the training period"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 707
        },
        "id": "hh_sSB1WW7Cs",
        "outputId": "2d563e75-fa42-46fb-f2a3-e8d9dc7ec8d5"
      },
      "outputs": [],
      "source": [
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='c')\n",
        "M01AB_df['Month'] = M01AB_df.index.month\n",
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]\n",
        "M01AB_lm_season = sm.ols(formula='M01AB_Sales ~ C(Month)', data=train_df).fit()\n",
        "M01AB_lm_season.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "F8qgz3krXK5v"
      },
      "source": [
        "##### Regression model with seasonality applied to the M01AB (a) and its forecast errors (b)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "IX8eSsG-XLow",
        "outputId": "fa9afa24-e171-44df-8185-dfc6346f3f9d"
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(55, 20))\n",
        "M01AB_lm_season.predict(train_df).plot(ax=axes[0], color='C1')\n",
        "M01AB_lm_season.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
        "residual = train_df.M01AB_Sales - M01AB_lm_season.predict(train_df)\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_df.M01AB_Sales - M01AB_lm_season.predict(valid_df)\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_df, valid_df)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "djPcGMSWXOEn",
        "outputId": "fcc101de-1756-4619-e4b6-bea7135c82b5"
      },
      "outputs": [],
      "source": [
        "train_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "qWo82ItoXPin",
        "outputId": "a3c52a3d-167b-4d42-b326-fa10e8423f21"
      },
      "outputs": [],
      "source": [
        "M01AB_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "vm5Ve4wpXQeb",
        "outputId": "1018cc90-5a04-4cf3-fd85-ec3605bb5b94"
      },
      "outputs": [],
      "source": [
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ct')\n",
        "M01AB_df['Month'] = M01AB_df.index.month\n",
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mvhJFGrfXRwL"
      },
      "outputs": [],
      "source": [
        "formula = 'M01AB_Sales ~ trend + np.square(trend) + C(Month)'\n",
        "M01AB_lm_trendseason = sm.ols(formula=formula, data=train_df).fit()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-OLCZUOOXUKf"
      },
      "source": [
        "##### Regression model with trend and seasonality applied to M01AB (a) and its forecast errors (b)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 242
        },
        "id": "g26_MtEBXTeG",
        "outputId": "9ac9f193-7fe0-4fa9-f9e2-4a94a7b91ef3"
      },
      "outputs": [],
      "source": [
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(55, 20))\n",
        "M01AB_lm_trendseason.predict(train_df).plot(ax=axes[0], color='C1')\n",
        "M01AB_lm_trendseason.predict(valid_df).plot(ax=axes[0], color='C1', linestyle='dashed')\n",
        "residual = train_df.M01AB_Sales - M01AB_lm_trendseason.predict(train_df)\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_df.M01AB_Sales - M01AB_lm_trendseason.predict(valid_df)\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_df, valid_df)\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vtMIW1jvXYgw"
      },
      "source": [
        "##### Summary of output from fitting trend and seasonality to M01AB in the training period"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 729
        },
        "id": "yV2tWaw6Xaii",
        "outputId": "5cfdd04d-3f6d-4f66-d3c7-fb57db552f0f"
      },
      "outputs": [],
      "source": [
        "M01AB_lm_trendseason.summary()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "N3oAEAThXbX2"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ToLixRvbXbVM"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EP6TaeUOc0e5"
      },
      "source": [
        "## Smoothing Methods\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SC1dN8bgc3dS"
      },
      "source": [
        "### Centered Moving Average for Visualization\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "oAnLYbwQfn5Q",
        "outputId": "cdb69df3-5538-4bee-f86c-5260a90ed882"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 468
        },
        "id": "1Min4XNegIzo",
        "outputId": "f1f07951-aab2-46e8-8a60-cb4fb9773cd9"
      },
      "outputs": [],
      "source": [
        "df_sales_daily"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "OaYKG_F6fjPz",
        "outputId": "69da852d-2d1b-40bf-ff8a-bf69d3455b0f"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum = df_sales_daily[['datum', 'M01AB']]\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 490
        },
        "id": "UV-qzMkPfqMI",
        "outputId": "08a1ee39-3e30-4513-edc4-c9a2e82a2c7c"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum['Date'] = pd.to_datetime(df_sales_daily_M01AB_datum.datum, format='%m/%d/%Y')\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "4dZ514P1mUsA",
        "outputId": "7d239765-b1ee-4d1e-c504-433db21ea8f5"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum.tail(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "DeJXL88rmceM",
        "outputId": "8c07caad-62d8-4f15-992c-d9f3813f471a"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum.dtypes"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "JVAbKEwJlOeT",
        "outputId": "70547812-046f-4abe-8434-71f921dbb578"
      },
      "outputs": [],
      "source": [
        "df_sales_endmonth_M01AB = df_sales_daily_M01AB_datum.groupby([df_sales_daily_M01AB_datum.Date.dt.year, df_sales_daily_M01AB_datum.Date.dt.month]).last()\n",
        "df_sales_endmonth_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "HtMFHvs2nX59",
        "outputId": "a095ad2c-6f32-48a5-a069-845281f2b891"
      },
      "outputs": [],
      "source": [
        "df_sales_endmonth_M01AB = df_sales_endmonth_M01AB[['datum', 'M01AB', 'Date']]\n",
        "df_sales_endmonth_M01AB"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yhqVmQd5csyE"
      },
      "outputs": [],
      "source": [
        "# Load data and convert to time series\n",
        "M01AB_ts = pd.Series(df_sales_endmonth_M01AB.M01AB.values, index=df_sales_endmonth_M01AB.Date, name='M01AB_Sales')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Texcrjf8c6IO"
      },
      "outputs": [],
      "source": [
        "M01AB_ts.index = pd.DatetimeIndex(M01AB_ts.index, freq=M01AB_ts.index.inferred_freq)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3ocKTTjVc6E7"
      },
      "outputs": [],
      "source": [
        "# Centered moving average with window size = 12\n",
        "ma_centered = M01AB_ts.rolling(12, center=True).mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pjXOL9Omc6Cc"
      },
      "outputs": [],
      "source": [
        "# Trailing moving average with window size = 12\n",
        "ma_trailing = M01AB_ts.rolling(12).mean()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tZ74MS3CdMqJ"
      },
      "source": [
        "### Centered moving average (smooth blue line) and trailing moving average (broken orange line) with window w = 12, overlaid on M01AB series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 251
        },
        "id": "QSo0MMl4dJw2",
        "outputId": "dc806e8c-a7a6-4dc8-8868-5e4dccb60912"
      },
      "outputs": [],
      "source": [
        "# Shift the average by one time unit to get the next day predictions\n",
        "ma_centered = pd.Series(ma_centered[:-1].values, index=ma_centered.index[1:])\n",
        "ma_trailing = pd.Series(ma_trailing[:-1].values, index=ma_trailing.index[1:])\n",
        "fig, ax = plt.subplots(figsize=(45, 15))\n",
        "ax = M01AB_ts.plot(ax=ax, color='black', linewidth=0.25)\n",
        "ma_centered.plot(ax=ax, linewidth=2)\n",
        "ma_trailing.plot(ax=ax, style='--', linewidth=2)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.legend(['M01AB', 'Centered Moving Average', 'Trailing Moving Average'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ZSqF0L7rdf9q"
      },
      "source": [
        "### Trailing Moving Average for Forecasting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jPVc1QaIqxlv",
        "outputId": "81883c4f-a376-4b92-ef17-bb9a74e63b55"
      },
      "outputs": [],
      "source": [
        "len(M01AB_ts)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "AA5Mntk9q071",
        "outputId": "7cb94c0a-e895-4240-aa3f-ab086ce69213"
      },
      "outputs": [],
      "source": [
        "len(M01AB_ts)*0.8"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "9QKmExtqdWID"
      },
      "outputs": [],
      "source": [
        "# partition the data\n",
        "nValid = 20\n",
        "nTrain = len(M01AB_ts) - nValid\n",
        "train_ts = M01AB_ts[:nTrain]\n",
        "valid_ts = M01AB_ts[nTrain:]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pr8XSIUudjgF"
      },
      "outputs": [],
      "source": [
        "# moving average on training\n",
        "ma_trailing = train_ts.rolling(12).mean()\n",
        "last_ma = ma_trailing[-1]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-sS25Z95dz-K"
      },
      "source": [
        "### Trailing moving average forecaster with w = 12 applied to M01AB series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 503
        },
        "id": "3RIZrir_dlPE",
        "outputId": "fefd4894-c536-4397-e05b-f874d1472971"
      },
      "outputs": [],
      "source": [
        "# create forecast based on last moving average in the training period\n",
        "ma_trailing_pred = pd.Series(last_ma, index=valid_ts.index)\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(45, 30))\n",
        "ma_trailing.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "ma_trailing_pred.plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - ma_trailing\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - ma_trailing_pred\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fgOVW8k3ev94"
      },
      "source": [
        "### Applying Moving Average to the residuals from the regression model (which lack trend and seasonality), to forecast the 2018-19 residual"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "zo3PKdhseH10",
        "outputId": "e68bd686-72ff-4130-92f8-3601f708231d"
      },
      "outputs": [],
      "source": [
        "# Build a model with seasonality, trend, and quadratic trend\n",
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ct')\n",
        "M01AB_df['Month'] = M01AB_df.index.month"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 455
        },
        "id": "lUFttDsme9FF",
        "outputId": "613a11c9-0486-4170-ef5a-cccc4f26b4ee"
      },
      "outputs": [],
      "source": [
        "M01AB_df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6EpEDtChe_FP"
      },
      "outputs": [],
      "source": [
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]\n",
        "formula = 'M01AB_Sales ~ trend + np.square(trend) + C(Month)'\n",
        "M01AB_lm_trendseason = sm.ols(formula=formula, data=train_df).fit()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gCgTRhSziCwn"
      },
      "outputs": [],
      "source": [
        "# create single-point forecast\n",
        "M01AB_prediction = M01AB_lm_trendseason.predict(valid_df.iloc[0, :])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "_vS6xzqPiGGu",
        "outputId": "034eee8e-37e5-4070-96d6-beb248622651"
      },
      "outputs": [],
      "source": [
        "# apply Moving Average to residuals\n",
        "ma_trailing = M01AB_lm_trendseason.resid.rolling(12).mean()\n",
        "print('Prediction', M01AB_prediction[0])\n",
        "print('Moving_Average_trailing', ma_trailing[-1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0WTD7jKpilox"
      },
      "source": [
        "### Simple Exponential Smoothing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Mr8ZLx8SipOw"
      },
      "source": [
        "### Output for simple exponential smoothing forecaster with = 0.2, applied to the series of residuals from the regression model (which lack trend and seasonality). The forecast value is -0.577178."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 267
        },
        "id": "MfR2N3bPiLds",
        "outputId": "8e78af9a-c51a-46e0-f3f3-0821a3cd3c2e"
      },
      "outputs": [],
      "source": [
        "residuals_ts = M01AB_lm_trendseason.resid\n",
        "residuals_pred = valid_df.M01AB_Sales - M01AB_lm_trendseason.predict(valid_df)\n",
        "fig, ax = plt.subplots(figsize=(45,15))\n",
        "M01AB_lm_trendseason.resid.plot(ax=ax, color='black', linewidth=0.5)\n",
        "residuals_pred.plot(ax=ax, color='black', linewidth=0.5)\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.set_xlabel('Time')\n",
        "ax.axhline(y=0, xmin=0, xmax=1, color='grey', linewidth=0.5)\n",
        "\n",
        "# Run Exponential Smoothing\n",
        "# with smoothing level alpha = 0.2\n",
        "expSmooth = ExponentialSmoothing(residuals_ts, freq='M')\n",
        "expSmoothFit = expSmooth.fit(smoothing_level=0.2)\n",
        "expSmoothFit.fittedvalues.plot(ax=ax)\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=ax, style='--', linewidth=2, color='C0')\n",
        "singleGraphLayout(ax, [-10, 10], train_df, valid_df)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "M7-4KxQkqUOJ"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.fittedvalues"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Yh58SDmppXIn"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.forecast(len(valid_ts))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4l30e94rorJJ",
        "outputId": "18227717-b804-4f7a-8c90-b0acb086c7a4"
      },
      "outputs": [],
      "source": [
        "residuals_ts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "I2uayzXOsHXo",
        "outputId": "51a6a91c-891f-4115-c7a8-8022b67bea7d"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.forecast(len(valid_ts))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gUE4d8enoaLs"
      },
      "source": [
        "### Advanced Exponential Smoothing"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 608
        },
        "id": "pKYEJdzgiziK",
        "outputId": "1fe26358-2762-4630-d692-c804714546b8"
      },
      "outputs": [],
      "source": [
        "# Run exponential smoothing with additive trend and additive seasonal\n",
        "expSmooth = ExponentialSmoothing(train_ts, trend='additive', seasonal='additive', seasonal_periods=12, freq='M')\n",
        "expSmoothFit = expSmooth.fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))\n",
        "expSmoothFit.fittedvalues.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - expSmoothFit.fittedvalues\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - expSmoothFit.forecast(len(valid_ts))\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "AE0aHBhDpuJ6"
      },
      "source": [
        "### Summary of a Holt-Winter’s exponential smoothing model applied to the M01AB data. Included are the initial and final states\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "iSkbBmveocqZ",
        "outputId": "f5d11715-3888-407d-bdab-89174321825e"
      },
      "outputs": [],
      "source": [
        "print(expSmoothFit.params)\n",
        "print('AIC: ', expSmoothFit.aic)\n",
        "print('AICc: ', expSmoothFit.aicc)\n",
        "print('BIC: ', expSmoothFit.bic)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pEXJmOb9HMem"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Lp5ouxxVyVKB"
      },
      "source": [
        "## Smoothing on all data points"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "YLWF10Jd4gWd",
        "outputId": "f7898c84-42f5-4820-a3b6-d702966faade"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "6EPgYCGm3XRW",
        "outputId": "e287fcf5-9a54-4c2e-a82f-5d793cf02281"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.tail()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oIithMUo3Dvj"
      },
      "source": [
        "### Single Graph Layout"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MYIYJN_i3Efv"
      },
      "outputs": [],
      "source": [
        "def singleGraphLayout(ax, ylim, train_df, valid_df):\n",
        "    ax.set_xlim('2014', '2019-12')\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",
        "    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",
        "    ypos = 0.925 * ylim[1] + 0.075 * ylim[0]\n",
        "    ax.text('2015', ypos, 'Training')\n",
        "    ax.text('2019-4', ypos, 'Validation')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hquA2vLa3Mh1"
      },
      "source": [
        "### Graph Layout"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OeCca-5E3OB7"
      },
      "outputs": [],
      "source": [
        "def graphLayout(axes, train_df, valid_df):\n",
        "    singleGraphLayout(axes[0], [-2, 20], train_df, valid_df)\n",
        "    singleGraphLayout(axes[1], [-7, 13], train_df, valid_df)\n",
        "    train_df.plot(y='M01AB_Sales', ax=axes[0], color='C0', linewidth=0.75)\n",
        "    valid_df.plot(y='M01AB_Sales', ax=axes[0], color='C0', linestyle='dashed',\n",
        "    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('M01AB_Sales')\n",
        "    axes[1].set_ylabel('Forecast Errors')\n",
        "    if axes[0].get_legend():\n",
        "        axes[0].get_legend().remove()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "y7BeTJIIye5g"
      },
      "source": [
        "### Centered Moving Average for Visualization\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "_tXDFDYxylGD",
        "outputId": "70816eda-5b18-4e9f-8df2-28fe8d59ae70"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "OKexdefVyYrH",
        "outputId": "7ee441a9-4280-4faa-8b9e-fe0e11900491"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum = df_sales_daily[['datum', 'M01AB']]\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "1erxTlSYyYmr",
        "outputId": "92f46fd6-64b9-42f8-c5fb-52a4aaddffbc"
      },
      "outputs": [],
      "source": [
        "df_sales_daily_M01AB_datum['Date'] = pd.to_datetime(df_sales_daily_M01AB_datum.datum, format='%m/%d/%Y')\n",
        "df_sales_daily_M01AB_datum.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "tulOAXAeyYhn"
      },
      "outputs": [],
      "source": [
        "# Load data and convert to time series\n",
        "M01AB_ts = pd.Series(df_sales_daily_M01AB_datum.M01AB.values, index=df_sales_daily_M01AB_datum.Date, name='M01AB_Sales')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "klVydQGEyYcb"
      },
      "outputs": [],
      "source": [
        "M01AB_ts.index = pd.DatetimeIndex(M01AB_ts.index, freq=M01AB_ts.index.inferred_freq)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GThVYxo90Z0Q"
      },
      "source": [
        "### Centered moving average (smooth blue line) and trailing moving average (broken orange line) with window w = 12, overlaid on M01AB series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "IuHUZ3190TDS"
      },
      "outputs": [],
      "source": [
        "# Centered moving average with window size = 12\n",
        "ma_centered = M01AB_ts.rolling(12, center=True).mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eFK0hdyt0S-R"
      },
      "outputs": [],
      "source": [
        "# Trailing moving average with window size = 12\n",
        "ma_trailing = M01AB_ts.rolling(12).mean()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 464
        },
        "id": "eqbIhWTq0S0W",
        "outputId": "fa6a33e9-87a2-4ec9-b68c-0acd87478f92"
      },
      "outputs": [],
      "source": [
        "# Shift the average by one time unit to get the next day predictions\n",
        "ma_centered = pd.Series(ma_centered[:-1].values, index=ma_centered.index[1:])\n",
        "ma_trailing = pd.Series(ma_trailing[:-1].values, index=ma_trailing.index[1:])\n",
        "fig, ax = plt.subplots(figsize=(45, 15))\n",
        "ax = M01AB_ts.plot(ax=ax, color='black', linewidth=0.25)\n",
        "ma_centered.plot(ax=ax, linewidth=2)\n",
        "ma_trailing.plot(ax=ax, style='--', linewidth=2)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.legend(['M01AB', 'Centered Moving Average', 'Trailing Moving Average'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "P78VJtuL0emp"
      },
      "source": [
        "### Centered moving average (smooth blue line) and trailing moving average (broken orange line) with window w = 24, overlaid on M01AB series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 464
        },
        "id": "we2qYKFy0gP7",
        "outputId": "fb89738b-1dea-4eb1-8b99-6fd34d024e4d"
      },
      "outputs": [],
      "source": [
        "# Centered moving average with window size = 24\n",
        "ma_centered = M01AB_ts.rolling(24, center=True).mean()\n",
        "\n",
        "# Trailing moving average with window size = 24\n",
        "ma_trailing = M01AB_ts.rolling(24).mean()\n",
        "\n",
        "# Shift the average by one time unit to get the next day predictions\n",
        "ma_centered = pd.Series(ma_centered[:-1].values, index=ma_centered.index[1:])\n",
        "ma_trailing = pd.Series(ma_trailing[:-1].values, index=ma_trailing.index[1:])\n",
        "fig, ax = plt.subplots(figsize=(45, 15))\n",
        "ax = M01AB_ts.plot(ax=ax, color='black', linewidth=0.25)\n",
        "ma_centered.plot(ax=ax, linewidth=2)\n",
        "ma_trailing.plot(ax=ax, style='--', linewidth=2)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.legend(['M01AB', 'Centered Moving Average', 'Trailing Moving Average'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Xs549Zlk2J2c"
      },
      "source": [
        "### Centered moving average (smooth blue line) and trailing moving average (broken orange line) with window w = 36, overlaid on M01AB series"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 464
        },
        "id": "qs30xiFo2Ftk",
        "outputId": "313548fb-d59c-4498-e159-cb8fda74e2c4"
      },
      "outputs": [],
      "source": [
        "# Centered moving average with window size = 36\n",
        "ma_centered = M01AB_ts.rolling(36, center=True).mean()\n",
        "\n",
        "# Trailing moving average with window size = 36\n",
        "ma_trailing = M01AB_ts.rolling(36).mean()\n",
        "\n",
        "# Shift the average by one time unit to get the next day predictions\n",
        "ma_centered = pd.Series(ma_centered[:-1].values, index=ma_centered.index[1:])\n",
        "ma_trailing = pd.Series(ma_trailing[:-1].values, index=ma_trailing.index[1:])\n",
        "fig, ax = plt.subplots(figsize=(45, 15))\n",
        "ax = M01AB_ts.plot(ax=ax, color='black', linewidth=0.25)\n",
        "ma_centered.plot(ax=ax, linewidth=2)\n",
        "ma_trailing.plot(ax=ax, style='--', linewidth=2)\n",
        "ax.set_xlabel('Time')\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.legend(['M01AB', 'Centered Moving Average', 'Trailing Moving Average'])\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bQccWTYg0c5J"
      },
      "source": [
        "### Trailing Moving Average for Forecasting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "X23ZaY-q0SwO",
        "outputId": "020f8a23-3f9b-4d9b-a859-161cadeddecc"
      },
      "outputs": [],
      "source": [
        "len(M01AB_ts)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "gFHrCS070l_6",
        "outputId": "15994b6c-13b9-4910-ef09-dceea0c3c51f"
      },
      "outputs": [],
      "source": [
        "len(M01AB_ts)*0.8"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OL5AKhBI48s5",
        "outputId": "f97a60f7-9203-4538-92c5-7455ae6de8ab"
      },
      "outputs": [],
      "source": [
        "(len(M01AB_ts)-len(M01AB_ts)*0.9)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "RCpmHZXE0l8O"
      },
      "outputs": [],
      "source": [
        "# partition the data\n",
        "nValid = 210\n",
        "nTrain = len(M01AB_ts) - nValid\n",
        "train_ts = M01AB_ts[:nTrain]\n",
        "valid_ts = M01AB_ts[nTrain:]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "tU8lJQBn5VTU",
        "outputId": "1b1e6fe2-601f-4905-ceae-8f18c72880e7"
      },
      "outputs": [],
      "source": [
        "valid_ts.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2TUJ-24U0u1j"
      },
      "source": [
        "### Rolling 12"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4bBKDlyc0l42"
      },
      "outputs": [],
      "source": [
        "# moving average on training\n",
        "ma_trailing = train_ts.rolling(12).mean()\n",
        "last_ma = ma_trailing[-1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 889
        },
        "id": "Wl-up9qA0l1c",
        "outputId": "4d3353d0-6f4c-4d89-c3fd-55863e000aea"
      },
      "outputs": [],
      "source": [
        "# create forecast based on last moving average in the training period\n",
        "ma_trailing_pred = pd.Series(last_ma, index=valid_ts.index)\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(45, 30))\n",
        "ma_trailing.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "ma_trailing_pred.plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - ma_trailing\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - ma_trailing_pred\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ezsGj0IB0w1g"
      },
      "source": [
        "### Rolling 24"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 889
        },
        "id": "oXbnY89b0zg2",
        "outputId": "375c2e04-46be-4d16-bee8-c0cfdce33834"
      },
      "outputs": [],
      "source": [
        "# moving average on training\n",
        "ma_trailing = train_ts.rolling(24).mean()\n",
        "last_ma = ma_trailing[-1]\n",
        "\n",
        "# create forecast based on last moving average in the training period\n",
        "ma_trailing_pred = pd.Series(last_ma, index=valid_ts.index)\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(45, 30))\n",
        "ma_trailing.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "ma_trailing_pred.plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - ma_trailing\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - ma_trailing_pred\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BT0_fOzM0yOS"
      },
      "source": [
        "### Rolling 36"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 889
        },
        "id": "gHbusCid0lxu",
        "outputId": "5c416673-6e4b-4331-ad72-5feb4e3084a2"
      },
      "outputs": [],
      "source": [
        "# moving average on training\n",
        "ma_trailing = train_ts.rolling(36).mean()\n",
        "last_ma = ma_trailing[-1]\n",
        "\n",
        "# create forecast based on last moving average in the training period\n",
        "ma_trailing_pred = pd.Series(last_ma, index=valid_ts.index)\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(45, 30))\n",
        "ma_trailing.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "ma_trailing_pred.plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - ma_trailing\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - ma_trailing_pred\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wfJSxKon05S-"
      },
      "source": [
        "### Applying Moving Average to the residuals from the regression model (which lack trend and seasonality), to forecast the 2019 Mar residual"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "6v5ypVxz0lsj"
      },
      "outputs": [],
      "source": [
        "# Build a model with seasonality, trend, and quadratic trend\n",
        "M01AB_df = tsatools.add_trend(M01AB_ts, trend='ct')\n",
        "M01AB_df['Month'] = M01AB_df.index.month"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jAwthKr21Rqu"
      },
      "outputs": [],
      "source": [
        "# partition the data\n",
        "train_df = M01AB_df[:nTrain]\n",
        "valid_df = M01AB_df[nTrain:]\n",
        "formula = 'M01AB_Sales ~ trend + np.square(trend) + C(Month)'\n",
        "M01AB_lm_trendseason = sm.ols(formula=formula, data=train_df).fit()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ZY5ia9FK1Rlv"
      },
      "outputs": [],
      "source": [
        "# create single-point forecast\n",
        "M01AB_prediction = M01AB_lm_trendseason.predict(valid_df.iloc[0, :])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "98WvU99U1Rhe",
        "outputId": "e1d11da1-cf15-4bc9-f1c2-b807d588d999"
      },
      "outputs": [],
      "source": [
        "# apply Moving Average to residuals\n",
        "ma_trailing = M01AB_lm_trendseason.resid.rolling(12).mean()\n",
        "print('Prediction', M01AB_prediction[0])\n",
        "print('Moving_Average_trailing', ma_trailing[-1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "o6v3ZPEg1Vvn"
      },
      "source": [
        "### Simple Exponential Smoothing"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nGTJPFjG1WzZ"
      },
      "source": [
        "### Output for simple exponential smoothing forecaster with = 0.2, applied to the series of residuals from the regression model (which lack trend and seasonality). The forecast value is 1.62731."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "SCkaeeGb2njD",
        "outputId": "aa1244ed-e510-4f0f-c177-224ba05ca921"
      },
      "outputs": [],
      "source": [
        "residuals_ts.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "r39i48-j2o5M",
        "outputId": "8edfca03-d864-4c59-a673-cf90f6112cc9"
      },
      "outputs": [],
      "source": [
        "residuals_pred.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 463
        },
        "id": "hiJX5Jdc1Raz",
        "outputId": "38dde233-4330-4b72-be11-e24d1273980c"
      },
      "outputs": [],
      "source": [
        "residuals_ts = M01AB_lm_trendseason.resid\n",
        "residuals_pred = valid_df.M01AB_Sales - M01AB_lm_trendseason.predict(valid_df)\n",
        "fig, ax = plt.subplots(figsize=(45,15))\n",
        "M01AB_lm_trendseason.resid.plot(ax=ax, color='black', linewidth=0.5)\n",
        "residuals_pred.plot(ax=ax, color='black', linewidth=0.5)\n",
        "ax.set_ylabel('M01AB')\n",
        "ax.set_xlabel('Time')\n",
        "ax.axhline(y=0, xmin=0, xmax=1, color='grey', linewidth=0.5)\n",
        "\n",
        "# Run Exponential Smoothing\n",
        "# with smoothing level alpha = 0.2\n",
        "expSmooth = ExponentialSmoothing(residuals_ts, freq=None)\n",
        "expSmoothFit = expSmooth.fit(smoothing_level=0.2)\n",
        "expSmoothFit.fittedvalues.plot(ax=ax)\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=ax, style='--', linewidth=2, color='C0')\n",
        "singleGraphLayout(ax, [-7, 13], train_df, valid_df)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5G9OX9et1Y5Y",
        "outputId": "17bc71a1-cfe0-4f4e-822b-5b0b97a41379"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.forecast(len(valid_ts))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KK_rQptx1cWy"
      },
      "source": [
        "### Advanced Exponential Smoothing"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "pVdFfpak_AC6",
        "outputId": "95802632-a14e-44af-d80e-743acc2058a1"
      },
      "outputs": [],
      "source": [
        "train_ts.max()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MYo1AvTCCJ5J"
      },
      "source": [
        "#### Quarterly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 881
        },
        "id": "gAEJBipv1Y1_",
        "outputId": "df89ea0f-b516-4b07-b0ce-d7699a4f7381"
      },
      "outputs": [],
      "source": [
        "# Run exponential smoothing with additive trend and additive seasonal\n",
        "expSmooth = ExponentialSmoothing(train_ts, trend='additive', seasonal='additive', seasonal_periods=4, freq=None)\n",
        "expSmoothFit = expSmooth.fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))\n",
        "expSmoothFit.fittedvalues.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - expSmoothFit.fittedvalues\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - expSmoothFit.forecast(len(valid_ts))\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "H4mnYFltCNBR"
      },
      "source": [
        "#### Yearly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 881
        },
        "id": "RDHica6JCQ-a",
        "outputId": "b0ea353f-4d34-4fba-cc6d-e55743bfa3dd"
      },
      "outputs": [],
      "source": [
        "# Run exponential smoothing with additive trend and additive seasonal\n",
        "expSmooth = ExponentialSmoothing(train_ts, trend='additive', seasonal='additive', seasonal_periods=12, freq=None)\n",
        "expSmoothFit = expSmooth.fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))\n",
        "expSmoothFit.fittedvalues.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - expSmoothFit.fittedvalues\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - expSmoothFit.forecast(len(valid_ts))\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-tDcAvNQCUi7"
      },
      "source": [
        "#### Daily"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 881
        },
        "id": "vB4cklsHCQ5Z",
        "outputId": "541b640e-8f79-4cd3-f2a3-b6ce88b69eef"
      },
      "outputs": [],
      "source": [
        "# Run exponential smoothing with additive trend and additive seasonal\n",
        "expSmooth = ExponentialSmoothing(train_ts, trend='additive', seasonal='additive', seasonal_periods=7, freq=None)\n",
        "expSmoothFit = expSmooth.fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))\n",
        "expSmoothFit.fittedvalues.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - expSmoothFit.fittedvalues\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - expSmoothFit.forecast(len(valid_ts))\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yNdUCUEn_NYa"
      },
      "source": [
        "### Data points must be strictly positive when using multiplicative trend or seasonal components"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2tDngJk6_Rxn"
      },
      "source": [
        "As data points is composed of values = 0 (in reality, it might be so, too), multiplicative trend can't be used in this case. So, while the max value of all pharma products is roughly 17.0, the min value of them all can be adjusted to a small figure, such as equal to or less than 0.01, to prepare data for this method using multiplcative trend, as long as they are not 0 (must be positive)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 365
        },
        "id": "SvHQvUfv_H65",
        "outputId": "6c61da16-dd99-4fa9-d20b-eadb068a9b02"
      },
      "outputs": [],
      "source": [
        "# Run exponential smoothing with multiplicative trend and multiplicative seasonal\n",
        "expSmooth = ExponentialSmoothing(train_ts, trend='multiplicative', seasonal='multiplicative', seasonal_periods=12, freq=None)\n",
        "expSmoothFit = expSmooth.fit()\n",
        "fig, axes = plt.subplots(nrows=2, ncols=1, figsize=(20, 15))\n",
        "expSmoothFit.fittedvalues.plot(ax=axes[0], linewidth=2, color='C1')\n",
        "expSmoothFit.forecast(len(valid_ts)).plot(ax=axes[0], linewidth=2, color='C1', linestyle='dashed')\n",
        "residual = train_ts - expSmoothFit.fittedvalues\n",
        "residual.plot(ax=axes[1], color='C1')\n",
        "residual = valid_ts - expSmoothFit.forecast(len(valid_ts))\n",
        "residual.plot(ax=axes[1], color='C1', linestyle='dashed')\n",
        "graphLayout(axes, train_ts, valid_ts)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "pzHTU6Lh_j3q",
        "outputId": "5c31e107-8ba6-461f-f158-70630c67ecff"
      },
      "outputs": [],
      "source": [
        "df_sales_daily.head()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_cFKtI3u1e18"
      },
      "source": [
        "### Summary of a Holt-Winter’s exponential smoothing model applied to the M01AB data. Included are the initial and final states\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "TBjoK83X1YyT",
        "outputId": "e0d9210b-a68e-4d00-ada2-ae013c84ebe3"
      },
      "outputs": [],
      "source": [
        "print(expSmoothFit.params)\n",
        "print('AIC: ', expSmoothFit.aic)\n",
        "print('AICc: ', expSmoothFit.aicc)\n",
        "print('BIC: ', expSmoothFit.bic)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 482
        },
        "id": "Kr7SwhfmHEMJ",
        "outputId": "52f291f1-ffb6-4d01-9c0d-0c188caf2a91"
      },
      "outputs": [],
      "source": [
        "expSmoothFit.summary()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XUg9t2t2fBoT"
      },
      "source": [
        "# Time Series Forecasting with Prophet"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "oPrrbvd_ok4h",
        "outputId": "97fb41c4-49b8-41db-88e1-fbb26dd56937"
      },
      "outputs": [],
      "source": [
        "df_sales_daily"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "vRaeM-zefLXS",
        "outputId": "391cb4ef-752a-40c2-821a-ffdc8d371ded"
      },
      "outputs": [],
      "source": [
        "df = df_sales_daily[[\"datum\", \"M01AB\"]]\n",
        "df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 363
        },
        "id": "C-Ezd0ksfu3P",
        "outputId": "4207c1ba-c63b-4e2b-92fc-ed986da29839"
      },
      "outputs": [],
      "source": [
        "df['ds'] = pd.to_datetime(df.datum, format='%m/%d/%Y')\n",
        "df.head(10)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "LXQ-eOTfgMXR",
        "outputId": "ad93a8b5-5961-40e0-9f0f-e659a094f6e5"
      },
      "outputs": [],
      "source": [
        "df.rename({'M01AB': 'y'}, axis=1, inplace=True)\n",
        "df"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "0zvCeBQFglg6",
        "outputId": "72f1904b-8685-4154-c771-87a976d26c06"
      },
      "outputs": [],
      "source": [
        "df = df[['ds', 'y']]\n",
        "df"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v8hmZhXZo47e"
      },
      "source": [
        "## Sub-daily data\n",
        "\n",
        "Prophet can make forecasts for time series with sub-daily observations by passing in a dataframe with timestamps in the ds column. The format of the timestamps should be YYYY-MM-DD HH:MM:SS. When sub-daily data are used, daily seasonality will automatically be fit. "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 620
        },
        "id": "VEEmiFltok1i",
        "outputId": "b803034d-49ff-4865-e13a-3fabc6099284"
      },
      "outputs": [],
      "source": [
        "m = Prophet(changepoint_prior_scale=0.01).fit(df)\n",
        "future = m.make_future_dataframe(periods=300, freq='H')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fUpERIp4pEoq"
      },
      "source": [
        "The daily seasonality will show up in the components plot:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 657
        },
        "id": "QkQPmOdXgdCh",
        "outputId": "60639cdd-d2d0-4829-f562-861eec8ed46e"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "YqhyzWKmDp1W",
        "outputId": "6bded549-900d-473d-8839-635ca1e50939"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 620
        },
        "id": "cP3iUBvjKy_V",
        "outputId": "3134ca66-b985-40b7-b721-cad6ba6219d2"
      },
      "outputs": [],
      "source": [
        "m = Prophet(changepoint_prior_scale=0.01).fit(df)\n",
        "future = m.make_future_dataframe(periods=150, freq='M')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 657
        },
        "id": "5N6N5otlLhLb",
        "outputId": "d7b6cd85-af64-49e8-de3e-29a6fa10b909"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 620
        },
        "id": "_KxZwX4GLNTV",
        "outputId": "c7a30467-f056-4a06-da60-81958bfff879"
      },
      "outputs": [],
      "source": [
        "m = Prophet(changepoint_prior_scale=0.1).fit(df)\n",
        "future = m.make_future_dataframe(periods=150, freq='M')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 657
        },
        "id": "0JnVahoNLj-h",
        "outputId": "0cb1eb3e-8055-4bdf-d969-5b955b00bcef"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 620
        },
        "id": "aTXgLsrjLZcw",
        "outputId": "8adcb4e4-50f0-42bd-a122-fb53916ac23a"
      },
      "outputs": [],
      "source": [
        "m = Prophet(changepoint_prior_scale=0.001).fit(df)\n",
        "future = m.make_future_dataframe(periods=150, freq='M')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 657
        },
        "id": "gMkGwR_EK6Gx",
        "outputId": "eee64857-0e13-4976-c395-06da4dedfca4"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mZfGGL6FpJ0A"
      },
      "source": [
        "## Data with regular gaps"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 427
        },
        "id": "fYvIy_LqgyIK",
        "outputId": "d2b596b7-6231-40ae-d50d-eaf455f510a3"
      },
      "outputs": [],
      "source": [
        "df2 = df.copy()\n",
        "df2['ds'] = pd.to_datetime(df2['ds'])\n",
        "df2 = df2[df2['ds'].dt.hour < 6]\n",
        "m = Prophet().fit(df2)\n",
        "future = m.make_future_dataframe(periods=300, freq='H')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CLg08By5pSgx"
      },
      "source": [
        "The forecast seems quite poor, with much larger fluctuations in the future than were seen in the history. The issue here is that we have fit a daily cycle to a time series that only has data for part of the day (12a to 6a). The daily seasonality is thus unconstrained for the remainder of the day and is not estimated well. The solution is to only make predictions for the time windows for which there are historical data. Here, that means to limit the future dataframe to have times from 12a to 6a:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "zv6H-cw3DhpF",
        "outputId": "d2c6715b-4ed1-4ed6-ea49-3a01a0d4a5d2"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 392
        },
        "id": "PBt6WFS2g34D",
        "outputId": "02cf9496-23d6-4956-979c-4cbff6a4943a"
      },
      "outputs": [],
      "source": [
        "future2 = future.copy()\n",
        "future2 = future2[future2['ds'].dt.hour < 6]\n",
        "fcst = m.predict(future2)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ywT0va7npjzz"
      },
      "source": [
        "The same principle applies to other datasets with regular gaps in the data. For example, if the history contains only weekdays, then predictions should only be made for weekdays since the weekly seasonality will not be well estimated for the weekends."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "BE6YVCmDDz1C",
        "outputId": "970402b2-2526-47b9-b850-c43aed2c3bea"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "RXuaLkFbpmiy"
      },
      "source": [
        "## Monthly data\n",
        "\n",
        "I can use Prophet to fit monthly data. However, the underlying model is continuous-time, which means that I can get strange results if I fit the model to monthly data and then ask for daily forecasts. Here I forecast Pharma sales volume for the next 10 years:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "NUSYNluuKchg"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='multiplicative').fit(df)\n",
        "future = m.make_future_dataframe(periods=3652)\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 427
        },
        "id": "R4dAcxhig7kx",
        "outputId": "a600bc44-7d94-4b26-a689-2333e688c851"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='multiplicative').fit(df)\n",
        "future = m.make_future_dataframe(periods=3652)\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "T79PEZu1pu_n"
      },
      "source": [
        "This is the same issue from above where the dataset has regular gaps. When I fit the yearly seasonality, it only has data for the first of each month and the seasonality components for the remaining days are unidentifiable and overfit. This can be clearly seen by doing MCMC to see uncertainty in the seasonality:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "TrD77TeUD1x4",
        "outputId": "83c6f7db-6ce6-4581-ebad-95f111a73e1e"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "id": "3v0Z5LbShC_s",
        "outputId": "a8615ee7-6892-4302-c07e-66d83f7bdf92"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='multiplicative', mcmc_samples=300).fit(df, show_progress=False)\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot_components(fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "QfCTEGl46iyR"
      },
      "outputs": [],
      "source": [
        "m"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SLZBalAup6gv"
      },
      "source": [
        "The seasonality has low uncertainty at the start of each month where there are data points, but has very high posterior variance in between. When fitting Prophet to monthly data, only make monthly forecasts, which can be done by passing the frequency into make_future_dataframe:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "xZw5qA17D4Br",
        "outputId": "49cf179d-78c9-4c06-e3fb-7ad913a2951c"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "iDlvKaX2hJqq",
        "outputId": "c158c9ea-a953-4eaa-e68d-214e3274c1e5"
      },
      "outputs": [],
      "source": [
        "future = m.make_future_dataframe(periods=120, freq='MS')\n",
        "fcst = m.predict(future)\n",
        "fig = m.plot(fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "adjwJDkYqFgX"
      },
      "source": [
        "In Python, the frequency can be anything from the pandas list of frequency strings here: https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timeseries-offset-aliases . Note that MS used here is month-start, meaning the data point is placed on the start of each month.\n",
        "\n",
        "In monthly data, yearly seasonality can also be modeled with binary extra regressors. In particular, the model can use 12 extra regressors like is_jan, is_feb, etc. where is_jan is 1 if the date is in Jan and 0 otherwise. This approach would avoid the within-month unidentifiability seen above. Be sure to use yearly_seasonality=False if monthly extra regressors are being added."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "-qTkXSAcD49S",
        "outputId": "c52871e8-6251-4102-ab83-a54232e637ff"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=fcst)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "_671kDaDqHJL"
      },
      "source": [
        "## Holidays with aggregated data\n",
        "Holiday effects are applied to the particular date on which the holiday was specified. With data that has been aggregated to weekly or monthly frequency, holidays that don't fall on the particular date used in the data will be ignored: for example, a Monday holiday in a weekly time series where each data point is on a Sunday. To include holiday effects in the model, the holiday will need to be moved to the date in the history dataframe for which the effect is desired. Note that with weekly or monthly aggregated data, many holiday effects will be well-captured by the yearly seasonality, so added holidays may only be necessary for holidays that occur in different weeks throughout the time series."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L7ww7fvpAflt"
      },
      "source": [
        "## Quick Start"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5jmIAzlHAhkX",
        "outputId": "5339145a-5648-4d1f-9a9b-d6b60b031d8b"
      },
      "outputs": [],
      "source": [
        "m = Prophet()\n",
        "m.fit(df)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "QZyBs9fNAhh7",
        "outputId": "309dea0b-5afd-4959-c448-bf3ada0ac592"
      },
      "outputs": [],
      "source": [
        "future = m.make_future_dataframe(periods=365)\n",
        "future.tail()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 206
        },
        "id": "ZnI9JUAmAhfX",
        "outputId": "9362350e-14d9-4009-c04e-810823cdb0bf"
      },
      "outputs": [],
      "source": [
        "forecast = m.predict(future)\n",
        "forecast[['ds', 'yhat', 'yhat_lower', 'yhat_upper']].tail()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 391
        },
        "id": "1aSFVor3Ahco",
        "outputId": "52578538-a78f-4f8d-f688-606ced9d8569"
      },
      "outputs": [],
      "source": [
        "fig1 = m.plot(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 433
        },
        "id": "ABHQi4zUAt57",
        "outputId": "d8f85ec8-47fb-403c-f049-5f9f22d0876d"
      },
      "outputs": [],
      "source": [
        "fig2 = m.plot_components(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WhjB3-zPBmur",
        "outputId": "0f88e02f-c10e-47de-c31f-47ecce307964"
      },
      "outputs": [],
      "source": [
        "!pip install chart-studio"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xNREfyiMBv1_"
      },
      "outputs": [],
      "source": [
        "from chart_studio.plotly import plot, iplot"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 637
        },
        "id": "BhUaErF3At2q",
        "outputId": "849d6d12-9735-4c31-87d1-f2187ee1501f"
      },
      "outputs": [],
      "source": [
        "from prophet.plot import plot_plotly, plot_components_plotly\n",
        "\n",
        "plot_plotly(m, forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "QVsdBwm6CcQH",
        "outputId": "8bc6cbfe-1b97-492d-bc6a-139803a000e8"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 437
        },
        "id": "nj2DQT5kAtzd",
        "outputId": "6a64d956-daa2-4ea4-fd88-7903f84254d0"
      },
      "outputs": [],
      "source": [
        "plot_components_plotly(m, forecast)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Sdk4Bw9w-7ly"
      },
      "source": [
        "## Uncertainty Intervals"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 424
        },
        "id": "5LZh9QQt_JzT",
        "outputId": "4b9357e8-b58a-436b-f73b-0cc363f6e9bc"
      },
      "outputs": [],
      "source": [
        "df_sales_monthly"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-9QzWIkThmec",
        "outputId": "0d5b22ec-1aed-40da-8e97-ac6dbf4fce93"
      },
      "outputs": [],
      "source": [
        "df = df.loc[:180,]  # Limit to first six months\n",
        "m = Prophet()\n",
        "m.fit(df)\n",
        "future = m.make_future_dataframe(periods=60)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "PkQzDjmc_Sj9",
        "outputId": "fb763e88-c83c-4f87-d7c1-bff1f5451090"
      },
      "outputs": [],
      "source": [
        "forecast = Prophet(interval_width=0.95).fit(df).predict(future)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "0ACP8qi__pUt",
        "outputId": "0f126b5c-06c0-4238-bcbe-f655fb9a156e"
      },
      "outputs": [],
      "source": [
        "m = Prophet(mcmc_samples=300)\n",
        "forecast = m.fit(df, show_progress=False).predict(future)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 429
        },
        "id": "uPrIHpbh_sYD",
        "outputId": "791c2732-6290-4e07-b5de-93fc3b679f53"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 661
        },
        "id": "JINc4onCDEaN",
        "outputId": "66afea3b-90ab-46e2-a8c5-34ec0f4b7329"
      },
      "outputs": [],
      "source": [
        "forecast"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 542
        },
        "id": "Vm-tgB_tC7oH",
        "outputId": "e1bedb5e-a947-4c94-bf59-a1d091b799bf"
      },
      "outputs": [],
      "source": [
        "import plotly.graph_objects as go\n",
        "from plotly.subplots import make_subplots\n",
        "fig = make_subplots(rows=1, cols=1, specs=[[{\"type\": \"surface\"}]])\n",
        "fig.add_surface(z=forecast)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2VxlbC2aALdh"
      },
      "source": [
        "## Multiplicative Seasonality"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 427
        },
        "id": "L4ndLTZu_v0L",
        "outputId": "00edd8ae-e9df-46da-b125-655842111c5b"
      },
      "outputs": [],
      "source": [
        "m = Prophet()\n",
        "m.fit(df)\n",
        "future = m.make_future_dataframe(150, freq='MS')\n",
        "forecast = m.predict(future)\n",
        "fig = m.plot(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 427
        },
        "id": "Bzzx0nY2EUH6",
        "outputId": "478c4a62-658c-4212-ad1d-c85b282146aa"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='multiplicative')\n",
        "m.fit(df)\n",
        "forecast = m.predict(future)\n",
        "fig = m.plot(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 645
        },
        "id": "a7fj0e1wEWTR",
        "outputId": "1e8166f3-1b8e-4d6d-b4fe-af3147f2933a"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "eigFibh5EgoQ",
        "outputId": "0c3ca13e-1f89-4580-c26e-c7eb5c5bd6e0"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='multiplicative')\n",
        "m.add_seasonality('quarterly', period=91.25, fourier_order=8, mode='additive')\n",
        "m.add_regressor('regressor', mode='additive')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 620
        },
        "id": "fOnomKgYEqPp",
        "outputId": "210d3e5e-b4f0-4368-c4ed-d8f4202c8195"
      },
      "outputs": [],
      "source": [
        "m = Prophet(seasonality_mode='additive')\n",
        "m.fit(df)\n",
        "forecast = m.predict(future)\n",
        "fig = m.plot(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 657
        },
        "id": "58aAYw7HMKLN",
        "outputId": "723ad19c-2099-46c5-fcf8-861c5c20cc3e"
      },
      "outputs": [],
      "source": [
        "fig = m.plot_components(forecast)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "S5SI8zFnMgD5"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "collapsed_sections": [],
      "provenance": [],
      "toc_visible": true
    },
    "gpuClass": "standard",
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}