Forecasting des séries temporelles avec la librairie fbprophet

En 2017, les équipes de recherche de Facebook publiaient ce papier qui introduira la librairie fbprophet, disponible en R et en Python (pas de jaloux !). Cet outil peut être rangé dans la catégorie des modèles additifs généraux car il décompose une série temporelle de la sorte :

y(t) = g(t) + s(t) + h(t) + e(t)

avec respectivement :

  • g(t) : la tendance (linéaire ou logistique)
  • s(s) : une ou plusieurs composantes saisonnières (annuelle, hebdomadaire ou quotidienne)
  • h(t) : l’effet des vacances ou de jours spécifiques qui pourront être paramétrés
  • e(t) : l’erreur, bruit aléatoire qui mesure l’écart entre le modèle et les données réelles

Si l’on retient le terme “additif”, il est bien sûr possible de modifier le modèle pour le rendre “multiplicatif” (observez les crêtes de votre série temporelles, si celles-ci forment un cône, le modèle est sans doute multiplicatif).

La grande force de ce type de modèle tient dans sa capacité à être interprété ainsi que dans la clarté des représentations graphiques de la décomposition. Il sera particulièrement adapté à des phénomènes comme… la fréquentation d’un réseau social mais aussi des mesures économiques fortement soumises à des saisonnalités et aux périodes de vacances d’un pays. Il est si simple à mettre en œuvre qu’il gagnera à être comparé à des modèles plus classiques comme SARIMA (nous en reparlerons en toute fin d’article).

Installation de la librairie

Pour ne pas “polluer” notre installation locale avec de nouveaux packages et leurs dépendances, nous créons au préalable un environnement virtuel, à l’aide du package pyenv pour Windows (voir ce GitHub). Ne pas oublier d’ajouter les variables d’environnement et de redémarrer votre terminal ou IDE pour terminer l’installation.

pyenv install 3.8.10

pyenv shell 3.8.10

Nous pouvons alors installer les librairies suivantes :

pip install pystan==2.19.1.1

pip install fbprophet

Pystan est une librairie pour l’inférence bayésienne. Si cette installation ne fonctionne pas (les dépendances sont parfois capricieuses…), vous pouvez utiliser un prompt Anaconda et tenter la commande suivante :

conda install -c conda-forge fbprophet

Vérifiez enfin la bonne installation en choisissant l’interpréteur voulu (ici, conda) dans Visual Studio Code.

Lancez ensuite Python dans le terminal et testez un import de la librairie.

Mise en pratique

Déroulons maintenant un exemple simple, de bout en bout, en réalisant quelques tentatives d’optimisation du modèle. RTE France met à disposition les données énergétiques “eco2mix” dont la profondeur d’historique et le niveau de détail sera intéressant pour évaluer notre outil. J’ai pris connaissance de ce jeu de données dans cet excellent article du blog de Publicis Sapient.

Une seule contrainte dans la façon dont les données doivent être soumises à Prophet : les colonnes de temps et de mesure quantitative doivent être respectivement nommées “ds” et “y”.

Une granularité plus fine que le jour peut être utilisée dans le champ datetime. C’est alors qu’il sera pertinent d’activer le daily effect qui sera visualisable graphiquement.

La régularité est un élément fondamental pour la modélisation d’une série temporelle, que viennent perturber les années bissextiles. Nous pouvons décider de supprimer les 29 février, par exemple avec le code ci-dessous :

df= df.loc[~(df['ds'].dt.month.eq(2) & df['ds'].dt.day.eq(29))]

Attention, nous dégradons alors inévitablement la régularité des semaines…

Création du modèle et évaluation

Lançons tout d’abord un modèle simple, sur l’intégralité de l’historique. Dans l’extrait de code ci-dessous, df représente un pandas dataframe regroupant les deux colonnes nommées ds et y.

De nombreux paramètres sont définis par défaut : additivité, pas de saisonnalité journalière, détection automatique des change points (changement de tendance), etc.

m = Prophet(daily_seasonality=False)
m.add_country_holidays(country_name='FR')
m.fit(df)

Ici, nous ajoutons au modèle les jours fériés français qui deviendront autant d’indicateurs dans notre modèle additif.

m.train_holiday_names

Comme l’indique cette discussion, il ne semble pas possible d’ajouter plusieurs pays à l’aide de cette méthode.

