Time Series Forecasting
Josue Afouda
2023-07-15
library(TSstudio)
library(forecast)Description du flux de travail en prévision de séries temporelles
La prévision traditionnelle des séries temporelles suit le même flux de travail que la plupart des domaines de l’analyse prédictive, tels que la régression ou la classification, et comprend généralement les étapes suivantes :
Préparation des données : Ici, nous préparons les données pour le processus de formation et de test du modèle. Cette étape comprend la division de la série en partitions d’apprentissage (dans l’échantillon) et de test (hors échantillon), la création de nouvelles fonctionnalités (le cas échéant) et l’application d’une transformation si nécessaire (par exemple, transformation de journal, mise à l’échelle, etc.).
Entraîner le modèle : ici, nous avons utilisé la partition d’entraînement pour entraîner un modèle statistique. L’objectif principal de cette étape est d’utiliser l’ensemble d’apprentissage pour former, régler et estimer les coefficients du modèle qui minimisent les critères d’erreur sélectionnés (plus loin dans ce chapitre, nous discuterons en détail des métriques d’erreur courantes). Les valeurs ajustées et l’estimation du modèle des observations de la partition d’apprentissage seront utilisées ultérieurement pour évaluer la performance globale du modèle.
Tester le modèle : Ici, nous utilisons le modèle entraîné pour prévoir les observations correspondantes de la partition de test. L’idée principale ici est d’évaluer les performances du modèle avec un nouvel ensemble de données (que le modèle n’a pas vu pendant le processus de formation).
Évaluation du modèle : Enfin et surtout, une fois le modèle formé et testé, il est temps d’évaluer les performances globales du modèle sur les partitions de formation et de test.
Sur la base du processus d’évaluation du modèle, si le modèle répond à un certain seuil ou à certains critères, nous conservons le modèle en utilisant la série complète afin de générer la prévision finale ou sélectionnons un nouveau paramètre de formation/modèle différent et répétons le processus de formation. .
D’un autre côté, ce flux de travail a ses propres caractéristiques qui la distinguent des autres domaines de l’analyse prédictive :
Les partitions de formation et de test doivent être classées par ordre chronologique, par opposition à un échantillonnage aléatoire.
En règle générale, une fois que nous avons entraîné et testé le modèle à l’aide des partitions d’entraînement et de test, nous réentraînons le modèle avec toutes les données (ou au moins l’observation la plus récente dans l’ordre chronologique). À première vue, cela peut être choquant et terrifiant pour les personnes ayant une formation en apprentissage automatique traditionnel ou en modélisation d’apprentissage supervisé, car cela entraîne généralement un surajustement et d’autres problèmes. Nous discuterons de la raison derrière cela et comment éviter le surajustement plus tard.
Pour illustrer le flux de travail en prévision de série temporelle, nous utiliserons la série USgas
Modélisation de la série temporelle USgas
data("USgas")
ts_plot(USgas)la principale caractéristique de la série USgas est la
forte tendance saisonnière.
Partion de la série
La partition en données d’entraînement et de test peut s’effectuer
avec la fonction ts_split() du package TSstudio
:
# par défaut la taille de la série Test est égale à 30% de la série totale
USgas_partitions <- ts_split(USgas)
train <- USgas_partitions$train
test <- USgas_partitions$test
ts_info(train)## The train series is a ts object with 1 variable and 167 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2013 11ts_info(test)## The test series is a ts object with 1 variable and 71 observations
## Frequency: 12
## Start time: 2013 12
## End time: 2019 10N.B : L’argument sample.out de la
fonction ts_split permet de spécifier explicitement la
taille de la série Test (et donc indirectement celle de la série
Train).
Modélisation automatique d’une série temporelle
La fonction “auto.arima” du package forecast est utilisée pour effectuer une modélisation automatique de séries chronologiques à l’aide de la méthode ARIMA (AutoRegressive Integrated Moving Average). Elle présente plusieurs avantages :
Sélection automatique du modèle : La fonction “auto.arima” utilise des techniques de sélection de modèle automatiques pour déterminer le meilleur modèle ARIMA pour une série chronologique donnée. Elle explore différents ordres de modèles ARIMA, tels que les ordres d’autorégression (p), d’intégration (d) et de moyenne mobile (q), et utilise des critères d’information tels que le critère d’information d’Akaike (AIC) ou le critère d’information bayésien (BIC) pour évaluer la qualité de chaque modèle. Cela permet de trouver rapidement et efficacement un modèle ARIMA approprié, sans avoir besoin d’effectuer manuellement des tests et des évaluations pour différents ordres de modèles.
Gain de temps : La sélection manuelle des ordres de modèle ARIMA peut être une tâche complexe et chronophage, nécessitant des connaissances spécialisées et des itérations multiples. La fonction “auto.arima” automatise ce processus, ce qui permet de gagner du temps dans le choix du meilleur modèle ARIMA.
Robustesse : La fonction “auto.arima” utilise des méthodes robustes de sélection de modèle, ce qui réduit les risques de sélection d’un mauvais modèle. Elle prend en compte différents aspects, tels que la tendance, la saisonnalité et la présence de valeurs aberrantes, lors de la sélection du modèle ARIMA.
Flexibilité : La fonction “auto.arima” est hautement configurable et permet de spécifier divers paramètres, tels que les ordres maximaux du modèle ARIMA à considérer, les restrictions sur les ordres de modèle (par exemple, en fixant une valeur spécifique pour certains ordres), et d’autres options avancées. Cela permet aux utilisateurs d’adapter la sélection automatique du modèle en fonction de leurs besoins spécifiques.
En résumé, la fonction “auto.arima” simplifie le processus de modélisation de séries chronologiques en effectuant automatiquement la sélection du modèle ARIMA le plus approprié. Elle permet de gagner du temps, d’obtenir des résultats robustes et de faciliter l’analyse et la prévision des séries chronologiques.
# Entraînement d'un modèle ARIMA
md <- auto.arima(train)
md## Series: train
## ARIMA(1,0,0)(1,1,2)[12]
##
## Coefficients:
## ar1 sar1 sma1 sma2
## 0.6430 -0.5680 -0.0912 -0.5406
## s.e. 0.0808 0.3297 0.3138 0.2236
##
## sigma^2 = 10628: log likelihood = -941.84
## AIC=1893.68 AICc=1894.08 BIC=1908.89Le résultat fourni est une représentation du modèle ARIMA sélectionné pour la série chronologique “train” à l’aide de la fonction “auto.arima”. Voici une explication détaillée des différents éléments du résultat :
“Series: train”: Indique le nom de la série chronologique sur laquelle le modèle a été ajusté, dans ce cas, “train”.
“ARIMA(1,0,0)(1,1,2)[12]”: C’est la spécification du modèle ARIMA sélectionné. Les nombres entre parenthèses représentent les ordres (p, d, q) de la partie ARIMA, tandis que les nombres entre crochets représentent les ordres saisonniers (P, D, Q) et la période saisonnière. Dans cet exemple, le modèle est spécifié comme ARIMA(1,0,0)(1,1,2)[12]. Cela signifie qu’il y a un terme d’autorégression d’ordre 1 (AR(1)), aucun terme de différenciation (d = 0) et aucun terme de moyenne mobile (MA = 0) dans la partie non saisonnière. Dans la partie saisonnière, il y a un terme d’autorégression saisonnier d’ordre 1 (SAR(1)), deux termes de moyenne mobile saisonniers d’ordre 1 et 2 (SMA(1) et SMA(2)), avec une période saisonnière de 12 (pour des données mensuelles par exemple).
“Coefficients”: Cette section affiche les estimations des coefficients du modèle ARIMA. Dans cet exemple, il y a quatre coefficients estimés : ar1, sar1, sma1 et sma2, correspondant respectivement aux coefficients du terme d’autorégression non saisonnier (AR(1)), du terme d’autorégression saisonnier (SAR(1)), du premier terme de moyenne mobile saisonnier (SMA(1)) et du deuxième terme de moyenne mobile saisonnier (SMA(2)). Les coefficients estimés sont respectivement 0.6430, -0.5680, -0.0912 et -0.5406.
“s.e.”: Cette section affiche les erreurs standard des coefficients estimés.
“sigma^2”: C’est l’estimation de la variance du bruit (ou de l’erreur) dans le modèle ARIMA. Dans cet exemple, sigma^2 est estimé à 10628.
“log likelihood”: Il s’agit de la valeur de la log-vraisemblance maximisée du modèle ARIMA. Dans cet exemple, la log-vraisemblance maximisée est -941.84.
“AIC”, “AICc” et “BIC”: Ce sont des critères d’information utilisés pour comparer différents modèles. AIC (Akaike Information Criterion), AICc (AIC corrigé) et BIC (Bayesian Information Criterion) sont souvent utilisés pour sélectionner le meilleur modèle. Dans cet exemple, les valeurs des critères sont respectivement AIC=1893.68, AICc=1894.08 et BIC=1908.89. Un modèle avec des valeurs de critères plus faibles est généralement préféré car il indique un meilleur ajustement aux données avec une complexité minimale.
L’équation représentant le modèle ARIMA trouvé à partir des résultats donnés est :
Y(t) = 0.6430 * Y(t-1) - 0.5680 * Y(t-12) - 0.0912 * e(t-1) - 0.5406 * e(t-2) + e(t)
où :
Y(t) est la valeur de la série chronologique à l’instant t.
Y(t-1) est la valeur de la série chronologique à l’instant t-1 (terme d’autorégression non saisonnier).
Y(t-12) est la valeur de la série chronologique à l’instant t-12 (terme d’autorégression saisonnier avec une période de 12).
e(t-1) est l’erreur (bruit) à l’instant t-1 (premier terme de moyenne mobile saisonnier).
e(t-2) est l’erreur (bruit) à l’instant t-2 (deuxième terme de moyenne mobile saisonnier).
e(t) est l’erreur (bruit) à l’instant t (bruit blanc gaussien).
Cela indique que la valeur de la série chronologique à un instant donné dépend de ses propres valeurs passées, des valeurs saisonnières passées, ainsi que des erreurs (bruit) passées.
# Performance du modèle sur les données d'entraînement
checkresiduals(md)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(1,1,2)[12]
## Q* = 29.796, df = 20, p-value = 0.07322
##
## Model df: 4. Total lags used: 24Voici une explication détaillée du code ci-dessus :
La fonction “checkresiduals” est appelée avec le modèle ARIMA “md” en tant qu’argument. Cette fonction est spécifique à la bibliothèque “forecast” et permet d’effectuer une analyse des résidus du modèle.
L’analyse des résidus est une étape importante pour évaluer la qualité du modèle ajusté. Elle permet de vérifier si les résidus du modèle présentent des comportements indésirables tels que des tendances, des corrélations, des valeurs aberrantes ou des non-stationnarités résiduelles.
La fonction “checkresiduals” génère plusieurs diagnostics et graphiques pour évaluer les résidus du modèle ARIMA :
Graphique des résidus : Il affiche le tracé des résidus du modèle au fil du temps. Cela permet de détecter des structures temporelles ou des schémas résiduels inappropriés.
Histogramme des résidus : Il représente la distribution des résidus. Une distribution normale des résidus est souhaitable.
Graphique de densité des résidus : Il montre la densité des résidus. Encore une fois, une distribution normale est souhaitable.
Graphique des autocorrélations partielles des résidus (PACF) : Il indique s’il reste des informations corrélées dans les résidus du modèle.
Test de Ljung-Box : Il effectue un test statistique pour vérifier si les résidus du modèle sont indépendants.
Test de normalité des résidus : Il teste l’hypothèse que les résidus suivent une distribution normale.
Test de stationnarité des résidus : Il vérifie si les résidus du modèle sont stationnaires.
L’objectif de cette analyse des résidus est de s’assurer que le modèle ARIMA ajusté capture correctement les motifs et la structure de la série chronologique, en minimisant les erreurs résiduelles et en obtenant des résidus qui sont aléatoires, indépendants et normalement distribués.
En utilisant la fonction “checkresiduals(md)”, on obtient ainsi une évaluation détaillée de la performance du modèle ARIMA sur les données d’entraînement, en se concentrant sur la qualité des résidus du modèle.
Q* = 29.796, df = 20, p-value = 0.07322”: Les résultats du test de Ljung-Box sont affichés ici. “Q*” est la statistique du test de Ljung-Box, qui mesure l’autocorrélation des résidus. “df” est le nombre de degrés de liberté du test, qui dépend de la longueur des résidus et du nombre de paramètres estimés dans le modèle. “p-value” est la valeur p associée au test, qui permet de décider si l’autocorrélation des résidus est significative. Dans cet exemple, la statistique du test est de 29.796 avec 20 degrés de liberté, et la valeur p est de 0.07322.
“Model df: 4. Total lags used: 24”: Cette partie fournit des informations supplémentaires sur le modèle ARIMA. “Model df” indique le nombre de degrés de liberté du modèle, qui dépend du nombre de paramètres estimés dans le modèle. “Total lags used” indique le nombre total de retards (lags) utilisés dans le test de Ljung-Box pour évaluer l’autocorrélation des résidus.
En interprétant les résultats, on peut dire que le test de Ljung-Box n’a pas trouvé suffisamment de preuves pour rejeter l’hypothèse nulle selon laquelle les résidus du modèle sont indépendants. Cela est indiqué par la valeur p de 0.07322, qui est supérieure à un seuil de signification de 0.05. Cependant, il est important de noter que l’interprétation des résultats du test de Ljung-Box doit être faite avec prudence et en considérant d’autres diagnostics et évaluations des résidus pour évaluer la qualité globale du modèle.
# Prédictions sur les données de test
fc <- forecast(md, h = length(test))
fc## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Dec 2013 2700.442 2568.291 2832.594 2498.334 2902.551
## Jan 2014 2911.024 2753.921 3068.127 2670.756 3151.292
## Feb 2014 2561.656 2395.326 2727.985 2307.277 2816.034
## Mar 2014 2357.730 2187.732 2527.728 2097.741 2617.719
## Apr 2014 1891.073 1719.582 2062.565 1628.800 2153.347
## May 2014 1684.252 1512.147 1856.357 1421.040 1947.464
## Jun 2014 1667.643 1495.285 1840.002 1404.044 1931.243
## Jul 2014 1837.714 1665.251 2010.176 1573.954 2101.473
## Aug 2014 1859.785 1687.279 2032.291 1595.960 2123.610
## Sep 2014 1656.029 1483.505 1828.552 1392.177 1919.881
## Oct 2014 1745.862 1573.332 1918.392 1482.000 2009.724
## Nov 2014 2048.286 1875.754 2220.818 1784.421 2312.151
## Dec 2014 2544.201 2365.708 2722.694 2271.219 2817.183
## Jan 2015 2783.789 2602.901 2964.677 2507.145 3060.433
## Feb 2015 2491.455 2309.582 2673.328 2213.304 2769.606
## Mar 2015 2251.138 2068.858 2433.417 1972.366 2529.910
## Apr 2015 1852.776 1670.329 2035.223 1573.748 2131.804
## May 2015 1682.902 1500.386 1865.418 1403.768 1962.036
## Jun 2015 1675.565 1493.020 1858.109 1396.387 1954.743
## Jul 2015 1846.786 1664.229 2029.342 1567.590 2125.982
## Aug 2015 1858.347 1675.786 2040.908 1579.144 2137.550
## Sep 2015 1648.344 1465.781 1830.907 1369.138 1927.550
## Oct 2015 1739.506 1556.943 1922.069 1460.300 2018.712
## Nov 2015 2016.149 1833.587 2198.711 1736.945 2295.354
## Dec 2015 2514.180 2330.132 2698.228 2232.702 2795.657
## Jan 2016 2779.692 2595.043 2964.341 2497.295 3062.088
## Feb 2016 2482.226 2297.325 2667.127 2199.444 2765.008
## Mar 2016 2280.107 2095.102 2465.113 1997.166 2563.049
## Apr 2016 1854.227 1669.179 2039.276 1571.220 2137.235
## May 2016 1670.616 1485.550 1855.682 1387.581 1953.650
## Jun 2016 1662.673 1477.599 1847.747 1379.627 1945.719
## Jul 2016 1836.236 1651.160 2021.313 1553.186 2119.287
## Aug 2016 1855.694 1670.616 2040.772 1572.642 2138.746
## Sep 2016 1650.477 1465.399 1835.556 1367.425 1933.530
## Oct 2016 1741.681 1556.604 1926.759 1458.629 2024.733
## Nov 2016 2033.479 1848.403 2218.556 1750.429 2316.529
## Dec 2016 2530.638 2342.147 2719.128 2242.367 2818.909
## Jan 2017 2781.638 2591.765 2971.510 2491.252 3072.023
## Feb 2017 2487.223 2296.777 2677.668 2195.961 2778.484
## Mar 2017 2263.496 2072.814 2454.178 1971.873 2555.119
## Apr 2017 1853.302 1662.522 2044.081 1561.529 2145.074
## May 2017 1677.529 1486.709 1868.349 1385.695 1969.363
## Jun 2017 1669.953 1479.116 1860.790 1378.093 1961.813
## Jul 2017 1842.201 1651.357 2033.045 1550.331 2134.071
## Aug 2017 1857.183 1666.337 2048.030 1565.309 2149.058
## Sep 2017 1649.254 1458.407 1840.102 1357.379 1941.130
## Oct 2017 1740.439 1549.591 1931.286 1448.563 2032.314
## Nov 2017 2023.632 1832.786 2214.478 1731.758 2315.506
## Dec 2017 2521.287 2328.297 2714.277 2226.135 2816.440
## Jan 2018 2780.530 2586.671 2974.390 2484.049 3077.012
## Feb 2018 2484.384 2290.162 2678.605 2187.347 2781.420
## Mar 2018 2272.930 2078.559 2467.301 1975.665 2570.195
## Apr 2018 1853.827 1659.394 2048.260 1556.467 2151.187
## May 2018 1673.602 1479.144 1868.061 1376.203 1971.001
## Jun 2018 1665.818 1471.349 1860.287 1368.403 1963.233
## Jul 2018 1838.813 1644.340 2033.287 1541.392 2136.235
## Aug 2018 1856.337 1661.862 2050.812 1558.913 2153.761
## Sep 2018 1649.949 1455.473 1844.425 1352.524 1947.374
## Oct 2018 1741.145 1546.669 1935.620 1443.720 2038.569
## Nov 2018 2029.225 1834.751 2223.699 1731.802 2326.648
## Dec 2018 2526.598 2329.398 2723.798 2225.006 2828.190
## Jan 2019 2781.159 2582.853 2979.465 2477.876 3084.442
## Feb 2019 2485.996 2287.230 2684.762 2182.010 2789.982
## Mar 2019 2267.572 2068.616 2466.527 1963.295 2571.848
## Apr 2019 1853.529 1654.495 2052.562 1549.133 2157.925
## May 2019 1675.833 1476.766 1874.899 1371.387 1980.278
## Jun 2019 1668.167 1469.087 1867.246 1363.701 1972.633
## Jul 2019 1840.737 1641.652 2039.823 1536.263 2145.212
## Aug 2019 1856.818 1657.730 2055.905 1552.340 2161.296
## Sep 2019 1649.555 1450.466 1848.643 1345.076 1954.033
## Oct 2019 1740.744 1541.656 1939.832 1436.265 2045.222accuracy(fc, test)## ME RMSE MAE MPE MAPE MASE
## Training set 17.42392 98.02719 71.97977 0.74708 3.621318 0.6947724
## Test set 308.83827 355.18763 312.96849 12.98995 13.166406 3.0208749
## ACF1 Theil's U
## Training set -0.08616356 NA
## Test set 0.66346316 1.203532La fonction “accuracy” est appelée avec deux arguments : “fc” et “test”. “fc” représente les prédictions du modèle sur les données de test, et “test” correspond aux valeurs réelles de la série chronologique sur lesquelles les prédictions ont été faites.
La fonction “accuracy” calcule différentes mesures d’évaluation pour évaluer la précision des prédictions par rapport aux valeurs réelles. Ces mesures comprennent généralement des métriques telles que l’erreur moyenne absolue (MAE), l’erreur quadratique moyenne (RMSE), la racine carrée de l’erreur quadratique moyenne (RMSE), ainsi que des mesures relatives telles que le coefficient de corrélation de Pearson et le coefficient de détermination (R²).
En résumé, le code “accuracy(fc, test)” évalue la précision des prédictions du modèle sur les données de test en calculant différentes mesures d’évaluation. Cela permet de quantifier la qualité des prédictions par rapport aux valeurs réelles et de fournir une évaluation de la performance du modèle.
test_forecast(actual = USgas,
forecast.obj = fc,
test = test) Pour effectuer des prévisions au-delà de la période de test, vous pouvez utiliser la fonction “forecast” sur le modèle “md” en spécifiant la longueur des prévisions souhaitées. Voici le code pour cela :
# Longueur des prévisions souhaitées au-delà de la période de test
horizon <- 12 # Prévisions sur les 12 prochains mois suivant la série
# Effectuer les prévisions au-delà de la période de test
future_predictions <- forecast(md, h = horizon)
future_predictions## Point Forecast Lo 80 Hi 80 Lo 95 Hi 95
## Dec 2013 2700.442 2568.291 2832.594 2498.334 2902.551
## Jan 2014 2911.024 2753.921 3068.127 2670.756 3151.292
## Feb 2014 2561.656 2395.326 2727.985 2307.277 2816.034
## Mar 2014 2357.730 2187.732 2527.728 2097.741 2617.719
## Apr 2014 1891.073 1719.582 2062.565 1628.800 2153.347
## May 2014 1684.252 1512.147 1856.357 1421.040 1947.464
## Jun 2014 1667.643 1495.285 1840.002 1404.044 1931.243
## Jul 2014 1837.714 1665.251 2010.176 1573.954 2101.473
## Aug 2014 1859.785 1687.279 2032.291 1595.960 2123.610
## Sep 2014 1656.029 1483.505 1828.552 1392.177 1919.881
## Oct 2014 1745.862 1573.332 1918.392 1482.000 2009.724
## Nov 2014 2048.286 1875.754 2220.818 1784.421 2312.151Voici comment, vous pouvez afficher les prévisions dans un graphique :
plot_forecast(future_predictions,
title = "The US Natural Gas Consumption Forecast",
Xtitle = "Year",
Ytitle = "Billion Cubic Feet")