Tarea 2

library(tidyverse)
library(fpp3)

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.

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

-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.

-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.

Si requieren realizar alguna transformación matemática a los datos para estabilizar la varianza, adelante.

Se puede realizar la descomposición de la serie en estacionalidad y serie desestacionalizada para intentar obtener mejores pronósticos.

Primer serie de tiempo

gas <- aus_production
gas

gas |> 
  autoplot(Gas)

Modelos de referencia

gas <- aus_production |>  filter(year(Quarter) >= 1990)
gas_train <- gas |>  filter(year(Quarter) <= 2007)

gas_fit <- gas_train |> 
  model(
    mean = MEAN(Gas),
    naive = NAIVE(Gas),
    snaive = SNAIVE(Gas),
    drift = RW(Gas ~ drift())
  )

gas_fc <- gas_fit |> 
  forecast(h = 10)

gas_fc |> 
  autoplot(gas, level = NULL)

accuracy(gas_fit) |> 
  arrange(MAPE)

gas_fc |> 
  accuracy(gas) |> 
  arrange(MAPE)

EL mejor modelo de los base fue el SNAIVE.

Descomposición de los datos de training

dcmp <- gas_train |> 
  model(
    STL(log(Gas), robust = TRUE)
  ) 

dcmp |> 
  components() |> 
  autoplot()

dcmp |> 
  components() |> 
  ggplot(aes(x = Quarter, y = season_adjust)) +
  geom_line()

gas_train |> 
  model(
    STL(Gas, robust = TRUE)
  ) |> 
  components() |> 
  autoplot()

Ajuste de los modelos

gas_fit <- gas_train |> 
  model(
    snaive = SNAIVE(Gas),
    ets_AAA = ETS(Gas ~ error("A") + trend("A") + season("A")),
    ets_MAM = ETS(Gas ~ error("M") + trend("A") + season("M")),
    ets_AAdA = ETS(Gas ~ error("A") + trend("Ad") + season("A")),
    ets_MAdM = ETS(Gas ~ error("M") + trend("Ad") + season("M")),
    ets_AAA_l = ETS(log(Gas) ~ error("A") + trend("A") + season("A")),
    ets_MAM_l = ETS(log(Gas) ~ error("M") + trend("A") + season("M")),
    ets_AAdA_l = ETS(log(Gas) ~ error("A") + trend("Ad") + season("A")),
    ets_MAdM_l = ETS(log(Gas) ~ error("M") + trend("Ad") + season("M")),
    stl_ets_A = decomposition_model(
      STL(log(Gas), robust = TRUE),
      ETS(season_year ~ error('A') + trend('N')+ season('A')),
      ETS(season_adjust~ error('A') + trend('A')+ season('N'))
    ),
    stl_ets_M = decomposition_model(
      STL(log(Gas), robust = TRUE),
      ETS(season_year ~ error('M') + trend('N')+ season('M')),
      ETS(season_adjust~ error('A') + trend('A')+ season('N'))
    )
  )
gas_fit

accuracy(gas_fit) |> 
  arrange(MAPE)

Análisis de residuos

gas_fit |> 
  augment() |> 
  features(.innov, ljung_box, lag = 8)

Pronóstico para test

gas_fc <- gas_fit |> 
  forecast(h = "2 years 6 months")

gas_fc

gas_fc |> 
  autoplot(gas, level = NULL)

gas_fc |> 
  filter(.model != "stl_ets_M") |> 
  autoplot(gas, level = NULL)

gas_fc |> 
  filter(.model != "stl_ets_M") |> 
  autoplot(gas |> filter_index("2005 Q1" ~ .)) + 
  facet_wrap(~ .model, ncol = 2) + theme(legend.position = NULL)

gas_fc |> 
  accuracy(gas) |> 
  arrange(MAPE)

Pronóstico

gas_fut <- gas |> 
  model(
    ets_AAdA = ETS(Gas ~ error("A") + trend("Ad") + season("A"))
  ) |> 
  forecast(h = "2 years 6 months")

gas_fut

gas_fut |> 
  autoplot(gas) + ggtitle("Producción de gas en Australia")

Para esta serie de tiempo sobre la producción de gas en Australia, los métodos de suavización expoenencial fueron mejores que los de referencia.

El mejor método de referencia fue el seasonal naive el cual en el test tuvo MAPE (error porcentual) de 2.94% y el mejor método de suavización exponencial tuvo un error de 1.21%. Este método considero error y estacionalidad aditiva y una tendencia amortiguada.

Segunda serie de tiempo

monaco <- global_economy |> 
  filter(Country=='Monaco') |> 
  select(-c(Country, Code))
monaco

monaco |> 
  autoplot(Population)

lambda <- monaco |> 
  features(Population, features = guerrero) |> 
  pull(lambda_guerrero)

Modelos de referencia

monaco_train <- monaco |> filter_index(.~"2006")

monaco_fit <- monaco_train |> 
  model(
    mean2 = MEAN(Population),
    naive2 = NAIVE(Population),
    drift2 = RW(Population ~ drift())
  )

monaco_fc <- monaco_fit |> 
  forecast(h = 11)

monaco_fc |> 
  autoplot(monaco, level = NULL)

accuracy(monaco_fit) |> 
  arrange(MAPE)

monaco_fc |> 
  accuracy(monaco) |> 
  arrange(MAPE)

Transformación

lambda <- monaco_train |> 
  features(Population, features = guerrero) |> 
  pull(lambda_guerrero)
lambda

[1] -0.06764032

monaco_train |> 
  autoplot(box_cox(Population, lambda))

Ajuste de los modelos