Nous pouvons également ajouter nos propres dates si nous considérons que des événements (répétitifs sur la saisonnalité attendue) ont un impact sur le phénomène observé. Nous pourrions par exemple identifier les journées les plus froides de l’année qui engendrent certainement une hausse de consommation électrique. Mais serons-nous en capacité de les prédire par la suite ? Mieux vaut se limiter à des événements connus à l’avance tels que des compétitions sportives (lors de la finale de la Coupe du Monde, tout le monde allume la télé !).

Les séries temporelles seraient si simples si toutes les tendances étaient linéaires ! Mais dans la vraie vie (et encore plus en période de pandémie…), les tendances fluctuent et il est indispensable que notre modèle les comprenne. Prophet identifie automatiquement les “change points” et on les visualise ainsi, en rouge, à l’aide du code ci-dessous.

from prophet.plot import add_changepoints_to_plot

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

Nous pouvons aussi les définir arbitrairement. Pourquoi pas en donnant les dates des changements de saison (été et hiver) ?

m = Prophet(daily_seasonality=False,
changepoints=['2013-12-21', '2014-06-21', '2014-12-21', '2015-06-21', '2015-12-21', '2016-06-21', '2016-12-21', '2017-06-21'],
changepoint_prior_scale=1)

Le paramètre changepoint_prior_scale indique à quel point le modèle doit respecter notre liste (ou sa détection automatique). C’est une valeur entre 0 et 1. Avec 1, nos change points sont bien retenus, comme l’atteste ce graphique.

Forecast

Pour calculer une prévision, nous allons avoir besoin d’un nouveau dataframe contenant une colonne de type datetime, toujours nommée “ds”, que nous soumettrons à la méthode .predict(), appliquée au modèle.

future = pd.date_range(start="2020-01-01",end="2021-12-31")
future = pd.DataFrame(future, columns=['ds'])

Tout comme pour la méthode .fit(), la syntaxe est ainsi tout à fait similaire à celle du package scikit-learn.

forecast = m.predict(future)

Une autre méthode pour créer la plage de forecast consiste à utiliser la fonction ci-dessous.

future = m.make_future_dataframe(periods=365)

Cette fonction génère automatiquement une plage de dates couvrant l’historique complet auquel s’ajoute la période définie en paramètre. Ceci a pour avantage de nous permettre de comparer les prévisions avec les données réelles qui ont été utilisées pour le modèle, ce que nous allons observer dans les sorties graphiques.

Sorties graphiques

Pour afficher la superposition des historiques et des prévisions, un simple plot suffit !

m.plot(forecast)

Nous pouvons ensuite obtenir les différents graphiques de décomposition.

m.plot_components(forecast)

D’autres graphiques interactifs sont disponibles conjointement avec la librairie plotly.

Cross validation

Nous allons maintenant évaluer la qualité de la prévision à l’aide des métriques d’évaluation traditionnelles que sont MSE, RMSE, MAE, MAPE, etc. et d’une méthode de validation croisée. Ici encore, tout est intégré dans une fonction ! Analysons le code ci-dessous.

df_cv = cross_validation(m, initial='730 days', period='365 days', horizon='365 days')
df_p = performance_metrics(df_cv, rolling_window=1)

Nous retrouvons le modèle préalablement entraîné (m). Un seul paramètre est obligatoire, celui de l’horizon de prévision. Mais nous pouvons également préciser la période initiale d’entrainement (par défaut, 3 horizons, et en effet avec ces méthodes dites “à court termes”, ne vous aventurez pas à prédire au delà d’un tiers de votre historique !). Enfin, le paramètre period indique la taille des “découpes” faites dans le jeu de données pour établir de nouvelles prévisions (par défaut, un demi horizon). Voici un exemple de sortie de la commande lancée sur un cluster Databricks.

A partir d’un entrainement sur les années 2013 et 2014, la validation croisée a effectué 4 prévisions sur les années 2015 à 2018.


Nous obtenons alors les métriques d’évaluation, en moyenne sur le nombre de “folds“.

L’argument rolling_window indique le pourcentage de données considérées dans la prévision (1 équivaut à 100%).

Nous pouvons enfin rechercher les meilleurs hyperparamètres pour notre modèle.

Utiliser une approche par hyperparameter tuning implique de lancer un nombre important d’entrainements et d’évaluations. C’est ici que nous tirerons profit d’un cluster de machines, en précisant le paramètre parallel=’processes’ dans la fonction cross_validation(). L’exemple de code donné sur le site officiel sera simple à adapter.

Et maintenant, Kats !

La R&D de Facebook ne s’est pas arrêtée là ! En 2021, une nouvelle librairie est mise à disposition : Kats. Celle-ci a pour vocation de simplifier les tâches des Data Scientists autour de l’analyse et de la modélisation des séries temporelles.

