Tarea - Suavización exponencial

Tarea

Seleccionar dos series de tiempo vistas en clase y realizar lo siguiente:

De preferencia, que una serie de tiempo tenga un patrón estacional y la otra sin estacionalidad.

library(tidyverse)
library(tsibble)
library(feasts)
library(fable)
library(tsibbledata)
library(fpp3)
library(plotly)
library(dplyr)

Selección de series de tiempo

tabaco <- aus_production |>
  select(Tobacco) |>
  filter_index(. ~ "2004 Q2")

tabaco

# A tsibble: 194 x 2 [1Q]
   Tobacco Quarter
     <dbl>   <qtr>
 1    5225 1956 Q1
 2    5178 1956 Q2
 3    5297 1956 Q3
 4    5681 1956 Q4
 5    5577 1957 Q1
 6    5651 1957 Q2
 7    5317 1957 Q3
 8    6152 1957 Q4
 9    5758 1958 Q1
10    5641 1958 Q2
# ℹ 184 more rows

autoplot(tabaco)

Plot variable not specified, automatically selected `.vars = Tobacco`

gdp_españa <- global_economy |>
  filter(Country == "Spain") |>
  select(GDP)

gdp_españa

# A tsibble: 58 x 2 [1Y]
            GDP  Year
          <dbl> <dbl>
 1 12072126075.  1960
 2 13834300571.  1961
 3 16138545209.  1962
 4 19074913948.  1963
 5 21343844644.  1964
 6 24756958695.  1965
 7 28721062242.  1966
 8 31647119228.  1967
 9 31475548481.  1968
10 36038711600.  1969
# ℹ 48 more rows

autoplot(gdp_españa)

Plot variable not specified, automatically selected `.vars = GDP`

Selección de datos de entrenamiento

tabaco_train <- tabaco |>
  filter_index("1975 Q1" ~ "2001 Q4")

tabaco_train

# A tsibble: 108 x 2 [1Q]
   Tobacco Quarter
     <dbl>   <qtr>
 1    6971 1975 Q1
 2    8021 1975 Q2
 3    8133 1975 Q3
 4    8349 1975 Q4
 5    6687 1976 Q1
 6    7991 1976 Q2
 7    6900 1976 Q3
 8    7119 1976 Q4
 9    7156 1977 Q1
10    8671 1977 Q2
# ℹ 98 more rows

autoplot(tabaco_train)

Plot variable not specified, automatically selected `.vars = Tobacco`

españa_train <- gdp_españa |>
  filter_index("1960" ~ "2012")

españa_train

# A tsibble: 53 x 2 [1Y]
            GDP  Year
          <dbl> <dbl>
 1 12072126075.  1960
 2 13834300571.  1961
 3 16138545209.  1962
 4 19074913948.  1963
 5 21343844644.  1964
 6 24756958695.  1965
 7 28721062242.  1966
 8 31647119228.  1967
 9 31475548481.  1968
10 36038711600.  1969
# ℹ 43 more rows

autoplot(españa_train)

Plot variable not specified, automatically selected `.vars = GDP`

Gráfica de las series de entrenamiento

t <- tabaco_train |> 
  autoplot(Tobacco) +
  labs(
    title = "Serie de tiempo del Tabaco",
    y = "Tobacco"
  )
 
ggplotly(t, dynamicTicks = TRUE) |> 
  rangeslider()

tabaco_train |> 
  gg_season(Tobacco) |> 
  ggplotly() |>
  layout(title = "Estacionalidad de la serie de tiempo del Tabaco")

e <- españa_train |>
  autoplot(GDP) +
  labs(
    title = "Serie de tiempo del GDP de España",
    y = "GDP"
  )

ggplotly(e, dynamicTicks = TRUE) |> 
  rangeslider()

Descomposición de las series de tiempo

#Gráfica con logaritmo
españa_train |>
  autoplot(log(GDP)) +
  ggtitle("Serie de tiempo del GDP de España con transformación logarítmica")

#Descomposición con logaritmo de españa
españa_train |> 
  model(stl = STL(log(GDP), robust = TRUE)) |> 
  components() |> 
  autoplot() |> 
  ggplotly() |>
  layout(title = "Descomposición STL del GDP de españa con transformación log")

#Descomposición del tabaco
tabaco_train |>
  model(stl = STL(Tobacco, robust = TRUE)) |>
  components() |>
  autoplot() |>
  ggplotly() |>
  layout(title = "Descomposición STL del Tabaco")

Seleccionar el método de referencia que mejores resultados haya dado para cada serie.

Creación de modelos (fit)

españa_fit_l <- españa_train |> 
  model(
    mean = MEAN(log(GDP)),
    naive = NAIVE(log(GDP)),
    drift  = RW(log(GDP) ~ drift())
  )

