4 - Lissage exponentiel

Analyse des séries temporelles avec

Auteur

Alain Quartier-la-Tente

L’objectif de ce TP est d’apprendre à appliquer les modèles ETS

Les packages suivants seront utilisés :

packages_to_install <- c("ggplot2", "patchwork", "forecast", "RJDemetra", "ggdemetra")

packages <- packages_to_install[! packages_to_install %in% installed.packages()[,"Package"]]
if (length(packages) > 0) {
    install.packages(packages)
}
library(forecast)
library(ggplot2)
library(patchwork)

1 Modèles exponentiels

Le but des deux premiers exercices est d’apprendre à manipuler les modèles exponentiels et à analyser la qualité des prévisions.

Exercice

Etudier les séries co2 et UKgas : quel modèle parait le plus adapté ? Faut-il transformer la série ? Comparer les prévisions en utilisant des schémas additifs et multiplicatifs et en transformant ou non la série.

Solution

autoplot(co2) + autoplot(UKgas)

Les deux séries ont une tendance et une saisonnalité. La saisonnalité parait additive pour co2 et multiplicative pour UKGas. Pas de raison de transformer les séries.

ets_co2_add <- ets(co2, model = "ZZA")
ets_co2_mult <- ets(co2, model = "ZZM")
ets_co2_add

ETS(M,A,A) 

Call:
 ets(y = co2, model = "ZZA") 

  Smoothing parameters:
    alpha = 0.5995 
    beta  = 0.0065 
    gamma = 0.129 

  Initial states:
    l = 315.2927 
    b = 0.0772 
    s = -0.8309 -1.8609 -3.0483 -2.782 -1.2615 0.7793
           2.1909 2.7066 2.1724 1.2282 0.6624 0.0439

  sigma:  9e-04

     AIC     AICc      BIC 
1748.989 1750.349 1819.513

ets_co2_mult

ETS(M,Ad,M) 

Call:
 ets(y = co2, model = "ZZM") 

  Smoothing parameters:
    alpha = 0.6824 
    beta  = 0.0415 
    gamma = 3e-04 
    phi   = 0.9773 

  Initial states:
    l = 315.3348 
    b = 0.1077 
    s = 0.9972 0.9939 0.9904 0.9909 0.9963 1.0024
           1.0069 1.0088 1.0074 1.004 1.0018 0.9999

  sigma:  8e-04

     AIC     AICc      BIC 
1721.105 1722.628 1795.777

autoplot(window(co2, start = 1993), y = "Millions de thermies") +
    autolayer(forecast(ets_co2_add, h = 60, PI = FALSE), "ETS(M,A,A)") +
    autolayer(forecast(ets_co2_mult, h = 60, PI = FALSE), "ETS(M,Ad,M)")

Les prévisions très proches sur le court terme mais s’éloignent sur le long terme. Cela vient notamment la tendance est amortie dans le modèle multiplicatif.

ets_gas_add <- ets(UKgas, model = "ZZA")
ets_gas_mult <- ets(UKgas, model = "ZZM")
ets_gas_add_log <- ets(UKgas, model = "ZZA", lambda = 0)
ets_gas_add_log_unb <- ets(UKgas, model = "ZZA", lambda = 0, biasadj = TRUE)

autoplot(window(UKgas, start = 1970), y = "co2") +
    autolayer(forecast(ets_gas_add, h = 24, PI = FALSE), "ETS(A,A,A)") +
    autolayer(forecast(ets_gas_mult, h = 24, PI = FALSE), "ETS(M,A,M)") +
    autolayer(forecast(ets_gas_add_log, h = 24, PI = FALSE), "ETS(A,A,A) sur log(UKgas)")

Ici la différence entre multiplicatif et additif est plus nette : Dans les deux modèles on prévoit une hausse de la tendance mais les amplitudes croissent de manière exponentielle dans le modèle multiplicatif alors qu’elles restent constantes dans le cas additif (ce qui est logique). Passer au log ne semble pas avoir beaucoup d’impact.

