{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "block_hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The rpy2.ipython extension is already loaded. To reload it, use:\n",
      "  %reload_ext rpy2.ipython\n"
     ]
    }
   ],
   "source": [
    "%load_ext rpy2.ipython\n",
    "%matplotlib inline\n",
    "\n",
    "from prophet import Prophet\n",
    "from matplotlib import pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import logging\n",
    "import warnings\n",
    "\n",
    "logging.getLogger('prophet').setLevel(logging.ERROR)\n",
    "logging.getLogger('cmdstanpy').setLevel(logging.ERROR)\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "block_hidden": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv')\n",
    "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": 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]: Disabling yearly seasonality. Run prophet with yearly.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"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "library(prophet)\n",
    "df <- read.csv('https://raw.githubusercontent.com/facebook/prophet/main/examples/example_wp_log_peyton_manning.csv')\n",
    "df <- df[1:180,]\n",
    "m <- prophet(df)\n",
    "future <- make_future_dataframe(m, periods=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default Prophet will return uncertainty intervals for the forecast `yhat`. There are several important assumptions behind these uncertainty intervals.\n",
    "\n",
    "There are three sources of uncertainty in the forecast: uncertainty in the trend, uncertainty in the seasonality estimates, and additional observation noise.\n",
    "\n",
    "### Uncertainty in the trend\n",
    "The biggest source of uncertainty in the forecast is the potential for future trend changes. The time series we have seen already in this documentation show clear trend changes in the history. Prophet is able to detect and fit these, but what trend changes should we expect moving forward? It's impossible to know for sure, so we do the most reasonable thing we can, and we assume that the *future will see similar trend changes as the history*. In particular, we assume that the average frequency and magnitude of trend changes in the future will be the same as that which we observe in the history. We project these trend changes forward and by computing their distribution we obtain uncertainty intervals.\n",
    "\n",
    "One property of this way of measuring uncertainty is that allowing higher flexibility in the rate, by increasing `changepoint_prior_scale`, will increase the forecast uncertainty. This is because if we model more rate changes in the history then we will expect more in the future, and makes the uncertainty intervals a useful indicator of overfitting.\n",
    "\n",
    "The width of the uncertainty intervals (by default 80%) can be set using the parameter `interval_width`:"
   ]
  },
  {
   "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",
      "R[write to console]: Disabling daily seasonality. Run prophet with daily.seasonality=TRUE to override this.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "m <- prophet(df, interval.width = 0.95)\n",
    "forecast <- predict(m, future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "forecast = Prophet(interval_width=0.95).fit(df).predict(future)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, these intervals assume that the future will see the same frequency and magnitude of rate changes as the past. This assumption is probably not true, so you should not expect to get accurate coverage on these uncertainty intervals.\n",
    "\n",
    "### Uncertainty in seasonality\n",
    "By default Prophet will only return uncertainty in the trend and observation noise. To get uncertainty in seasonality, you must do full Bayesian sampling. This is done using the parameter `mcmc.samples` (which defaults to 0). We do this here for the first six months of the Peyton Manning data from the Quickstart:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "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",
      "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 8.6e-05 seconds\n",
      "Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.86 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: 2.10541 seconds (Warm-up)\n",
      "Chain 1:                2.37589 seconds (Sampling)\n",
      "Chain 1:                4.4813 seconds (Total)\n",
      "Chain 1: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 2).\n",
      "Chain 2: \n",
      "Chain 2: Gradient evaluation took 7.4e-05 seconds\n",
      "Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.74 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: 2.04288 seconds (Warm-up)\n",
      "Chain 2:                2.12795 seconds (Sampling)\n",
      "Chain 2:                4.17083 seconds (Total)\n",
      "Chain 2: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 3).\n",
      "Chain 3: \n",
      "Chain 3: Gradient evaluation took 7.1e-05 seconds\n",
      "Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.71 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: 2.24488 seconds (Warm-up)\n",
      "Chain 3:                2.5196 seconds (Sampling)\n",
      "Chain 3:                4.76448 seconds (Total)\n",
      "Chain 3: \n",
      "\n",
      "SAMPLING FOR MODEL 'prophet' NOW (CHAIN 4).\n",
      "Chain 4: \n",
      "Chain 4: Gradient evaluation took 6.5e-05 seconds\n",
      "Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.65 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: 2.2803 seconds (Warm-up)\n",
      "Chain 4:                2.73488 seconds (Sampling)\n",
      "Chain 4:                5.01518 seconds (Total)\n",
      "Chain 4: \n"
     ]
    }
   ],
   "source": [
    "%%R\n",
    "m <- prophet(df, mcmc.samples = 300)\n",
    "forecast <- predict(m, future)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Prophet(mcmc_samples=300)\n",
    "forecast = m.fit(df, show_progress=False).predict(future)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This replaces the typical MAP estimation with MCMC sampling, and can take much longer depending on how many observations there are - expect several minutes instead of several seconds. If you do full sampling, then you will see the uncertainty in seasonal components when you plot them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "output_hidden": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 32644 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%R -w 9 -h 6 -u in\n",
    "prophet_plot_components(m, forecast)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 60112 characters]\n",
      "text/plain": [
       "<Figure size 648x432 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = m.plot_components(forecast)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can access the raw posterior predictive samples in Python using the method `m.predictive_samples(future)`, or in R using the function `predictive_samples(m, future)`."
   ]
  }
 ],
 "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
}