tabaco_fit <- tabaco_train |> 
  model(
    mean = MEAN(Tobacco),
    naive = NAIVE(Tobacco),
    snaive = SNAIVE(Tobacco),
    drift  = RW(Tobacco ~ drift())
  )

españa_fit_l

# A mable: 1 x 3
     mean   naive         drift
  <model> <model>       <model>
1  <MEAN> <NAIVE> <RW w/ drift>

tabaco_fit

# A mable: 1 x 4
     mean   naive   snaive         drift
  <model> <model>  <model>       <model>
1  <MEAN> <NAIVE> <SNAIVE> <RW w/ drift>

Análisis de residuos del log(GDP) de España

españa_fit_l |> 
  select(drift) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Drift Model")

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

españa_fit_l |> 
  select(mean) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Mean Model")

españa_fit_l |> 
  select(naive) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Naive Model")

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

españa_accuracy_l <- accuracy(españa_fit_l) |>
  arrange(MAPE)

españa_accuracy_l

# A tibble: 3 × 10
  .model .type              ME    RMSE     MAE      MPE   MAPE  MASE RMSSE  ACF1
  <chr>  <chr>           <dbl>   <dbl>   <dbl>    <dbl>  <dbl> <dbl> <dbl> <dbl>
1 drift  Training     -1.74e10 7.98e10 4.42e10   -0.669   9.12 0.901  1.09 0.393
2 naive  Training      2.55e10 7.31e10 4.90e10    8.04   11.7  1      1    0.361
3 mean   Training      2.49e11 5.36e11 3.73e11 -182.    249.   7.62   7.33 0.948

#El modelo con mejor error MAPE es el DRIFT con 9.11

Análisis de residuos del Tabaco

tabaco_fit |> 
  select(drift) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Drift Model")

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

tabaco_fit |> 
  select(mean) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Mean Model")

tabaco_fit |> 
  select(naive) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Naive Model")

Warning: Removed 1 row containing missing values (`geom_line()`).

Warning: Removed 1 rows containing missing values (`geom_point()`).

Warning: Removed 1 rows containing non-finite values (`stat_bin()`).

tabaco_fit |> 
  select(snaive) |> 
  gg_tsresiduals() +
  ggtitle("Residuals Diagnostics for the Seasonal Naive Model")

Warning: Removed 4 rows containing missing values (`geom_line()`).

Warning: Removed 4 rows containing missing values (`geom_point()`).

Warning: Removed 4 rows containing non-finite values (`stat_bin()`).

tabaco_accuracy <- accuracy(tabaco_fit) |>
  arrange(MAPE)

tabaco_accuracy

# A tibble: 4 × 10
  .model .type           ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
  <chr>  <chr>        <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
1 snaive Training -1.21e+ 2  587.  453. -2.44   7.40  1     1     0.306
2 drift  Training -3.31e-13  815.  655. -0.883 10.4   1.45  1.39 -0.321
3 naive  Training -2.29e+ 1  815.  656. -1.25  10.4   1.45  1.39 -0.321
4 mean   Training  8.00e-14 1102.  924. -3.20  15.3   2.04  1.88  0.713

#El modelo con mejor error MAPE es el SNAIVE con 7.39

Análisis de correlación de Portmanteau

españa_fit_l |> 
  augment() |> 
  features(.innov, ljung_box, lag = 5, dof = 0)

# A tibble: 3 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 drift     15.1    0.0101
2 mean     195.     0     
3 naive     15.1    0.0101

tabaco_fit |> 
  augment() |>
  features(.innov, ljung_box, lag = 10, dof = 0)

# A tibble: 4 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 drift    129.  0        
2 mean     447.  0        
3 naive    129.  0        
4 snaive    36.7 0.0000644

Modelado usando descomposición

españa_fit_desc <- españa_train |> 
      model(
        stlf = decomposition_model(
          STL(log(GDP) ~ season(window = "periodic"), robust = TRUE),
          RW(season_adjust ~ drift())
        )
      )

españa_fit_desc

# A mable: 1 x 1
                       stlf
                    <model>
1 <STL decomposition model>

Análisis de residuos del modelo de descomposición del GDP de España

españa_fit_desc |> 
  accuracy() |> 
  arrange(MAPE)

# A tibble: 1 × 10
  .model .type               ME      RMSE     MAE    MPE  MAPE  MASE RMSSE  ACF1
  <chr>  <chr>            <dbl>     <dbl>   <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
1 stlf   Training -17440494957.   7.98e10 4.42e10 -0.669  9.12 0.901  1.09 0.393

#Tiene el mismo MAPE que el drift con 9.11

españa_fit_desc |> 
  augment() |> 
  features(.innov, ljung_box)