Sur le plus long terme les résultats sont en revanche différents. Passer au log conduit à des estimations plus importantes (notamment lorsque l’on corrige du bais. Rappels : sans correction du biais cela revient à avoir une estimation de la médiane plutôt que moyen, ce qui n’est pas forcément incohérent). C’est logique car lorsque l’on passe au log cela revient à supposer que la tendance est également multiplicative.

autoplot(window(UKgas, start = 1985), y = "co2")  +
    autolayer(forecast(ets_gas_mult, h = 60, PI = FALSE), "ETS(M,A,M)") +
    autolayer(forecast(ets_gas_add_log, h = 60, PI = FALSE), "ETS(A,A,A) sur log(UKgas)") +
    autolayer(forecast(ets_gas_add_log_unb, h = 60, PI = FALSE), "ETS(A,A,A) sur log(UKgas) corrigé biais")

Appliquer le modèle $ETS(A,A,A)$ sur la série en logarithme et utiliser le modèle $ETS(M,M,M)$ donnent des résultats similaires.

autoplot(window(UKgas, start = 1985), y = "co2")  +
    autolayer(forecast(ets_gas_add_log, h = 60, PI = FALSE), "ETS(A,A,A) sur log(UKgas)") +
    autolayer(forecast(ets(UKgas, model = "MMM"), h = 60, PI = FALSE), "ETS(M,M,M)")

Exercice

Sur les modèles précédents, comparer les qualités des modèles :

En termes d’AIC.
En termes de BIC.
En termes de qualité de prévision dans l’échantillon (MAE).
En termes de qualité de prévision hors échantillon (MAE). On pourra s’aider de la fonction forecast::tsCV().

Solution

Les AIC et les BIC peuvent être extraits avec les fonctions AIC() et BIC(). L’AICc peut être extrait à partir du modèle estimé, par exemple ets_co2_add$aicc.

Pour simplifier les fonctions, nous construisons une liste pour rassembler tous les modèles :

modeles_co2 <- list(
    "ETS(M,A,A)" = ets_co2_add,
    "ETS(M,Ad,M)" = ets_co2_mult
)
modeles_gas <- list(
    "ETS(A,A,A)" = ets_gas_add,
    "ETS(M,A,M)" = ets_gas_mult,
    "ETS(A,A,A) sur log(UKgas)" = ets_gas_add_log
)

Ce sont toujours les mêmes modèles qui minimisent l’AIC, l’AICc, le BIC et les erreurs de prévision dans l’échantillon :

sapply(modeles_co2, AIC)

 ETS(M,A,A) ETS(M,Ad,M) 
   1748.989    1721.105

sapply(modeles_co2, `[[`, "aicc")

 ETS(M,A,A) ETS(M,Ad,M) 
   1750.349    1722.628

sapply(modeles_co2, BIC)

 ETS(M,A,A) ETS(M,Ad,M) 
   1819.513    1795.777

sapply(modeles_co2, function(x) mean(abs(resid(x))))

  ETS(M,A,A)  ETS(M,Ad,M) 
0.0006932697 0.0006665409

sapply(modeles_gas, AIC)

               ETS(A,A,A)                ETS(M,A,M) ETS(A,A,A) sur log(UKgas) 
               1310.37144                1254.72156                  29.18537

sapply(modeles_gas, `[[`, "aicc")

               ETS(A,A,A)                ETS(M,A,M) ETS(A,A,A) sur log(UKgas) 
               1312.20817                1256.55829                  31.02211

sapply(modeles_gas, BIC)

               ETS(A,A,A)                ETS(M,A,M) ETS(A,A,A) sur log(UKgas) 
               1334.51062                1278.86074                  53.32455

sapply(modeles_gas, function(x) mean(abs(resid(x))))

               ETS(A,A,A)                ETS(M,A,M) ETS(A,A,A) sur log(UKgas) 
              26.26259517                0.06945649                0.06602586

Les statistiques sur les résidus peuvent aussi être calculés avec la fonction forecast::accuracy() :

accuracy_co2 <- do.call(rbind, lapply(modeles_co2, accuracy))
rownames(accuracy_co2) <- names(modeles_co2)
accuracy_gas <- do.call(rbind, lapply(modeles_gas, accuracy))
rownames(accuracy_gas) <- names(modeles_gas)
accuracy_co2

                    ME      RMSE       MAE         MPE       MAPE      MASE