Ici, pas de nommage particulier des colonnes, mais nous transformerons le dataframe en objet spécifique à Kats : TimeSeriesData().

Une première fonctionnalité consiste à déterminer automatiquement 65 features de la série temporelle, c’est-à-dire des caractéristiques de cette série (moyenne, variance, entropie, etc.) qui pourront être par la suite intégrées à des modèles de Machine Learning ou à une approche par régression de la modélisation de la série temporelle.

De nombreuses fonctionnalités s’orientent autour de la détection : seasonalities, outlier, change point, and slow trend changes

Enfin, Kats intègre un grand nombre de modèles (ARIMA, HW, stlf… et l’inévitable Prophet !). Bref, c’est le couteau suisse rêvé de tout.e Data Scientist qui s’attaque à un sujet de séries temporelles.

Forecasting par automated ML (UI)

Le portail Azure Machine Learning propose une interface graphique (UI) pour réaliser de premières modèles d’apprentissage dans trois cas de figure :

  • régression
  • classification
  • forecasting

Nous allons nous intéresser ici au forecasting, c’est-à-dire à la prévision sur des données issues d’une série temporelle (ou dite encore série chronologique).

Une série temporelle est un jeu de données composé tout simplement de deux colonnes : une mesure numérique et une variable de temps indiquant le moment de cette mesure. Il est indispensable que les intervalles de temps soient réguliers. Par exemple, nous pouvons disposer d’une mesure par jour, par semaine ou bien par mois, ou encore à des granularités plus fines comme l’heure, la minute, voire la seconde.

S’il existent des valeurs manquantes, il sera nécessaire de les “imputer” avant d’utiliser le jeu de données pour l’apprentissage automatisé.

Depuis le menu latéral, nous lançons une nouvelle exécution (“run“) d’automated ML.

Nous sélectionnons ensuite le jeu de données (“dataset“) répondant à la définition donnée ci-dessous d’une série temporelle.

Pour un premier exemple d’utilisation, nous choisirons le jeu de données du nombre de filles nées par jour en Californie en 1959, disponible sur ce lien.

Il faut ensuite donner un nom pour l’expérience (la trace de l’exécution dans le portail) et choisir une ressource de calcul de type compute cluster.

C’est à l’écran suivant que l’on choisira explicitement une tâche de type “time series forecasting“.

Avant de lancer l’exécution (bouton Finish), nous pouvons réaliser quelques paramétrages (attention, aucun retour arrière ne sera ensuite possible !).

Si les données ne sont pas agrégées pour ne disposer que d’une seule ligne par échelon de temps, il sera nécessaire de lister les “time series identifiers“, c’est-à-dire les colonnes qui permettraient de traiter plusieurs séries temporelles au travers de la même expérience. <Cela se traduira par la présence d’autres paramètres en entrée du service web prédictif.>

Deux options sont par défaut positionnées sur “autodetect” :

  • Frequency : il s’agit de la granularité temporelle du jeu de données. Les valeurs possibles sont présentées ci-dessous. Utilisez explicitement ce paramètre si un doute est possible ou qu’une agrégation est nécessaire (sum, min, max, mean).
  • Forecast horizon : il s’agit du nombre de périodes (dans l’unité du paramètre Frequency) qui seront prédites. La documentation officielle ne mentionne pas la règle appliquée dans le cas de l’auto-détection. Je vous recommande de le spécifier explicitement, en étant “raisonnable”, c’est-à-dire en ne dépassant pas le tiers de votre historique (ex.: pour 3 années d’entrainement, se limiter à une année de prévision).

Nous pouvons enfin définir des éléments de configuration additionnelle (cliquer sur “View additional configuration settings”).

Pour bloquer des algorithmes, il suffit de cliquer dans la liste déroulante.

Je vous conseille de conserver les méthodes dites “naïves” qui donneront un premier seuil pour les métriques d’évaluation. Vous pouvez également choisir de ne pas utiliser les méthodes non linéaires comme les RandomForest ou KNN. Les méthodes traditionnelles du domaine des séries temporelles sont présentes (autoARIMA, exponential smoothing) ainsi que la librairie issue de Facebook : Prophet.

Les méthodes d’apprentissage profond seront également évaluées en cochant la case “enable Deep Learning”. Ces méthodes ont montré leur efficacité mais rappelons qu’elles sont coûteuses en temps d’entrainement puis d’inférence et sur des problématiques relativement simples, elles ne donneront sans doute pas un gain significatif de performance.

