Séries M4

Neste post vou implementar algumas rotinas simples que modelam e preveem séries de tempo univariadas. Para avaliar modelos em série de tempo é bastante comum separar a série em duas partes: uma primeira é usada para alimentar o modelo: serve para estimar os parâmetros; a segunda parte serve para testar a sua capacidade preditiva. A terminologia usual é de train (treino) e test (teste). Usamos o train para estimar o modelo e depois testamos as suas previsões contra as observações reservadas no test. Vale enfatizar que o objetivo de um modelo depende muito da pergunta que se quer responder. Em alguns casos pode ser de maior interesse conseguir uma estimativa precisa dos parâmetros, com clara interpretação, do que um modelo que prevê melhor.

Para baixar os dados da competição pelo R é preciso instalar o pacote M4comp2018 que está disponível no github usando o comando install_github("carlanetto/M4comp2018") do pacote devtools. Além disso, vou usar o pacote forecast que fornece vários métodos estatísticos para previsão de séries de tempo.

library(forecast)
library(M4comp2018)
# Carregar os dados
data(M4)
M4[[2]]$period
## [1] Daily
## Levels: Daily Hourly Monthly Quarterly Weekly Yearly

Os dados estão organizados numa lista, em que cada elemento é um conjunto que inclui:

  • Nome
  • Uma série de tempo separada em train e test
  • A categoria da série (há dados de finanças, macroeconomia, microeconomia, indústria além de dados demográficos e outros).
  • A periodicidade da série (anual, trimestral, mensal, semanal, diária e horária)
str(M4[[2]])
## List of 7
##  $ st    : chr "D2"
##  $ x     : Time-Series [1:1006] from 9132 to 10137: 2794 2794 2804 2806 2802 ...
##  $ n     : int 1006
##  $ type  : Factor w/ 6 levels "Demographic",..: 4
##  $ h     : int 14
##  $ period: Factor w/ 6 levels "Daily","Hourly",..: 1
##  $ xx    : Time-Series [1:14] from 10138 to 10151: 2986 3001 2976 2996 2982 ...

Para construir uma série espécifica pode-se usar o código abaixo

serie <- ts(M4[[2]]$x,
            start = start(M4[[2]]$x),
            frequency = frequency(M4[[2]]$x))
# Gráfico usando autoplot
autoplot(serie) +
  theme_vini

Neste post vou focar apenas nas séries de macroeconomia. A ideia, como exposto acima, é modelar os dados usando alguns modelos automáticos e testar a capacidade preditiva. Vou primeiro trabalhar com uma única série.

# séries de macro
macro_m4 <- Filter(function(l) l$type == "Macro", M4)

serie <- ts(macro_m4[[1]]$x, start = start(macro_m4[[1]]$x,
            frequency = frequency(macro_m4[[1]]$x)))
autoplot(serie) +
  theme_vini

A série, dividida em train e test é plotada abaixo.

test <- ts(macro_m4[[1]]$xx, start = start(macro_m4[[1]]$xx),
           frequency = frequency(macro_m4[[1]]$xx))

autoplot(serie) + 
  autolayer(test) +
  theme_vini

#plot(ts(c(serie, test), start = start(macro_m4[[1]]$x)))
#lines(test, col = "red")

Os modelos

Serão usados os seguintes modelos

  • SARIMA: a ordem do modelo é escolhida automaticamente usando a função auto.arima. Além de minimizar critérios de informação a função tem várias checagens para verificar a qualidade dos resíduos e também para evitar sobre-parametrização.
  • ETS: modelo de suavização exponencial com sazonalidade e tendência; a forma específica do modelo é escolhida automaticamente pela função ets()
  • Holt-Winters: moodelo Holt-Winters com
  • SES: modelo de suavização exponencial simples
  • Naive: modelo ingênuo (random-walk com tendência)
  • H1: média simples entre o SARIMA e o ETS
  • H2: média simples entre o SARIMA, ETS e Holt-Winters
  • H3: média simples entre todos os modelos exceto o Naive
# Modelos
h <- length(test)
(m1 <- auto.arima(serie))
## Series: serie 
## ARIMA(0,1,0) with drift 
## 
## Coefficients:
##        drift
##       1.0076
## s.e.  0.2285
## 
## sigma^2 estimated as 52.53:  log likelihood=-3416.09
## AIC=6836.18   AICc=6836.19   BIC=6846
(m2 <- ets(serie))
## ETS(A,A,N) 
## 
## Call:
##  ets(y = serie) 
## 
##   Smoothing parameters:
##     alpha = 0.9697 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 1020.3793 
##     b = 1.0011 
## 
##   sigma:  7.2531
## 
##      AIC     AICc      BIC 
## 10947.86 10947.92 10972.43
m3 <- holt(serie, h)
m4 <- holt(serie, damped = TRUE, h)
m5 <- ses(serie, h)
m6 <- naive(serie, h)

Construo um data.frame com as previsões de todos os modelos

previsao <- data.frame(p1 = forecast(m1, h)$mean,
                       p2 = forecast(m2, h)$mean,
                       p3 = m3$mean,
                       p4 = m4$mean,
                       p5 = m5$mean,
                       p6 = m6$mean
                       )