ETS(M,A,A)  0.01617189 0.2890652 0.2334053 0.004603584 0.06932078 0.1843375
ETS(M,Ad,M) 0.04391736 0.2795208 0.2242251 0.012834152 0.06663769 0.1770872
                   ACF1
ETS(M,A,A)   0.06310175
ETS(M,Ad,M) -0.01244155

accuracy_gas

                                ME     RMSE      MAE         MPE      MAPE
ETS(A,A,A)                3.561575 38.17327 26.26260 -0.36141773 10.168961
ETS(M,A,M)                4.706001 32.16226 21.63151  0.32062812  6.581003
ETS(A,A,A) sur log(UKgas) 2.735393 31.94691 21.40806  0.06739374  6.475240
                               MASE        ACF1
ETS(A,A,A)                0.9328563 -0.06195286
ETS(M,A,M)                0.7683585 -0.09735967
ETS(A,A,A) sur log(UKgas) 0.7604214 -0.08546581

Pour la prévision hors échantillon, nous allons construire une fonction par modèle prenant comme paramètre la série brute et l’horizon de prévision et qui renvoie la série prévue :

f_co2 <- list(
    "ETS(M,A,A)" = function(x, h) {
    forecast(ets(x, model = "MAA"), h = h)
},
    "ETS(M,Ad,M)" = function(x, h) {
    forecast(ets(x, model = "MAM", damped = TRUE), h = h)
}
)
f_gas <- list(
    "ETS(A,A,A)" = function(x, h) {
    forecast(ets(x, model = "AAA"), h = h)
},
    "ETS(M,A,M)" = function(x, h) {
    forecast(ets(x, model = "MAM"), h = h)
},
    "ETS(A,A,A) sur log(UKgas)" = function(x, h) {
    forecast(ets(x, model = "AAA", lambda = 0), h = h)
}
)
# Exemple : 
f_co2[[1]](co2, 1)

         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Jan 1998       365.1537 364.7453 365.5622 364.5291 365.7784

Pour estimer les modèles hors échantillon, nous allons utiliser au moins 10 ans de données pour la base co2 (sur 39 années) et pour la base UKgas (sur 27 ans). C’est toujours les mêmes modèles qui minimisent ces erreurs :

e_co2 <- lapply(f_co2, tsCV, y = co2, h = 1, initial = 10 * 12)
e_gas <- lapply(f_gas, tsCV, y = UKgas, h = 1, initial = 10 * 4)
sapply(e_co2, function(x) mean(abs(x), na.rm = TRUE))

 ETS(M,A,A) ETS(M,Ad,M) 
  0.2429932   0.2338688

sapply(e_gas, function(x) mean(abs(x), na.rm = TRUE))

               ETS(A,A,A)                ETS(M,A,M) ETS(A,A,A) sur log(UKgas) 
                 35.90959                  34.17049                  33.38889

2 Décomposition et modèles exponentiels

Le but du prochain exercice est de voir comment coupler une méthode de décomposition (STL ou X-13ARIMA) avec une méthode de prévision (ici ETS mais cela peut être généralisé). Pour faire cet exercice il faut donc avoir vu le cours sur la décomposition d’une série temporelle et le TP associé.

Exercice

L’objectif de cette exercice est d’étudier la série AirPassengers en utilisant plusieurs méthodes :

Enlever les 12 derniers mois de AirPassengers.
Désaisonnaliser la série en utilisant stl(., s.window = "periodic") (après transformation de la série) et RJDemetra::x13(., spec = RJDemetra::x13_spec(easter.enabled = FALSE, transform.function = "Log")).
Appliquer un modèle ETS sur la série désaisonnalisée.
Prévoir la série désaisonnalisée sur 12 mois puis la série brute réintroduisant :

la saisonnalité sur la dernière année pour la méthode STL (peut se faire en une étape avec la fonction forecast::stlf()) ;
les prévisions de la saisonnalité pour la méthode X-13-ARIMA.

