Combinação de preditores de séries temporais

1 Introdução

A previsão de séries temporais fornece informações para a tomada de decisão baseada em dados em diversas áreas do conhecimento. Nesse contexto, a obtenção de estimativas precisas pelos modelos de séries temporais é importante. Para melhorar o desempenho dos modelos de previsão, estratégias de modelagem combinada são comumente usadas. A combinação pode aumentar a qualidade dos modelos e reduzir a variância das previsões. Uma alternativa de combinação consiste em se combinar as previsões de modelos individuais de uma mesma série temporal [1].

Neste texto, diferentes técnicas de combinação de preditores de séries temporais serão avaliadas. Como caso de estudo, uma série temporal conhecida será utilizada. Para esta aplicação, os modelos individuais e as técnicas de combinação usadas serão apresentados brevemente.

2 Métodos

2.1 Modelos ARIMA

Os modelos autorregressivos (AR) usam os valores passados de uma série temporal e os modelos de médias móveis (MA) usam os erros de previsões passadas para realizar previsões. Pode-se combinar modelos AR e MA, dando origem aos modelos autorregressivos de médias móveis (ARMA). Para a aplicação dos modelos ARMA, a série temporal deve ser estacionária, isto é, as propriedades estatísticas que descrevem seu comportamento devem ser constantes ao longo do tempo. Para atingir a estacionaridade da série, pode-se usar a diferenciação, ou seja, tomar diferenças entre os valores da série, por exemplo \(y_t-y_{t-1}\). Obtida uma série estacionária, pode-se usar os modelos autorregressivos integrados de médias móveis (ARIMA). O modelo ARIMA não sazonal é mostrado abaixo:

\[y'_{t} = c + \phi_1 y'_{t-1} + ... + \phi_p y'_{t-p} + \theta_1 \varepsilon_{t-1} + ... + \theta_q \varepsilon_{t-q} + \varepsilon_t\]

em que \(y'_t\) é a série diferenciada, \(\varepsilon_t\), \(\varepsilon_{t-1}\), …, \(\varepsilon_{t-q}\) são ruídos brancos, \(\theta_i\) e \(\phi_i\) são coeficientes e \(c\) é a constante da equação de regressão.

O modelo ARIMA pode ser indicado \(ARIMA(p,d,q)\), em que \(p\) é a ordem da parcela autorregressiva, \(d\) é o número de diferenciações, \(q\) é a ordem da parcela de médias móveis. Os modelos ARIMA também podem ser usados para séries temporais que possuam sazonalidade, incluindo-se termos sazonais adicionais ao modelo. O modelo ARIMA sazonal pode ser indicado por \(ARIMA(p,d,q)(P,D,Q)[m]\), em que \(P\), \(D\) e \(Q\) são as ordens da parcela sazonal, e \(m\) é a sazonalidade da série [2].

2.2 Modelos ETS

Os modelos ETS são modelos de previsão referentes aos métodos de suavização exponencial simples, de Holt e Holt-Winters. A tríade (E,T,S) representa as componentes de erro, tendência e sazonalidade, respectivamente. A notação ETS(\(\cdot,\cdot,\cdot\)) indica a manifestação de cada componente, podendo ser aditiva (letra A), multiplicativa (letra M), ou sem determinada componente (letra N). Existem diversas opções de modelos ETS e pode-se representá-los usando o conceito de Modelos de Estados de Espaço:

\[\begin{align} y_{t} & = \psi \dot{\textbf{x}}_{t-1} + \varepsilon_t\label{eq:state-space-model-a}\\ \dot{\textbf{x}}_t & = F \dot{\textbf{x}}_{t-1} + g \varepsilon_t\label{eq:state-space-model-b} \end{align}\]

em que \(y_t\) é o valor da série no instante de tempo \(t\), \({\varepsilon_t}\) é um ruído branco e \(F\), \(g\) e \(\psi\) são coeficientes do modelo. \(\dot{\textbf{x}}_t\) é chamado vetor de estado e contém as componentes do modelo A primeira equação descreve a relação entre os valores da série temporal e as componentes, e a segunda equação expressa a evolução temporal das componentes [2].

2.3 Combinação de preditores

A combinação leva em consideração as estimativas de \(k\) modelos individuais \(y_{k,t}\) para um determinado valor a que se deseja prever, \(y_t\), no instante de tempo \(t\). No geral, a combinação pode ser expressa da seguinte maneira:

\[\hat{y}_t=\varphi(y_{1,t},y_{2,t},\cdots,y_{i,t},\cdots,y_{k,t})\] em que \(y_{i,t}\) é a previsão do \(iº\) modelo para o valor da série no instante de tempo \(t\) e \(\varphi\) é uma função de combinação [3].

2.3.1 Média

Neste método, a média aritmética simples (SA) dos preditores para cada observação é calculada. Desse modo, todos os preditores recebem o mesmo peso [4].

2.3.2 Mediana

A mediana (MED) dos preditores para cada observação é calculada neste método. Os pesos atribuídos variam com o tempo, já que a cada instante de tempo, o resultado da combinação será a mediana dos preditores individuais [3].

2.3.3 Método dos Mínimos Quadrados

Para a combinação de preditores, o método dos mínimos quadrados (MMQ) é usado para estimar um modelo de regressão linear múltipla. O conjunto de treinamento é usado para construir o modelo, de modo que a série temporal original é a variável dependente e os valores ajustados pelos modelos são as variáveis independentes. Assim, os pesos dos preditores são estimados pelo modelo de regressão construído, para se obter a combinação no conjunto de teste [4].

2.4 Avaliação do desempenho

Existem várias medidas para avaliar o desempenho de modelos de previsão. O objetivo é medir o erro cometido pelo método ao realizar previsões. Neste texto, serão usados a Raiz do Erro Quadrático Médio (RMSE), ou Root Mean Squared Error, em inglês, e o Erro Absoluto Médio (MAE), ou Mean Absolute Error, em inglês. O RMSE avalia a raiz da média dos erros quadráticos e o MAE calcula a média dos erros absolutos. Quanto menores os seus valores, melhor o desempenho do modelo. O RMSE e o MAE podem ser calculados conforme mostrado a seguir.

\[RMSE=\sqrt{\frac{\sum^{n}_{t=1}{(\hat{y}_t-y_t)^2}}{n}}\] \[MAE=\frac{\sum^{n}_{t=1}{|\hat{y}_t-y_t|}}{n}\]

em que \(\hat{y}_t\) é o valor previsto por um modelo, \(y_t\) é o respectivo valor original da série no instante de tempo \(t\) e \(n\) é o número de valores previstos.

3 Caso de estudo

3.1 Conjunto de dados

Neste texto, foi usado uma série temporal conhecida e disponível no próprio R, a série ldeaths. São dados mensais do número de mortes decorrentes de doenças no pulmão registradas no Reino Unido, de 1974 a 1980. A série possui 72 observações. Para usá-la, basta carregar o conjunto de dados ao qual ela pertence UKLungDeaths usando a função data.

# Carrega o conjunto de dados da série ldeaths
data(UKLungDeaths)

# Visualiza a série
plot(ldeaths, type = 'o', pch = 20)

3.2 Pré-processamento dos dados

A série ldeaths foi separada em conjunto de treinamento e teste. Os modelos foram construídos usando o conjunto de treinamento e destinou-se o conjunto de teste para realizar previsões e avaliar o desempenho dos modelos. O conjunto de treino correspondeu aos 70% iniciais da série e o conjunto de teste foi os 30% finais da série.

# Tamanho da série temporal
n <- length(ldeaths)

# Proporção dos conjuntos de treinamento e teste
p.treino <- 0.7
p.teste <- 0.3

# Tamanho dos conjuntos de treinamento e teste
n.treino <- round(n * p.treino)
n.teste <- round(n * p.teste)

# Separação dos dados
ldeaths.treino <- head(ldeaths, n.treino)
ldeaths.teste <- tail(ldeaths, n.teste)

A tabela abaixo mostra o tamanho dos conjuntos de treinamento e teste para a série ldeaths.

##  ldeaths treino teste
##       72     50    22

Após a separação, os dados de treinamento e teste foram convertidos em objetos da classe ts, para possuírem atributos de séries temporais. Em R, é importante se trabalhar com objetos ts, pois os pacotes muitas vezes recorrem à atributos desta classe, como frequência e sazonalidade, para a construção dos modelos.

# Salva como série temporal
ldeaths.treino <- ts(ldeaths.treino, frequency = frequency(ldeaths), start = start(ldeaths))
ldeaths.teste <- ts(ldeaths.teste, frequency = frequency(ldeaths), end = end(ldeaths))

3.3 Construção dos modelos

Os modelos foram construídos usando o pacote forecast [5]. Primeiramente, deve-se instalar e carregar o pacote no R.

install.packages(forecast)
library(forecast)

Em seguida, os modelos ARIMA e ETS foram construídos usando respectivamente as funções auto.arima e ets.

# Constrói os modelos ARIMA e ETS
mod.arima <- auto.arima(ldeaths.treino)
mod.ets <- ets(ldeaths.treino)

Para verificar os modelos construídos, geralmente pode-se usar a função summary, ou, no caso dos modelos do pacote forecast, a função print também funciona, trazendo uma descrição mais sucinta dos modelos.

print(mod.arima)

## Series: ldeaths.treino 
## ARIMA(0,0,2)(1,1,0)[12] with drift 
## 
## Coefficients:
##         ma1     ma2    sar1   drift
##       0.195  -0.340  -0.569  -6.505
## s.e.  0.178   0.195   0.123   2.647
## 
## sigma^2 estimated as 105618:  log likelihood=-274.1
## AIC=558.2   AICc=560.1   BIC=566.4

print(mod.ets)

## ETS(M,N,M) 
## 
## Call:
##  ets(y = ldeaths.treino) 
## 
##   Smoothing parameters:
##     alpha = 0.1374 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 2196.0248 
##     s = 1.244 0.9391 0.8361 0.6724 0.678 0.7328
##            0.7553 0.8705 1.132 1.306 1.445 1.389
## 
##   sigma:  0.1063
## 
##   AIC  AICc   BIC 
## 749.7 763.8 778.4

Foi selecionado um modelo ARIMA sazonal com drift, usando um termo autorregressivo sazonal, uma diferenciação sazonal e dois termos de médias móveis não sazonal. O modelo ETS usou componentes de erro e sazonalidade multiplicativas, e sem componente de tendência.

3.4 Previsão

Os modelos construídos foram usados para realizar previsões no conjunto de teste. O objetivo é avaliar a capacidade dos modelos de realizar previsões em um intervalo desconhecido, que não tinha sido usado ainda. Para os modelos do pacote forecast, foi usada a função de mesmo nome, passando o argumento h para indicar o horizonte de previsões, isto é, o número de valores a serem previstos. O horizonte é igual ao tamanho do conjunto de teste, n.teste.

# Realiza previsões dos modelos ARIMA e ETS
fct.arima <- forecast(mod.arima, h = n.teste)
fct.ets <- forecast(mod.ets, h = n.teste)

A figura abaixo mostra o gráfico das previsões dos modelos ARIMA e ETS no conjunto de teste.

Por fim, as previsões dos modelos foram armazenadas em um data.frame, para facilitar a combinação em R.

# Tabela com as previsões no conjunto de teste
fct.teste <- data.frame(cbind(arima = fct.arima$mean, ets = fct.ets$mean))

3.5 Combinação

Após serem realizadas as previsões, foi possível combinar os modelos. A média artimética simples, mediana e o método dos mínimos quadrados foram implementados através das seguintes funções:

# Função para a combinação pela média simples
SA <- function(forecasts) {
  return(apply(forecasts, 1, mean))
}

# Função para a combinação pela mediana
MED <- function(forecasts) {
  return(apply(forecasts, 1, median))
}

# Funções para a combinação pelo Método dos Mínimos Quadrados
MMQ <- function(target, fitted) {
  col.names <- colnames(fitted)
  df <- cbind(target, fitted)
  colnames(df) <- c("target", col.names)
  mod <- lm(formula = target~., data = df)
  return(mod)
}
MMQ.forecast <- function(model, forecasts) {
  return(predict(model, newdata = forecasts))
}

Em seguida, as funções foram aplicadas para realizar as combinações.

# Média simples
fct.mean <- SA(fct.teste)

# Mediana
fct.median <- MED(fct.teste)

# Método dos Mínimos Quadrados
mod.mmq <- MMQ(ldeaths.treino, cbind(arima = mod.arima$fitted, ets = mod.ets$fitted))
fct.mmq <- MMQ.forecast(mod.mmq, fct.teste)

3.6 Avaliação do desempenho

A partir dos valores originais da série ldeaths no conjunto de testes e dos valores previstos por cada modelo de previsão, foi possível avaliar o desempenho dos modelos usando o RMSE e o MAE. Para calcular o RMSE e MAE, foram usadas as seguintes funções:

# Função para calcular o RMSE
rmse <- function(previsto, original) {
  return(sqrt(mean((as.numeric(previsto) - as.numeric(original))^2)))
}

# Função para calcular o MAE
mae <- function(previsto, original) {
  return(mean(abs((as.numeric(previsto) - as.numeric(original)))))
}

Em seguida, as funções foram usadas para avaliar o desempenho dos modelos individuais e combinados.

# Calcula o RMSE
rmse.arima <- rmse(fct.arima$mean, ldeaths.teste)
rmse.ets <- rmse(fct.ets$mean, ldeaths.teste)
rmse.mean <- rmse(fct.mean, ldeaths.teste)
rmse.median <- rmse(fct.median, ldeaths.teste)
rmse.mmq <- rmse(fct.mmq, ldeaths.teste)

# Calcula o MAE
mae.arima <- mae(fct.arima$mean, ldeaths.teste)
mae.ets <- mae(fct.ets$mean, ldeaths.teste)
mae.mean <- mae(fct.mean, ldeaths.teste)
mae.median <- mae(fct.median, ldeaths.teste)
mae.mmq <- mae(fct.mmq, ldeaths.teste)

Os resultados do RMSE e MAE são mostrados na tabela abaixo.

##      arima   ets  mean median   mmq
## RMSE 169.4 190.5 157.2  157.2 162.8
## MAE  149.1 126.7 118.9  118.9 120.4

Comparando os dois modelos individuais, o ARIMA obteve RMSE inferior em comparação ao ETS. Por outro lado, o ETS atingiu MAE inferior em comparação ao ARIMA. Segundo as medidas de erro usadas, não houve diferenças significativas entre os dois modelos. Em relação aos modelos combinados, a média e a mediana atingiram resultados iguais, o que era esperado, visto que foram combinados apenas dois modelos individuais. Dito isto, a média e a mediana obtiveram valores de RMSE e MAE inferiores em relação ao MMQ. Portanto, a média e mediana foram melhores com comparação ao método dos mínimos quadrados. Ao se comparar os resultados dos modelos individuais e combinados, as três alternativas de combinação usadas obtiveram valores de RMSE e MAE inferiores aos modelos individuais. Segundo estes resultados, para a série estudada e o horizonte de previsões adotado, as combinações foram superiores em relação modelos individuais.

4 Considerações finais

Este texto avaliou a aplicação de modelos esemble para previsão de séries temporais em linguagem R. Uma série temporal conhecida foi usada como caso de estudo. A combinação mostrou desempenho superior em relação aos modelos individuais usados para o caso de estudo adotado.

5 Referências

[1]

J. M. Bates and C. W. Granger, “The combination of forecasts,” Journal of the Operational Research Society, vol. 20, no. 4, pp. 451–468, 1969, doi: https://doi.org/10.1057/jors.1969.103.

[2]

R. J. Hyndman and G. Athanasopoulos, Forecasting: Principles and practice. Melbourne, Australia: OTexts, 2018. Available: https://www.OTexts.com/fpp2

[3]

Z. Hajirahimi and M. Khashei, “Hybrid structures in time series modeling and forecasting: A review,” Engineering Applications of Artificial Intelligence, vol. 86, pp. 83–106, 2019, doi: https://doi.org/10.1016/j.engappai.2019.08.018.

[4]

C. Aksu and S. I. Gunter, “An empirical analysis of the accuracy of SA, OLS, ERLS and NRLS combination forecasts,” International Journal of Forecasting, vol. 8, no. 1, pp. 27–43, 1992, doi: https://doi.org/10.1016/0169-2070(92)90005-T.

[5]

R. J. Hyndman and Y. Khandakar, “Automatic time series forecasting: The forecast package for R,” Journal of Statistical Software, vol. 26, no. 3, pp. 1–22, 2008, Available: https://www.jstatsoft.org/article/view/v027i03