Forecasting Weekly Data with Prophet

Prepare Notebook

import numpy as np
import pandas as pd
from fbprophet import Prophet

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('darkgrid', {'axes.facecolor': '.9'})
sns.set_palette(palette='deep')
sns_c = sns.color_palette(palette='deep')
%matplotlib inline

from pandas.plotting import register_matplotlib_converters
register_matplotlib_converters()

Warning: At the time of writing this post I encounter a problem when importing fbprophet, see the issue on GitHub. The (partial) solution is given on the thread. The problem as in the structure change of the holidays package.

Brief Model Description

The Prophet model has the form \(y(t) = g(t) + s(t) + h(t) + \varepsilon_t\), where:

• \(g(t)\) is the trend function.
• \(s(t)\) is the periodic component (seasonalities)
• \(h(t)\) represents holidays/events which occur on potentially irregular basis.
• \(\varepsilon_t\) is the error term (which is often assumed to be normally distributed)

Let us describe these components in more detail:

Trend:

The basic model for this term is

\[ g(t) = \frac{C}{1 + e^{-k(t-m)}} \]

where \(C\) is known as the capacity, \(k\) the growth rate and \(m\) the offset parameter. Observe that

\[ \lim_{t\rightarrow \infty} g(t) = C \]

Let us plot \(g(t)\) for various values of \(k\).

plt.rcParams['figure.figsize'] = [12, 7]

def g(t, k, C=1, m=0):
    """Trend model."""
    return C/(1 + np.exp(-k*(t - m)))

# grid of k values. 
k_grid = [0.5, 1, 3, 5, 10, 30, 80]
# Time range. 
t = np.linspace(start=0, stop=1, num=50)

fig, ax = plt.subplots()

for k in k_grid:
    
    sns.lineplot(
        x=t, 
        y=np.apply_along_axis(lambda t: g(t, k, C=2, m=0), axis=0, arr=t),
        label=f'k = {k}',
        ax=ax
    )

ax.axvline(x=0.0, color='gray', linestyle='--', label='m = 0')
ax.axhline(y=2.0, color='black', linestyle='--', label='C = 2')
ax.legend(bbox_to_anchor=(1.02, 1.0))
ax.set(title='Trend Model', ylim=(0.9, 2.1));

Hence, \(C\) can be understood as the saturation point.

In the actual implementation there are some extensions:

• The capacity is a function of time \(C=C(t)\).
• The growth rate is not constant. There is a change-point-grid on which the growth rate is allowed to change. To avoid overfitting, the vector of rate growth adjustments \(\bf{\delta}\in \mathbb{R}^{s}\) has a Laplace prior (related to an \(L^1\) regularization). For the details please refer to the original paper.²

Seasonality

The seasonal components are approximated by Fourier modes:

\[ s(t) = \sum_{n=1}^{N} \left( a_n\cos\left(\frac{2\pi nt}{P}\right) + b_n\sin\left(\frac{2\pi nt}{P}\right) \right) \]

where \(P\) is the period.

Let us plot some of these Fourier modes:

P = 1.0
n_grid = [1, 2 , 5, 20]
t = np.linspace(start=0, stop=P, num=100)

fig, ax = plt.subplots(2, 1, figsize=(12, 9))

for n in n_grid:
    
    sns.lineplot(
        x=t, 
        y=np.apply_along_axis(lambda t: np.cos(2*np.pi*n*t/P), axis=0, arr=t), 
        label=f'n = {n}', 
        ax=ax[0]
    )
    
    sns.lineplot(
        x=t, 
        y=np.apply_along_axis(lambda t: np.sin(2*np.pi*n*t/P), axis=0, arr=t), 
        label=f'n = {n}', 
        ax=ax[1]
    )
    
ax[0].legend(bbox_to_anchor=(1.12, 1.01))
ax[0].set(title=r'$\cos(2\pi n t/P)$')

ax[1].legend(bbox_to_anchor=(1.12, 1.01))
ax[1].set(title=r'$\sin(2\pi n t/P)$')

fig.suptitle('Fourier Modes');

Generate Data

We generate a time series by including the following components:

• We include a non-linear trend of the form \(\text{trend}(w) = (w + 1)^{2/5} + \log(w + 3)\), where \(w\) denotes the week index.