Comparer les prévisions des précédentes méthodes avec un ETS directement calculé sur la série désaisonnalisée (avec ou sans transformation de Box-Cox).
Qu’en est-il de la prévision en temps-réel (en utilisant la fonction tsCV()) ?

Indice

Pour récupérer la composante saisonnière avec X-13 on pourra utiliser le code suivant :

library(RJDemetra)
y <- window(ipi_c_eu[,"FR"], start = 2010)
x13_spec <- x13_spec(easter.enabled = FALSE, transform.function = "Log")
mod_x13 <- x13(y, x13_spec)
ggdemetra::seasonal(mod_x13)

           Jan       Feb       Mar       Apr       May       Jun       Jul
2010 0.9479753 0.9734184 1.1314166 1.0194466 0.9149364 1.1040334 1.0253971
2011 0.9486036 0.9738929 1.1146702 0.9990949 0.9435734 1.0943653 0.9971672
2012 0.9744944 0.9952811 1.0894809 0.9994940 0.9609052 1.0773442 1.0308806
2013 0.9881143 0.9734166 1.0586928 1.0414531 0.9661151 1.0343284 1.0633799
2014 0.9887445 0.9708317 1.0618924 1.0340629 0.9456969 1.0716857 1.0449424
2015 0.9611904 0.9664643 1.0947446 1.0268935 0.9248260 1.1112972 1.0452371
2016 0.9336952 0.9939049 1.1148269 1.0091687 0.9628352 1.0902251 0.9872514
2017 0.9643035 0.9598990 1.1181909 0.9641950 1.0015167 1.0999538 0.9910340
2018 0.9940038 0.9589041 1.0890126 0.9941445 0.9936340 1.0782580 1.0244067
2019 0.9762525 0.9595422 1.0576157 1.0265143 1.0011358 1.0398603 1.0629867
2020 0.9768492 0.9563783 1.0931460 1.0048644 0.9500369 1.1224476 1.0549500
           Aug       Sep       Oct       Nov       Dec
2010 0.7473449 1.0778667 1.0402228 1.0450781 0.9908451
2011 0.7725295 1.0855555 1.0436457 1.0336021 0.9645023
2012 0.7668908 1.0190248 1.1114204 1.0328655 0.9392665
2013 0.7511679 1.0558980 1.0943162 1.0101039 0.9672872
2014 0.7347217 1.0923548 1.0977123 0.9715733 0.9995163
2015 0.7410440 1.0774877 1.0685360 1.0097237 0.9855990
2016 0.7929649 1.0745879 1.0421123 1.0381299 0.9631477
2017 0.7846862 1.0498446 1.0744776 1.0308174 0.9322877
2018 0.7891563 1.0084274 1.1117543 1.0382887 0.9294119
2019 0.7697163 1.0458578 1.0967583 1.0142887 0.9512178
2020 0.7506828 1.0740572 1.0727587 1.0092022 0.9628865

ggdemetra::seasonaladj(mod_x13)

           Jan       Feb       Mar       Apr       May       Jun       Jul
2010  95.25565  95.64233  96.78133  98.48481 104.37885 101.26505  98.30338
2011 104.36394 104.32359 103.43866 101.69204 116.68408  99.14422 101.28692
2012 101.89900 100.37365 101.24088  99.85052 100.00987 100.71062 100.69062
2013  97.45836  98.31350  98.13990  99.28435  99.57405 101.99855  98.92984
2014  97.90194 100.01734  98.78591  99.51039  97.70572  97.78987  98.85712
2015  98.41963  99.43461  99.01853  99.71823  98.39689 100.60315  96.72447
2016 101.85337 100.21080  97.14512 102.36148 102.40590 101.17176  97.03709
2017 102.56108 101.15647 102.39754 101.63919 102.54447 101.27698 100.40019
2018 102.41410 102.82571 103.67189 103.60666  99.33235 104.61318 104.06023
2019 106.32495 106.30070 105.52037 104.43109 105.08065 101.93677 103.29386
2020 103.39364 104.66569  83.97780  66.37711  77.57594  87.48738  92.32665
           Aug       Sep       Oct       Nov       Dec