# A tibble: 1 × 3
  .model lb_stat lb_pvalue
  <chr>    <dbl>     <dbl>
1 stlf      11.7  0.000626

Modelar cada serie de tiempo mediante suavización exponencial.

Pueden utilizar dos o tres variantes del modelo, según sea el caso. Llevar a cabo el proceso completo de pronóstico.

Modelo nuevo GDP España

españa_fit_exp <- españa_train |> 
  model(
    ets_AAN = ETS(GDP ~ error("A") + trend("A") + season("N")),
    ets_AAdN = ETS(GDP ~ error("A") + trend("Ad") + season("N")),
    ets_MAdN = ETS(GDP ~ error("M") + trend("Ad") + season("N")),
  )

españa_fit_exp

# A mable: 1 x 3
       ets_AAN      ets_AAdN      ets_MAdN
       <model>       <model>       <model>
1 <ETS(A,A,N)> <ETS(A,Ad,N)> <ETS(M,Ad,N)>

Modelos nuevos Tabaco

tabaco_fit_exp <- tabaco_train |> 
  model(
    snaive = SNAIVE(Tobacco),
    ets_ANA = ETS(Tobacco ~ error("A") + trend("N") + season("A")),
    ets_AAdA = ETS(Tobacco ~ error("A") + trend("Ad") + season("A")),
    ets_MAdM = ETS(Tobacco ~ error("M") + trend("Ad") + season("M")),
    ets_ANA_l = ETS(log(Tobacco) ~ error("A") + trend("N") + season("A")),
    ets_AAdA_l = ETS(log(Tobacco) ~ error("A") + trend("Ad") + season("A")),
    ets_MAdM_l = ETS(log(Tobacco) ~ error("M") + trend("Ad") + season("M")),
    stl_ets_A = decomposition_model(
      STL(log(Tobacco), robust = TRUE),
      ETS(season_year ~ error("A") + trend("N") + season("A")),
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
    ),
    stl_ets_M = decomposition_model(
      STL(log(Tobacco), robust = TRUE),
      ETS(season_year ~ error("M") + trend("N") + season("M")),
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
    )
  )

tabaco_fit_exp

# A mable: 1 x 9
    snaive      ets_ANA      ets_AAdA      ets_MAdM    ets_ANA_l    ets_AAdA_l
   <model>      <model>       <model>       <model>      <model>       <model>
1 <SNAIVE> <ETS(A,N,A)> <ETS(A,Ad,A)> <ETS(M,Ad,M)> <ETS(A,N,A)> <ETS(A,Ad,A)>
# ℹ 3 more variables: ets_MAdM_l <model>, stl_ets_A <model>, stl_ets_M <model>

Analizar si los modelos de suavización exponencial superan a los métodos de referencia. En caso de que el método de referencia sea mejor, intenten cambiar las características de la suavización exponencial, logrando un mejor pronóstico.

Métricas de error de modelos nuevos

accuracy(españa_fit_exp) |>
  arrange(MAPE)

# A tibble: 3 × 10
  .model   .type              ME    RMSE     MAE   MPE  MAPE  MASE RMSSE    ACF1
  <chr>    <chr>           <dbl>   <dbl>   <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
1 ets_AAN  Training -2661528789. 6.77e10 3.89e10 0.621  8.33 0.794 0.925  0.0965
2 ets_MAdN Training  1571614796. 6.66e10 3.89e10 1.81   8.44 0.793 0.911 -0.0770
3 ets_AAdN Training  6550867499. 6.44e10 3.88e10 3.10   8.60 0.792 0.880  0.0546

# Tenemos un nuevo mejor modelo, el ets_AAN con un MAPE de 8.3264

accuracy(tabaco_fit_exp) |>
  arrange(MAPE)

# A tibble: 9 × 10
  .model     .type        ME  RMSE   MAE   MPE  MAPE  MASE RMSSE    ACF1
  <chr>      <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
1 stl_ets_A  Training  -31.1  445.  347. -1.02  5.62 0.766 0.759  0.103 
2 ets_MAdM   Training  -75.6  470.  362. -1.75  5.92 0.799 0.802  0.0423
3 ets_MAdM_l Training  -42.7  471.  365. -1.26  5.92 0.806 0.803  0.0268
4 ets_AAdA_l Training  -40.4  470.  365. -1.23  5.93 0.807 0.801  0.0265
5 ets_AAdA   Training  -66.1  472.  367. -1.57  6.00 0.810 0.804  0.0306
6 ets_ANA    Training  -98.0  479.  368. -2.16  6.06 0.813 0.816 -0.0349
7 ets_ANA_l  Training  -82.6  478.  371. -1.88  6.07 0.820 0.815  0.0220
8 snaive     Training -121.   587.  453. -2.44  7.40 1     1      0.306 
9 stl_ets_M  Training  484.  1898.  956.  5.85 14.0  2.11  3.24   0.197 

