Tijdreeksanalyse met Python: Van Pandas tot Prophet
Tijdreeksdata is overal om ons heen. Aandelenkoersen, weersvoorspellingen, websiteverkeer, energieverbruik — zodra je data hebt die over de tijd verandert, heb je te maken met een tijdreeks. En eerlijk gezegd, als data scientist wil je niet alleen snappen wat er is gebeurd, maar ook voorspellen wat er gaat komen. Dat is waar het pas echt interessant wordt.
Python biedt een ongelooflijk rijke set aan tools hiervoor. Van basismanipulatie met pandas tot geavanceerde voorspelmodellen met statsmodels en Prophet — het hele spectrum is afgedekt. In deze handleiding neem ik je stap voor stap mee: data inladen en verkennen, decompositie en stationariteitstoetsen, en uiteindelijk het bouwen van ARIMA- en Prophet-modellen die echte voorspellingen genereren.
Of je nu een ervaren data scientist bent die z'n tijdreeksvaardigheden wil opfrissen, of een beginner die voor het eerst met temporele data werkt — pak een kop koffie en lees verder.
Wat Maakt Tijdreeksdata Bijzonder?
Voordat we in de code duiken, even dit: tijdreeksdata is fundamenteel anders dan andere vormen van data. Het belangrijkste verschil? De volgorde doet ertoe.
Bij een gewone dataset kun je rijen door elkaar husselen zonder dat je informatie verliest. Bij een tijdreeks is dat onmogelijk — de temporele structuur bevat cruciale informatie.
Dit heeft drie grote implicaties:
- Autocorrelatie: Waarden zijn vaak gecorreleerd met hun eigen verleden. De temperatuur van vandaag hangt sterk samen met die van gisteren. Dit is eigenlijk de kern van alle tijdreeksmodellering.
- Trend en seizoensgebondenheid: Data kan een langetermijntrend vertonen (stijgend of dalend) en seizoenspatronen hebben die zich herhalen — dagelijks, wekelijks, jaarlijks.
- Niet-willekeurig splitsen: Bij het trainen en testen van modellen kun je niet zomaar willekeurig splitsen. Je moet altijd het verleden gebruiken om de toekomst te voorspellen, nooit andersom. (Dit is een veelgemaakte beginnersfout, trouwens.)
Goed, laten we beginnen met de basis: tijdreeksdata manipuleren met pandas.
Tijdreeksdata Laden en Bewerken met Pandas
Pandas is de ruggengraat van elke tijdreeksanalyse in Python. Met het krachtige DatetimeIndex en methoden als resample() en rolling() kun je vrijwel alles doen wat je nodig hebt.
DatetimeIndex: De Basis
Het allereerste wat je moet doen bij tijdreeksdata? Zorgen dat je datumkolom wordt herkend als een datetime-type en ingesteld als index:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# Data laden en datetime-index instellen
df = pd.read_csv("energieverbruik.csv", parse_dates=["datum"])
df = df.set_index("datum")
df = df.sort_index()
# Controleer het type van de index
print(f"Index type: {type(df.index)}")
print(f"Frequentie: {df.index.freq}")
print(f"Bereik: {df.index.min()} tot {df.index.max()}")
print(df.head())Heb je geen bestaand CSV-bestand bij de hand? Geen probleem. Je kunt eenvoudig een voorbeeldtijdreeks genereren om mee te oefenen:
# Voorbeelddata genereren: dagelijks energieverbruik (kWh)
np.random.seed(42)
datums = pd.date_range(start="2022-01-01", end="2025-12-31", freq="D")
# Realistisch patroon: trend + seizoensgebondenheid + ruis
trend = np.linspace(100, 130, len(datums))
seizoen = 20 * np.sin(2 * np.pi * np.arange(len(datums)) / 365.25)
week_effect = 5 * np.sin(2 * np.pi * np.arange(len(datums)) / 7)
ruis = np.random.normal(0, 5, len(datums))
verbruik = trend + seizoen + week_effect + ruis
df = pd.DataFrame({"verbruik_kwh": verbruik}, index=datums)
df.index.name = "datum"
print(df.describe())
print(f"\nAantal datapunten: {len(df)}")
print(f"Periode: {df.index[0].strftime('%Y-%m-%d')} tot {df.index[-1].strftime('%Y-%m-%d')}")Resamplen: Frequentie Aanpassen
Met resample() kun je de frequentie van je tijdreeks veranderen. Superhandig wanneer je dagelijkse data wilt aggregeren naar wekelijks of maandelijks niveau:
# Dagelijks naar maandelijks gemiddelde (downsampling)
df_maandelijks = df["verbruik_kwh"].resample("ME").mean()
print("Maandelijks gemiddelde:")
print(df_maandelijks.head(12))
# Dagelijks naar wekelijks totaal
df_wekelijks = df["verbruik_kwh"].resample("W").sum()
# Dagelijks naar kwartaal met meerdere aggregaties
df_kwartaal = df["verbruik_kwh"].resample("QE").agg(
["mean", "std", "min", "max"]
)
print("\nKwartaaloverzicht:")
print(df_kwartaal.head(8))
# Visualisatie van verschillende frequenties
fig, axes = plt.subplots(3, 1, figsize=(14, 10), sharex=False)
df["verbruik_kwh"].plot(ax=axes[0], alpha=0.5, title="Dagelijks")
df_wekelijks.plot(ax=axes[1], title="Wekelijks totaal")
df_maandelijks.plot(ax=axes[2], title="Maandelijks gemiddelde", marker="o")
plt.tight_layout()
plt.savefig("frequentie_vergelijking.png", dpi=150)
plt.show()Rolling Windows: Voortschrijdende Statistieken
Rolling windows zijn onmisbaar voor het gladstrijken van ruis en het ontdekken van trends. Ik gebruik ze eigenlijk bij elk project als eerste verkennende stap. Met rolling() bereken je statistieken over een verschuivend venster:
# Voortschrijdend gemiddelde en standaarddeviatie
df["rolling_7d"] = df["verbruik_kwh"].rolling(window=7).mean()
df["rolling_30d"] = df["verbruik_kwh"].rolling(window=30).mean()
df["rolling_std_30d"] = df["verbruik_kwh"].rolling(window=30).std()
# Exponentieel gewogen voortschrijdend gemiddelde (EWM)
# Geeft meer gewicht aan recente observaties
df["ewm_30d"] = df["verbruik_kwh"].ewm(span=30).mean()
# Visualisatie
fig, ax = plt.subplots(figsize=(14, 6))
df["verbruik_kwh"].plot(ax=ax, alpha=0.3, label="Dagelijks")
df["rolling_7d"].plot(ax=ax, label="7-daags gemiddelde")
df["rolling_30d"].plot(ax=ax, label="30-daags gemiddelde", linewidth=2)
df["ewm_30d"].plot(ax=ax, label="EWM 30 dagen", linestyle="--")
ax.legend()
ax.set_title("Energieverbruik met Voortschrijdende Gemiddelden")
ax.set_ylabel("kWh")
plt.tight_layout()
plt.show()Het verschil tussen een gewoon voortschrijdend gemiddelde en een exponentieel gewogen gemiddelde (EWM)? EWM legt meer nadruk op recente data. Dat maakt het responsiever op veranderingen, maar ook gevoeliger voor uitschieters. Het hangt er dus maar net van af wat je nodig hebt.
Tijdreeksdecompositie: Trend, Seizoen en Residu
Een van de krachtigste technieken in tijdreeksanalyse is decompositie: het opsplitsen van je tijdreeks in z'n samenstellende componenten. We onderscheiden er typisch drie:
- Trend: De langetermijnrichting van de data — stijgend, dalend of vlak.
- Seizoensgebondenheid: Herhalende patronen met een vaste periode (dagelijks, wekelijks, jaarlijks).
- Residu (ruis): Wat overblijft nadat trend en seizoen zijn verwijderd. Eigenlijk het deel dat je model niet kan verklaren.
Klassieke Decompositie met statsmodels
from statsmodels.tsa.seasonal import seasonal_decompose
# Maandelijkse data gebruiken voor duidelijkere decompositie
df_maand = df["verbruik_kwh"].resample("ME").mean()
# Additieve decompositie (als seizoenseffect constant is)
decomp_add = seasonal_decompose(df_maand, model="additive", period=12)
# Multiplicatieve decompositie (als seizoenseffect groeit met trend)
decomp_mult = seasonal_decompose(df_maand, model="multiplicative", period=12)
# Visualisatie van additieve decompositie
fig = decomp_add.plot()
fig.set_size_inches(14, 10)
fig.suptitle("Additieve Tijdreeksdecompositie", fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
# Componenten apart bekijken
print("Trend (eerste 5 waarden):")
print(decomp_add.trend.dropna().head())
print("\nSeizoenscomponent (eerste 12 waarden):")
print(decomp_add.seasonal.head(12))STL Decompositie: De Moderne Aanpak
De klassieke decompositie heeft z'n beperkingen. Met name: de trend kan niet van richting veranderen en het seizoenspatroon moet strikt constant zijn. STL (Seasonal and Trend decomposition using Loess) lost deze problemen op:
from statsmodels.tsa.seasonal import STL
# STL decompositie — flexibeler en robuuster
stl = STL(df_maand, period=12, robust=True)
result = stl.fit()
fig = result.plot()
fig.set_size_inches(14, 10)
fig.suptitle("STL Decompositie (Robuust)", fontsize=14, y=1.02)
plt.tight_layout()
plt.show()
# De sterkte van trend en seizoensgebondenheid kwantificeren
resid_var = result.resid.var()
trend_sterkte = max(0, 1 - resid_var / (result.trend + result.resid).var())
seizoen_sterkte = max(0, 1 - resid_var / (result.seasonal + result.resid).var())
print(f"Sterkte van de trend: {trend_sterkte:.3f}")
print(f"Sterkte van seizoensgebondenheid: {seizoen_sterkte:.3f}")Mijn advies? Gebruik STL in plaats van klassieke decompositie, tenzij je een specifieke reden hebt om dat niet te doen. Het is robuuster tegen uitschieters (zeker met robust=True) en staat toe dat seizoenspatronen geleidelijk veranderen over de tijd.
Stationariteit: De Sleutel tot Goede Voorspellingen
Oké, nu wordt het even technisch — maar blijf erbij, want dit is echt belangrijk. Voordat je een voorspelmodel als ARIMA kunt bouwen, moet je weten of je tijdreeks stationair is. Dat betekent: constante statistische eigenschappen over de tijd. Het gemiddelde, de variantie en de autocorrelatie veranderen niet.
Waarom dit ertoe doet? De meeste klassieke voorspelmodellen gaan ervan uit dat de data stationair is. Als dat niet zo is, moet je de data eerst transformeren.
De Augmented Dickey-Fuller Test
from statsmodels.tsa.stattools import adfuller
def test_stationariteit(tijdreeks, naam=""):
"""Voer de Augmented Dickey-Fuller test uit."""
resultaat = adfuller(tijdreeks.dropna(), autolag="AIC")
print(f"=== ADF Test: {naam} ===")
print(f"Test Statistiek: {resultaat[0]:.4f}")
print(f"P-waarde: {resultaat[1]:.4f}")
print(f"Aantal lags: {resultaat[2]}")
print(f"Observaties: {resultaat[3]}")
print("Kritieke waarden:")
for sleutel, waarde in resultaat[4].items():
print(f" {sleutel}: {waarde:.4f}")
if resultaat[1] < 0.05:
print(">>> RESULTAAT: Tijdreeks is STATIONAIR (p < 0.05)")
else:
print(">>> RESULTAAT: Tijdreeks is NIET stationair (p >= 0.05)")
print()
# Test op originele data
test_stationariteit(df_maand, "Origineel maandelijks verbruik")
# Eerste differentiatie toepassen
df_diff = df_maand.diff().dropna()
test_stationariteit(df_diff, "Na eerste differentiatie")Differentiatie: Niet-stationariteit Oplossen
Is je tijdreeks niet stationair? Dan is differentiatie je beste vriend. Het idee is simpel: je trekt de vorige waarde af van de huidige.
# Eerste orde differentiatie
df_diff1 = df_maand.diff(1).dropna()
# Seizoensdifferentiatie (verschil met dezelfde maand vorig jaar)
df_diff_seizoen = df_maand.diff(12).dropna()
# Combinatie: eerst seizoensdifferentiatie, dan eerste differentiatie
df_diff_beide = df_maand.diff(12).diff(1).dropna()
# Visualisatie
fig, axes = plt.subplots(4, 1, figsize=(14, 12))
df_maand.plot(ax=axes[0], title="Origineel")
df_diff1.plot(ax=axes[1], title="Eerste differentiatie (d=1)")
df_diff_seizoen.plot(ax=axes[2], title="Seizoensdifferentiatie (D=1, m=12)")
df_diff_beide.plot(ax=axes[3], title="Beide differentiaties")
for ax in axes:
ax.axhline(y=0, color="r", linestyle="--", alpha=0.3)
plt.tight_layout()
plt.show()ACF en PACF: De Parameters van ARIMA Bepalen
Nu komen we bij een stuk dat ik persoonlijk erg leuk vind. De AutoCorrelation Function (ACF) en Partial AutoCorrelation Function (PACF) zijn je kompas voor het kiezen van de juiste ARIMA-parameters. Ze laten zien hoe sterk de relatie is tussen een observatie en haar verleden.
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
# ACF en PACF van originele data
plot_acf(df_maand.dropna(), ax=axes[0, 0], lags=36,
title="ACF - Origineel")
plot_pacf(df_maand.dropna(), ax=axes[0, 1], lags=36,
title="PACF - Origineel", method="ywm")
# ACF en PACF na differentiatie
plot_acf(df_diff1, ax=axes[1, 0], lags=36,
title="ACF - Na differentiatie")
plot_pacf(df_diff1, ax=axes[1, 1], lags=36,
title="PACF - Na differentiatie", method="ywm")
plt.tight_layout()
plt.show()En hoe lees je deze plots dan? Hier zijn de vuistregels:
- ACF daalt geleidelijk, PACF knipt scherp af na lag p: Gebruik een AR(p)-model.
- ACF knipt scherp af na lag q, PACF daalt geleidelijk: Gebruik een MA(q)-model.
- Beide dalen geleidelijk: Dan heb je een ARMA(p,q)-model nodig.
- Significante pieken bij seizoenslags (12, 24, 36): Voeg seizoenscomponenten toe — je hebt SARIMA nodig.
ARIMA en SARIMA: Klassieke Voorspelmodellen
ARIMA staat voor AutoRegressive Integrated Moving Average en is, laten we eerlijk zijn, het werkpaard van de tijdreeksvoorspelling. Het combineert drie componenten:
- AR (AutoRegressief): Gebruikt de relatie tussen een observatie en een aantal vertraagde observaties (parameter p).
- I (Integrated): Het aantal keer dat de data gedifferentieerd is om stationariteit te bereiken (parameter d).
- MA (Moving Average): Gebruikt de relatie tussen een observatie en de restfouten van vertraagde voorspellingen (parameter q).
Een ARIMA-model Bouwen
from statsmodels.tsa.arima.model import ARIMA
from sklearn.metrics import mean_absolute_error, mean_squared_error
# Train-test split (altijd chronologisch!)
train = df_maand[:"2024-12"]
test = df_maand["2025-01":]
print(f"Training: {train.index[0]} tot {train.index[-1]} ({len(train)} punten)")
print(f"Test: {test.index[0]} tot {test.index[-1]} ({len(test)} punten)")
# ARIMA(1,1,1) model fitten
model = ARIMA(train, order=(1, 1, 1))
resultaat = model.fit()
# Modelsamenvatting bekijken
print(resultaat.summary())
# Voorspellingen maken
voorspelling = resultaat.forecast(steps=len(test))
# Evaluatie
mae = mean_absolute_error(test, voorspelling)
rmse = np.sqrt(mean_squared_error(test, voorspelling))
mape = np.mean(np.abs((test - voorspelling) / test)) * 100
print(f"\nModel Evaluatie:")
print(f"MAE: {mae:.2f} kWh")
print(f"RMSE: {rmse:.2f} kWh")
print(f"MAPE: {mape:.2f}%")SARIMA: Seizoensgebonden ARIMA
Als je data seizoenspatronen vertoont (en dat is eerlijk gezegd vaker wel dan niet), heb je SARIMA nodig. Dit voegt seizoensparameters (P, D, Q, m) toe aan het standaard ARIMA-model:
from statsmodels.tsa.statespace.sarimax import SARIMAX
# SARIMA model met seizoenscomponent
# order = (p, d, q) voor het niet-seizoensdeel
# seasonal_order = (P, D, Q, m) voor het seizoensdeel
model_sarima = SARIMAX(
train,
order=(1, 1, 1),
seasonal_order=(1, 1, 1, 12),
enforce_stationarity=False,
enforce_invertibility=False
)
resultaat_sarima = model_sarima.fit(disp=False)
print(resultaat_sarima.summary())
# Voorspellingen met betrouwbaarheidsinterval
voorspelling_sarima = resultaat_sarima.get_forecast(steps=len(test))
gemiddelde = voorspelling_sarima.predicted_mean
ci = voorspelling_sarima.conf_int(alpha=0.05)
# Visualisatie
fig, ax = plt.subplots(figsize=(14, 6))
train[-60:].plot(ax=ax, label="Training data")
test.plot(ax=ax, label="Werkelijke waarden", color="green")
gemiddelde.plot(ax=ax, label="SARIMA voorspelling", color="red")
ax.fill_between(ci.index, ci.iloc[:, 0], ci.iloc[:, 1],
color="red", alpha=0.15, label="95% betrouwbaarheidsinterval")
ax.legend()
ax.set_title("SARIMA Voorspelling vs. Werkelijkheid")
ax.set_ylabel("kWh")
plt.tight_layout()
plt.show()
# Evaluatie SARIMA
mae_s = mean_absolute_error(test, gemiddelde)
rmse_s = np.sqrt(mean_squared_error(test, gemiddelde))
print(f"\nSARIMA Evaluatie:")
print(f"MAE: {mae_s:.2f} kWh")
print(f"RMSE: {rmse_s:.2f} kWh")Auto ARIMA: Automatisch de Beste Parameters Vinden
Handmatig de juiste (p, d, q) en (P, D, Q, m) kiezen is... eerlijk gezegd best tijdrovend. Gelukkig automatiseert de bibliotheek pmdarima dit hele proces:
import pmdarima as pm
# Auto ARIMA: zoekt automatisch de optimale parameters
auto_model = pm.auto_arima(
train,
seasonal=True,
m=12, # seizoensperiode
start_p=0, max_p=3, # bereik voor p
start_q=0, max_q=3, # bereik voor q
d=None, # laat auto_arima d bepalen
D=None, # laat auto_arima D bepalen
trace=True, # toon zoekproces
error_action="ignore",
suppress_warnings=True,
stepwise=True # snellere zoekstrategie
)
print(f"\nBeste model: {auto_model.summary()}")
print(f"AIC: {auto_model.aic():.2f}")
# Voorspellingen met het automatisch gekozen model
voorspelling_auto, ci_auto = auto_model.predict(
n_periods=len(test),
return_conf_int=True,
alpha=0.05
)Die stepwise=True parameter is trouwens belangrijk: het gebruikt een slimme zoekstrategie die veel sneller is dan brute-force alle combinaties doorzoeken. Voor de meeste datasets vindt het een goed of optimaal model.
Prophet: Voorspellen op Schaal
Prophet, ontwikkeld door het data science team van Meta (voorheen Facebook), is een ander populair framework voor tijdreeksvoorspellingen. Het is specifiek ontworpen voor zakelijke tijdreeksen met sterke seizoenseffecten en meerdere seizoenen aan historische data.
Waarom Prophet?
Prophet heeft een paar unieke voordelen ten opzichte van ARIMA:
- Gebruiksvriendelijkheid: Minder parameterkeuzes nodig. Je kunt sneller aan de slag.
- Meerdere seizoenspatronen: Kan tegelijkertijd dagelijkse, wekelijkse en jaarlijkse seizoensgebondenheid modelleren.
- Feestdagen en speciale events: Ingebouwde ondersteuning voor feestdagen — super handig voor retail en e-commerce data.
- Robuust tegen ontbrekende data: Gaat prima om met gaten in je tijdreeks.
- Interpreteerbaarheid: De componenten zijn eenvoudig te visualiseren en te begrijpen.
Prophet Installeren en Gebruiken
# Installatie
# pip install prophet
from prophet import Prophet
# Prophet vereist specifieke kolomnamen: ds (datum) en y (waarde)
df_prophet = df_maand.reset_index()
df_prophet.columns = ["ds", "y"]
# Train-test split
train_prophet = df_prophet[df_prophet["ds"] <= "2024-12-31"]
test_prophet = df_prophet[df_prophet["ds"] > "2024-12-31"]
# Model initialiseren en trainen
model_prophet = Prophet(
yearly_seasonality=True,
weekly_seasonality=False, # maandelijkse data, geen weekpatroon
daily_seasonality=False,
changepoint_prior_scale=0.05, # flexibiliteit van trendveranderingen
seasonality_prior_scale=10.0, # sterkte van seizoenscomponent
interval_width=0.95 # 95% betrouwbaarheidsinterval
)
model_prophet.fit(train_prophet)
# Toekomstige datums genereren voor voorspelling
toekomst = model_prophet.make_future_dataframe(periods=len(test_prophet), freq="ME")
voorspelling_prophet = model_prophet.predict(toekomst)
# Visualisatie
fig = model_prophet.plot(voorspelling_prophet)
plt.title("Prophet Voorspelling: Energieverbruik")
plt.ylabel("kWh")
plt.tight_layout()
plt.show()
# Componentenplot
fig2 = model_prophet.plot_components(voorspelling_prophet)
plt.tight_layout()
plt.show()Feestdagen Toevoegen aan Prophet
Een van de sterkste features van Prophet is het toevoegen van feestdagen en speciale events. In mijn ervaring maakt dit voor zakelijke data echt een verschil:
# Nederlandse feestdagen definiëren
feestdagen_nl = pd.DataFrame({
"holiday": "koningsdag",
"ds": pd.to_datetime(["2022-04-27", "2023-04-27",
"2024-04-27", "2025-04-27"]),
"lower_window": -1, # 1 dag voor het event
"upper_window": 1 # 1 dag na het event
})
kerst = pd.DataFrame({
"holiday": "kerst",
"ds": pd.to_datetime(["2022-12-25", "2023-12-25",
"2024-12-25", "2025-12-25"]),
"lower_window": -2,
"upper_window": 2
})
alle_feestdagen = pd.concat([feestdagen_nl, kerst], ignore_index=True)
# Model met feestdagen
model_feest = Prophet(
yearly_seasonality=True,
holidays=alle_feestdagen
)
model_feest.fit(train_prophet)Cross-validatie met Prophet
Prophet biedt ingebouwde cross-validatie die speciaal is ontworpen voor tijdreeksdata. Dit is echt een van de handiger features:
from prophet.diagnostics import cross_validation, performance_metrics
# Cross-validatie uitvoeren
# initial: trainingsperiode
# period: afstand tussen cutoff-datums
# horizon: voorspelhorizon
cv_resultaten = cross_validation(
model_prophet,
initial="730 days", # 2 jaar initiële training
period="90 days", # elke 90 dagen een nieuw cutoff
horizon="180 days" # 6 maanden vooruit voorspellen
)
# Prestatiemetrieken berekenen
metrieken = performance_metrics(cv_resultaten)
print(metrieken[["horizon", "mae", "rmse", "mape"]].tail(10))
# Visualisatie van prestatie over horizon
from prophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(cv_resultaten, metric="mape")
plt.title("MAPE per Voorspelhorizon")
plt.tight_layout()
plt.show()Moderne Alternatieven: Het Nixtla Ecosysteem
Naast de klassieke tools is er de afgelopen jaren iets interessants ontstaan: het Nixtla-ecosysteem. Dit open-source project biedt een suite van bibliotheken die specifiek zijn gebouwd voor snelheid en schaalbaarheid. Als je met grote hoeveelheden tijdreeksen werkt, is dit het verkennen waard.
StatsForecast: Razendsnel Statistisch Voorspellen
StatsForecast van Nixtla bevat geoptimaliseerde implementaties van alle klassieke modellen — maar dan véél sneller dan statsmodels. We hebben het over AutoARIMA, AutoETS, AutoCES en meer:
# pip install statsforecast
from statsforecast import StatsForecast
from statsforecast.models import AutoARIMA, AutoETS, MSTL
# Data voorbereiden in het juiste formaat
df_sf = df_maand.reset_index()
df_sf.columns = ["ds", "y"]
df_sf["unique_id"] = "energieverbruik"
# Meerdere modellen tegelijkertijd fitten
sf = StatsForecast(
models=[
AutoARIMA(season_length=12),
AutoETS(season_length=12),
MSTL(season_length=12),
],
freq="ME",
n_jobs=-1 # gebruik alle CPU-cores
)
# Fitten en voorspellen
sf.fit(df_sf[df_sf["ds"] <= "2024-12-31"])
voorspellingen = sf.predict(h=12)
print(voorspellingen)Het grote voordeel van StatsForecast is niet alleen de snelheid. Het kan ook moeiteloos miljoenen tijdreeksen tegelijkertijd verwerken. Heb je voorspellingen nodig voor duizenden producten in een webshop? Dit is de tool die je wilt.
MLForecast: Machine Learning voor Tijdreeksen
Wil je machine learning modellen zoals LightGBM of XGBoost inzetten voor tijdreeksvoorspellingen? MLForecast maakt dit verrassend eenvoudig door automatisch lag-features en rolling statistics te genereren:
# pip install mlforecast lightgbm
from mlforecast import MLForecast
from lightgbm import LGBMRegressor
# MLForecast automatiseert feature engineering voor tijdreeksen
mlf = MLForecast(
models=[LGBMRegressor(n_estimators=100, learning_rate=0.1)],
freq="ME",
lags=[1, 2, 3, 6, 12], # lag-features
lag_transforms={
1: ["expanding_mean"], # cumulatief gemiddelde
12: ["rolling_mean_3", "rolling_mean_6"] # voortschrijdende gemiddelden
},
date_features=["month", "quarter"] # kalenderfeatures
)
mlf.fit(df_sf[df_sf["ds"] <= "2024-12-31"])
ml_voorspellingen = mlf.predict(h=12)
print(ml_voorspellingen)Modellen Vergelijken: Welk Model Past Bij Jouw Data?
Met zoveel opties is de grote vraag natuurlijk: welk model moet je kiezen? Er is helaas geen universeel antwoord (dat zou te makkelijk zijn), maar hier zijn wat praktische richtlijnen:
- ARIMA/SARIMA: Goed voor univariate tijdreeksen met duidelijke autocorrelatie. Betrouwbaar, goed begrepen, en ideaal wanneer interpreteerbaarheid belangrijk is.
- Prophet: Uitstekend voor zakelijke tijdreeksen met meerdere seizoenspatronen, feestdagen en ontbrekende data. Minder geschikt voor heel korte tijdreeksen.
- StatsForecast: Wanneer je snelheid en schaal nodig hebt. Perfect voor duizenden tot miljoenen tijdreeksen tegelijk.
- MLForecast/LightGBM: Wanneer je veel exogene variabelen hebt (weer, prijs, promoties). Let op: minder sterk in het extrapoleren van trends.
Een Eerlijke Vergelijking Opzetten
import warnings
warnings.filterwarnings("ignore")
# Alle voorspellingen verzamelen
resultaten = {}
# 1. SARIMA
model_s = SARIMAX(train, order=(1,1,1), seasonal_order=(1,1,1,12))
res_s = model_s.fit(disp=False)
resultaten["SARIMA"] = res_s.forecast(steps=len(test))
# 2. Auto ARIMA
auto = pm.auto_arima(train, seasonal=True, m=12,
suppress_warnings=True, stepwise=True)
resultaten["Auto ARIMA"] = pd.Series(
auto.predict(n_periods=len(test)), index=test.index
)
# Evaluatiemetrieken berekenen
print(f"{'Model':<15} {'MAE':<10} {'RMSE':<10} {'MAPE (%)':<10}")
print("-" * 45)
for naam, pred in resultaten.items():
mae = mean_absolute_error(test, pred)
rmse = np.sqrt(mean_squared_error(test, pred))
mape = np.mean(np.abs((test - pred) / test)) * 100
print(f"{naam:<15} {mae:<10.2f} {rmse:<10.2f} {mape:<10.2f}")
# Visualisatie van alle modellen
fig, ax = plt.subplots(figsize=(14, 6))
train[-36:].plot(ax=ax, label="Training", color="gray", alpha=0.5)
test.plot(ax=ax, label="Werkelijk", color="black", linewidth=2)
for naam, pred in resultaten.items():
pred.plot(ax=ax, label=naam, linestyle="--")
ax.legend()
ax.set_title("Modelvergelijking: Tijdreeksvoorspellingen")
ax.set_ylabel("kWh")
plt.tight_layout()
plt.show()Praktische Tips en Veelgemaakte Fouten
Na het werken met talloze tijdreeksprojecten zijn er een aantal lessen die steeds terugkomen. Hier zijn de belangrijkste — bewaar ze ergens, want je hebt ze nodig.
1. Pas Op voor Data-lekkage
De meest voorkomende fout bij tijdreeksmodellering is het per ongeluk gebruiken van toekomstige informatie. En het is subtieler dan je denkt:
# FOUT: Scaler fitten op de hele dataset (inclusief test)
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
df_scaled = scaler.fit_transform(df) # Lekkage!
# GOED: Scaler alleen fitten op trainingsdata
scaler = StandardScaler()
train_scaled = scaler.fit_transform(train.values.reshape(-1, 1))
test_scaled = scaler.transform(test.values.reshape(-1, 1))2. Gebruik Altijd Meerdere Evaluatiemetrieken
Vertrouw nooit op slechts één metriek. MAE, RMSE en MAPE vertellen elk een ander verhaal:
- MAE (Mean Absolute Error): Gemiddelde absolute fout. Makkelijk te interpreteren in de eenheid van je data.
- RMSE (Root Mean Squared Error): Straft grote fouten harder af. Nuttig als grote afwijkingen kostbaar zijn.
- MAPE (Mean Absolute Percentage Error): Procentuele fout, handig voor vergelijking tussen datasets. Maar pas op bij waarden dicht bij nul!
3. Controleer de Residuen
Na het fitten van je model zijn de residuen je beste diagnose-instrument. Als het model goed werkt, moeten ze eruitzien als witte ruis — geen patronen, geen autocorrelatie:
from statsmodels.stats.diagnostic import acorr_ljungbox
# Residuen van het SARIMA-model
residuen = resultaat_sarima.resid
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
# Tijdplot van residuen
residuen.plot(ax=axes[0, 0], title="Residuen over Tijd")
axes[0, 0].axhline(y=0, color="r", linestyle="--")
# Histogram
residuen.hist(ax=axes[0, 1], bins=30, edgecolor="black")
axes[0, 1].set_title("Verdeling van Residuen")
# ACF van residuen
plot_acf(residuen.dropna(), ax=axes[1, 0], lags=24,
title="ACF van Residuen")
# Q-Q plot
from scipy import stats
stats.probplot(residuen.dropna(), dist="norm", plot=axes[1, 1])
axes[1, 1].set_title("Q-Q Plot")
plt.tight_layout()
plt.show()
# Ljung-Box test voor autocorrelatie in residuen
lb_test = acorr_ljungbox(residuen.dropna(), lags=20, return_df=True)
print("Ljung-Box Test (p-waarden > 0.05 = geen significante autocorrelatie):")
print(lb_test)4. Denk Na Over Je Voorspelhorizon
Hoe verder je vooruit voorspelt, hoe breder je betrouwbaarheidsintervallen worden. Dat is geen bug in je model — dat is gewoon hoe de werkelijkheid werkt. Wees realistisch:
- Korte termijn (1-7 dagen): De meeste modellen doen het hier prima.
- Middellange termijn (1-3 maanden): Seizoensmodellen als SARIMA en Prophet zijn hier sterk.
- Lange termijn (6+ maanden): Betrouwbaarheidsintervallen worden breed. Gebruik voorspellingen als richtlijn, niet als waarheid.
Van Analyse naar Productie: Volgende Stappen
Je hebt nu een solide basis in tijdreeksanalyse met Python. Van data-exploratie met pandas tot het bouwen en evalueren van voorspelmodellen — de tools zijn er.
Maar het echte werk begint eigenlijk pas wanneer je deze modellen gaat inzetten in productie. Hier zijn wat suggesties voor als je verder wilt:
- Automatiseer je pipeline: Gebruik tools als Apache Airflow of Prefect om je voorspelpipeline op schema te draaien.
- Monitor je modellen: Modelprestaties verslechteren over tijd (concept drift). Implementeer monitoring en automatische hertraining.
- Experimenteer met deep learning: Voor complexe patronen en grote datasets zijn NeuralForecast (Nixtla) en Darts het verkennen waard.
- Verken multivariate tijdreeksen: VAR-modellen en ML-benaderingen kunnen relaties tussen meerdere tijdreeksen modelleren.
Tijdreeksanalyse is eerlijk gezegd een vakgebied waar je nooit uitgeleerd raakt. Elke dataset brengt nieuwe uitdagingen. Maar met de technieken uit deze gids heb je een stevig fundament om op te bouwen. Succes met je volgende project!