2010  99.68623 101.12568 100.93991  98.27017 102.84151
2011 101.35535 101.33060 101.95031 102.84422 103.88777
2012 102.75257 100.97889  96.81305  98.65757  97.41645
2013  97.71450  97.83142  99.87972  98.00971  97.79929
2014  97.31577  98.13661  98.11314  98.29418  98.44762
2015 102.69296 101.25406 100.79212 100.81966 101.35969
2016 100.50886 100.41059  99.22155 100.85443 103.51476
2017 103.73573 103.06287 105.35353 107.68154 106.72671
2018 104.66875 103.82503 104.33960 105.55831 105.01264
2019 102.37538 104.22067 106.22213 102.53491 102.81557
2020  95.51304  97.48084  99.46319 100.67358 100.32336

ggdemetra::seasonal(mod_x13, forecast = TRUE)

           Jan       Feb       Mar       Apr       May       Jun       Jul
2021 0.9240691 0.9625840 1.1283879 1.0116767 0.9542506 1.1122949 1.0301236
           Aug       Sep       Oct       Nov       Dec
2021 0.7736290 1.0661433 1.0421204 1.0454394 0.9628614

Si l’on veut une version plus rapide du code on peut également utiliser cette option :

mod_jx13 <- jx13(y, x13_spec)
saisonnalite <- get_indicators(mod_jx13, c("s", "sa", "s_f", "y_f"))
saisonnalite[["s"]]

           Jan       Feb       Mar       Apr       May       Jun       Jul
2010 0.9479753 0.9734184 1.1314166 1.0194466 0.9149364 1.1040334 1.0253971
2011 0.9486036 0.9738929 1.1146702 0.9990949 0.9435734 1.0943653 0.9971672
2012 0.9744944 0.9952811 1.0894809 0.9994940 0.9609052 1.0773442 1.0308806
2013 0.9881143 0.9734166 1.0586928 1.0414531 0.9661151 1.0343284 1.0633799
2014 0.9887445 0.9708317 1.0618924 1.0340629 0.9456969 1.0716857 1.0449424
2015 0.9611904 0.9664643 1.0947446 1.0268935 0.9248260 1.1112972 1.0452371
2016 0.9336952 0.9939049 1.1148269 1.0091687 0.9628352 1.0902251 0.9872514
2017 0.9643035 0.9598990 1.1181909 0.9641950 1.0015167 1.0999538 0.9910340
2018 0.9940038 0.9589041 1.0890126 0.9941445 0.9936340 1.0782580 1.0244067
2019 0.9762525 0.9595422 1.0576157 1.0265143 1.0011358 1.0398603 1.0629867
2020 0.9768492 0.9563783 1.0931460 1.0048644 0.9500369 1.1224476 1.0549500
           Aug       Sep       Oct       Nov       Dec
2010 0.7473449 1.0778667 1.0402228 1.0450781 0.9908451
2011 0.7725295 1.0855555 1.0436457 1.0336021 0.9645023
2012 0.7668908 1.0190248 1.1114204 1.0328655 0.9392665
2013 0.7511679 1.0558980 1.0943162 1.0101039 0.9672872
2014 0.7347217 1.0923548 1.0977123 0.9715733 0.9995163
2015 0.7410440 1.0774877 1.0685360 1.0097237 0.9855990
2016 0.7929649 1.0745879 1.0421123 1.0381299 0.9631477
2017 0.7846862 1.0498446 1.0744776 1.0308174 0.9322877
2018 0.7891563 1.0084274 1.1117543 1.0382887 0.9294119
2019 0.7697163 1.0458578 1.0967583 1.0142887 0.9512178
2020 0.7506828 1.0740572 1.0727587 1.0092022 0.9628865

saisonnalite[["sa"]]

           Jan       Feb       Mar       Apr       May       Jun       Jul
