{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "block_hidden": true
   },
   "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": 2,
   "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]: \n",
      "Attaching package: ‘dplyr’\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:stats’:\n",
      "\n",
      "    filter, lag\n",
      "\n",
      "\n",
      "R[write to console]: The following objects are masked from ‘package:base’:\n",
      "\n",
      "    intersect, setdiff, setequal, union\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "library(prophet)\n",
    "library(dplyr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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 - see the example csv [here](https://github.com/facebook/prophet/blob/main/examples/example_yosemite_temps.csv). When sub-daily data are used, daily seasonality will automatically be fit. Here we fit Prophet to data with 5-minute resolution (daily temperatures at Yosemite):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 210016 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "df <- read.csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_yosemite_temps.csv')\n",
    "m <- prophet(df, changepoint.prior.scale=0.01)\n",
    "future <- make_future_dataframe(m, periods = 300, freq = 60 * 60)\n",
    "fcst <- predict(m, future)\n",
    "plot(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 263452 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_yosemite_temps.csv')\n",
    "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": {},
   "source": [
    "The daily seasonality will show up in the components plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 56284 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 9 -h 9 -u in\n",
    "prophet_plot_components(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 73092 characters]",
      "text/plain": [
       "<Figure size 648x648 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = m.plot_components(fcst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data with regular gaps\n",
    "\n",
    "Suppose the dataset above only had observations from 12a to 6a:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Disabling yearly seasonality. Run prophet with yearly.seasonality=TRUE to override this.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 102404 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "df2 <- df %>%\n",
    "  mutate(ds = as.POSIXct(ds, tz=\"GMT\")) %>%\n",
    "  filter(as.numeric(format(ds, \"%H\")) < 6)\n",
    "m <- prophet(df2)\n",
    "future <- make_future_dataframe(m, periods = 300, freq = 60 * 60)\n",
    "fcst <- predict(m, future)\n",
    "plot(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 120364 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": {},
   "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": 6,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 140476 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "future2 <- future %>% \n",
    "  filter(as.numeric(format(ds, \"%H\")) < 6)\n",
    "fcst <- predict(m, future2)\n",
    "plot(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 171760 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": {},
   "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.\n",
    "\n",
    "## Monthly data\n",
    "\n",
    "You can use Prophet to fit monthly data. However, the underlying model is continuous-time, which means that you can get strange results if you fit the model to monthly data and then ask for daily forecasts. Here we forecast US retail sales volume for the next 10 years:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.\n",
      "\n",
      "R[write to console]: Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 70136 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "df <- read.csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_retail_sales.csv')\n",
    "m <- prophet(df, seasonality.mode = 'multiplicative')\n",
    "future <- make_future_dataframe(m, periods = 3652)\n",
    "fcst <- predict(m, future)\n",
    "plot(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 87300 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_retail_sales.csv')\n",
    "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": {},
   "source": [
    "This is the same issue from above where the dataset has regular gaps. When we 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": 8,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: Disabling weekly seasonality. Run prophet with weekly.seasonality=TRUE to override this.\n",
      "\n",
      "R[write to console]: Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 1).\n",
      "Chain 1: \n",
      "Chain 1: Gradient evaluation took 0.000184 seconds\n",
      "Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.84 seconds.\n",
      "Chain 1: Adjust your expectations accordingly!\n",
      "Chain 1: \n",
      "Chain 1: \n",
      "Chain 1: Iteration:   1 / 300 [  0%]  (Warmup)\n",
      "Chain 1: Iteration:  30 / 300 [ 10%]  (Warmup)\n",
      "Chain 1: Iteration:  60 / 300 [ 20%]  (Warmup)\n",
      "Chain 1: Iteration:  90 / 300 [ 30%]  (Warmup)\n",
      "Chain 1: Iteration: 120 / 300 [ 40%]  (Warmup)\n",
      "Chain 1: Iteration: 150 / 300 [ 50%]  (Warmup)\n",
      "Chain 1: Iteration: 151 / 300 [ 50%]  (Sampling)\n",
      "Chain 1: Iteration: 180 / 300 [ 60%]  (Sampling)\n",
      "Chain 1: Iteration: 210 / 300 [ 70%]  (Sampling)\n",
      "Chain 1: Iteration: 240 / 300 [ 80%]  (Sampling)\n",
      "Chain 1: Iteration: 270 / 300 [ 90%]  (Sampling)\n",
      "Chain 1: Iteration: 300 / 300 [100%]  (Sampling)\n",
      "Chain 1: \n",
      "Chain 1:  Elapsed Time: 11.649 seconds (Warm-up)\n",
      "Chain 1:                17.3646 seconds (Sampling)\n",
      "Chain 1:                29.0136 seconds (Total)\n",
      "Chain 1: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 2).\n",
      "Chain 2: \n",
      "Chain 2: Gradient evaluation took 0.000176 seconds\n",
      "Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.76 seconds.\n",
      "Chain 2: Adjust your expectations accordingly!\n",
      "Chain 2: \n",
      "Chain 2: \n",
      "Chain 2: Iteration:   1 / 300 [  0%]  (Warmup)\n",
      "Chain 2: Iteration:  30 / 300 [ 10%]  (Warmup)\n",
      "Chain 2: Iteration:  60 / 300 [ 20%]  (Warmup)\n",
      "Chain 2: Iteration:  90 / 300 [ 30%]  (Warmup)\n",
      "Chain 2: Iteration: 120 / 300 [ 40%]  (Warmup)\n",
      "Chain 2: Iteration: 150 / 300 [ 50%]  (Warmup)\n",
      "Chain 2: Iteration: 151 / 300 [ 50%]  (Sampling)\n",
      "Chain 2: Iteration: 180 / 300 [ 60%]  (Sampling)\n",
      "Chain 2: Iteration: 210 / 300 [ 70%]  (Sampling)\n",
      "Chain 2: Iteration: 240 / 300 [ 80%]  (Sampling)\n",
      "Chain 2: Iteration: 270 / 300 [ 90%]  (Sampling)\n",
      "Chain 2: Iteration: 300 / 300 [100%]  (Sampling)\n",
      "Chain 2: \n",
      "Chain 2:  Elapsed Time: 10.7148 seconds (Warm-up)\n",
      "Chain 2:                14.0044 seconds (Sampling)\n",
      "Chain 2:                24.7192 seconds (Total)\n",
      "Chain 2: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 3).\n",
      "Chain 3: \n",
      "Chain 3: Gradient evaluation took 0.000114 seconds\n",
      "Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.14 seconds.\n",
      "Chain 3: Adjust your expectations accordingly!\n",
      "Chain 3: \n",
      "Chain 3: \n",
      "Chain 3: Iteration:   1 / 300 [  0%]  (Warmup)\n",
      "Chain 3: Iteration:  30 / 300 [ 10%]  (Warmup)\n",
      "Chain 3: Iteration:  60 / 300 [ 20%]  (Warmup)\n",
      "Chain 3: Iteration:  90 / 300 [ 30%]  (Warmup)\n",
      "Chain 3: Iteration: 120 / 300 [ 40%]  (Warmup)\n",
      "Chain 3: Iteration: 150 / 300 [ 50%]  (Warmup)\n",
      "Chain 3: Iteration: 151 / 300 [ 50%]  (Sampling)\n",
      "Chain 3: Iteration: 180 / 300 [ 60%]  (Sampling)\n",
      "Chain 3: Iteration: 210 / 300 [ 70%]  (Sampling)\n",
      "Chain 3: Iteration: 240 / 300 [ 80%]  (Sampling)\n",
      "Chain 3: Iteration: 270 / 300 [ 90%]  (Sampling)\n",
      "Chain 3: Iteration: 300 / 300 [100%]  (Sampling)\n",
      "Chain 3: \n",
      "Chain 3:  Elapsed Time: 10.1486 seconds (Warm-up)\n",
      "Chain 3:                14.4732 seconds (Sampling)\n",
      "Chain 3:                24.6218 seconds (Total)\n",
      "Chain 3: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 4).\n",
      "Chain 4: \n",
      "Chain 4: Gradient evaluation took 0.00012 seconds\n",
      "Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.2 seconds.\n",
      "Chain 4: Adjust your expectations accordingly!\n",
      "Chain 4: \n",
      "Chain 4: \n",
      "Chain 4: Iteration:   1 / 300 [  0%]  (Warmup)\n",
      "Chain 4: Iteration:  30 / 300 [ 10%]  (Warmup)\n",
      "Chain 4: Iteration:  60 / 300 [ 20%]  (Warmup)\n",
      "Chain 4: Iteration:  90 / 300 [ 30%]  (Warmup)\n",
      "Chain 4: Iteration: 120 / 300 [ 40%]  (Warmup)\n",
      "Chain 4: Iteration: 150 / 300 [ 50%]  (Warmup)\n",
      "Chain 4: Iteration: 151 / 300 [ 50%]  (Sampling)\n",
      "Chain 4: Iteration: 180 / 300 [ 60%]  (Sampling)\n",
      "Chain 4: Iteration: 210 / 300 [ 70%]  (Sampling)\n",
      "Chain 4: Iteration: 240 / 300 [ 80%]  (Sampling)\n",
      "Chain 4: Iteration: 270 / 300 [ 90%]  (Sampling)\n",
      "Chain 4: Iteration: 300 / 300 [100%]  (Sampling)\n",
      "Chain 4: \n",
      "Chain 4:  Elapsed Time: 11.4876 seconds (Warm-up)\n",
      "Chain 4:                16.2248 seconds (Sampling)\n",
      "Chain 4:                27.7124 seconds (Total)\n",
      "Chain 4: \n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 40172 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 9 -h 6 -u in\n",
    "m <- prophet(df, seasonality.mode = 'multiplicative', mcmc.samples = 300)\n",
    "fcst <- predict(m, future)\n",
    "prophet_plot_components(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING:pystan:481 of 600 iterations saturated the maximum tree depth of 10 (80.2 %)\n",
      "WARNING:pystan:Run again with max_treedepth larger than 10 to avoid saturation\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 68468 characters]",
      "text/plain": [
       "<Figure size 648x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "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": "markdown",
   "metadata": {},
   "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": 9,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 60764 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 10 -h 6 -u in\n",
    "future <- make_future_dataframe(m, periods = 120, freq = 'month')\n",
    "fcst <- predict(m, future)\n",
    "plot(m, fcst)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 84392 characters]",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "future = m.make_future_dataframe(periods=120, freq='MS')\n",
    "fcst = m.predict(future)\n",
    "fig = m.plot(fcst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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": "markdown",
   "metadata": {},
   "source": [
    "## Holidays with aggregated data\n",
    "\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."
   ]
  }
 ],
 "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
}
