{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "0893f116",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from prophet import Prophet\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import logging\n",
    "logging.getLogger('prophet').setLevel(logging.ERROR)\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "plt.rcParams['figure.figsize'] = 9, 6"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f84512d1",
   "metadata": {},
   "source": [
    "As a result of the lockdowns caused by the COVID-19 pandemic, many time series experienced \"shocks\" during 2020, e.g. spikes in media consumption (Netflix, YouTube), e-commerce transactions (Amazon, eBay), whilst attendance to in-person events declined dramatically.\n",
    "\n",
    "Most of these time series would also maintain their new level for a period of time, subject to fluctuations driven by easing of lockdowns and/or vaccines.\n",
    "\n",
    "Seasonal patterns could have also changed: for example, people may have consumed less media (in total hours) on weekdays compared to weekends before the COVID lockdowns, but during lockdowns weekday consumption could be much closer to weekend consumption.\n",
    "\n",
    "In this page we'll explore some strategies for capturing these effects using Prophet's functionality:\n",
    "\n",
    "1. Marking step changes / spikes due to COVID events as once-off.\n",
    "2. Sustained changes in behaviour leading to trend and seasonality changes."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36216eb9",
   "metadata": {},
   "source": [
    "#### Case Study - Pedestrian Activity\n",
    "\n",
    "For this case study we'll use [Pedestrian Sensor data from the City of Melbourne](https://data.melbourne.vic.gov.au/Transport/Pedestrian-Counting-System-Monthly-counts-per-hour/b2ak-trbp). This data measures foot traffic from sensors in various places in the central business district, and we've chosen one sensor (`Sensor_ID = 4`) and aggregated the values to a daily grain. \n",
    "\n",
    "The aggregated dataset can be found in the examples folder [here](https://github.com/facebook/prophet/tree/master/examples/example_pedestrians_covid.csv)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ac3aa527",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_pedestrians_covid.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "284f936b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 47760 characters]",
      "text/plain": [
       "<Figure size 648x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.set_index('ds').plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6ed9e13",
   "metadata": {},
   "source": [
    "We can see two key events in the time series:\n",
    "\n",
    "* The initial drop in foot traffic around 21 March 2020, which started to recover around 6 June 2020. This corresponds to the declaration of a pandemic by WHO and subsequent lockdowns mandated by the Victorian government.\n",
    "* After some slow recovery, a second drop in foot traffic around July 9 2020, which began to recover around 27 October 2020. This corresponds to the \"second wave\" of the pandemic in metropolitan Melbourne.\n",
    "\n",
    "There are also shorter periods of strict lockdown that lead to sudden tips in the time series: 5 days in February 2021, and 14 days in early June 2021."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea888620",
   "metadata": {},
   "source": [
    "#### Default model without any adjustments\n",
    "\n",
    "First we'll fit a model with the default Prophet settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f94487ef",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02:53:41 - cmdstanpy - INFO - Chain [1] start processing\n",
      "02:53:41 - cmdstanpy - INFO - Chain [1] done processing\n"
     ]
    }
   ],
   "source": [
    "m = Prophet()\n",
    "m = m.fit(df)\n",
    "future = m.make_future_dataframe(periods=366)\n",
    "forecast = m.predict(future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ef22e986",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 165916 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m.plot(forecast)\n",
    "plt.axhline(y=0, color='red')\n",
    "plt.title('Default Prophet');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "bd0a2b76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 72696 characters]",
      "text/plain": [
       "<Figure size 648x648 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m.plot_components(forecast);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "709a0e93",
   "metadata": {},
   "source": [
    "The model seems to fit reasonably to past data, but notice how we're capturing the dips, and the spikes after the dips, as a part of the trend component. \n",
    "\n",
    "By default, the model assumes that these large spikes are possible in the future, even though we realistically won't see something of the same magnitude within our forecast horizon (1 year in this case). This leads to a fairly optimistic forecast of the recovery of foot traffic in 2022."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff5f2833",
   "metadata": {},
   "source": [
    "### Treating COVID-19 lockdowns as a one-off holidays\n",
    "\n",
    "To prevent large dips and spikes from being captured by the trend component, we can treat the days impacted by COVID-19 as holidays that will not repeat again in the future. Adding custom holidays is explained in more detail [here](https://facebook.github.io/prophet/docs/seasonality,_holiday_effects,_and_regressors.html#modeling-holidays-and-special-events). We set up a DataFrame like so to describe the periods affected by lockdowns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e566cf71",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>holiday</th>\n",
       "      <th>ds</th>\n",
       "      <th>lower_window</th>\n",
       "      <th>ds_upper</th>\n",
       "      <th>upper_window</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>lockdown_1</td>\n",
       "      <td>2020-03-21</td>\n",
       "      <td>0</td>\n",
       "      <td>2020-06-06</td>\n",
       "      <td>77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>lockdown_2</td>\n",
       "      <td>2020-07-09</td>\n",
       "      <td>0</td>\n",
       "      <td>2020-10-27</td>\n",
       "      <td>110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>lockdown_3</td>\n",
       "      <td>2021-02-13</td>\n",
       "      <td>0</td>\n",
       "      <td>2021-02-17</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>lockdown_4</td>\n",
       "      <td>2021-05-28</td>\n",
       "      <td>0</td>\n",
       "      <td>2021-06-10</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      holiday         ds  lower_window   ds_upper  upper_window\n",
       "0  lockdown_1 2020-03-21             0 2020-06-06            77\n",
       "1  lockdown_2 2020-07-09             0 2020-10-27           110\n",
       "2  lockdown_3 2021-02-13             0 2021-02-17             4\n",
       "3  lockdown_4 2021-05-28             0 2021-06-10            13"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lockdowns = pd.DataFrame([\n",
    "    {'holiday': 'lockdown_1', 'ds': '2020-03-21', 'lower_window': 0, 'ds_upper': '2020-06-06'},\n",
    "    {'holiday': 'lockdown_2', 'ds': '2020-07-09', 'lower_window': 0, 'ds_upper': '2020-10-27'},\n",
    "    {'holiday': 'lockdown_3', 'ds': '2021-02-13', 'lower_window': 0, 'ds_upper': '2021-02-17'},\n",
    "    {'holiday': 'lockdown_4', 'ds': '2021-05-28', 'lower_window': 0, 'ds_upper': '2021-06-10'},\n",
    "])\n",
    "for t_col in ['ds', 'ds_upper']:\n",
    "    lockdowns[t_col] = pd.to_datetime(lockdowns[t_col])\n",
    "lockdowns['upper_window'] = (lockdowns['ds_upper'] - lockdowns['ds']).dt.days\n",
    "lockdowns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "338f8e50",
   "metadata": {},
   "source": [
    "* We have an entry for each lockdown period, with `ds` specifying the start of the lockdown. `ds_upper` is not used by Prophet, but it's a convenient way for us to calculate `upper_window`.\n",
    "* `upper_window` tells Prophet that the lockdown spans for x days after the start of the lockdown. Note that the holidays regression is inclusive of the upper bound.\n",
    "\n",
    "Note that since we don't specify any future dates, Prophet will assume that these holidays will _not_ occur again when creating the future dataframe (and hence they won't affect our projections). This is different to how we would specify a recurring holiday."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "61ac656c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02:53:44 - cmdstanpy - INFO - Chain [1] start processing\n",
      "02:53:45 - cmdstanpy - INFO - Chain [1] done processing\n"
     ]
    }
   ],
   "source": [
    "m2 = Prophet(holidays=lockdowns)\n",
    "m2 = m2.fit(df)\n",
    "future2 = m2.make_future_dataframe(periods=366)\n",
    "forecast2 = m2.predict(future2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6ea726f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 148728 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m2.plot(forecast2)\n",
    "plt.axhline(y=0, color='red')\n",
    "plt.title('Lockdowns as one-off holidays');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "27f83233",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 94688 characters]",
      "text/plain": [
       "<Figure size 648x864 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m2.plot_components(forecast2);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "227a7af2",
   "metadata": {},
   "source": [
    "* Prophet is sensibly assigning a large negative effect to the days within the lockdown periods.\n",
    "* The forecast for the trend is not as strong / optimistic and seems fairly reasonable.\n",
    "\n",
    "### Sense checking the trend\n",
    "\n",
    "In an environment when behaviours are constantly changing, it's important to ensure that the trend component of the model is capturing to emerging patterns without overfitting to them.\n",
    "\n",
    "The [trend changepoints](https://facebook.github.io/prophet/docs/trend_changepoints.html) documentation explains two things we could tweak with the trend component:\n",
    "\n",
    "* The changepoint locations, which by default are evenly spaced across 80% of the history. We should be mindful of where this range ends, and extend the range (either by increasing the % or adding changepoints manually) if we believe the most recent data better reflects future behaviour.\n",
    "* The strength of regularisation (`changepoint_prior_scale`), which determines how flexible the trend is; the default value is `0.05` and increasing this will allow the trend to fit more closely to the observed data.\n",
    "\n",
    "We plot the trend component and changepoints detected by our current model below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "f3675c69",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 158560 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from prophet.plot import add_changepoints_to_plot\n",
    "fig = m2.plot(forecast2)\n",
    "a = add_changepoints_to_plot(fig.gca(), m2, forecast2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6cd90ff",
   "metadata": {},
   "source": [
    "The detected changepoints look reasonable, and the future trend tracks the latest upwards trend in activity, but not to the extent of late 2020. This seems suitable for a best guess of future activity. \n",
    "\n",
    "We can see what the forecast would look like if we wanted to emphasize COVID patterns more in model training; we can do this by adding more potential changepoints after 2020 and making the trend more flexible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "85ccc556",
   "metadata": {},
   "outputs": [],
   "source": [
    "m3_changepoints = (\n",
    "    # 10 potential changepoints in 2.5 years\n",
    "    pd.date_range('2017-06-02', '2020-01-01', periods=10).date.tolist() + \n",
    "    # 15 potential changepoints in 1 year 2 months\n",
    "    pd.date_range('2020-02-01', '2021-04-01', periods=15).date.tolist()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6df44e0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02:53:49 - cmdstanpy - INFO - Chain [1] start processing\n",
      "02:53:52 - cmdstanpy - INFO - Chain [1] done processing\n"
     ]
    }
   ],
   "source": [
    "# Default changepoint_prior_scale is 0.05, so 1.0 will lead to much more flexibility in comparison.\n",
    "m3 = Prophet(holidays=lockdowns, changepoints=m3_changepoints, changepoint_prior_scale=1.0)\n",
    "m3 = m3.fit(df)\n",
    "forecast3 = m3.predict(future2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "effcdd7e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 166392 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from prophet.plot import add_changepoints_to_plot\n",
    "fig = m3.plot(forecast3)\n",
    "a = add_changepoints_to_plot(fig.gca(), m3, forecast3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab8ce52d",
   "metadata": {},
   "source": [
    "We're seeing many changepoints detected post-COVID, matching the various fluctuations from loosening / tightening lockdowns. Overall the trend curve and forecasted trend are quite similar to our previous model, but we're seeing a lot more uncertainty because of the higher number of trend changes we picked up in the history.\n",
    "\n",
    "We probably wouldn't pick this model over the model with default parameters as a best estimate, but it's a good demonstration of how we can incorporate our beliefs into the model about which patterns are important to capture."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dd2dc42",
   "metadata": {},
   "source": [
    "### Changes in seasonality between pre- and post-COVID\n",
    "\n",
    "The seasonal component plots in the previous sections show a peak of activity on Friday compared to other days of the week. If we're not sure whether this will still hold post-lockdown, we can add _conditional seasonalities_ to the model. Conditional seasonalities are explained in more detail [here](https://facebook.github.io/prophet/docs/seasonality,_holiday_effects,_and_regressors.html#seasonalities-that-depend-on-other-factors).\n",
    "\n",
    "First we define boolean columns in the history dataframe to flag \"pre covid\" and \"post covid\" periods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "75559257",
   "metadata": {},
   "outputs": [],
   "source": [
    "df2 = df.copy()\n",
    "df2['pre_covid'] = pd.to_datetime(df2['ds']) < pd.to_datetime('2020-03-21')\n",
    "df2['post_covid'] = ~df2['pre_covid']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b75a9e13",
   "metadata": {},
   "source": [
    "The conditional seasonality we're interested in modelling here is the day-of-week (\"weekly\") seasonality. To do this, we firstly turn off the default `weekly_seasonality` when we create the Prophet model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "cd65ad72",
   "metadata": {},
   "outputs": [],
   "source": [
    "m4 = Prophet(holidays=lockdowns, weekly_seasonality=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b63e33",
   "metadata": {},
   "source": [
    "We then add this weekly seasonality manually, as two different model components - one for pre-covid, one for post-covid. Note that `fourier_order=3` is the default setting for weekly seasonality. After this we can run `.fit()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "538a2914",
   "metadata": {},
   "outputs": [],
   "source": [
    "m4.add_seasonality(\n",
    "    name='weekly_pre_covid',\n",
    "    period=7,\n",
    "    fourier_order=3,\n",
    "    condition_name='pre_covid',\n",
    ")\n",
    "m4.add_seasonality(\n",
    "    name='weekly_post_covid',\n",
    "    period=7,\n",
    "    fourier_order=3,\n",
    "    condition_name='post_covid',\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "10614e70",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "02:53:55 - cmdstanpy - INFO - Chain [1] start processing\n",
      "02:53:56 - cmdstanpy - INFO - Chain [1] done processing\n"
     ]
    }
   ],
   "source": [
    "m4 = m4.fit(df2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee8bed18",
   "metadata": {},
   "source": [
    "We also need to create the `pre_covid` and `post_covid` flags in the future dataframe. This is so that Prophet can apply the correct weekly seasonality parameters to each future date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b65e2ec7",
   "metadata": {},
   "outputs": [],
   "source": [
    "future4 = m4.make_future_dataframe(periods=366)\n",
    "future4['pre_covid'] = pd.to_datetime(future4['ds']) < pd.to_datetime('2020-03-21')\n",
    "future4['post_covid'] = ~future4['pre_covid']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b7dc89f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "forecast4 = m4.predict(future4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "4dc6c2d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 147452 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m4.plot(forecast4)\n",
    "plt.axhline(y=0, color='red')\n",
    "plt.title('Lockdowns as one-off holidays + Conditional weekly seasonality');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "47921893",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 119520 characters]",
      "text/plain": [
       "<Figure size 648x1080 with 5 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m4.plot_components(forecast4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1234edf",
   "metadata": {},
   "source": [
    "Interestingly, the model with conditional seasonalities suggests that, post-COVID, pedestrian activity peaks on Saturdays, instead of Fridays. This could make sense if most people are still working from home and are hence less likely to go out on Friday nights. From a prediction perspective this would only be important if we care about predicting weekdays vs. weekends accurately, but overall this kind of exploration helps us gain insight into how COVID has changed behaviours."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b44dbd7",
   "metadata": {},
   "source": [
    "#### Further reading\n",
    "\n",
    "A lot of the content in this page was inspired by this [GitHub discussion](https://github.com/facebook/prophet/issues/1416). We've covered a few low hanging fruits for tweaking Prophet models when faced with shocks such as COVID, but there are many other possible approaches as well, such as:\n",
    "\n",
    "* Using external regressors (e.g. the lockdown stringency index). This would only be fruitful if we a) have regressor data that aligns well (in terms of location) with the series we're forecasting, and b) have control over or can predict the regressor much more accurately than the time series alone.\n",
    "* Detecting and removing outlier data from the training period, or throwing away older training data completely. This might be a better approach for sub-daily time series that don't have yearly seasonal patterns.\n",
    "\n",
    "Overall though it's difficult to be confident in our forecasts in these environments when rules are constantly changing and outbreaks occur randomly. In this scenario it's more important to constantly re-train / re-evaluate our models and clearly communicate the increased uncertainty in forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa9f4101",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.6 ('python-P8zs_nH7')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "b89c51048a83bff423ddf9419a4921628c4e382d20a75037c2900b180e09c6be"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