2010  95.25565  95.64233  96.78133  98.48481 104.37885 101.26505  98.30338
2011 104.36394 104.32359 103.43866 101.69204 116.68408  99.14422 101.28692
2012 101.89900 100.37365 101.24088  99.85052 100.00987 100.71062 100.69062
2013  97.45836  98.31350  98.13990  99.28435  99.57405 101.99855  98.92984
2014  97.90194 100.01734  98.78591  99.51039  97.70572  97.78987  98.85712
2015  98.41963  99.43461  99.01853  99.71823  98.39689 100.60315  96.72447
2016 101.85337 100.21080  97.14512 102.36148 102.40590 101.17176  97.03709
2017 102.56108 101.15647 102.39754 101.63919 102.54447 101.27698 100.40019
2018 102.41410 102.82571 103.67189 103.60666  99.33235 104.61318 104.06023
2019 106.32495 106.30070 105.52037 104.43109 105.08065 101.93677 103.29386
2020 103.39364 104.66569  83.97780  66.37711  77.57594  87.48738  92.32665
           Aug       Sep       Oct       Nov       Dec
2010  99.68623 101.12568 100.93991  98.27017 102.84151
2011 101.35535 101.33060 101.95031 102.84422 103.88777
2012 102.75257 100.97889  96.81305  98.65757  97.41645
2013  97.71450  97.83142  99.87972  98.00971  97.79929
2014  97.31577  98.13661  98.11314  98.29418  98.44762
2015 102.69296 101.25406 100.79212 100.81966 101.35969
2016 100.50886 100.41059  99.22155 100.85443 103.51476
2017 103.73573 103.06287 105.35353 107.68154 106.72671
2018 104.66875 103.82503 104.33960 105.55831 105.01264
2019 102.37538 104.22067 106.22213 102.53491 102.81557
2020  95.51304  97.48084  99.46319 100.67358 100.32336

saisonnalite[["s_f"]]

           Jan       Feb       Mar       Apr       May       Jun       Jul
2021 0.9240691 0.9625840 1.1283879 1.0116767 0.9542506 1.1122949 1.0301236
           Aug       Sep       Oct       Nov       Dec
2021 0.7736290 1.0661433 1.0421204 1.0454394 0.9628614

# On a d'ailleurs directement une prévision de la série brute qui est faite par modèle ARIMA
saisonnalite[["y_f"]]

           Jan       Feb       Mar       Apr       May       Jun       Jul
2021  93.80530  98.83526 116.70635 105.64151  99.21541 117.25615 107.92506
           Aug       Sep       Oct       Nov       Dec
2021  81.39014 112.06948 110.18417 110.21518 101.77148

Pour la fonction tsCV() on pourra utiliser la fonction suivante pour X-13 :

fx13_ets <- function(x, h){
    mod_jx13 <- jx13(x, x13_spec)
    mod_jx13 <- get_indicators(mod_jx13, c("sa", "s_f"))
    ets_x13 <- ets(mod_jx13$sa, model = "AAN") # modèle fixé pour gagner du temps
    ets_x13_f <- forecast(ets_x13, h = h)
    ets_x13_f$mean <- ets_x13_f$mean * mod_jx13$s_f[1:h] # on ajoute 1:h pour éviter quelques bugs
    ets_x13_f$model <- "X-13 + ETS"
    # Pas la peine d'actualiser autres paramètres
    ets_x13_f
}

Solution

library(RJDemetra)
y <- window(AirPassengers, end = end(AirPassengers) - c(1, 0))
autoplot(y)

# Saisonnalité présente qui dépend de la tendance : passage au log nécessaire pour STL
mod_stl <- stl(log(y), s.window = "periodic")

x13_spec <- x13_spec(easter.enabled = FALSE)
mod_jx13 <- jx13(y, x13_spec)
mod_jx13 <- get_indicators(mod_jx13, c("s", "sa", "s_f", "y_f"))
autoplot(exp(seasadj(mod_stl))) + 
    autolayer(mod_jx13$sa)

ets_x13 <- ets(mod_jx13$sa)
ets_stl <- ets(seasadj(mod_stl))
ets_x13_f <- forecast(ets_x13, h = 12)
ets_stl_f <- forecast(ets_stl, h = 12)
x13_f <- ets_x13_f$mean * mod_jx13$s_f # Il faut multiplier car schéma multiplicatif
ets_f <- ets_stl_f$mean + lag(seasonal(mod_stl), -12)
ets_f <- exp(ets_f)

On aurait directement pu obtenir les résultats avec fonction forecast::stlf() :

