Forecasting in Python met Prophet

Vraag aan kinderen de dab te doen en ze zullen waarschijnlijk allemaal net even een andere houding aannemen. Zo is het ook onwaarschijnlijk dat een voorspelling perfect samenvalt met observaties in de toekomst. Het leuke aan kinderen is dat ze nog onbezonnen in het leven staan. Met deze gedachte gaan we Prophet uitproberen als forecasting algoritme. Dit artikel volgt de volgende verhaallijn:

Wie is Prophet?

Een profeet is iemand met een visioen over de toekomst. Op haar website wordt Prophet als volgt omschreven:

Prophet is een algoritme voor het voorspellen van tijdreeksen met een additief model waarbij niet-lineaire trends passen bij seizoensinvloeden plus vakantie-effecten. Het werkt het beste voor tijdreeksen met sterke seizoenseffecten. Robuust tegen ontbrekende gegevens, wijzigingen in de trend en gaat doorgaans goed om met uitschieters. Prophet is open source en uitgebracht door Facebook.

Een grondige uitleg kan je lezen in de publicatie Forecasting at Scale.

Onderzoeksopzet

De kwaliteit van de forecast gaan we beoordelen met steekproeven op historische data. Elke steekproef gebruikt de volledige tijdsreeks tot en met de cutoff datum om het algoritme te trainen.

mijn richtingsgevoel zegt dat het alle kanten op kan

Loesje

Vervolgens wordt een forecast berekent met een horizon van 6 maanden na de cutoff datum. Door de forecast te vergelijken met daadwerkelijke observaties krijgen we inzicht in de kwaliteit van het voorspelmodel. In vakjargon noemen we deze methode cross-validation.

Wat gaan we voorspellen

In het voorgaande artikel heb ik beschreven hoe je CBS open data kan downloaden en samenvoegen tot een kerncijfers dataset. Vervolgens hebben we gezien dat er een sterke correlatie is tussen het aantal lopende bijstandsuitkering en het werkloosheidspercentage van 24 maanden terug (de lag).

Voorspelmodel bijstandsuitkeringen

Het doel van dit artikel is het voorspellen van het aantal lopende bijstandsuitkeringen met een horizon van 6 maanden. Het is aannemelijk dat het werkloosheidspercentage invloed heeft op het aantal bijstandsuitkeringen. We nemen het daarom op als extra variabele.

Welke Python packages hebben we nodig?

from copy import copy 
from datetime import datetime
from dateutil.relativedelta import relativedelta
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from prophet import Prophet
sns.set_theme()
sns.set(font_scale=1.2)

Voorspelmodel initiëren

De eerste keer probeerde het voorspelmodel met standaard instellingen. Bij het visualiseren van de steekproeven leek het erop dat de seizoensinvloeden worden overdreven. Ook mag het model naar mijn smaak iets flexibeler omgaan met trendwijzigingen en daarvoor best meer datapunten gebruiken uit de voorliggende historie.

model = Prophet(interval_width=0.8 
            ,uncertainty_samples=200  
            ,daily_seasonality=False 
            ,weekly_seasonality=False
            ,yearly_seasonality=True 
            ,seasonality_mode='additive' 
            ,seasonality_prior_scale = 5 
            ,changepoint_range=0.95 
            ,changepoint_prior_scale=0.07 
           )

model.add_regressor('ww_lag',
                standardize=True
               );

Samenvattend heb ik gekozen voor de volgende instellingen:

• uncertainty samples: standaard worden 1.000 datapunten berekend voor het bepalen van de onzekerheidsintervallen. Kost systeemtijd en daarom verlaagd naar 200. Heeft geen effect op de kwaliteit van de voorspelling
• daily seasonality & weekly seasonality: maandcijfers kennen geen seizoenspatronen op dag- of weekniveau. Expliciet uitgeschakeld. Signaleert Prophet overigens zelf ook
• yearly seasonality: ingeschakeld
• seasonality mode: dit artikel heeft een wat onbezonnen karakter, ingesteld op ‘additive’ zonder verder na te denken over ‘multiplicative’
• seasonality prior scale: bepaald flexibiliteit van de seizoenseffecten. Volgens de documentatie kan verlagen zinvol zijn bij patronen op jaarbasis. Ik kies ervoor om te verlagen van 10 naar 5
• changepoint prior scale: bepaalt de flexibiliteit van de trend, en in het bijzonder hoeveel de trend verandert. Verhoogd van 0.5 naar 0.7
• changepoint range: deel van de train data (historie) waarin de trend kan veranderen. Verhoogd van 0.8 naar 0.95

ook voor Windows supporters draait open source beter op Linux

Loepje