• Seasonality:

  • We use the formula \(\sin(2\pi \times \text{day_of_month}/\text{daysinmonth})\) to generate a monthly seasonality.
    
  • We use the formula \(\sin(2\pi \times\text{month}/3) + \cos(2\pi \times \text{month}/4)\) to generate the yearly seasonality. Note that as \(3\) and \(4\) are relative primes, which implies that the period is the least common multiple \(\text{lcm}(3, 4)=12\).
  • We model the End of the year season as a bump function.

• Gaussian Noise

np.random.seed(seed=42)

def generate_time_series_df(start_date, end_date, freq):
    """Generate time series sample data."""
    date_range = pd.date_range(start=start_date, end=end_date, freq=freq)

    df = pd.DataFrame(data={'ds': date_range})
    
    # Get date variables. 
    df['day_of_month'] = df['ds'].dt.day
    df['month'] = df['ds'].dt.month
    df['daysinmonth'] = df['ds'].dt.daysinmonth
    df['week'] = df['ds'].dt.week
    
    # Time Series Components 
    ## Trend
    df['trend'] = np.power(df.index.values + 1, 2/5) + np.log(df.index.values + 3)
    ## Seasonal
    df['monthly_seas'] = np.cos(2*np.pi*df['day_of_month']/df['daysinmonth'])
    df['yearly_seas'] = 1.2*(np.sin(np.pi*df['month']/3) + np.cos(2*np.pi*df['month']/4))
    df['end_of_year']= - 8.5*np.exp(- ((df['week'] - 51.5)/1.0)**2) \
    ## Gaussian noise
    df['noise'] = np.random.normal(loc=0.0, scale=0.3, size=df.shape[0])
    
    # Target variable.
    df['y'] = df['trend'] \
        + df['monthly_seas'] \
        + df['yearly_seas'] \
        + df['end_of_year'] \
        + df['noise']
    
    return df

df = generate_time_series_df(
    start_date='2016-06-30', 
    end_date='2020-10-31', 
    freq='W'
)

df.head()

Let us plot the resulting time series:

fig, ax = plt.subplots()
sns.lineplot(x='ds', y='y', label='y', data=df, ax=ax)
ax.legend(loc='upper left')
ax.set(title='Raw Data', xlabel='date', ylabel='');

Let us now plot the individual components:

fig, ax = plt.subplots()
sns.lineplot(x='ds', y='y', label='y', data=df, ax=ax)
sns.lineplot(x='ds', y='trend', label='trend', data=df, alpha=0.7, ax=ax)
sns.lineplot(x='ds', y='monthly_seas', label='monthly_seas', data=df, alpha=0.7, ax=ax)
sns.lineplot(x='ds', y='yearly_seas', label='yearly_seas', data=df, alpha=0.7, ax=ax)
sns.lineplot(x='ds', y='end_of_year', label='end_of_year', data=df, alpha=0.7, ax=ax)
sns.lineplot(x='ds', y='noise', label='noise', data=df, alpha=0.7, ax=ax)
ax.legend(loc='upper left')
ax.set(title='Raw Data -  Components', xlabel='date', ylabel='');

Training - Test Split

Let us split the data into a training and test set in order to evaluate our model.

# Define threshold date.
threshold_date = pd.to_datetime('2019-11-01')
mask = df['ds'] < threshold_date
# Split the data and select `ds` and `y` columns.
df_train = df[mask][['ds', 'y']]
df_test = df[~ mask][['ds', 'y']]

Warning: The input to Prophet is always a data frame with two columns: ds and y.³

fig, ax = plt.subplots()
sns.lineplot(x='ds', y='y', label='y_train', data=df_train, ax=ax)
sns.lineplot(x='ds', y='y', label='y_test', data=df_test, ax=ax)
ax.axvline(threshold_date, color=sns_c[3], linestyle='--', label='train test split')
ax.legend(loc='upper left')
ax.set(title='Dependent Variable', ylabel='');

Time Series Decomposition

We begin the analysis by decomposing the (training) time series into the trend, seasonal and residual components.

from statsmodels.tsa.seasonal import seasonal_decompose

decomposition_obj = seasonal_decompose(
    x=df_train.set_index('ds'), 
    model='additive'
)

