Introdução

Você certamente já leu um pouco sobre Machine Learning nas suas redes sociais ou escutou no trabalho, certo? Ela tem impactado diretamente diversas áreas, com pesquisadores tentando melhorar desempenho de algumas tarefas e conquistar novas descobertas. A área de séries tempo não ficou impune, sobretudo quando o objetivo é fazer previsões.

Vamos apresentar aqui alguns comandos de pacotes disponíveis no R, sem entrar na teoria e na formalidade dos modelos, mas deixando uma referência apropriada para tal. Nosso objetivo aqui é bem direto e foca na apresentação dessas novas possibilidades, sobretudo os códigos.

Como série candidata a cobaia, vamos trabalhar com a série de Produção Industrial Brasileira (PIM-PF), divulgada todo mês pelo IBGE e amplamente acompanhada pelo mercado financeiro.

Vamos iniciar baixando diretamente do IPEA e transformar em série de tempo, indo até setembro de 2018 enquanto escrevo. Destaque para os pacotes forecast e nnfor.

library(ipeaData)
library(forecast)
library(nnfor)

h<-search_serie('produção') # para procurar o código da indústria
f<-ipeadata('PIMPFN12_QIIGN12')
## [1] "Instituto Brasileiro de Geografia e Estatística, Pesquisa Industrial Mensal - Produção Física (IBGE/PIM-PF). 01/2002 a 09/2018. Acesso em: 29/11/2018"
g<-f$VALVALOR[1:201]
pim<-ts(g,start=2002,frequency = 12)
plot(pim,xlab='Anos',ylab='',main='Produção Industrial Brasileira (PIM-PF)',sub='Fonte: IBGE, Elaboração: Arthur Lula Mota',col='blue')

Critério de Seleção de Modelo

Antes de entrar nos modelos, precisamos estabelecer um critério para saber se o modelo é bom ou não, ou seja, se ele acerta bem as projeções. Existe diversas formas de compara, sobretudo com o termo de erro entre a série original e a projetada. Vamos escolher o critério chamado MASE (mean absolute scaled error), proposto por Hyndman e Koehler (2006), que compararam diversos critérios de acurácia de modelos e sugeriram esse como o melhor.

\[MASE=\frac{1}{T}\displaystyle\sum_{t=1}^{T} \Bigg(\frac{|e_t|}{\sum_{i=2}^{T}|Y_t-Y_{t-1}|}\Bigg)\]

calculateMASE <- function(f,y) { 
  if(length(f)!=length(y)){ stop("Vector length is not equal") }
  n <- length(f)
  return(mean(abs((y - f) / ((1/(n-1)) * sum(abs(y[2:n]-y[1:n-1]))))))
}

Multilayer Perceptron (MLP)

Vamos começar a modelar com o Multilayer Perceptron (MLP), ou Perceptron multicamadas em português. A parte teórica do modelo pode ser encontrada em Gardner e Dorling (1998), trabalhando com redes neurais do tipo Feed-forward (que diferem daquelas do tipo Feedback).

Vamos utilizar a função mlp do pacakge nnfor para fazer duas coisas: a projeção de outubro de 2017 até setembro de 2018, ou seja, 12 meses, e comparar com o que de fato ocorreu. Além disso, vamos fazer a projeção de outubro de 2018 em diante, algo que ainda não aconteceu de fato, apenas para mostrar a ferramenta e como vocês podem elaborar projeções e não apenas checar se o modelo acerta o passado.

As linhas azuis são as projeções.

par(mfrow=c(2,1))
#projeção entre out17 e set 18
mlp.fit <- mlp(ts(pim[-c(190:201)],start=2002,frequency = 12),hd=NULL,
               hd.auto.type=TRUE,lags=12,reps=20,difforder=c(1,12),sel.lag=TRUE)
mlp.frc <- forecast(mlp.fit,h=12)
plot(mlp.frc)
lines(pim)

## modelo de out18 em diante
mlp.fit <- mlp(ts(pim,start=2002,frequency = 12),hd=NULL,
               hd.auto.type=TRUE,lags=12,reps=20,difforder=c(1,12),sel.lag=TRUE)