Meer weten over Prophet? Lees dan de documentie.

Steekproeven voor beoordelen forecast

De code hieronder helpt bij de volgende taken: opsplitsen dataset in een train dataset waarop het algoritme kan leren en een test dataset voor de voorspelling. De test dataset heeft het algoritme nog niet eerder gezien. Daarmee sluiten we voorkennis uit.

def train_test_split(df, cutoff_datum, forecast_horizon = 6):
    """
    """
    date_format = '%Y-%m-%d'
    train_tm_datum = datetime.strptime(cutoff_datum, date_format)
    
    forecast_vanaf_datum = train_tm_datum + relativedelta(days=1)
    forecast_tm_datum = forecast_vanaf_datum  + relativedelta(months=forecast_horizon)
    
    forecast_vanaf_datum = forecast_vanaf_datum.strftime(date_format)
    forecast_tm_datum = forecast_tm_datum.strftime(date_format)
    
    # Train inclusief cutoff
    train_df = df.loc['2000-01-01': train_tm_datum,:].reset_index()
    # Voorspel vanaf eerst volgende datum
    future_df = df.loc[forecast_vanaf_datum:forecast_tm_datum,:].reset_index()
    return train_df, future_df

Daarna hebben we functie nodig die het voorspelmodel start voor het berekenen van de forecast. Voorafgaande aan elke voorspelling (dus per steekproef) wordt het algoritme opnieuw getraind op de train dataset.

def forecast_prophet(df, cutoff_datum, model, forecast_horizon = 6):
    train_df, future_df = train_test_split(df=df, cutoff_datum=cutoff_datum)
    model = copy(model)
    model.fit(train_df)
    forecast_df =  model.predict(future_df.drop('y', axis=1))
    return forecast_df.tail(forecast_horizon)

Het laatste stukje code bevat een lijst met de datums waarop de steekproeven worden uitgevoerd. Daarna wordt per datum de forecast berekend en opgeslagen in de forecast dataset.

cutoff_datums = ['2015-04-30','2015-12-31',
                 '2016-06-30',
                 '2017-01-31','2017-12-31',
                 '2018-07-31','2019-05-31','2019-12-31',
                 '2020-03-31','2020-05-31','2020-11-30',
                 '2021-03-31','2021-07-31']

with warnings.catch_warnings():
    # FutureWarning: The frame.append method is deprecated and will be removed 
    warnings.simplefilter("ignore", category=FutureWarning)
    forecasts = pd.DataFrame()
    for d in cutoff_datums:
        f = forecast_prophet(df=df, cutoff_datum=d, model=model)
        f['cutoff'] = d
        forecasts = pd.concat([forecasts,f])
forecasts.set_index('cutoff', inplace=True)

Steekproeven visualiseren

Data visualisatie doen we met Matplotlip. De code gebruikt de dataset met forecasts om de voorspellingen in het lijndiagram te tekenen inclusief de betrouwbaarheidsintervallen.

# Canvas
fig_dims = (16, 6)
fig, ax = plt.subplots(figsize=fig_dims)
fig.canvas.draw()
# Grafiek
plt.ioff() 
ax.plot(df.reset_index()['ds'], df['y'], 'o', c='black', label='Aantal lopende bijstandsuitkeringen per maand')
ax.set_ylabel('')
ax.set_xlabel('')
ax.yaxis.grid(False)
for d in cutoff_datums:    
    ax.plot(forecasts.loc[d]['ds'], forecasts.loc[d]['yhat'], ls='-', c='#0072B2')
    ax.fill_between(forecasts.loc[d]['ds'], forecasts.loc[d]['yhat_lower'], forecasts.loc[d]['yhat_upper'], color='#0072B2', alpha=0.2)
plt.title("FORECASTING MET PROPHET, RESULTAAT STEEKPROEVEN", y=1.0);

Conclusie

Kijken we naar de uitkomst van de steekproeven, de cross validation, dan lijkt Prophet goed te passen bij de data. Wat vooral indruk maakt is het herstel van het voorspellend vermogen na de eerste lockdown. Aan het begin van de coronacrisis zit het model er vooral erg naast.

Prophet is relatief eenvoudig in gebruik, werkt vrijwel ‘out-of-the-box’ en draait duidelijk het beste op Linux. Ben je een Windows gebruiker? Overweeg dan Ubuntu via WSL.

What’s next

Rob Hyndman is positief over Prophet in zijn boek Forecasting: Principles and Practice, maar komt ook tot de conclusie dat het algoritme nauwelijks beter presenteert dan ETS en Arima.

In komende artikelen zoom ik daarom in op klassieke methoden voor forecasting om deze te vergelijken met de uitkomsten in dit artikel.