{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "block_hidden": true
   },
   "outputs": [],
   "source": [
    "%load_ext rpy2.ipython\n",
    "%matplotlib inline\n",
    "from prophet import Prophet\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "import pandas as pd\n",
    "import logging\n",
    "import warnings\n",
    "\n",
    "logging.getLogger('prophet').setLevel(logging.ERROR)\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "block_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO:numexpr.utils:NumExpr defaulting to 8 threads.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8395cd1a302349a6b766e256fff7aa02",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=3.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv')\n",
    "m = Prophet()\n",
    "m.fit(df)\n",
    "future = m.make_future_dataframe(periods=366)\n",
    "\n",
    "from prophet.diagnostics import cross_validation\n",
    "df_cv = cross_validation(\n",
    "    m, '365 days', initial='1825 days', period='365 days')\n",
    "cutoff = df_cv['cutoff'].unique()[0]\n",
    "df_cv = df_cv[df_cv['cutoff'].values == cutoff]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "block_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Loading required package: Rcpp\n",
      "\n",
      "R[write to console]: Loading required package: rlang\n",
      "\n",
      "R[write to console]: Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "library(prophet)\n",
    "df <- read.csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv')\n",
    "m <- prophet(df)\n",
    "future <- make_future_dataframe(m, periods=366)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross validation\n",
    "\n",
    "Prophet includes functionality for time series cross validation to measure forecast error using historical data. This is done by selecting cutoff points in the history, and for each of them fitting the model using data only up to that cutoff point. We can then compare the forecasted values to the actual values. This figure illustrates a simulated historical forecast on the Peyton Manning dataset, where the model was fit to an initial history of 5 years, and a forecast was made on a one year horizon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "input_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 70564 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(facecolor='w', figsize=(10, 6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(m.history['ds'].values, m.history['y'], 'k.')\n",
    "ax.plot(df_cv['ds'].values, df_cv['yhat'], ls='-', c='#0072B2')\n",
    "ax.fill_between(df_cv['ds'].values, df_cv['yhat_lower'],\n",
    "                df_cv['yhat_upper'], color='#0072B2',\n",
    "                alpha=0.2)\n",
    "ax.axvline(x=pd.to_datetime(cutoff), c='gray', lw=4, alpha=0.5)\n",
    "ax.set_ylabel('y')\n",
    "ax.set_xlabel('ds')\n",
    "ax.text(x=pd.to_datetime('2010-01-01'),y=12, s='Initial', color='black',\n",
    "       fontsize=16, fontweight='bold', alpha=0.8)\n",
    "ax.text(x=pd.to_datetime('2012-08-01'),y=12, s='Cutoff', color='black',\n",
    "       fontsize=16, fontweight='bold', alpha=0.8)\n",
    "ax.axvline(x=pd.to_datetime(cutoff) + pd.Timedelta('365 days'), c='gray', lw=4,\n",
    "           alpha=0.5, ls='--')\n",
    "ax.text(x=pd.to_datetime('2013-01-01'),y=6, s='Horizon', color='black',\n",
    "       fontsize=16, fontweight='bold', alpha=0.8);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[The Prophet paper](https://peerj.com/preprints/3190.pdf) gives further description of simulated historical forecasts.\n",
    "\n",
    "This cross validation procedure can be done automatically for a range of historical cutoffs using the `cross_validation` function. We specify the forecast horizon (`horizon`), and then optionally the size of the initial training period (`initial`) and the spacing between cutoff dates (`period`). By default, the initial training period is set to three times the horizon, and cutoffs are made every half a horizon.\n",
    "\n",
    "The output of `cross_validation` is a dataframe with the true values `y` and the out-of-sample forecast values `yhat`, at each simulated forecast date and for each cutoff date. In particular, a forecast is made for every observed point between `cutoff` and `cutoff + horizon`. This dataframe can then be used to compute error measures of `yhat` vs. `y`.\n",
    "\n",
    "Here we do cross-validation to assess prediction performance on a horizon of 365 days, starting with 730 days of training data in the first cutoff and then making predictions every 180 days. On this 8 year time series, this corresponds to 11 total forecasts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Making 11 forecasts with cutoffs between 2010-02-15 and 2015-01-20\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         y         ds     yhat yhat_lower yhat_upper     cutoff\n",
      "1 8.242493 2010-02-16 8.956252   8.451796   9.483227 2010-02-15\n",
      "2 8.008033 2010-02-17 8.722669   8.206799   9.174078 2010-02-15\n",
      "3 8.045268 2010-02-18 8.606339   8.136408   9.059811 2010-02-15\n",
      "4 7.928766 2010-02-19 8.528198   8.044159   9.012228 2010-02-15\n",
      "5 7.745003 2010-02-20 8.270095   7.787403   8.780680 2010-02-15\n",
      "6 7.866339 2010-02-21 8.601325   8.051539   9.112439 2010-02-15\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "df.cv <- cross_validation(m, initial = 730, period = 180, horizon = 365, units = 'days')\n",
    "head(df.cv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f4355ae2480f4f7cb37c613c84427237",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=11.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "from prophet.diagnostics import cross_validation\n",
    "df_cv = cross_validation(m, initial='730 days', period='180 days', horizon = '365 days')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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>ds</th>\n",
       "      <th>yhat</th>\n",
       "      <th>yhat_lower</th>\n",
       "      <th>yhat_upper</th>\n",
       "      <th>y</th>\n",
       "      <th>cutoff</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2010-02-16</td>\n",
       "      <td>8.959678</td>\n",
       "      <td>8.470035</td>\n",
       "      <td>9.451618</td>\n",
       "      <td>8.242493</td>\n",
       "      <td>2010-02-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2010-02-17</td>\n",
       "      <td>8.726195</td>\n",
       "      <td>8.236734</td>\n",
       "      <td>9.219616</td>\n",
       "      <td>8.008033</td>\n",
       "      <td>2010-02-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2010-02-18</td>\n",
       "      <td>8.610011</td>\n",
       "      <td>8.104834</td>\n",
       "      <td>9.125484</td>\n",
       "      <td>8.045268</td>\n",
       "      <td>2010-02-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2010-02-19</td>\n",
       "      <td>8.532004</td>\n",
       "      <td>7.985031</td>\n",
       "      <td>9.041575</td>\n",
       "      <td>7.928766</td>\n",
       "      <td>2010-02-15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2010-02-20</td>\n",
       "      <td>8.274090</td>\n",
       "      <td>7.779034</td>\n",
       "      <td>8.745627</td>\n",
       "      <td>7.745003</td>\n",
       "      <td>2010-02-15</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          ds      yhat  yhat_lower  yhat_upper         y     cutoff\n",
       "0 2010-02-16  8.959678    8.470035    9.451618  8.242493 2010-02-15\n",
       "1 2010-02-17  8.726195    8.236734    9.219616  8.008033 2010-02-15\n",
       "2 2010-02-18  8.610011    8.104834    9.125484  8.045268 2010-02-15\n",
       "3 2010-02-19  8.532004    7.985031    9.041575  7.928766 2010-02-15\n",
       "4 2010-02-20  8.274090    7.779034    8.745627  7.745003 2010-02-15"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_cv.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In R, the argument `units` must be a type accepted by `as.difftime`, which is weeks or shorter. In Python, the string for `initial`, `period`, and `horizon` should be in the format used by Pandas Timedelta, which accepts units of days or shorter.\n",
    "\n",
    "Custom cutoffs can also be supplied as a list of dates to the `cutoffs` keyword in the `cross_validation` function in Python and R. For example, three cutoffs six months apart, would need to be passed to the `cutoffs` argument in a date format like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%R\n",
    "cutoffs <- as.Date(c('2013-02-15', '2013-08-15', '2014-02-15'))\n",
    "df.cv2 <- cross_validation(m, cutoffs = cutoffs, horizon = 365, units = 'days')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bd72143956e3426abd28c01b3a0b75b0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=3.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "cutoffs = pd.to_datetime(['2013-02-15', '2013-08-15', '2014-02-15'])\n",
    "df_cv2 = cross_validation(m, cutoffs=cutoffs, horizon='365 days')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `performance_metrics` utility can be used to compute some useful statistics of the prediction performance (`yhat`, `yhat_lower`, and `yhat_upper` compared to `y`), as a function of the distance from the cutoff (how far into the future the prediction was). The statistics computed are mean squared error (MSE), root mean squared error (RMSE), mean absolute error (MAE), mean absolute percent error (MAPE), median absolute percent error (MDAPE) and coverage of the `yhat_lower` and `yhat_upper` estimates. These are computed on a rolling window of the predictions in `df_cv` after sorting by horizon (`ds` minus `cutoff`). By default 10% of the predictions will be included in each window, but this can be changed with the `rolling_window` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  horizon       mse      rmse       mae       mape      mdape      smape\n",
      "1 37 days 0.4945480 0.7032411 0.5055601 0.05858532 0.04934318 0.05887239\n",
      "2 38 days 0.5002964 0.7073163 0.5105297 0.05916131 0.04959087 0.05950813\n",
      "3 39 days 0.5225149 0.7228519 0.5166123 0.05975950 0.04934318 0.06023204\n",
      "4 40 days 0.5297905 0.7278671 0.5194979 0.06006795 0.04934318 0.06061026\n",
      "5 41 days 0.5372450 0.7329700 0.5204945 0.06014704 0.04923865 0.06075311\n",
      "6 42 days 0.5410136 0.7355363 0.5208659 0.06016509 0.04934318 0.06082444\n",
      "   coverage\n",
      "1 0.6815898\n",
      "2 0.6804477\n",
      "3 0.6777067\n",
      "4 0.6838739\n",
      "5 0.6914116\n",
      "6 0.6930105\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "df.p <- performance_metrics(df.cv)\n",
    "head(df.p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>horizon</th>\n",
       "      <th>mse</th>\n",
       "      <th>rmse</th>\n",
       "      <th>mae</th>\n",
       "      <th>mape</th>\n",
       "      <th>mdape</th>\n",
       "      <th>smape</th>\n",
       "      <th>coverage</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>37 days</td>\n",
       "      <td>0.493764</td>\n",
       "      <td>0.702683</td>\n",
       "      <td>0.504754</td>\n",
       "      <td>0.058485</td>\n",
       "      <td>0.049922</td>\n",
       "      <td>0.058774</td>\n",
       "      <td>0.674052</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>38 days</td>\n",
       "      <td>0.499522</td>\n",
       "      <td>0.706769</td>\n",
       "      <td>0.509723</td>\n",
       "      <td>0.059060</td>\n",
       "      <td>0.049389</td>\n",
       "      <td>0.059409</td>\n",
       "      <td>0.672910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>39 days</td>\n",
       "      <td>0.521614</td>\n",
       "      <td>0.722229</td>\n",
       "      <td>0.515793</td>\n",
       "      <td>0.059657</td>\n",
       "      <td>0.049540</td>\n",
       "      <td>0.060131</td>\n",
       "      <td>0.670169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>40 days</td>\n",
       "      <td>0.528760</td>\n",
       "      <td>0.727159</td>\n",
       "      <td>0.518634</td>\n",
       "      <td>0.059961</td>\n",
       "      <td>0.049232</td>\n",
       "      <td>0.060504</td>\n",
       "      <td>0.671311</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>41 days</td>\n",
       "      <td>0.536078</td>\n",
       "      <td>0.732174</td>\n",
       "      <td>0.519585</td>\n",
       "      <td>0.060036</td>\n",
       "      <td>0.049389</td>\n",
       "      <td>0.060641</td>\n",
       "      <td>0.678849</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  horizon       mse      rmse       mae      mape     mdape     smape  \\\n",
       "0 37 days  0.493764  0.702683  0.504754  0.058485  0.049922  0.058774   \n",
       "1 38 days  0.499522  0.706769  0.509723  0.059060  0.049389  0.059409   \n",
       "2 39 days  0.521614  0.722229  0.515793  0.059657  0.049540  0.060131   \n",
       "3 40 days  0.528760  0.727159  0.518634  0.059961  0.049232  0.060504   \n",
       "4 41 days  0.536078  0.732174  0.519585  0.060036  0.049389  0.060641   \n",
       "\n",
       "   coverage  \n",
       "0  0.674052  \n",
       "1  0.672910  \n",
       "2  0.670169  \n",
       "3  0.671311  \n",
       "4  0.678849  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from prophet.diagnostics import performance_metrics\n",
    "df_p = performance_metrics(df_cv)\n",
    "df_p.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cross validation performance metrics can be visualized with `plot_cross_validation_metric`, here shown for MAPE. Dots show the absolute percent error for each prediction in `df_cv`. The blue line shows the MAPE, where the mean is taken over a rolling window of the dots. We see for this forecast that errors around 5% are typical for predictions one month into the future, and that errors increase up to around 11% for predictions that are a year out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 90728 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "plot_cross_validation_metric(df.cv, metric = 'mape')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 78132 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from prophet.plot import plot_cross_validation_metric\n",
    "fig = plot_cross_validation_metric(df_cv, metric='mape')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size of the rolling window in the figure can be changed with the optional argument `rolling_window`, which specifies the proportion of forecasts to use in each rolling window. The default is 0.1, corresponding to 10% of rows from `df_cv` included in each window; increasing this will lead to a smoother average curve in the figure. The `initial` period should be long enough to capture all of the components of the model, in particular seasonalities and extra regressors: at least a year for yearly seasonality, at least a week for weekly seasonality, etc.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parallelizing cross validation\n",
    "\n",
    "Cross-validation can also be run in parallel mode in Python, by setting specifying the `parallel` keyword. Four modes are supported\n",
    "\n",
    "* `parallel=None` (Default, no parallelization)\n",
    "* `parallel=\"processes\"`\n",
    "* `parallel=\"threads\"`\n",
    "* `parallel=\"dask\"`\n",
    "\n",
    "For problems that aren't too big, we recommend using `parallel=\"processes\"`. It will achieve the highest performance when the parallel cross validation can be done on a single machine. For large problems, a [Dask](https://dask.org) cluster can be used to do the cross validation on many machines. You will need to [install Dask](https://docs.dask.org/en/latest/install.html) separately, as it will not be installed with `prophet`.\n",
    "\n",
    "\n",
    "```python\n",
    "from dask.distributed import Client\n",
    "\n",
    "client = Client()  # connect to the cluster\n",
    "df_cv = cross_validation(m, initial='730 days', period='180 days', horizon='365 days',\n",
    "                         parallel=\"dask\")\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hyperparameter tuning\n",
    "\n",
    "Cross-validation can be used for tuning hyperparameters of the model, such as `changepoint_prior_scale` and `seasonality_prior_scale`. A Python example is given below, with a 4x4 grid of those two parameters, with parallelization over cutoffs. Here parameters are evaluated on RMSE averaged over a 30-day horizon, but different performance metrics may be appropriate for different problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    changepoint_prior_scale  seasonality_prior_scale      rmse\n",
      "0                     0.001                     0.01  0.757694\n",
      "1                     0.001                     0.10  0.743399\n",
      "2                     0.001                     1.00  0.753387\n",
      "3                     0.001                    10.00  0.762890\n",
      "4                     0.010                     0.01  0.542315\n",
      "5                     0.010                     0.10  0.535546\n",
      "6                     0.010                     1.00  0.527008\n",
      "7                     0.010                    10.00  0.541544\n",
      "8                     0.100                     0.01  0.524835\n",
      "9                     0.100                     0.10  0.516061\n",
      "10                    0.100                     1.00  0.521406\n",
      "11                    0.100                    10.00  0.518580\n",
      "12                    0.500                     0.01  0.532140\n",
      "13                    0.500                     0.10  0.524668\n",
      "14                    0.500                     1.00  0.521130\n",
      "15                    0.500                    10.00  0.522980\n"
     ]
    }
   ],
   "source": [
    "import itertools\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "param_grid = {  \n",
    "    'changepoint_prior_scale': [0.001, 0.01, 0.1, 0.5],\n",
    "    'seasonality_prior_scale': [0.01, 0.1, 1.0, 10.0],\n",
    "}\n",
    "\n",
    "# Generate all combinations of parameters\n",
    "all_params = [dict(zip(param_grid.keys(), v)) for v in itertools.product(*param_grid.values())]\n",
    "rmses = []  # Store the RMSEs for each params here\n",
    "\n",
    "# Use cross validation to evaluate all parameters\n",
    "for params in all_params:\n",
    "    m = Prophet(**params).fit(df)  # Fit model with given params\n",
    "    df_cv = cross_validation(m, cutoffs=cutoffs, horizon='30 days', parallel=\"processes\")\n",
    "    df_p = performance_metrics(df_cv, rolling_window=1)\n",
    "    rmses.append(df_p['rmse'].values[0])\n",
    "\n",
    "# Find the best parameters\n",
    "tuning_results = pd.DataFrame(all_params)\n",
    "tuning_results['rmse'] = rmses\n",
    "print(tuning_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'changepoint_prior_scale': 0.1, 'seasonality_prior_scale': 0.1}\n"
     ]
    }
   ],
   "source": [
    "best_params = all_params[np.argmin(rmses)]\n",
    "print(best_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, parallelization could be done across parameter combinations by parallelizing the loop above.\n",
    "\n",
    "The Prophet model has a number of input parameters that one might consider tuning. Here are some general recommendations for hyperparameter tuning that may be a good starting place.\n",
    "\n",
    "**Parameters that can be tuned**\n",
    "- `changepoint_prior_scale`: This is probably the most impactful parameter. It determines the flexibility of the trend, and in particular how much the trend changes at the trend changepoints. As described in this documentation, if it is too small, the trend will be underfit and variance that should have been modeled with trend changes will instead end up being handled with the noise term. If it is too large, the trend will overfit and in the most extreme case you can end up with the trend capturing yearly seasonality. The default of 0.05 works for many time series, but this could be tuned; a range of [0.001, 0.5] would likely be about right. Parameters like this (regularization penalties; this is effectively a lasso penalty) are often tuned on a log scale.\n",
    "\n",
    "- `seasonality_prior_scale`: This parameter controls the flexibility of the seasonality. Similarly, a large value allows the seasonality to fit large fluctuations, a small value shrinks the magnitude of the seasonality. The default is 10., which applies basically no regularization. That is because we very rarely see overfitting here (there's inherent regularization with the fact that it is being modeled with a truncated Fourier series, so it's essentially low-pass filtered). A reasonable range for tuning it would probably be [0.01, 10]; when set to 0.01 you should find that the magnitude of seasonality is forced to be very small. This likely also makes sense on a log scale, since it is effectively an L2 penalty like in ridge regression.\n",
    "\n",
    "- `holidays_prior_scale`: This controls flexibility to fit holiday effects. Similar to seasonality_prior_scale, it defaults to 10.0 which applies basically no regularization, since we usually have multiple observations of holidays and can do a good job of estimating their effects. This could also be tuned on a range of [0.01, 10] as with seasonality_prior_scale.\n",
    "\n",
    "- `seasonality_mode`: Options are [`'additive'`, `'multiplicative'`]. Default is `'additive'`, but many business time series will have multiplicative seasonality. This is best identified just from looking at the time series and seeing if the magnitude of seasonal fluctuations grows with the magnitude of the time series (see the documentation here on multiplicative seasonality), but when that isn't possible, it could be tuned.\n",
    "\n",
    "**Maybe tune?**\n",
    "- `changepoint_range`: This is the proportion of the history in which the trend is allowed to change. This defaults to 0.8, 80% of the history, meaning the model will not fit any trend changes in the last 20% of the time series. This is fairly conservative, to avoid overfitting to trend changes at the very end of the time series where there isn't enough runway left to fit it well. With a human in the loop, this is something that can be identified pretty easily visually: one can pretty clearly see if the forecast is doing a bad job in the last 20%. In a fully-automated setting, it may be beneficial to be less conservative. It likely will not be possible to tune this parameter effectively with cross validation over cutoffs as described above. The ability of the model to generalize from a trend change in the last 10% of the time series will be hard to learn from looking at earlier cutoffs that may not have trend changes in the last 10%. So, this parameter is probably better not tuned, except perhaps over a large number of time series. In that setting, [0.8, 0.95] may be a reasonable range.\n",
    "\n",
    "**Parameters that would likely not be tuned**\n",
    "- `growth`: Options are 'linear' and 'logistic'. This likely will not be tuned; if there is a known saturating point and growth towards that point it will be included and the logistic trend will be used, otherwise it will be linear.\n",
    "\n",
    "- `changepoints`: This is for manually specifying the locations of changepoints. None by default, which automatically places them.\n",
    "\n",
    "- `n_changepoints`: This is the number of automatically placed changepoints. The default of 25 should be plenty to capture the trend changes in a typical time series (at least the type that Prophet would work well on anyway). Rather than increasing or decreasing the number of changepoints, it will likely be more effective to focus on increasing or decreasing the flexibility at those trend changes, which is done with `changepoint_prior_scale`.\n",
    "\n",
    "- `yearly_seasonality`: By default ('auto') this will turn yearly seasonality on if there is a year of data, and off otherwise. Options are ['auto', True, False]. If there is more than a year of data, rather than trying to turn this off during HPO, it will likely be more effective to leave it on and turn down seasonal effects by tuning `seasonality_prior_scale`.\n",
    "\n",
    "- `weekly_seasonality`: Same as for `yearly_seasonality`.\n",
    "\n",
    "- `daily_seasonality`: Same as for `yearly_seasonality`.\n",
    "\n",
    "- `holidays`: This is to pass in a dataframe of specified holidays. The holiday effects would be tuned with `holidays_prior_scale`.\n",
    "\n",
    "- `mcmc_samples`: Whether or not MCMC is used will likely be determined by factors like the length of the time series and the importance of parameter uncertainty (these considerations are described in the documentation).\n",
    "\n",
    "- `interval_width`: Prophet `predict` returns uncertainty intervals for each component, like `yhat_lower` and `yhat_upper` for the forecast `yhat`. These are computed as quantiles of the posterior predictive distribution, and `interval_width` specifies which quantiles to use. The default of 0.8 provides an 80% prediction interval. You could change that to 0.95 if you wanted a 95% interval. This will affect only the uncertainty interval, and will not change the forecast `yhat` at all and so does not need to be tuned.\n",
    "\n",
    "- `uncertainty_samples`: The uncertainty intervals are computed as quantiles from the posterior predictive interval, and the posterior predictive interval is estimated with Monte Carlo sampling. This parameter is the number of samples to use (defaults to 1000). The running time for predict will be linear in this number. Making it smaller will increase the variance (Monte Carlo error) of the uncertainty interval, and making it larger will reduce that variance. So, if the uncertainty estimates seem jagged this could be increased to further smooth them out, but it likely will not need to be changed. As with `interval_width`, this parameter only affects the uncertainty intervals and changing it will not affect in any way the forecast `yhat`; it does not need to be tuned."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit Metadata",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
