<jupyter_start><jupyter_text>Computer Lab 4: Forecasting Daily Data using Facebook Prophet**In this practical you will learn:*** How to wrangle time series data into `Prophet` format* How to fit a `Prophet` model to a time series* How to obtain a point forecast and prediction intervals using `Prophet`* How to model special calender events.* How to use `Prophet`'s built in diagnostic and cross validation tools---* Video to introduce the notebook: https://bit.ly/intro_prophet> Note that in this video I use a slightly older version of Prophet.  The only material difference is that prior to version 1.0 (that we are using) you imported from `fbprophet` instead of `prophet`!---  Standard Imports<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys<jupyter_output><empty_output><jupyter_text>FB Prophet ImportsIf you are using the provided conda environment `hds_forecast` you will already have `fbprophet` installed.  > At the time of writing a slightly older version of Prophet is also installed by default in Google Colab.  For this version you will need to use `import fbprophet` instead of `import prophet`<jupyter_code>from prophet import Prophet

import prophet

# should be version 1.0
prophet.__version__<jupyter_output><empty_output><jupyter_text>forecast-tools imports<jupyter_code># if running in Google Colab install forecast-tools
if 'google.colab' in sys.modules:
    !pip install forecast-tools
from forecast_tools.model_selection import auto_naive
from forecast_tools.metrics import mean_absolute_error<jupyter_output><empty_output><jupyter_text>ARIMA importsAt the end of the notebook we will compare Prophet's performance to ARIMA<jupyter_code># if running in Google Colab install pmdarima
if 'google.colab' in sys.modules:
    !pip install pmdarima

