Time Series Forecasting with Machine Learning
Josue Afouda
2023-07-16
# Importation des librairies nécessaires
library(TSstudio)
library(dplyr)
library(lubridate)
library(psych)
library(rpart)
library(rpart.plot)Dans cet article, nous nous concentrerons sur les prévisions des ventes mensuelles de véhicules aux États-Unis en utilisant les modèles suivants :
Régression linéaire ;
Arbre de décision.
Données
# Chargement des données
data(USVSales)
ts_info(USVSales)## The USVSales series is a ts object with 1 variable and 528 observations
## Frequency: 12
## Start time: 1976 1
## End time: 2019 12Analyse exploratoire
Visulaisation de la série originelle
La série des ventes mensuelles totales de véhicules aux États-Unis, comme on peut le voir dans le graphique suivant, est un exemple de série avec des modèles à la fois saisonniers et cycliques :
# Visualisation de la série
ts_plot(USVSales,
title = "US Total Monthly Vehicle Sales",
Ytitle = "Thousands of Units",
Xtitle = "Year")Comme vous pouvez le voir dans le graphique précédent, la série a des modèles de cycle, ce qui est courant pour un indicateur macroéconomique. Dans ce cas, il s’agit d’un indicateur macro de l’économie américaine.
Décomposition de la série
# Décomposition de la série
ts_decompose(USVSales)La principale caractéristique de la série USgas est la
forte tendance saisonnière.
# Analyse de la saisonnalité
USVSales_detrend <- USVSales - decompose(USVSales)$trend
ts_seasonal(USVSales_detrend, type = "box")Nous pouvons voir sur le graphique saisonnier précédent que, généralement, le pic de l’année s’est produit pendant les mois de mars, mai et juin. De plus, vous pouvez voir que les ventes diminuent à partir des mois d’été et culminent à nouveau en décembre pendant les périodes de vacances. D’autre part, le mois de janvier est généralement le mois le plus bas de l’année en termes de ventes.
Analyse de la corrélation entre une série temporelle et plusieurs de ses lags
Tout d’abord, qu’appelle t-on lag d’une série temporelle ?
Dans le domaine des séries temporelles, le terme “lag” fait référence à la différence entre une observation à un certain moment dans le temps et une observation précédente à un moment antérieur. Il mesure le décalage ou la retardation d’une variable par rapport à elle-même dans le temps.
Plus précisément, un “lag” est le nombre d’intervalles de temps entre deux observations successives dans une série temporelle. Par exemple, pour une série chronologique mensuelle, un lag de 1 représente le décalage d’une observation d’un mois par rapport à l’observation précédente, tandis qu’un lag de 12 correspond à un décalage d’une année.
Le concept de “lag” est largement utilisé dans l’analyse et la modélisation des séries temporelles. Il est utilisé pour capturer les dépendances et les structures temporelles dans les données. Par exemple, en utilisant les lags d’une variable, on peut construire des modèles autorégressifs (AR) qui prennent en compte les valeurs précédentes de la série pour prédire les valeurs futures.
Les lags peuvent être utilisés pour calculer les autocorrélations, les fonctions d’autocorrélation partielle (PACF) et d’autres métriques pour identifier les relations de dépendance dans une série temporelle. Ils peuvent également être utilisés pour créer des variables retardées (lagged variables) qui servent d’entrées aux modèles de prévision ou de régression temporelle.
En résumé, dans le contexte des séries temporelles, un “lag” mesure la différence de temps entre deux observations successives et joue un rôle important dans l’analyse, la modélisation et la prévision des séries chronologiques.
Passons maintenant à l’analyse :
En raison de la nature continue et chronologique des données des séries temporelles, il est probable qu’il y aura un certain degré de corrélation entre les observations des séries. Par exemple, la température de l’heure suivante n’est pas un événement aléatoire puisque, dans la plupart des cas, elle a une forte relation avec la température actuelle ou les températures qui se sont produites au cours des dernières 24 heures. Dans de nombreux cas, la série d’observations passées contient des informations prédictives sur des événements futurs, qui peuvent être utilisées pour prévoir les observations futures de la série.
Le but de l’analyse des lags est d’identifier et de quantifier la relation entre une série et ses retards. Cette relation est généralement mesurée en calculant la corrélation entre les deux et en utilisant des outils de visualisation de données. Le niveau de corrélation entre une série et ses retards est dérivé des caractéristiques de la série.
ts_lags(USgas) En regardant les tracés de décalage de la série USgas, vous pouvez voir que, du premier décalage au sixième décalage, la relation entre la série et ses décalages devient moins linéaire. Ce processus commence à s’inverser à partir du septième décalage à mesure que la relation devient progressivement plus linéaire, où le décalage saisonnier (ou décalage 12) a la relation la plus forte avec la série.
Par défaut, la fonction ts_lags() trace les 12 premiers
décalages de la série, où l’argument lags vous permet de
définir le nombre de décalages.
Par exemple, nous pouvons tracer les décalages saisonniers les plus récents (c’est-à-dire les décalages 12, 24, 36 et 48) en définissant le nombre de décalages avec l’argument décalages sur les décalages correspondants :
ts_lags(USgas, lags = c(12, 24, 36, 48))La relation de la série avec les premier et deuxième décalages saisonniers a une forte relation linéaire, comme le montre le graphique des décalages précédent.
Nous pouvons conclure notre courte analyse exploratoire de la série USVSales par les observations suivantes :
La série USVSales est une série mensuelle avec une saisonnalité mensuelle claire ;
La tendance de la série a une forme cyclique, et donc la série a une composante de cycle intégrée dans la tendance ;
Le cycle le plus récent de la série commence juste après la fin de la crise économique de 2008, entre 2009 et 2010 ;
Il semble que le cycle actuel ait atteint son apogée alors que la tendance commence à s’aplatir ;
La série a une forte corrélation avec son premier décalage saisonnier.
Passons à présent à la phase de preprocessing et Feature Engineering afin de mettre les données dans un format convenable aux algorithmes de Machine Learning.
Feature Engineering et Préparation des données
Comme nous avons l’intention d’avoir une prévision à court terme (de 12 mois), il ne sert à rien d’utiliser la série complète, car elle peut introduire du bruit dans le modèle en raison du changement de direction de la tendance tous les deux ans. Par conséquent, nous utiliserons les observations de formation du modèle à partir de 2010 et au-delà.
# Mettre l'objet ts sous format dataframe
df <- ts_to_prophet(window(USVSales, start = c(2010,1)))
names(df) <- c("date", "y")
class(df)## [1] "data.frame"head(df)## date y
## 1 2010-01-01 712.469
## 2 2010-02-01 793.362
## 3 2010-03-01 1083.953
## 4 2010-04-01 997.334
## 5 2010-05-01 1117.570
## 6 2010-06-01 1000.455ts_plot(df,
title = "US Total Monthly Vehicle Sales (Subset)",
Ytitle = "Thousands of Units",
Xtitle = "Year")L’ingénierie des fonctionnalités (Feature Engineering) joue un rôle central lors de la modélisation avec des algorithmes ML. Notre prochaine étape, basée sur les observations précédentes, consiste à créer de nouvelles fonctionnalités qui peuvent être utilisées comme entrée informative pour le modèle. Dans le contexte de la prévision de séries chronologiques, voici quelques exemples de nouvelles fonctionnalités possibles qui peuvent être créées à partir de la série elle-même :
La tendance de la série : elle utilise un index numérique. De plus, comme la tendance de la série n’est pas linéaire, nous utiliserons un deuxième polynôme de l’indice pour capturer la courbure globale de la tendance de la série.
Composant saisonnier : cela crée une variable catégorielle pour le mois de l’année afin de capturer la saisonnalité de la série.
Corrélation des séries : utilise la forte corrélation de la série avec son décalage saisonnier et utilise le décalage saisonnier (lag12) comme entrée du modèle.
# Création de nouvelles colonnes
df_eng <- df %>%
mutate(
#month = factor(month(date, label = TRUE), ordered = FALSE),
lag12 = lag(y, n = 12),
month = as.integer(format(date, "%m")), # Extraction du mois
Weekday = format(date, "%A"), # Jour de la semaine
IsWeekend = Weekday %in% c("samedi", "dimanche")
) %>%
filter(!is.na(lag12))
head(df_eng)## date y lag12 month Weekday IsWeekend
## 1 2011-01-01 836.366 712.469 1 samedi TRUE
## 2 2011-02-01 1007.082 793.362 2 mardi FALSE
## 3 2011-03-01 1276.843 1083.953 3 mardi FALSE
## 4 2011-04-01 1173.520 997.334 4 vendredi FALSE
## 5 2011-05-01 1081.272 1117.570 5 dimanche TRUE
## 6 2011-06-01 1071.229 1000.455 6 mercredi FALSE# Structure de df_eng
glimpse(df_eng)## Rows: 108
## Columns: 6
## $ date <date> 2011-01-01, 2011-02-01, 2011-03-01, 2011-04-01, 2011-05-01,…
## $ y <dbl> 836.366, 1007.082, 1276.843, 1173.520, 1081.272, 1071.229, 1…
## $ lag12 <dbl> 712.469, 793.362, 1083.953, 997.334, 1117.570, 1000.455, 106…
## $ month <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6, 7, …
## $ Weekday <chr> "samedi", "mardi", "mardi", "vendredi", "dimanche", "mercred…
## $ IsWeekend <lgl> TRUE, FALSE, FALSE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE,…# Données d'entraînement et de test
h <- 12
train_df <- df_eng[1:(nrow(df_eng) - h), ]
test_df <- df_eng[(nrow(df_eng) - h + 1):nrow(df_eng), ]# Matrice de corrélation
cor(train_df[c("y", "lag12", "month", "IsWeekend")])## y lag12 month IsWeekend
## y 1.0000000 0.93254721 0.22589347 -0.14118169
## lag12 0.9325472 1.00000000 0.19748978 0.01906031
## month 0.2258935 0.19748978 1.00000000 0.01642957
## IsWeekend -0.1411817 0.01906031 0.01642957 1.00000000# Visualisation des corrélations
pairs.panels(train_df[c("y", "lag12", "month", "IsWeekend")])
Nous devons aussi créer un dataframe avec les dates des 12 mois suivants (dates de prévision) et construire le reste des fonctionnalités :
forecast_df <- data.frame(
date = seq.Date(from = max(df$date) + month(1), length.out = h, by = "month")
)
forecast_df$lag12 <- tail(df$y, 12)
forecast_df <- forecast_df %>%
mutate(
month = as.integer(format(date, "%m")),
Weekday = format(date, "%A"),
IsWeekend = Weekday %in% c("samedi", "dimanche")
)
forecast_df## date lag12 month Weekday IsWeekend
## 1 2019-12-02 1171.503 12 lundi FALSE
## 2 2020-01-02 1288.278 1 jeudi FALSE
## 3 2020-02-02 1642.750 2 dimanche TRUE
## 4 2020-03-02 1372.659 3 lundi FALSE
## 5 2020-04-02 1628.074 4 jeudi FALSE
## 6 2020-05-02 1554.748 5 samedi TRUE
## 7 2020-06-02 1443.947 6 mardi FALSE
## 8 2020-07-02 1685.342 7 jeudi FALSE
## 9 2020-08-02 1315.632 8 dimanche TRUE
## 10 2020-09-02 1380.180 9 mercredi FALSE
## 11 2020-10-02 1446.483 10 vendredi FALSE
## 12 2020-11-02 1565.023 11 lundi FALSEModélisation
Modèle de base
La performance d’un modèle de prévision doit être mesurée par le taux d’erreur, principalement sur la partition de test, mais aussi sur la partition d’apprentissage. Vous devez évaluer les performances du modèle par rapport à un modèle de base. Puisque nous utilisons une famille de modèles de régression ML, il est plus logique d’utiliser un modèle de régression comme référence pour les modèles d’Apprentissage Automatique. Par conséquent, nous formerons un modèle de régression linéaire de séries temporelles.
En utilisant les partitions d’entraînement et de test que nous avons créées précédemment, entraînons le modèle de régression linéaire et évaluons ses performances avec les partitions de test :
# Régression linéaire multiple
lr <- lm(
y ~ lag12 + month + Weekday + IsWeekend,
data = train_df
)
# Résultats du modèle lr
summary(lr)##
## Call:
## lm(formula = y ~ lag12 + month + Weekday + IsWeekend, data = train_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -121.910 -48.509 1.453 46.949 132.512
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 246.16457 40.68845 6.050 3.55e-08 ***
## lag12 0.80358 0.02925 27.470 < 2e-16 ***
## month 1.87748 1.91644 0.980 0.329964
## Weekdayjeudi 92.58316 23.33064 3.968 0.000148 ***
## Weekdaylundi 62.21509 24.22913 2.568 0.011942 *
## Weekdaymardi 112.54238 23.27660 4.835 5.69e-06 ***
## Weekdaymercredi 60.43498 23.62451 2.558 0.012254 *
## Weekdaysamedi 27.99248 23.36427 1.198 0.234136
## Weekdayvendredi 71.52795 23.21831 3.081 0.002765 **
## IsWeekendTRUE NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 62.32 on 87 degrees of freedom
## Multiple R-squared: 0.9058, Adjusted R-squared: 0.8971
## F-statistic: 104.6 on 8 and 87 DF, p-value: < 2.2e-16# Prédictions sur les données de test
test_df$yhat <- predict(lr, newdata = test_df)
head(test_df)## date y lag12 month Weekday IsWeekend yhat
## 97 2019-01-01 1171.503 1181.715 1 mardi FALSE 1310.184
## 98 2019-02-01 1288.278 1328.140 2 vendredi FALSE 1388.711
## 99 2019-03-01 1642.750 1687.609 3 vendredi FALSE 1679.450
## 100 2019-04-01 1372.659 1391.226 4 lundi FALSE 1433.848
## 101 2019-05-01 1628.074 1626.484 5 mercredi FALSE 1622.993
## 102 2019-06-01 1554.748 1586.664 6 samedi TRUE 1560.430ts_plot(
test_df[c("date", "y", "yhat")],
title = "Total Vehicle Sales - Actual vs. Prediction (Linear Regression)"
)# Evaluation du modèle lr
mape_lr <- mean(abs(test_df$y - test_df$yhat) / test_df$y)
mape_lr## [1] 0.04828065Nous utiliserons cette valeur pour comparer les performances des futurs modèles de Machine Learning que nous allons construits.
Modèle d’arbre de décision
# Arbre de décision
rpart_model <- rpart(
y ~ .,
data = train_df %>% select(-date)
)
print(rpart_model)## n= 96
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 96 3587259.00 1370.152
## 2) lag12< 1251.997 38 509047.90 1169.959
## 4) lag12< 1068.514 16 125666.30 1069.467 *
## 5) lag12>=1068.514 22 104291.70 1243.044
## 10) Weekday=dimanche,samedi 8 24570.15 1185.924 *
## 11) Weekday=jeudi,lundi,mardi,mercredi,vendredi 14 38703.48 1275.685 *
## 3) lag12>=1251.997 58 557491.40 1501.313
## 6) lag12< 1450.765 26 81433.11 1422.361 *
## 7) lag12>=1450.765 32 182303.90 1565.462
## 14) lag12< 1571.188 22 96352.48 1537.631 *
## 15) lag12>=1571.188 10 31420.75 1626.691 *print(summary(rpart_model))## Call:
## rpart(formula = y ~ ., data = train_df %>% select(-date))
## n= 96
##
## CP nsplit rel error xerror xstd
## 1 0.70268686 0 1.0000000 1.0110190 0.12294400
## 2 0.08188827 1 0.2973131 0.3172365 0.04455734
## 3 0.07780031 2 0.2154249 0.2663600 0.04227104
## 4 0.01520119 3 0.1376246 0.1902166 0.02915301
## 5 0.01143438 4 0.1224234 0.1926968 0.02974052
## 6 0.01000000 5 0.1109890 0.1808070 0.02942341
##
## Variable importance
## lag12 month Weekday IsWeekend
## 77 19 3 1
##
## Node number 1: 96 observations, complexity param=0.7026869
## mean=1370.152, MSE=37367.29
## left son=2 (38 obs) right son=3 (58 obs)
## Primary splits:
## lag12 < 1251.997 to the left, improve=0.70268690, (0 missing)
## month < 2.5 to the left, improve=0.23668160, (0 missing)
## Weekday splits as LRLLLLR, improve=0.02732364, (0 missing)
## IsWeekend < 0.5 to the right, improve=0.01993227, (0 missing)
## Surrogate splits:
## month < 2.5 to the left, agree=0.708, adj=0.263, (0 split)
##
## Node number 2: 38 observations, complexity param=0.07780031
## mean=1169.959, MSE=13396
## left son=4 (16 obs) right son=5 (22 obs)
## Primary splits:
## lag12 < 1068.514 to the left, improve=0.54825860, (0 missing)
## month < 1.5 to the left, improve=0.18136110, (0 missing)
## IsWeekend < 0.5 to the right, improve=0.07864444, (0 missing)
## Weekday splits as LRRRRLR, improve=0.07864444, (0 missing)
## Surrogate splits:
## Weekday splits as RLRRLRR, agree=0.658, adj=0.188, (0 split)
## month < 1.5 to the left, agree=0.632, adj=0.125, (0 split)
##
## Node number 3: 58 observations, complexity param=0.08188827
## mean=1501.313, MSE=9611.921
## left son=6 (26 obs) right son=7 (32 obs)
## Primary splits:
## lag12 < 1450.766 to the left, improve=0.52692190, (0 missing)
## Weekday splits as LRLRLLR, improve=0.11638900, (0 missing)
## IsWeekend < 0.5 to the right, improve=0.09303089, (0 missing)
## month < 11.5 to the left, improve=0.06792437, (0 missing)
## Surrogate splits:
## month < 9.5 to the right, agree=0.655, adj=0.231, (0 split)
## Weekday splits as LRRRRRR, agree=0.621, adj=0.154, (0 split)
##
## Node number 4: 16 observations
## mean=1069.467, MSE=7854.143
##
## Node number 5: 22 observations, complexity param=0.01143438
## mean=1243.044, MSE=4740.533
## left son=10 (8 obs) right son=11 (14 obs)
## Primary splits:
## Weekday splits as LRRRRLR, improve=0.39330160, (0 missing)
## IsWeekend < 0.5 to the right, improve=0.39330160, (0 missing)
## month < 3.5 to the left, improve=0.05488139, (0 missing)
## lag12 < 1172.664 to the right, improve=0.02887205, (0 missing)
## Surrogate splits:
## IsWeekend < 0.5 to the right, agree=1.000, adj=1.00, (0 split)
## lag12 < 1183.343 to the right, agree=0.727, adj=0.25, (0 split)
##
## Node number 6: 26 observations
## mean=1422.361, MSE=3132.043
##
## Node number 7: 32 observations, complexity param=0.01520119
## mean=1565.462, MSE=5696.996
## left son=14 (22 obs) right son=15 (10 obs)
## Primary splits:
## lag12 < 1571.188 to the left, improve=0.29911940, (0 missing)
## Weekday splits as RRLRLLR, improve=0.18740060, (0 missing)
## IsWeekend < 0.5 to the right, improve=0.11432110, (0 missing)
## month < 5.5 to the right, improve=0.04601182, (0 missing)
## Surrogate splits:
## month < 11 to the left, agree=0.75, adj=0.2, (0 split)
##
## Node number 10: 8 observations
## mean=1185.924, MSE=3071.269
##
## Node number 11: 14 observations
## mean=1275.685, MSE=2764.534
##
## Node number 14: 22 observations
## mean=1537.631, MSE=4379.658
##
## Node number 15: 10 observations
## mean=1626.691, MSE=3142.075
##
## n= 96
##
## node), split, n, deviance, yval
## * denotes terminal node
##
## 1) root 96 3587259.00 1370.152
## 2) lag12< 1251.997 38 509047.90 1169.959
## 4) lag12< 1068.514 16 125666.30 1069.467 *
## 5) lag12>=1068.514 22 104291.70 1243.044
## 10) Weekday=dimanche,samedi 8 24570.15 1185.924 *
## 11) Weekday=jeudi,lundi,mardi,mercredi,vendredi 14 38703.48 1275.685 *
## 3) lag12>=1251.997 58 557491.40 1501.313
## 6) lag12< 1450.765 26 81433.11 1422.361 *
## 7) lag12>=1450.765 32 182303.90 1565.462
## 14) lag12< 1571.188 22 96352.48 1537.631 *
## 15) lag12>=1571.188 10 31420.75 1626.691 *# Visualisation du modèle d'abre de décision
rpart.plot(rpart_model, digits = 3)
# Evaluation de la performance du modèle
test_df$y_preds_tree <- predict(
rpart_model,
test_df %>% select(-date)
)
test_df$y_preds_tree## [1] 1275.685 1422.361 1626.691 1422.361 1626.691 1626.691 1422.361 1537.631
## [9] 1537.631 1422.361 1422.361 1626.691On peut mesurer la corrélation entre les valeurs prédictes et les valeurs réelles de l’ensemble de test :
cor(test_df$y_preds_tree, test_df$y)## [1] 0.8168067# Caclul du MAPE
MAPE <- function(actual, predicted) {
mean(abs(actual - predicted) / actual)
}
mape_tree <- MAPE(test_df$y, test_df$y_preds_tree)
mape_tree## [1] 0.05367428En comparant cette erreur avec celle du modèle de base, on conclut que le modèle de base est plus performant que le modèle d’arbre de décision.
# Visualisation des prédictions
ts_plot(
test_df[c("date", "y", "y_preds_tree")],
title = "Total Vehicle Sales - Actual vs. Prediction (Decision Tree)"
)Prévisions
# Prévisions sur un an
futures_ventes <- predict(
lr,
forecast_df %>% select(-date)
)
forecast_df$Previsions <- futures_ventes
forecast_df## date lag12 month Weekday IsWeekend Previsions
## 1 2019-12-02 1171.503 12 lundi FALSE 1272.303
## 2 2020-01-02 1288.278 1 jeudi FALSE 1375.857
## 3 2020-02-02 1642.750 2 dimanche TRUE 1569.997
## 4 2020-03-02 1372.659 3 lundi FALSE 1417.050
## 5 2020-04-02 1628.074 4 jeudi FALSE 1654.542
## 6 2020-05-02 1554.748 5 samedi TRUE 1532.906
## 7 2020-06-02 1443.947 6 mardi FALSE 1530.296
## 8 2020-07-02 1685.342 7 jeudi FALSE 1706.194
## 9 2020-08-02 1315.632 8 dimanche TRUE 1318.397
## 10 2020-09-02 1380.180 9 mercredi FALSE 1432.579
## 11 2020-10-02 1446.483 10 vendredi FALSE 1498.829
## 12 2020-11-02 1565.023 11 lundi FALSE 1586.650