We plot and compare the results.

fig, ax = plt.subplots(4, 1, figsize=(12, 12))

# Observed time series.
decomposition_obj.observed.plot(ax=ax[0])
ax[0].set(title='observed')
# Trend component. 
decomposition_obj.trend.plot(label='fit', ax=ax[1])
df[mask][['ds', 'trend']].set_index('ds').plot(c=sns_c[1], ax=ax[1])
ax[1].legend(loc='lower right')
ax[1].set(title='trend')
# Seasonal component. 
decomposition_obj.seasonal.plot(label='fit', ax=ax[2])
df.assign(seasonal = lambda x: x['yearly_seas'] + x['monthly_seas'] + x['end_of_year']) \
    [mask][['ds', 'seasonal']] \
    .set_index('ds')\
    .plot(c=sns_c[2], ax=ax[2])
ax[2].legend(loc='lower right')
ax[2].set(title='seasonal')
# Residual.
decomposition_obj.resid.plot(label='fit', ax=ax[3])
df[mask][['ds', 'noise']].set_index('ds').plot(c=sns_c[3], ax=ax[3])
ax[3].legend(loc='lower right')
ax[3].set(title='residual')

fig.suptitle('Time Series Decomposition', y=1.01)
plt.tight_layout()

Remarks:

• The trend component is overall correct. Nevertheless, it fails to capture it for the initial (non-linear) period.
• From the seasonal component obtained we clearly see the yearly and monthly contributions.

Define Model

Based on the analysis above, we now define the forecasting model structure.

# End of the year season.
mask_eoy = (df['month']==12) & (df['day_of_month'] > 21)

df[mask_eoy][['ds', 'end_of_year', 'y']]

def create_end_of_year_holydays_df():
    """Create holidays data frame for the end of the year season."""
    holidays = pd.DataFrame({
      'holiday': 'end_of_year',
      'ds': pd.to_datetime(
          ['2016-12-25', '2017-12-24', '2018-12-23', '2019-12-22']
      ),
      'lower_window': -7,
      'upper_window': 7,
    })
    return holidays

Build Model

Now we define the forecasting model object.

def build_model():
    """Define forecasting model."""
    # Create holidays data frame. 
    holidays = create_end_of_year_holydays_df()
    
    model = Prophet(
        yearly_seasonality=True,
        weekly_seasonality=False,
        daily_seasonality=False, 
        holidays = holidays, 
        interval_width=0.95, 
        mcmc_samples = 500
    )

    model.add_seasonality(
        name='monthly', 
        period=30.5, 
        fourier_order=5
    )
    
    return model
    
model = build_model()

Remarks:

• We specify to fit the yearly_seasonality with the auto option for the Fourier modes.
• We ask for the 0.95 interval_with instead of the default (0.8).
• We include the mcmc_samples option to get uncertainty in seasonality (via Bayesian sampling).
• We add monthly seasonality by specifying the period and fourier_order. This is the general strategy for adding any type of seasonality.⁵

# We train the model with the training data. 
model.fit(df_train)

Generate Predictions

Let us get the model predictions. First we extend the dates from the training data.

# Extend dates and features. 
future = model.make_future_dataframe(periods=df_test.shape[0], freq='W')
# Generate predictions. 
forecast = model.predict(df=future)

Let us see the columns of the forecast data frame.

for c in forecast.columns.sort_values():
    print(c)

    additive_terms
    additive_terms_lower
    additive_terms_upper
    ds
    end_of_year
    end_of_year_lower
    end_of_year_upper
    holidays
    holidays_lower
    holidays_upper
    monthly
    monthly_lower
    monthly_upper
    multiplicative_terms
    multiplicative_terms_lower
    multiplicative_terms_upper
    trend
    trend_lower
    trend_upper
    yearly
    yearly_lower
    yearly_upper
    yhat
    yhat_lower
    yhat_upper

• The variable yhat represents the model predictions.
• For all the components we have the _lower and _upper bounds.

Let us verify the encoding of the end_of_year season in the forecast data frame.

forecast[forecast['end_of_year'].abs()>0][['ds', 'end_of_year']]