ets_f - stlf(y, lambda = 0, h = 12, s.window = "periodic")$mean

     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
1960   0   0   0   0   0   0   0   0   0   0   0   0

est_direct <- ets(y)
est_direct_bc <- ets(y, lambda = 0)
autoplot(window(AirPassengers, start = end(AirPassengers) - c(1, 0)),
         y = "AirPassengers") +
    autolayer(x13_f, series = "X-13 + ETS") +
    autolayer(mod_jx13$y_f, series = "X-13")+ 
    autolayer(ets_f, series = "STL + ETS")+ 
    autolayer(forecast(est_direct, PI=FALSE, h =12), series = "ETS(M,Ad,M)")+ 
    autolayer(forecast(est_direct_bc, PI=FALSE, h =12), series = "ETS(A,N,A) sur log(AirPassengers)")

Les séries qui semblent avoir les meilleurs prévisions sont X-13 (modèle ARIMA, voir TP 5) et X-13 + ETS. Le moins bon semble l’ETS directement calculé.

fx13_ets <- function(x, h = 12){
    mod_jx13 <- jx13(x, x13_spec)
    mod_jx13 <- get_indicators(mod_jx13, c("sa", "s_f"))
    ets_x13 <- ets(mod_jx13$sa, model = "AAN") # modèle fixé pour gagner du temps
    ets_x13_f <- forecast(ets_x13, h = h)
    ets_x13_f$mean <- ets_x13_f$mean * mod_jx13$s_f[1:h] # on ajoute 1:h pour éviter quelques bugs
    ets_x13_f$model <- "X-13 + ETS"
    # Pas la peine d'actualiser autres paramètres
    # ets_x13_f$upper = ets_x13_f$upper * mod_jx13$s_f
    # ets_x13_f$lower = ets_x13_f$lower * mod_jx13$s_f
    ets_x13_f
}
fx13 <- function(x, h = 12){
    mod_jx13 <- jx13(x, x13_spec)
    x13_f <- get_indicators(mod_jx13, c("y_f"))$y_f
    ## Ici il y a des calculs en trop car la prévision est uniquement faite avec le pré-ajustement, on pourrait encore optimiser en ne faisant que le pré-ajustement
    # mod_jx13 <- jregarima(y, x13_spec$regarima)
    # x13_f <- get_indicators(mod_jx13, c("model.y_f"))[[1]]
    x13_f <- window(x13_f, end = time(x13_f)[h])
    return(structure(list(method = "X-13ARIMA", 
                          model = mod_jx13, 
                          mean = x13_f,
                          level = NULL, 
                          lower = NULL, upper = NULL, x = NULL, series = NULL, 
                          fitted = NULL, residuals =NULL), class = "forecast"))
}
fstl_ets <- function(x, h = 12){
    stlf(x, lambda = 0, h = h, s.window = "periodic", etsmodel = "AAN")
}
fets <- function(x, h){
    forecast(ets(x, model = "MAM", damped = TRUE), h = h)
}
fets_bc <- function(x, h){
    forecast(ets(x, lambda = 0, model = "AAA"), h = h)
}
# On enlève les 3 premières années pour X-13
e_x13_ets <- tsCV(AirPassengers, fx13_ets, h = 1, initial = 3*12)
e_x13 <- tsCV(AirPassengers, fx13, h = 1, initial = 3*12)
e_stl_est <- tsCV(AirPassengers, fstl_ets, h = 1, initial = 3*12)
e_ets <- tsCV(AirPassengers, fets, h = 1, initial = 3*12)
e_ets_bc <- tsCV(AirPassengers, fets_bc, h = 1, initial = 3*12)

erreur <- ts.union(e_x13_ets, e_x13, e_stl_est,
                  e_ets, e_ets_bc)
colnames(erreur) <- c("X-13+ETS", "X-13", "STL+ETS", "ETS", "log(ETS)")

Le modèle ARIMA (issu de X-13) semble de meilleure qualité :

colMeans(erreur^2, na.rm = TRUE)

X-13+ETS     X-13  STL+ETS      ETS log(ETS) 
146.4192 140.2659 202.0902 206.2573 210.5955