Rolling-window forecast exercise

Exemplo de previsão de séries temporais

Previsão de Séries Temporais - Rolling-window Forecast

Exemplo: h = 1

Exemplo: h = 3

• Aqui ilustramos a opção de se manter o test set fixo. O training set precisa ir mais para o passado, à medida que o horizonte de previsão aumenta.

Previsão da Inflação Brasileira

• Vamos implementar um exercício empírico de previsão da série de taxa de inflação brasileira, medida pelo IPCA.
• Vamos coletar os dados no IPEADATA.

API Wrapper for Ipeadata

• Site para mais informações: ipeadatar

# install.packages("ipeadatar")

library(ipeadatar)
library(dplyr)

# Baixar a série do IPCA

# ipeadatar::available_series() %>% print(n = 1000)

ipca <- ipeadata("PRECOS12_IPCA12") %>%
  select(date, value) %>%
  as.data.frame() %>%
  mutate(date = as.Date(date)) %>%
  filter(date >= "1995-12-01")


data_inf <- ipca %>%
  mutate(inf = c(NA, value) %>% log() %>% diff()) %>%
  na.omit()

inf_series <- data_inf$inf

plot(
  x = data_inf$date,
  y = inf_series,
  main = "Taxa de Inflação Mensal - IPCA",
  xlab = "Data",
  ylab = "Inflação Mensal (%)",
  type = "l"
)

# Salva em CSV - opcional

write.csv(data_inf, "ipca_data.csv")

Modelos candidatos

• Vamos considerar os seguintes modelos:
  • RW (passeio aleatório);
  • AR(1) iterado;
  • Modelo escolhido pela função auto.arima, disponível no pacote forecasting.

Colocando a mão na massa

Preparando os dados

dados_y <- data.frame(y = inf_series)

y_vec <- dados_y$y

# forecasting exercise

## basic parameters

n_available <- nrow(dados_y) # tamanho da amostra

n_1st_for <- floor(0.75 * n_available) # índice do primeiro ponto de previsão

for_set <- seq( 
  from = n_1st_for,
  to = n_available
) # índices dos pontos de previsão

Exemplo

# horizonte de previsão

for_hor <- 1 # horizonte de previsão

# tamanho da amostra de treino (in-sample)

n_train <- for_set[1] - for_hor

# valores observados "não conhecidos"

y_for_real <- y_vec[for_set]

# Modelos RW

y_for_rw <- y_vec[for_set - for_hor] # última informação disponível

Exercício de previsão

• A ideia é implementar o experimento de previsão do tipo rolling-window.

library(forecast)

f_forecast <- function(for_hor) {
  # tamanho da amostra de treino (in-sample)

  n_train <- for_set[1] - for_hor

  # valores observados "não conhecidos"

  y_for_real <- y_vec[for_set]

  # Modelos RW

  y_for_rw <- y_vec[for_set - for_hor]

  ## primeira janela - training set

  train_set <- seq(from = 1, to = n_train)

  # objeto para guardar as previsões

  y_for_ar <-
    y_for_auto <-
    rep(NA, length(for_set))

  # Loop rolling-window

  for (i in seq_along(for_set)) {

    y_train <- y_vec[train_set]

    ## AR

    reg_ar_1 <- arima(
      x = y_train,
      order = c(1, 0, 0)
    )

    y_for_ar[i] <- predict(
      reg_ar_1,
      for_hor
    )$pred %>% tail(1) %>% as.numeric()

    ## Auto.arima

    reg_arima <- auto.arima(
      y = y_train,
      stepwise = FALSE,
      approximation = FALSE
    )

    for_arima_aux <- forecast(
      object = reg_arima,
      h = for_hor
    )

    y_for_auto[i] <- tail(for_arima_aux$mean, 1)

    ## atualizando os índices das janelas

    train_set <- train_set + 1

    # print(i / length(for_set))
  }

  ### Salvando os resultados

  ## results

  file_name <- paste(
    "inf_forecast_horizon_",
    for_hor,
    ".RData",
    sep = ""
  )

  ## salvando

  save(
    y_for_real,
    y_for_rw,
    y_for_ar,
    y_for_auto,
    for_hor,
    file = file_name
  )
}