fig, ax = plt.subplots()
sns.lineplot(x='ds', y='end_of_year', data=forecast, ax=ax)
ax.axvline(threshold_date, color=sns_c[3], linestyle='--', label='train test split')
ax.set(title='end_of_year', ylabel='');

The seasonality works as expected (see Warning above).

Let us split the predictions into training and test set.

mask2 = forecast['ds'] < threshold_date

forecast_train = forecast[mask2]
forecast_test = forecast[~ mask2]

fig, ax = plt.subplots()

ax.fill_between(
    x=forecast['ds'],
    y1=forecast['yhat_lower'],
    y2=forecast['yhat_upper'],
    color=sns_c[2], 
    alpha=0.25,
    label=r'0.95 credible_interval'
)

sns.lineplot(x='ds', y='y', label='y_train', data=df_train, ax=ax)
sns.lineplot(x='ds', y='y', label='y_test', data=df_test, ax=ax)
sns.lineplot(x='ds', y='yhat', label='y_hat', data=forecast, ax=ax)
ax.axvline(threshold_date, color=sns_c[3], linestyle='--', label='train test split')
ax.legend(loc='upper left')
ax.set(title='Dependent Variable', ylabel='');

Zooming in:

fig, ax = plt.subplots()

ax.fill_between(
    x=forecast_test['ds'],
    y1=forecast_test['yhat_lower'],
    y2=forecast_test['yhat_upper'],
    color=sns_c[2], 
    alpha=0.25,
    label=r'0.95 credible_interval'
)

sns.lineplot(x='ds', y='y', label='y_test', data=df_test, ax=ax)
sns.lineplot(x='ds', y='yhat', label='y_hat', data=forecast_test, ax=ax)
ax.legend(loc='lower right')
ax.set(title='Dependent Variable', ylabel='');

fig, ax = plt.subplots(figsize=(8,8))

# Generate diagonal line to plot. 
d_x = np.linspace(start=df_train['y'].min() - 1, stop=df_train['y'].max() + 1, num=100)

sns.regplot(x=df_train['y'], y=forecast_train['yhat'], color=sns_c[0], label='train', ax=ax)
sns.regplot(x=df_test['y'], y=forecast_test['yhat'], color=sns_c[1], label='test', ax=ax)
sns.lineplot(x=d_x, y=d_x, dashes={'linestyle': ''}, color=sns_c[3], ax=ax)
ax.lines[2].set_linestyle('--')
ax.set(title='Test Data vs Predictions');

Let us compute the r2_score and mean_absolute_error on the training and test set respectively:

from sklearn.metrics import r2_score, mean_absolute_error

print('r2 train: {}'.format(r2_score(y_true=df_train['y'], y_pred=forecast_train['yhat'])))
print('r2 test: {}'.format(r2_score(y_true=df_test['y'], y_pred=forecast_test['yhat'])))
print('---'*10)
print('mae train: {}'.format(mean_absolute_error(y_true=df_train['y'], y_pred=forecast_train['yhat'])))
print('mae test: {}'.format(mean_absolute_error(y_true=df_test['y'], y_pred=forecast_test['yhat'])))

    r2 train: 0.9830369776415686
    r2 test: 0.9320449074303058
    ------------------------------
    mae train: 0.3113940440897326
    mae test: 0.38975341010250913

This might indicate a potential overfit. In a second iteration one could modify the prior_scale in the model definition to add more regularization.⁶

Error Analysis

Let us study the forecast errors.

• Distribution

forecast_test.loc[:, 'errors'] = forecast_test.loc[:, 'yhat'] - df_test.loc[:, 'y']

errors_mean = forecast_test['errors'].mean()
errors_std = forecast_test['errors'].std()

fig, ax = plt.subplots()

sns.distplot(a=forecast_test['errors'], ax=ax, bins=15, rug=True)
ax.axvline(x=errors_mean, color=sns_c[2], linestyle='--', label=r'$\mu$')
ax.axvline(x=errors_mean + 2*errors_std, color=sns_c[3], linestyle='--', label=r'$\mu \pm 2\sigma$')
ax.axvline(x=errors_mean - 2*errors_std, color=sns_c[3], linestyle='--')
ax.legend()
ax.set(title='Model Errors (Test Set)');

• Autocorrelation

fig, ax = plt.subplots()