mlp.frc <- forecast(mlp.fit,h=12)
plot(mlp.frc)
lines(pim)

Vamos calcular o erro agora da parte projetada até setembro de 2018, comparando o que efetivamente ocorreu.

calculateMASE(pim[190:201],mlp.frc$mean)
## [1] 0.6348726323

O ajuste ficou muito bom e o erro foi relativamente pequeno, ficando inclusive difícil de ver no primeiro gráfico os erros (diferença entre linha azul e preta). Vamos agorar olhar outros modelos.

TBATS (Exponential smoothing state space model with Box-Cox transformation, ARMA errors, Trend and Seasonal components)

Esse modelo é uma mistura de muito outros modelos clássicos já conhecidos, além de transformações Box-Cox, com o objetivo de projetar séries de tempo com sazonalidade. Para uma melhor descrição veja Livera et al.(2011).

Vamos elaborar a projeção para os dois períodos, conforme feito no modelo MLP e comparar o indicador de acurácia.

par(mfrow=c(2,1))
#projeção entre out17 e set 18
pimt<-tbats(ts(pim[-c(190:201)],start=2002,frequency = 12),use.box.cox=TRUE,
            use.trend=TRUE,use.arma.errors = TRUE,
            seasonal.periods = 12,biasadj = FALSE)
f1<-forecast(pimt,h=12)
plot(f1)
lines(pim)

## modelo de out18 em diante
pimt<-tbats(ts(pim,start=2002,frequency = 12),use.box.cox=TRUE,
            use.trend=TRUE,use.arma.errors = TRUE,
            seasonal.periods = 12,biasadj = FALSE)
f1<-forecast(pimt,h=12)
plot(f1)
lines(pim)

Vamos verificar nossa medida de acurácia:

calculateMASE(pim[190:201],f1$mean)
## [1] 0.4882311979

Percebam que para essa série, o resultado foi melhor (até agora). Vamos para o último modelo.

Autoregressive Neural Network (AR-NN, ou NNAR)

Esse modelo foi aperfeiçoado por Thielbar e Dickey (2011), e trabalha com funções lineares e não lineares, conjuntamente com o ferramental de redes neurais, para treinar e testar as séries de tempo desejadas. Vamos seguir com nossos modelos:

par(mfrow=c(2,1))
#projeção entre out17 e set 18
fit <- nnetar(ts(pim[-c(190:201)],start=2002,frequency = 12),lambda=0.6,maxit=150) #com a série menor, se for maior lambda0.3
fcast <- forecast(fit,PI=TRUE,h=12)
plot(fcast)
lines(pim)


## modelo de out18 em diante
fit <- nnetar(ts(pim,start=2002,frequency = 12),lambda=0.6,maxit=150)
fcast <- forecast(fit,PI=TRUE,h=12)
plot(fcast)
lines(pim)

e seguindo com o cálculo do MASE

calculateMASE(pim[190:201],fcast$mean)
## [1] 1.254111534

Conclusão

Bom, notávelmente podemos ver que para esse série de produção industrial a metodologia que exibiu o menor erro, ou o melhor indicador de acurácia, foi O TBATS, seguIdo por MLP e ARNN.

Certamente vocês estão curiosos pelo o que está por trás de cada metodologia, então recomendo olhar esses artigos, todos disponíveis na internet. No futuro pretendo explica cada um deles.

REFERÊNCIAS

DE LIVERA, A.M.; HYNDMAN, R.J.; SNYDER, R. D. , “Forecasting time series with complex seasonal patterns using exponential smoothing”, Journal of the American Statistical Association, 2011

GARDNER, M.; DORLING, S.R. “Artificial neural networks (the multilayer perceptron)-a review of applications in the atmospheric sciences” Atmospheric Environment Vol. 32, No. 14/15, pp. 2627-2636, 1998

HYNDMAN, R.; KOEHLER, A. “Another look at measures of forecast accuracy”, International Journal of Forecasting, Volume 22, 2006.

THIELBAR, M.; DICKEY, D.; “Neural Networks for Time Series Forecasting: Practical Implications of Theoretical Results”, 2010.

’