from pmdarima import auto_arima, ARIMA
from pmdarima.model_selection import RollingForecastCV, cross_val_score<jupyter_output><empty_output><jupyter_text>Data - ED reattendancesWe will use an time series of patients that reattend an ED within 7 days.  The data are held at the **daily** level.The data is held in the file `ed_reattends_day.csv` Exercise 1: Read in and adjust**Task:*** Load `ed_reattends_day.csv` (url provided below)* Plot the data* Name the `pd.DataFrame` you create `y_train`**Hints*** The data is stored in UK day first format.* Remember to set the freq of the DataTimeIndex**Questions**:* How would you describe the time series?<jupyter_code># your code here ...
url = 'https://raw.githubusercontent.com/health-data-science-OR/' \
      + 'hpdm097-datasets/master/ed_reattends_day.csv'<jupyter_output><empty_output><jupyter_text>Exercise 2: Wrangle the data into Prophet FormatA `Prophet` model requires the data in a `pd.DataFrame` with two columns: 'ds' - the datetimestamp and `y` the observations used in training (the index is an arbitrary numeric value).  But your current (and standard form) dataframe has a datetimeindex **Task:*** Code a function that converts the training data into the correct format for Prophet.* The end product should be a dataframe with 2 columns (in the following order) * 'ds': the date * 'y': the number of reattends on that day.* Convert your data<jupyter_code># your code goes here...<jupyter_output><empty_output><jupyter_text>Exercise 3: Fit and predict using a basic Prophet modelFitting a basic Prophet model is relatively straightforward.  We need to create a `Prophet` object```pythonmodel = Prophet()```Unlike ARIMA where we specified the prediction interval width at prediction time, with Prophet we need to specify the interval width as we create the model.  We do this using the parameter `interval_width`.  By default this produces a 80\% prediction interval.  Note we specify a $1-\alpha$ interval width (unlike in ARIMA where we specified $\alpha$)```pythonmodel = Prophet(interval_width=0.95)```and call `.fit()` passing in the training data.```pythonmodel.fit(y_train)``` Exercise 3.1 Fitting**Task:*** Fit a basic `Prophet` model to the training data* Set the model up to produce 95\% prediction intervals<jupyter_code># your code here ...<jupyter_output><empty_output><jupyter_text>Exercise 3.2 ForecastingThere are two steps to makeing a forecast with a fitted `Prophet` model.  You must first create a future dataframe.  This is a `pd.DataFrame` with a single column 'ds' that runs from the start of the training period until the end of the training period + a horizon.  Prophet makes this easy by including a `make_future_dataframe()` method.```pythonfuture = model.make_future_dataframe(periods=28)```Once we have the future dataframe we can pass that to the predict method.```pythonprophet_forecast = model.predict(future)```This returns a `pd.DataFrame` that contains both the in-sample fitted values and the out of sample forecasts. It contains all of the components of the fitted model. For example the last 5 rows of a forecast are:|      | ds                  |   trend |   yhat_lower |   yhat_upper |   trend_lower |   trend_upper |   additive_terms |   additive_terms_lower |   additive_terms_upper |    weekly |   weekly_lower |   weekly_upper |   yearly |   yearly_lower |   yearly_upper |   multiplicative_terms |   multiplicative_terms_lower |   multiplicative_terms_upper |    yhat ||-----:|:--------------------|--------:|-------------:|-------------:|--------------:|--------------:|-----------------:|-----------------------:|-----------------------:|----------:|---------------:|---------------:|---------:|---------------:|---------------:|-----------------------:|-----------------------------:|-----------------------------:|--------:|| 1629 | 2018-09-18 00:00:00 | 187.353 |      163.694 |      236.271 |       186.22  |       188.383 |        11.0807   |              11.0807   |              11.0807   |  -2.83494 |       -2.83494 |       -2.83494 |  13.9156 |        13.9156 |        13.9156 |                      0 |                            0 |                            0 | 198.433 || 1630 | 2018-09-19 00:00:00 | 187.305 |      157.7   |      226.979 |       186.153 |       188.356 |         5.91039  |               5.91039  |               5.91039  |  -7.49286 |       -7.49286 |       -7.49286 |  13.4032 |        13.4032 |        13.4032 |                      0 |                            0 |                            0 | 193.215 || 1631 | 2018-09-20 00:00:00 | 187.257 |      150.343 |      224.271 |       186.088 |       188.328 |         0.527382 |               0.527382 |               0.527382 | -12.3179  |      -12.3179  |      -12.3179  |  12.8452 |        12.8452 |        12.8452 |                      0 |                            0 |                            0 | 187.785 || 1632 | 2018-09-21 00:00:00 | 187.209 |      156.873 |      226.588 |       186.023 |       188.3   |         3.90459  |               3.90459  |               3.90459  |  -8.34281 |       -8.34281 |       -8.34281 |  12.2474 |        12.2474 |        12.2474 |                      0 |                            0 |                            0 | 191.114 || 1633 | 2018-09-22 00:00:00 | 187.162 |      169.299 |      241.272 |       185.954 |       188.272 |        15.9777   |              15.9777   |              15.9777   |   4.36149 |        4.36149 |        4.36149 |  11.6162 |        11.6162 |        11.6162 |                      0 |                            0 |                            0 | 203.139 |We can the plot the forecast using Prophet's built in functionality```pythonmodel.plot(prophet_forecast);```In the plot you will see* Black dots - these represent the actual observations (training data)* Blue line - this represents the point forecast* Light blue shaded region - this is the 95\% prediction interval**Task:*** Use your Prophet model to make a 84 day prediction.* Plot the prediction<jupyter_code># your code here ...<jupyter_output><empty_output><jupyter_text>Exercise 4: Plot the model componentsTo decompose the Prophet model and see the trend and seasonal components you can call:```pythonmodel.plot_components(forecast)```**Task:*** Plot the Prophet model's trend and seasonal components.**Questions:*** What observations can you make about the weekly (i.e. days of the week) and yearly (i.e months of the year) seasonality?<jupyter_code># your code here...<jupyter_output><empty_output><jupyter_text>Exercise 5: Adding standard holidaysAs this is ED reattendance data at the daily level, it is likely that we will be seeing some calender day/holiday effects.  Prophet has been designed to deal with 'holidays'.  Effectively a series of binary variables are added for each holiday.Prophet has a number of standard holidays built in by country.  To add them in for 'England' you would use the following code:```pythonmodel = Prophet(interval_width=0.95)model.add_country_holidays(country_name='England')model.fit(y_train)```Note that you call the `add_country_holidays` **before** you fit the model.  This is because it is command telling Prophet what you would like to fit.To see what models have been fitted you can call:```pythonmodel.train_holiday_names.to_list()```When you plot the components of the model a new panel will be added for holidays.  It can sometimes be a bit difficult to which holidays are having which effect so you can plot individual holidays as follows:```pythonfrom prophet.plot import plot_forecast_componentplot_forecast_component(model, prophet_forecast, 'Christmas Day');```**Task:*** Modify your code to include a holidays for 'England'* Make an 84 day forecast* Plot the model components. * Explore the holidays**Questions:*** What do you notice about Christmas day, Easter Monday and Good Friday?  Are there any other days that appear important?**Hints:*** Take a look at the new forecast dataframe.  This includes new columns with the effect of each `holidays`.<jupyter_code># your code here ...<jupyter_output><empty_output><jupyter_text>Exercise 6: Cross ValidationFBProphet provides its own cross validation procedures in the `diagnostics` module.  Two key functions are the `cross_validation` and `performance_metrics` functions.```pythonfrom prophet.diagnostics import cross_validation, performance_metrics```Prophet uses a rolling origin forecast type procedure for cross validation (but rather than working forwards Prophet works backwards from the final fold).  The `cross_validation` function requires some parameters to be **time delta** compatable.  This means that we need to pass in a `str` such as `"7 days"`To use the `cross_validation` function with your model to do a CV where the initial length of the training set is 1095 days (3 years), the step size (period) is 28 days and the forecast horizon is 84 days we would use the following code:```pythondf_cv = cross_validation(model=model, initial='1095 days', period='28 days',                          horizon='84 days')```The variable `df_cv` is of type `pandas.DataFrame`.  It contains all of the predictions and actual values by date and cut off.  In our example above the first row will be day 366; the next 83 rows will be all part of the same cross-validation fold.  See below for example, note that all the select rows (0, 1 and 83) have the same **cut-off** - this the **forecast origin** of the cross validation fold.```pythondf_cv.iloc[[0, 1, 83]```|    | ds                  |    yhat |   yhat_lower |   yhat_upper |   y | cutoff              ||---:|:--------------------|--------:|-------------:|-------------:|----:|:--------------------||  0 | 2017-04-09 00:00:00 | 215.292 |      180.674 |      249.161 | 241 | 2017-04-08 00:00:00 ||  1 | 2017-04-10 00:00:00 | 222.806 |      187.367 |      255.534 | 259 | 2017-04-08 00:00:00 || 83 | 2017-07-01 00:00:00 | 234.83  |      200.091 |      270.619 | 221 | 2017-04-08 00:00:00 |If we then look at rows 84 and 85 we see that the cutoff (forecast origin) has rolled forward.  It is equal to the previous cut-off + 28 days (period in prophet terminology).```pythondf_cv.iloc[[84, 85]```|    | ds                  |    yhat |   yhat_lower |   yhat_upper |   y | cutoff              ||---:|:--------------------|--------:|-------------:|-------------:|----:|:--------------------|| 84 | 2017-05-07 00:00:00 | 226.971 |      192.175 |      260.083 | 217 | 2017-05-06 00:00:00 || 85 | 2017-05-08 00:00:00 | 234.145 |      198.177 |      269.963 | 213 | 2017-05-06 00:00:00 |So in essence `df_cv` is a record of all of the cross-validation folds stacked on top of each other.  This is useful, because it means have all of the data to calculate any forecast error measure we would like for over any sub-forecast horizon.  Prophet provides the function `performance_metrics` and the plotting function `fbprophet.plot.plot_cross_validation_metric` to automate this for you.  For example:**Note: mdape is median absolute percentage error**```pythondf_p = performance_metrics(df_cv)df_p.head()```|    | horizon          |     mse |    rmse |     mae |      mape |     mdape |   coverage ||---:|:-----------------|--------:|--------:|--------:|----------:|----------:|-----------:||  0 | 9 days 00:00:00  | 387.234 | 19.6783 | 15.6597 | 0.0798998 | 0.0632103 |   0.91453  ||  1 | 10 days 00:00:00 | 355.57  | 18.8566 | 14.9167 | 0.0775167 | 0.0632103 |   0.931013 ||  2 | 11 days 00:00:00 | 317.469 | 17.8176 | 13.8917 | 0.0726941 | 0.0535132 |   0.945665 ||  3 | 12 days 00:00:00 | 311.386 | 17.6461 | 13.8011 | 0.0719934 | 0.0535132 |   0.954212 ||  4 | 13 days 00:00:00 | 321.142 | 17.9204 | 13.908  | 0.0719082 | 0.0575246 |   0.948718 |---**Task:*** Run Prophet's built in cross validation* Use an initial training period of 1095 days, period=28 and horizon of '84 days'* Use Prophet's built in `cross_validation` and `performance_metrics` functions<jupyter_code># your code here...<jupyter_output><empty_output><jupyter_text>Exercise 7: Is Prophet better than ARIMA?**Task:*** Choose and fit an ARIMA model to the dataset.* Use cross-validation of the 84 day MAE compare the performance of ARIMA and Prophet.**Questions*** Which of the two methods would you choose in practice? * What else might you do?**Hints**:* See Introduction to ARIMA notebook for a reminder about how to use ARIMA models in Python* This is daily level data.  For ARIMA this is `m=7`* Remember that Prophet and ARIMA use data in different formats. * It might take a few minutes for `auto_arima` to return results.  Once it has write down the order of the model so that you can easily recreate it with a `ARIMA` object.* ARIMA fitting is not always stable for daily level data and procedures may output some warnings.  Types of warning you may encounter are:**Convergence warning**: Where the the Maximum likelihood optimisation has failed to converge. If you get a convergence warning when fitting `ARIMA` to your chosen specification then try passing increasing the parameter `maxiter` (default=50) to 100 or more in the `auto_arima` function and `ARIMA` constructor method.**Non-invertible starting MA parameters warning**.  Details of what is meant by an invertible MA process can be found [here](https://otexts.com/fpp2/MA.html).  For our practical purposes the chances are if you see this then when using `auto_arima` and it is unlikely that this model will be selected by the procedure. There is a good chance the warning will also be encountered in cross validation when the model is using less data.  Its worth knowing that `statsmodels` will still fit a model in these circumstances.* If you wish you can surpress the warnings in `auto_arima` by setting `supress_warnings=True`* It will also take a few minutes for ARIMA cross validation to run.<jupyter_code># your code here<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12294664459442252, 'fraction_numerical': 0.07322042531516618, 'mean_word_length': 2.8403422982885087, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 51, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '17882', 'license': 'MIT', 'n_tokens_mistral': 5223, 'n_tokens_neox': 4339, 'n_words': 1799, 'path': '03_prophet/03_practicals_notebooks/prophet_exercises.ipynb', 'repo_name': 'health-data-science-OR/forecasting', 'size': '15706'}