sns.scatterplot(x='index', y='errors', data=forecast_test.reset_index(), ax=ax)
ax.axhline(y=errors_mean, color=sns_c[2], linestyle='--', label=r'$\mu$ ')
ax.axhline(y=errors_mean + 2*errors_std, color=sns_c[3], linestyle='--', label=r'$\mu \pm 2\sigma$')
ax.axhline(y=errors_mean - 2*errors_std, color=sns_c[3], linestyle='--')
ax.legend()
ax.set(title='Model Errors (Test Set)');

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

fig, ax = plt.subplots(2, 1)
plot_acf(x=forecast_test['errors'], ax=ax[0])
plot_pacf(x=forecast_test['errors'], ax=ax[1]);

We do not see any (partial)-autocorrelation on the errors.

Model Deep Dive

• Model Components

# Plot model components.
fig = model.plot_components(forecast)

• Trend Fit

Let us plot where

from fbprophet.plot import add_changepoints_to_plot

fig = model.plot(forecast)
a = add_changepoints_to_plot(fig.gca(), model, forecast)

Let us compare the true and the fitted trend:

fig, ax = plt.subplots()

sns.lineplot(x='ds', y='trend', data=df, label='trend_true', ax=ax)

ax.fill_between(
    x=forecast['ds'],
    y1=forecast['trend_lower'],
    y2=forecast['trend_upper'],
    color=sns_c[1], 
    alpha=0.25,
    label=r'0.95 credible_interval'
)

sns.lineplot(x='ds', y='trend', data=forecast, label='trend_fit', ax=ax)
ax.legend(loc='upper left')
ax.set(title='Trend Fit', ylabel='');

The trend fit is quite good.

Diagnostics - Time-slice Cross-Validation

The cross_validation function from the fbprophet.diagnostics module allow us to run a time-slice cross-validation on the model by specifying⁷:

• initial: initial training period.
• period: spacing between cutoff dates.
• horizon: forecast horizon on each cross-validation step.

From the documentation:

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.⁸

from fbprophet.diagnostics import cross_validation

df_cv = cross_validation(
    model=model, 
    initial='730 days', 
    period='35 days', 
    horizon = '56 days'
)

df_cv.head()

Remark: For weekly data it might be convenient to choose horizon as a multiple of \(7\) to have the same number of observations on each week of the horizon.

It is then easy to compute some error metrics via the performance_metrics function.

Remark: The metric computations are done on a rolling window specified by the rolling_window parameter. There is a clear explanation of it on the function doc-strings:

Metrics are calculated over a rolling window of cross validation predictions, after sorting by horizon. Averaging is first done within each value of horizon, and then across horizons as needed to reach the window size. The size of that window (number of simulated forecast points) is determined by the rolling_window argument, which specifies a proportion of simulated forecast points to include in each window. rolling_window=0 will compute it separately for each horizon. The default of rolling_window=0.1 will use 10% of the rows in df in each window. rolling_window=1 will compute the metric across all simulated forecast points. The results are set to the right edge of the window.⁹

from fbprophet.diagnostics import performance_metrics

df_p = performance_metrics(df=df_cv, rolling_window=0.1)

df_p.head()

Let us see the how to compute the mae for the case rolling_window = 0.0

df_p = performance_metrics(df=df_cv, rolling_window=0.0)

df_p.head()

We can get the mae column explicitly as:

df_cv.assign(abs_error = lambda x: (x['y'] - x['yhat']).abs()) \
    .assign(horizon = lambda x: x['ds'] - x['cutoff']) \
    .assign(horizon = lambda x: x['horizon']) \
    .groupby('horizon', as_index=False) \
    .agg({'abs_error': np.mean}) \
    .rename(columns={'abs_error': 'mae'})

We can plot these metrics as a function of the horizon parameter.

from fbprophet.plot import plot_cross_validation_metric

fig = plot_cross_validation_metric(df_cv=df_cv, metric='mae', rolling_window=0.1)

Prophet makes the time-slice cross-validation metric computation procedure quite straightforward.

I strongly recommend going to Prophet Documentation and the source code to keep learning abut this forecasting framework (I definitely will!). Moreover, the notebook shows a concrete use case when external regressors are added to the model. Finally, I want to point out an article on how to scale this via Spark.