Loop que itera nos horizontes de previsão

• Vamos definir os horizontes de previsão.

for_hor_can <- c(1, 3, 6, 9, 12)

• Vamos iterar sobre os horizontes de previsão.

for(for_hor in for_hor_can) {
  f_forecast(for_hor)
}

• Aqui estamos obtendo e guardando os resultados.

Análise dos resultados

Medidas de performance

• A partir do arquivo com os resultados, podemos construir uma rotina para analisar os resultados.

• Vamos nos concentrar no RMSE.

library(ggplot2)

# função para o cálculo do RMSE

f_rmse <- function(x, y) {

  res <- sqrt(cumsum((x - y)^2))

  return(res)
}

# horizontes de previsão considerados

for_hor_sel <- for_hor_can

# matriz para salvar os resultados

rmse_table <- matrix(
  data = NA,
  nrow = 2,
  ncol = length(for_hor_sel)
)

rownames(rmse_table) <- c("AR", "auto")

# loop que itera nos horizontes de previsão

for (i in seq_along(for_hor_sel)) {

  # arquivo salvo anteriormente

  file_name <- paste(
    "inf_forecast_horizon_",
    for_hor_sel[i],
    ".RData",
    sep = ""
  )
    
  load(file_name)

  # Cálculo do RMSE

  rmse_rw <- f_rmse(y_for_rw, y_for_real)
  rmse_ar <- f_rmse(y_for_ar, y_for_real)
  rmse_auto <- f_rmse(y_for_auto, y_for_real)

  # RMSE table

  rmse_table["AR", i] <- tail(rmse_ar, 1) / tail(rmse_rw, 1)
  rmse_table["auto", i] <- tail(rmse_auto, 1) / tail(rmse_rw, 1)

  # Forecast plot

  data_plot_for <- data.frame(
    x = seq( # a data poderia estar automatizada
      from = data_inf$date[n_1st_for],
      to = data_inf$date[n_available],
      by = "months"
    ),
    y_for_real,
    y_for_rw,
    y_for_ar,
    y_for_auto
  )

  g1 <- ggplot(
    data = data_plot_for
  ) +
    geom_line(
      mapping = aes(x = x, y = y_for_real, colour = "Real")
    ) +
    geom_line(
      mapping = aes(x = x, y = y_for_rw, colour = "RW")
    ) +
    geom_line(
      mapping = aes(x = x, y = y_for_ar, colour = "AR")
    ) +
    geom_line(
      mapping = aes(x = x, y = y_for_auto, colour = "auto")
    ) +
    scale_color_manual(
      values = c(
        "red",
        "blue",
        "black",
        "green"
      )
    ) +
    labs(
      title = paste(
        "h = ",
        for_hor,
        sep = ""
      )
    )

    g1$labels$colour <- c(" ")
    
    print(g1)

    # RMSE plot

    data_plot_for_g2 <- data.frame(
      x = seq( # a data poderia estar automatizada
        from = data_inf$date[n_1st_for],
        to = data_inf$date[n_available],
        by = "months"
      ),
      rmse_rw,
      rmse_ar,
      rmse_auto
    )

    g2 <- ggplot(
      data = data_plot_for_g2
    ) +
    ylim(c(0.3, 2.0)) +
    geom_line(
      mapping = aes(x = x, y = rmse_ar / rmse_rw, colour = "AR / RW")
    ) +
    geom_line(
      mapping = aes(x = x, y = rmse_auto / rmse_rw, colour = "Auto / RW")
    ) +
    geom_abline(intercept = 1, slope = 0) +
    scale_color_manual(
      values = c("blue", "red")
    ) +
    labs(
      title = paste(
        "h = ",
        for_hor,
        sep = ""
      )
    )

    g2$labels$colour <- c(" ")
    print(g2)
}

save(rmse_table, file = "rmse_table.RData")

colnames(rmse_table) <- paste(
  "h = ",
  for_hor_sel,
  sep = ""
)

print(rmse_table)

         h = 1     h = 3     h = 6     h = 9    h = 12
AR   0.8767038 0.7664647 0.6927798 0.7746988 0.7650873
auto 0.8861851 0.7770749 0.7055156 0.8001324 0.7907458