Oficina de R

Exemplo de previsão de séries temporais

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.

Preparando os dados

rm(list = ls())

data_1 <- read.table(
    file = "ipeadata.csv",
    row.names = 1,
    sep = ",",
    header = TRUE
)

colnames(data_1) <- "ipca"

data_1b <- data.frame(
    tail(data_1, -1),
    inf = diff(log(data_1$ipca))
)

# Restringindo a amostra (período de inflação alta)

data_2 <- subset(
    x = data_1b,
    subset = rownames(data_1b) >= "1996.01"
) # Aqui não precisamos da estrutura "Date"

dados_inf <- data_2$inf

plot.ts(dados_inf)

Medidas de performance

• Vamos considerar o RMSE.

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 = dados_inf)

y_vec <- dados_y$y

# forecasting exercise

## basic parameters

n_available <- nrow(dados_y)

for_set <- seq(
    from = floor(0.75 * n_available),
    to = n_available
)

# horizonte de previsão

for_hor <- 12 

# 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]

Exercício de previsão

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

## 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
)

Análise dos resultados

• 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 <- c(1, 3, 6, 9, 12)

# 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
            as.Date("2016-09-01"),
            as.Date("2023-09-01"),
            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(
            as.Date("2016-09-01"),
            as.Date("2023-09-01"),
            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.8926101 0.7926106 0.7384853 0.7912964 0.7393729
auto 0.9184002 0.8349397 0.7742244 0.8195137 0.7775673