monaco_fit <- monaco_train |> 
  model(
    drift2 = RW(Population ~ drift()),
    ets_AAN = ETS(Population ~ error("A") + trend("A") + season("N")),
    ets_MAN = ETS(Population ~ error("M") + trend("A") + season("N")),
    ets_AAdN = ETS(Population ~ error("A") + trend("Ad") + season("N")),
    ets_MAdN = ETS(Population ~ error("M") + trend("Ad") + season("N")),
    
    ets_AAN_bc = ETS(box_cox(Population, lambda) ~ error("A") + trend("A") + season("N")),
    ets_MAN_bc = ETS(box_cox(Population, lambda) ~ error("M") + trend("A") + season("N")),
    ets_AAdN_bc = ETS(box_cox(Population, lambda) ~ error("A") + trend("Ad") + season("N")),
    ets_MAdN_bc = ETS(box_cox(Population, lambda) ~ error("M") + trend("Ad") + season("N")),
  )
monaco_fit

accuracy(monaco_fit) |> 
  arrange(MAPE)

Análisis de residuos

monaco_fit |> 
  augment() |> 
  features(.innov, ljung_box, lag = 12)

Pronóstico para test

monaco_fc <- monaco_fit |> 
  forecast(h = "11 years")

monaco_fc

monaco_fc |> 
  autoplot(monaco, level = NULL)

monaco_fc |> 
  autoplot(monaco |> filter_index("2005" ~ .), level = NULL) + 
  theme(legend.position = NULL)

monaco_fc |> 
  autoplot(monaco |> filter_index("2005" ~ .)) + 
  facet_wrap(~ .model, ncol = 3) + theme(legend.position = NULL)

monaco_fc |> 
  accuracy(monaco) |> 
  arrange(MAPE)

Pronóstico

monaco_fut <- monaco |> 
  model(
    ets_MAdN = ETS(Population ~ error("M") + trend("Ad") + season("N"))
  ) |> 
  forecast(h = "5 years")

monaco_fut

monaco_fut |> 
  autoplot(monaco) + ggtitle("Población de Monaco")

monaco_fut |> 
  autoplot(monaco |>filter_index("2005" ~ .)) + ggtitle("Población de Monaco")

Para esta serie de tiempo sobre la población de Monaco, los métodos de suavización exponencial también fueron mejores que los de referencia.

El mejor método de referencia fue el drift el cual en el test tuvo MAPE (error porcentual) de 3.82%. En cuanto a los métodos de suavización exponencial se escogió el tercero mejor ya que los dos mejores rechazaron la hipótesis nula, lo que significa que podrían tener error.

El método de suavización elegido tuvo un error de 1.36%. Este método considero un error multiplicativo, una tendencia amortiguada y sin estacionalidad.

Modelo ARIMA (Segunda serie de tiempo)

monaco_train |> 
  autoplot(Population)

Prueba de diferenciación

monaco_train |> 
  features(Population, unitroot_ndiffs)

Diferenciación

monaco_train2 <- monaco_train |> 
  mutate(
    "pop_diff" = difference(Population, 1)
  )
monaco_train2 |> 
  autoplot(pop_diff)

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

Análisis de grados

monaco_train2 |> 
  gg_tsdisplay(pop_diff, plot_type = "partial")

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

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

Los grados de acuerdo al análisis de la autocorrelación y autocorrelación parcial son \(p = 2\), \(d = 1\) y \(q = 3\)

Ajuste de los modelos

monaco_fit2 <- monaco_train2 |> 
  model(
    arima_213 = ARIMA(Population ~ pdq(2,1,3) + PDQ(0,0,0) + 1),
    arima_113 = ARIMA(Population ~ pdq(1,1,3) + PDQ(0,0,0) + 1),
    arima_110 = ARIMA(Population ~ pdq(1,1,0) + PDQ(0,0,0)),
    arima_011 = ARIMA(Population ~ pdq(0,1,1) + PDQ(0,0,0) + 1),
  )

monaco_fit2

monaco_fit2 |> 
  glance() |> 
  arrange(AICc)

accuracy(monaco_fit2) |> 
  arrange(MAPE)

Análisis de residuos

monaco_fit2 |> 
  select(arima_213) |> 
  gg_tsresiduals()

monaco_fit2 |> 
  select(arima_113) |> 
  gg_tsresiduals()

monaco_fit2 |> 
  augment() |> 
  features(.innov, ljung_box, lag = 10) |> 
  mutate(residuos_autocorrelacionados = if_else(lb_pvalue >=0.05, "No autocorrelación", "Autocorrelacionados"))

Pronóstico para test

monaco_fc2 <- monaco_fit2 |> 
  forecast(h = "11 years")

monaco_fc2

monaco_fc2 |> 
  autoplot(monaco, level = NULL)

monaco_fc2 |> 
  autoplot(monaco |> filter_index("2005" ~ .), level = NULL) + 
  theme(legend.position = NULL)

monaco_fc2 |> 
  autoplot(monaco |> filter_index("2005" ~ .)) + 
  facet_wrap(~ .model, ncol = 3) + theme(legend.position = NULL)

monaco_fc2 |> 
  accuracy(monaco) |> 
  arrange(MAPE)

Pronóstico

monaco_fut2 <- monaco |> 
  model(
    arima_213 = ARIMA(Population ~ pdq(2,1,3) + PDQ(0,0,0) + 1)
  ) |> 
  forecast(h = "5 years")

monaco_fut2

monaco_fut2 |> 
  autoplot(monaco) + ggtitle("Población de Monaco")

monaco_fut2 |> 
  autoplot(monaco |>filter_index("2005" ~ .)) + ggtitle("Población de Monaco")