previsao$p7 <- rowMeans(previsao[, 1:2])
previsao$p8 <- rowMeans(previsao[, 1:3])
previsao$p9 <- rowMeans(previsao[, 1:5])

Computo o erro de previsão segundo algumas medidas convencionais. Por simplicidade vou escolher o melhor modelo segundo a raiz do erro quadrado médio.

metodos <- 9
erros <- matrix(nrow = metodos, ncol = 5)
scale <- mean(abs(diff(macro_m4[[1]]$xx, lag = frequency(macro_m4[[1]]$xx))))
for(i in 1:metodos){
    erros[i, 1] <- mean(abs(previsao[, i] - test))
    erros[i, 2] <- sqrt(mean((previsao[, i] - test)^2))
    erros[i, 3] <- mean(abs( (test - previsao[, i]) / test) ) * 100
    erros[i, 4] <- mean(abs( (test - previsao[, i]) / (abs(test) + abs(previsao[, i]) )))
    erros[i, 5] <- mean(abs(test - previsao[, i]) / scale)
}
t(apply(erros, 2, which.min))
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    3    3    3    3    3

Aqui o modelo 3 (Holt-Winters) teve a melhor performance. Abaixo mostro a previsão contra o observado.

autoplot(forecast(m3, h), include = 50) +
  autolayer(test) +
  theme_vini

Desempenho dos modelos

MAE RMSE MAPE sMAPE MASE
SARIMA 24.47186 26.42446 1.183958 0.0059608 3.648327
ETS 24.43099 26.39245 1.181970 0.0059508 3.642235
Holt 24.42801 26.38925 1.181826 0.0059501 3.641791
Holt (damped) 27.87388 30.15279 1.348443 0.0067957 4.155510
SES 32.01123 34.64532 1.548563 0.0078135 4.772316
Naive 32.02857 34.66134 1.549404 0.0078178 4.774902
H1 24.45142 26.40845 1.182964 0.0059558 3.645281
H2 24.44362 26.40205 1.182585 0.0059539 3.644118
H3 26.64319 28.78481 1.288952 0.0064935 3.972036

Loop para várias séries

Abaixo proponho uma forma simples de repetir o procedimento acima para várias séries. Por simplicidade vou considerar apenas o MSQE o MASE. \[\begin{align} \text{MSQE} & = \sqrt(\frac{1}{T}\sum_{i = 1}^{h}(y_{T + 1} - \hat{y_{T + 1}})^2)\\ \text{MASE} & = \frac{1}{T}\sum_{i = 1}^{h}\frac{\|y_{T + 1} - \hat{y_{T + 1}}\|}{\sigma} \\ \sigma & = \frac{1}{T + m} \sum^{T}_{m + 1}\| y_{t} - y_{t - m}\| \end{align}\] em que \(m\) é a frequência sazonal da série (sendo igual a um na ausência de sazonalidade). Essencialmente, o MASE é uma medida de erro normalizada pela previsão de modelo ingênuo sazonal que prevê sempre a média amostral do processo. Neste sentido, um bom modelo deve ter valores de MASE menores do que um.

Note que o loop abaixo é bastante ineficiente, mas resole o trabalho. Num outro post discuto como melhorá-lo usando processamento paralelo e funções mais eficientes.

N <- length(macro_m4)
melhor_modelo <- matrix(ncol = 5, nrow = N)
merros <- matrix(ncol = 4, nrow = N * metodos)
modelos <- c("SARIMA", "ETS", "HW", "HW (damped)", "SES", "Naive",
             "SARIMA + ETS", "SARIMA + ETS + HW", "Todos (-Naive)")

for(n in 1:N) {
    train <- ts(macro_m4[[n]]$x,
                start = start(macro_m4[[n]]$x),
                frequency = frequency(macro_m4[[n]]$x)
                )
    test <- as.numeric(macro_m4[[n]]$xx)
    h <- length(test)
    p <- data.frame(p1 = forecast(auto.arima(train, max.d = 2, max.D = 1), h)$mean,
                    p2 = forecast(ets(train), h)$mean,
                    p3 = holt(train, h)$mean,
                    p4 = holt(train, h, damped = TRUE)$mean,
                    p5 = ses(train, h)$mean,
                    p6 = naive(train, h)$mean
                    )
    p$p7 <- rowMeans(p[, 1:2])
    p$p8 <- rowMeans(p[, 1:3])
    p$p9 <- rowMeans(p[, 1:5])

    scale <- mean(abs(diff(train, lag = frequency(train))))
    e <- apply(p, 2, function(x) x - test)
    erro1 <- apply(e, 2, function(x) mean(abs(x)))
    erro2 <- apply(e, 2, function(x) sqrt(mean(x^2)))
    erro3 <- apply(e, 2, function(x) mean(abs(100 * x / test)))
    erro5 <- apply(e, 2, function(x) mean(abs(x / scale)))

    erros_at <- cbind(n, modelos, erro2, erro5)
    merros[1:9 + 9 * (n - 1), 1:4] <- erros_at
}