# Tenemos un nuevo mejor modelo, el stl_ets_A con un MAPE de 5.61

Creación de pronósticos

españa_fc_l <- españa_fit_l |>
  forecast(h = "5 years")

españa_fc_e <- españa_fit_exp |>
  forecast(h = "5 years")

graf1 = españa_fc_l |>
  autoplot(gdp_españa |>
             filter_index("2000" ~ .), level = NULL) +
  facet_wrap(~.model, ncol=1) +
  xlab("Año") + ylab("GDP") +
  ggtitle("Test para el log(GDP) de España (Modelos tradicionales)") +
  guides(colour=guide_legend(title="Forecast"))



graf2 = españa_fc_e |>
  autoplot(gdp_españa |>
             filter_index("2000" ~ .), level = NULL) +
  facet_wrap(~.model, ncol=1) +
  xlab("Año") + ylab("GDP") +
  ggtitle("Test para el log(GDP) de España (Modelo nuevo)") +
  guides(colour=guide_legend(title="Forecast"))

graf1

graf2

tabaco_fc <- tabaco_fit |>
  forecast(h = 10)

tabaco_fc_e <- tabaco_fit_exp |>
  forecast(h = 10)

graf3 <- tabaco_fc |>
  autoplot(tabaco |>
             filter_index("1999" ~.), level = NULL) +
  facet_wrap(~.model, ncol=2) +
  xlab("Año") + ylab("GDP") +
  ggtitle("Test para el Tabaco") +
  guides(colour=guide_legend(title="Forecast"))

graf4 <- tabaco_fc_e |>
    autoplot(tabaco |>
             filter_index("1999" ~.), level = NULL) +
  facet_wrap(~.model, ncol=3) +
  xlab("Año") + ylab("GDP") +
  ggtitle("Test para el Tabaco") +
  guides(colour=guide_legend(title="Forecast"))


graf3

graf4

Pronóstico hacia el futuro con los mejores modelos

españa_fut <- gdp_españa |>
  model(
    ets_AAN = ETS(log(GDP) ~ error("A") + trend("A") + season("N"))
  ) |>
  forecast(h = "5 years")

tabaco_fut <- tabaco |>
  model(
    stl_ets_A = decomposition_model(
      STL(log(Tobacco), robust = TRUE),
      ETS(season_year ~ error("A") + trend("N") + season("A")),
      ETS(season_adjust ~ error("A") + trend("Ad") + season("N"))
    )
  ) |>
  forecast(h = 10)
    


españa_fut

# A fable: 5 x 4 [1Y]
# Key:     .model [1]
  .model   Year             GDP   .mean
  <chr>   <dbl>          <dist>   <dbl>
1 ets_AAN  2018 t(N(28, 0.013)) 1.35e12
2 ets_AAN  2019 t(N(28, 0.043)) 1.40e12
3 ets_AAN  2020 t(N(28, 0.099)) 1.47e12
4 ets_AAN  2021  t(N(28, 0.19)) 1.57e12
5 ets_AAN  2022  t(N(28, 0.31)) 1.70e12

tabaco_fut

# A fable: 10 x 4 [1Q]
# Key:     .model [1]
   .model    Quarter           Tobacco .mean
   <chr>       <qtr>            <dist> <dbl>
 1 stl_ets_A 2004 Q3 t(N(8.5, 0.0047)) 5090.
 2 stl_ets_A 2004 Q4 t(N(8.4, 0.0049)) 4615.
 3 stl_ets_A 2005 Q1 t(N(8.3, 0.0052)) 4108.
 4 stl_ets_A 2005 Q2 t(N(8.5, 0.0054)) 4845.
 5 stl_ets_A 2005 Q3 t(N(8.5, 0.0057)) 5039.
 6 stl_ets_A 2005 Q4  t(N(8.4, 0.006)) 4570.
 7 stl_ets_A 2006 Q1 t(N(8.3, 0.0063)) 4068.
 8 stl_ets_A 2006 Q2 t(N(8.5, 0.0067)) 4799.
 9 stl_ets_A 2006 Q3  t(N(8.5, 0.007)) 4993.
10 stl_ets_A 2006 Q4 t(N(8.4, 0.0074)) 4529.

Gráfica de los pronósticos

final_e = españa_fut |>
  autoplot(gdp_españa) +
  ggtitle("Pronóstico final para el GDP de España")

final_t = tabaco_fut |>
  autoplot(tabaco) +
  ggtitle("Pronóstico final para el Tabaco")



final_e

final_t