Les trois métriques de comparaison des modèles sont :

  • normalized RMSE
  • R2
  • normalized MAE

Nous retrouverons cette métrique comme critère d’arrêt de l’expérience. Il faut alors renseigner un seuil entre 0 et 1 (les métriques sont normalisées et le R2 est par définition compris entre 0 et 1).

Il est important de ne pas se limiter à un seul indicateur d’évaluation, tous seront accessibles quand les différents modèles seront entrainés.

Je vous recommande aussi de borner le temps d’entrainement, afin de limiter les coûts mais également parce que vous obtiendrez sans doute assez rapidement une idée des modèles adaptés à vos données. Il sera ensuite plus efficace de bloquer les algorithmes les moins efficaces.

Si les deux critères sont utilisés, c’est le premier atteint qui arrêtera l’expérience.

Trois paramètres additionnels de configuration sont ensuite disponibles.

Ils correspondent respectivement à :

  • target lags : il s’agit ici d’exploiter des valeurs précédentes pour prédire les nouvelles valeurs. Par exemple, un lag de 1 permettra d’utiliser la valeur à t pour prédire t+1. Je vous recommande cet article pour approfondir ce sujet.
  • rolling window size : permet de ne tenir compte que d’une partie des données d’apprentissage. Par défaut, l’intégralité du jeu de données est considéré. Plus de détails sur cette méthode ici.
  • season and trend : active la méthode de décomposition STL (“Seasonal and Trend decomposition using Loess”), permettant d’isoler saisonnalité, tendance et bruit.

Une dernière case permet de sélectionner un (et un seul) pays pour tenir compte des vacances dans les calculs de saisonnalité, ce que fait par exemple le modèle Prophet.

La seule méthode de validation du modèle est la méthode ce validation croisée “k-fold”, avec k = 5 par défaut.

Enfin, les algorithmes se paralléliseront en fonction du nombre de nœuds du compute cluster défini préalablement.

Nous pouvons maintenant analyser les sorties produites par l’expérience d’automated ML.

Le run principal est constitué de children runs, correspondant chacun à un algorithme, précédé éventuellement d’une préparation de données (PCA, MinMaxScaler, etc.) , et avec un choix d’ hyperparamètres. L’onglet Models présente ces différentes exécutions, triées selon la valeur décroissante de la métrique principale.

Notons également l’onglet Data guardails (garde-fous) qui nous avertit sur d’éventuels problèmes de constitution de notre jeu de données initial (historique trop court, valeurs manquantes, etc.).

Il est très vraisemblable que le modèle “VotingEnsemble” soit le meilleur. C’est en effet un “méta-modèle” qui utilise plusieurs modèles simples et les fait voter pour obtenir les meilleures prévisions. C’est également un modèle plus lourd à exposer.

Il est tout à fait possible de choisir un autre modèle pour l’exposer.

Pour un test simple, nous choisissons une ressource de type Azure Container Instance (ACI) et n’enclenchons pas l’authentification (un jeton serait alors demandé lors de l’appel au service web).

Dans les paramètres avancés, nous pouvons réduire le CPU et la RAM utilisés (et donc la facture associée… mais aussi les performances !).

Un “virtual” CPU peut être défini par un nombre décimal.

Un clic sur “Deploy” lance le déploiement et nous pouvons basculer sur la page dédiée aux endpoints.

La ressource ACI se retrouvera également dans le portail Azure, dans le même groupe de ressources que le service Azure Machine Learning.

Pour consommer le modèle à partir du service web, nous devons utiliser des dates postérieures à la dernière date ayant servi à l’entrainement. Ceci peut se faire au travers de l’interface de test ou à l’aide des exemples de code donnés dans les langages C#, Python et R.

Nous testons ici la journée du 1er janvier de l’année suivante (1960).

En conclusion, nous avons ici un outil dédié aux “citizen data scientists” puisqu’il n’est jamais nécessaire d’écrire de code mais une (très) bonne connaissance de la théorie des différentes méthodes algorithmes et de l’interprétation des métriques d’évaluation sera indispensable. Nous obtenons très rapidement un service web prédictif efficient qui pourra être exploité par d’autres applications. N’oublions pas enfin que pour “vivre en production”, une telle approche ne sera pas suffisante. Nous ne maîtrisons pas ici en particulier le réentrainement du modèle ni l’archivage de ces versions. Une approche par le code, exploitant les librairies du SDK azureml, sera alors une démarche plus pérenne et respectant les bonnes pratiques du MLOps. Prochain article à venir !