merros <- as.data.frame(merros)
names(merros) <- c("serie", "modelo", "reqm", "mase")
merros$reqm <- as.numeric(levels(x$reqm))[x$reqm]
merros$mase <- as.numeric(levels(x$mase))[x$mase]

Resultados

Os modelos mistos, como esperado, têm um desempenho bastante razoável.

library(extrafont)
library(ggplot2)
library(dplyr)

merros %>%
  ggplot(aes(mase)) +
  geom_histogram() +
  labs(x = "", y = "", title = "Erro absoluto médio escalado") +
  theme_vini

merros %>%
  ggplot(aes(reqm)) +
  geom_histogram() +
  labs(x = "", y = "", title = "Raiz do erro quadrado médio") +
  theme_vini

# reqm
merros %>%
  group_by(modelo) %>%
  summarise(mase_medio = mean(mase)) %>%
  mutate(modelo = reorder(modelo, mase_medio)) %>%
  ggplot(aes(x = mase_medio, y = modelo)) +
  geom_point() +
  geom_segment(aes(yend = modelo), xend = 0)

# mase
merros %>%
  group_by(modelo) %>%
  summarise(reqm_medio = mean(reqm)) %>%
  mutate(modelo = reorder(modelo, reqm_medio)) %>%
  ggplot(aes(x = reqm_medio, y = modelo)) +
  geom_point() +
  geom_segment(aes(yend = modelo), xend = 0)

Exemplo de bom desempenho

Em geral os modelos univaridos têm uma performance bastante razoável, especialmente a média das previsões de diferentes modelos. Séries relativamente contínuas e sem quebras estruturais são particularmente fáceis de prever. A série abaixo fornece um exemplo

n <- 16044
serie_mape <- ts(macro_m4[[n]]$x, start = start(macro_m4[[n]]$x),
                 frequency = frequency(macro_m4[[n]]$x))
teste_mape <- ts(macro_m4[[n]]$xx, start = start(macro_m4[[n]]$xx),
                 frequency = frequency(macro_m4[[n]]$xx))
autoplot(serie_mape) + autolayer(teste_mape) +
  theme_vini +
  labs(x = "", y = "", title = "Série 16044") +
  scale_colour_discrete(name = "", labels = "Test")

No gráfico abaixo mostro a previsão dos melhores modelos univariados considerados.

m1 <- auto.arima(serie_mape)
m2 <- ets(serie_mape)
p1 <- forecast(m1, h = length(teste_mape))$mean
p2 <- forecast(m2, h = length(teste_mape))$mean
p3 <- holt(serie_mape, damped = TRUE, h = length(teste_mape))$mean
p <- data.frame(p1, p2, p3)
p4 <- ts(rowMeans(p[, 1:2]),
         start = start(teste_mape),
         frequency = frequency(teste_mape))
p5 <- ts(rowMeans(p[, 1:3]),
         start = start(teste_mape),
         frequency = frequency(teste_mape))
nomes_modelos <- c("SARIMA", "ETS", "Holt", "SARIMA + ETS", "SARIMA + HW + ETS")
autoplot(teste_mape) +
  autolayer(p[, 1]) +
  autolayer(p[, 2]) +
  autolayer(p[, 3]) +
  autolayer(p4) +
  autolayer(p5) +
  autolayer(teste_mape) +
  theme_vini +
  scale_colour_discrete(name = "", labels = c(nomes_modelos, "Observado"))

Neste caso específico o modelo ARIMA (3, 2, 1) teve a melhor performance segundo o REQM.

autoplot(teste_mape) +
  autolayer(p[, 1]) +
  autolayer(teste_mape) +
  theme_vini +
  labs(title = "Previsão do modelo ARIMA (3, 2, 1)",
       x = "", y = "") + 
  scale_colour_manual(name = "", labels = c("Previsto", "Observado"),
                      values = c("indianred", "black"))

Exemplo de mau desempenho

Como era de se esperar estes modelos lineares não são apropriados para séries com quebras estruturais ou descontinuidades repentinas. Tome como exemplo a série abaixo que tem uma explosão nas últimas observações.

n <- 3133
serie_mape <- ts(macro_m4[[n]]$x, start = start(macro_m4[[n]]$x),
                 frequency = frequency(macro_m4[[n]]$x))
teste_mape <- ts(macro_m4[[n]]$xx, start = start(macro_m4[[n]]$xx),
                 frequency = frequency(macro_m4[[n]]$xx))
autoplot(serie_mape) +
  autolayer(teste_mape) +
  labs(title = "Série 3133", x = "", y = "") +
  theme_vini

A série fica estável durante um período de mais de cinco anos até explodir em janeiro de 2008. A previsão de um modelo autoregressivo, por exemplo, coloca um peso nas observações recentes (todas constantes) e prevê que a série vai continuar nesta tendência.

autoplot(ses(serie_mape, h = length(teste_mape)), include = 50) +
  autolayer(teste_mape) +
  labs(x = "", y = "") +
  theme_vini