Présentation et remerciements
Bonjour à tous et bienvenue sur cet article qui vise à développer avec vous une pipeline pour l’exécution d’un modèle VAR sous Python, de l’intégration de la data jusqu’au plotting des résultats. De formation initiale en management et finance de marché (ESCP BS), j’ai eu le plaisir de réaliser un cursus de formation en Data Science chez DataScientest, qui a aujourd’hui la gentillesse de publier cet article !
Je souhaite tout d’abord remercier l’entreprise « Les Mousquetaires », avec qui j’ai pu développer une version de ce modèle pour réaliser des prédictions sur le marché du gaz naturel. À cet égard, je salue particulièrement Hélène Laroche, responsable de l’achat du gaz, Quentin Pailloux, directeur des énergies nouvelles, Arnaud de Ligniville, Brett Chauchard & Ellias Ouidir de la direction financement et trésorerie, Christelle Schneider, Virginie Hyppolite de la direction des risques.
Pour le reste de l’article, chacune des parties suivantes sera découpée de la manière suivante :
- Explications mathématiques et statistiques de nos choix et contraintes
- Réalisation des fonctions (disponible sur le lien Git)
- Exécution des fonctions et output
Qu’est-ce que le modèle VAR ?
Le modèle VAR pour Vector Autoregression est un modèle de séries temporelles qui permet de prédire plusieurs variables dans le temps. Chacune des prédictions de chaque variable contribue à la prédiction des autres. Ainsi, si on l’écrivait sous un format linéaire, nous aurions ceci :
Soit 2 variables, Y et X, et (t) un moment dans le temps.
Y(t+1) = f(Y(t), X(t))
X(t+1) = f(X(t), Y(t))
Y(t+1) = α1 + β1Y(t) + β2X(t) + ε1
X(t+1) = α2 + β3X(t) + β4Y(t) + ε2
Avec α les constantes, β les coefficients associés aux variables, ε les erreurs
Le modèle calcule sur la base de l’historique les constantes, coefficients pour chaque variable
De plus, le modèle inclut des « lags » (λ), qui permettent d’inclure dans le calcul de X(t+1), Y(t+1) toutes les valeurs de X, Y jusqu’à X(t- λ), Y(t- λ)
On retrouve bien l’explication du nom : un vecteur auto-regressif.
Quel intérêt au modèle VAR ?
Le modèle VAR peut-être très intéressant pour concurrencer des modèles de machine learning tels que les RandomForestRegressor, le XGBoostRegressor, lorsque votre variable est une série temporelle. Les bases de sa littérature sont solides, bien que certains aspects soient encore sujet à discussion (nous y reviendrons).
Avantage de ce modèle : il est très facilement explicable auprès des métiers, car il s’appuie sur des bases statistiques aisées à expliquer et à comprendre.
Le modèle VAR peut être utilisé pour prédire un ensemble de variables, mais on peut aussi l’utiliser pour se focaliser sur une seule variable, ce qui est notre cas ici.
Définition du problème et de la fenêtre de prédiction
Avant de prédire, il nous faut un problème, et des variables !
Il m’est impossible de dévoiler mon Dataset, mais je vais vous guider sur la méthodologie de choix des variables pour que vous puissiez l’appliquer dans tous les cas.
Supposons que votre métier vous demande un modèle permettant de prédire un nombre (type « regression » en data science donc).
-> Cette variable est-elle dépendante ou s’observe-t-elle en fonction du temps ?
bh
-> Cette variable est-elle a priori liée à d’autres variables ?
mp,g
-> Ces autres variables sont-elles aussi observables au cours du temps ?
Si vous avez répondu oui à ces 3 questions, alors le modèle VAR pourrait bien correspondre à votre besoin.
Comment choisir la fenêtre de prédiction ?
Ceci est une question épineuse mais très intéressante : combien de jours dois-je (et puis-je !) prédire ? Le métier peut avoir un besoin précis comme 30 jours, 120, 360… Mais peut-être ne disposez vous pas de la volumétrie de data nécessaire ?
Ce que je vous recommande à ce sujet : démarrez avec la demande du métier, et réduisez la voilure si les tests statistiques ne vous permettent pas d’aller aussi loin. Inversement si tout fonctionne pour la durée que vous demande votre métier, pourquoi ne pas aller plus loin dans le temps ? De manière générale, je ne peux que recommander de tester votre pipeline sur différentes cibles de temps.
Les variables explicatives :
Premièrement, tentez d’identifier un maximum de variables qui pourraient avoir un impact sur votre variable cible.
Plus vous en aurez, plus vous aurez de chance de découvrir celles qui influencent votre variable cible. Nous raisonnons ici dans un mode design thinking : quality through quantity. Ne vous inquiétez pas, nous sélectionnerons ensuite les variables les plus intéressantes à retenir pour le modèle final.
Assurez-vous que vos datas possèdent la même unité temporelle : si votre variable cible est exprimée en unité journalière, faites attention à ce que le reste le soit aussi !
La gestion des NaN
En cas de quelques manques dans votre data set, vous pouvez évidemment vous appuyer sur les méthodes fillna() de Python. Pour une approche peut-être plus solide notamment si le nombre de NaN commence à être trop élevé, je vous recommande un peu de lecture et le modèle MICE.
Dans le cas où vous auriez de la data sur certaines séries pendant une longue période de temps, mais pas sur d’autres, je vous recommande de réfléchir à soit carrément les drop, soit à tronquer votre data set pour avoir, dès votre première ligne, toutes vos variables complètes.
Voici quelques fonctions pour intégrer la data, concaténée en un seul dataframe : def agreg_df (voir le lien GIT)
Et une fonction de remplacement : def replace_na (voir le lien GIT)
Et on exécute ensuite le code :
Besoins statistiques et transformations adéquates
Nos séries temporelles doivent répondre à un critère fondamental : la stationnarité.
Quel que soit le moment (t) dans le temps, l’espérance de notre variable est toujours la même. Il y a absence de tendance et sa variance est toujours la même et non infinie, et l’auto-corrélation entre deux moments n’est fonction que du décalage dans le temps
Je recommande ce feed sur researchgate pour plus de détails sur la nécessité de la stationnarité ça se passe ici.
Maintenant que le problème est posé, que faisons-nous ?
La transformation en utilisant Box-Cox
Il peut être très utile de transformer nos datas avant de tester leur stationnarité. En effet, les transformations ont de multiples avantages, comme le fait de supprimer (ou diminuer) l’aspect saisonnier, la trend. Une transformation extrêmement pratique est celle de Box-Cox : en fonction d’un certain paramètre Lambda compris entre [-1 ;1], la série sera transformée de différente manière.Voici la documentation de la fonction sous Python.
Ce qui va nous intéresser en particulier est l’argument Lambda et sa valeur « None » : dans ce cas, la fonction choisira d’elle-même le paramètre optimal ! Dans ce cas, le deuxième élément retourné par la fonction sera le « lambda » associé à la série transformée, ce qui va nous intéresser afin de dé-transformer la série.
Le test de stationnarité et la différenciation
Le test de Dickey-Fuller augmenté permet de tester la stationnarité d’une série temporelle.
L’hypothèse H0 du test est la suivante : « La data présente une racine unitaire, la série est donc non-stationnaire ».
Nous allons donc nous attacher à rejeter cette hypothèse, avec un niveau de confiance de 95%. Cela signifie que le test doit nous renvoyer un p-valeur inférieure à 0.05. Voici la documentation de la fonction.
Je vous propose la fonction suivante pour une mise en forme du test
def dickey_fuller_test (voir le lien GIT)
Dans le cas où au moins une des variables n’est pas stationnaire, nous allons procéder à une différenciation :
Data_diff(t) = observation(t) – observation(t-1)
Avec (t) un moment dans le temps.
Différencier est une forme de transformation, et cela marche très bien pour arriver à la stationnarité des séries.
Il existe une méthode Pandas sur les dataframe pour différencier : df.diff()
Dans le cas où une première différenciation ne suffirait pas pour atteindre la stationnarité sur nos séries, on peut différencier une deuxième fois. Si cela ne suffit toujours pas, il vaut mieux abandonner la variable ou chercher d’autres méthodes de stationarisation
Passons maintenant à l’exécution du code pour cette partie !
Il se peut que vous ayez des « 0 » dans vos datas : il faut alors ajouter ‘1’ à l’ensemble du DF pour que la transformation fonctionne :
Problème, toutes les variables sont non-stationnaires…
Deux de mes variables sont toujours non-stationnaires, alors je vais redifférencier et retester :
Cette fois-ci, tout est ok ! Mes séries sont stationnarisées.
CEPENDANT, attention : il faut noter qu’une seconde différenciation fait perdre de l’information à nos variables.
UNE CARRIÈRE DANS LA DATA VOUS TEND LES BRAS !
UNE CARRIÈRE DANS LA DATA VOUS TEND LES BRAS !
Une reconversion dans le big data vous intéresse, mais vous ne savez pas par où commencer ?
Ainsi, si jamais toutes vos séries ne sont pas stationnaires après une première différenciation, il faudra optimiser et comparer deux familles de modèles :
Une première famille sur la data différenciée UNE fois, SANS les variables non stationnaires
Une seconde famille sur la data différenciée DEUX fois, AVEC les variables qui sont devenues stationnaire.
Dans mon cas, je sais que le modèle différencié deux fois perd trop de signal par rapport à l’ajout des deux nouvelles variables stationnaires. Je vais donc travailler sur la data différenciée une fois, sans les deux variables qui ne le sont pas. On va donc les drop, puis créer le train et test set
Sélections des variables et du maxlags : itération par le test de Granger
Toutes les variables que nous avons sélectionnées ne sont pas forcément à retenir pour notre modèle. La question à se poser est la suivante : toutes les variables participent-elles efficacement à la prédiction des autres, et notamment de la variable cible ?
Un test permet d’y répondre : le test de Granger. Il permet de tester si les prédictions de la variables X influencent significativement les prédictions de la variable Y. L’hypothèse HO du test est : « Les prédictions de la variable X n’influencent pas les prédictions de la variable Y ». On souhaite donc, au niveau de confiance 95%, des p-valeur inférieures à 0.05. La page wikipédia du test l’explique très bien de mon point de vue. Et voici la documentation du test sous python.
Ce qui est particulièrement intéressant et à noter ici, est l’argument « maxlag » : il va déterminer combien de lag maximum notre test va prendre. Et c’est ici que nous allons devoir réaliser des tests itératifs, car le maxlag change les p-valeurs, potentiellement de manière drastique : des variables qui ne passent pas le test à maxlag = 10 peuvent complètement le passer au maxlag 20.
Je vous propose le code suivant pour une fonction du test de granger :
def grangers_test (voir le lien GIT)
Comment lire le DF qui nous est renvoyé ?
Les variables en colonne (les X) sont celles dont le poids des prédictions est testé sur les prédictions des variables Y (nous allons voir un exemple plus bas). Si la valeur est inférieure à 0.05, alors en effet on rejette H0 et donc les prédictions de X influencent les prédictions de Y.
Notre objectif est d’avoir un tableau où toutes les valeurs (sauf la diagonale) sont inférieures à 0.05 : cela signifie que tous les couples de variables rejettent l’hypothèse H0.
Mais avant de tester nos variables avec cette fonction, il faut déterminer le maxlag.
Pour cela, on va entraîner notre modèle VAR avec un maxlag choisi au hasard. Si l’entraînement fonctionne, on augmente le maxlag, jusqu’à ce que le code plante : un maxlag trop élevé entraîne un problème sur les équations. Une fois le maxlag déterminé, on passe à la sélection des variables.
Pour cela, nous allons procéder à des itérations de la démarche suivante :
- On teste avec notre fonction et le maxlag, et on affiche le DF
- On regarde la colonne OU la ligne qui contient le plus de valeurs >0.05. Dans le cas des colonnes cela signifie : ma variable X détermine trop peu de Y. Dans le cas des lignes : ma variable Y est trop peu déterminée par les X. Si on arrive à trouver une variable qui a la fois en ligne et en colonne contient beaucoup de 0.05, c’est parfait, puisque l’on supprimera d’un coup une variable qui explique peu les autres et qui est peu expliquée.
- On supprime la variable choisie de notre DF actuel, mais aussi de celui de base, car il va nous servir plus tard pour les transformations inverses, il faut les mêmes colonnes.
- On entraîne notre modèle pour trouver le nouveau « maxlag » qui ne fasse pas planter le modèle
- Retour à l’étape 1 jusqu’à ce que le DF ne contienne que des valeurs <0.05
Voici le code associé à ce passage, avec quelques itérations et des screenshots pour illustrer :
Voici un screenshot de ce que j’obtiens (attention, cette opération est la plus longue de notre modélisation).
Pour rappel, notre variable 1 est la variable cible du projet, on doit la garder.
En un coup d’œil, je pense que la variable 8 est un bon candidat : c’est là où le nombre de p-valeur >0,05 est le plus grand.
On la drop, on modifie les set :
On va chercher un nouveau maxlag optimal (le plus élevé possible), en refaisant tourner le bloc de code précédent (maxlag, model_var, res).
Dans mon cas, il passe à 19. Je reconduis le test de granger :
Maintenant, c’est la variable 6 qui est le meilleur candidat.
Je vais la drop, et réitérer les étapes précédentes jusqu’à ce que le DF ressemble à ça (aucune p-valeur >0.05) :
Itérer à travers le maxlag pour trouver le meilleur modèle
Une fois le df_granger nettoyé, on récupère notre maxlag le plus élevé (le dernier donc). Dans mon cas, maxlag = 67.
Nous allons créer un modèle par maxlag, et récupérer à chaque fois les éléments suivants :
La MAE (ou votre metric cible), la p-valeur du test de normalité des résidus, la p-valeur du test de blancheur des résidus.
Quelques spécifications concernant la normalité et la blancheur :
Des résidus qui réussissent le test de normalité signifient que les erreurs suivent une loi normale. Cela est une condition nécessaire afin de créer les intervalles de confiance classiques. Dans le cas où les résidus ne passent pas le test, il faut alors procéder par bootstraping, ce qui sera le sujet d’un futur article !
Un test de blancheur des résidus réussi signifie que les résidus ont la forme d’un bruit blanc, ce qui signifie qu’il n’y a plus (ou à la marge) d’information à tirer de nos variables. On peut toujours améliorer un modèle qui ne passe pas le test de blancheur des résidus, mais ce n’est pas parce qu’un modèle le passe qu’il n’est pas améliorable !
On peut améliorer un modèle aisément en ajoutant à toutes les prédictions la moyenne des biais, ce qui est réalisé dans notre pipeline.
Les points précédents sont une presque paraphrase de l’excellent travail de Rob J Hyndman et George Athanasopoulos, dans « Forecasting, principles and practice » que vous pouvez trouver gratuitement ici et que je vous recommande chaudement : https://otexts.com/fpp2/
Je vous propose la fonction suivante pour notre pipeline concrète :
def pipeline_var (voir le lien GIT)
Les plus attentifs d’entre vous auront noté que l’on utilise la fonction « inv_diff » à l’intérieur de « pipeline_var« , disponible aussi sur le git. Attention au paramètre « second_diff » ! False si on utilise de la data différenciée une fois, True sinon
On va maintenant exécuter le tout et stocker nos MAE et p-valeur des tests. Je print « i » pour un suivi de l’exécution. N’oubliez pas de tester différentes valeurs pour les paramètres « ic » et « trend » dans le fitting du modèle !
Maintenant que l’on a fait tourner tous nos modèles (de 1 jusqu’au nombre de maxlag, 67 pour moi), on peut regarder la meilleure MAE et le lag associé dans « array_mae_var ». Plusieurs choses à noter :
Si durant l’itération le code est en échec, c’est qu’à partir de ce lag, le modèle prédit des données extrêmes. Cela signifie que l’on n’ira pas plus haut en termes de lag
Le meilleur lag est : index de la plus faible MAE + 1 (nous sommes sur Python !)
Dans mon cas, l’index qui renvoie la meilleure MAE est de 43, avec une mae de 1.44. Donc, le meilleur lag est de 44
Or, 44 < 67, qui était mon maxlag.
Vérifier le test de granger pour le nouveau lag optimal
Comme ma meilleure MAE est obtenu pour un maxlag inférieur à 67, c’est-à-dire le lag le plus élevé avec lequel j’ai conduit mes tests de Granger, je dois re-tester ma matrice de Granger avec un lag de 44. En effet, pour rappel, les P-values du test renvoyées dépendent du lag, et ont tendance à augmenter en même temps que le maxlag diminue… On va donc retester notre df avec maxlag = 44, et voici ce que j’obtiens :
Aucune de mes P-Valeur ne met le test en échec, je peux donc conclure sereinement que mon meilleur modèle est obtenu avec lag = 44
Résultats, graphiques, métriques
J’obtiens donc sur 120 jours une prédiction avec une MAE de 1.44. Mon test de blancheur ne passe pas, ce qui signifie qu’il y a de l’information restante dans mes variables. Le test de normalité ne passe pas non plus, je ne peux donc pas utiliser les intervalles de confiance classique.
Merci de m’avoir lu jusqu’ici ! J’espère que cet article pourra vous aider. N’hésitez pas à commenter ou à me contacter pour une éventuelle question 😊 LinkedIn Perso.
Si tout comme Victor vous souhaitez vous former et maitriser les données, nous sommes là pour ça !
