LSTM

Ao falar em redes neurais artificiais, é comum pensarmos em problemas de classificação ou regressão, contudo, as redes neurais artificiais também podem ser usadas para gerar previsões de séries temporais. Dentre os algoritmos utilizados para esta finalidade são conhecidos os MLP (Multilayer Perceptron) e o LSTM (Long Shot Term Memory). Este relatório irá focar no modelo de tipo LSTM para previsão de séries temporais univariadas.

Arquitetura LSTM (Long Short Term Memory)

As redes LSTM são projetadas para evitar o problema de dependência de longo prazo. Afinal, recordar a informação por longos períodos de tempo é praticamente o seu comportamento padrão, algo bem previsível e fácil de aprender.

No lugar de ter uma única camada de rede neural, existem quatro, esquematizados no diagrama a seguir. Nele, cada linha carrega um vetor inteiro (desde a saída de um nó até as entradas dos outros); os círculos cor-de-rosa representam operações pontuais (como a adição de vetores); as caixas amarelas são camadas de rede neural aprendidas; as linhas que se fundem significam concatenação, enquanto que uma linha de bifurcação denota que o seu conteúdo é copiado e enviado para locais diferentes.

O primeiro passo é decidir qual informação vai sair do estado celular (função sigmoid “forget gate”) e depois qual deve ficar, o que é dividido em duas partes: decidir quais valores devem ser atualizados (função sigmoid “input gate”) e cria um vetor de novos valores candidatos que podem ser adicionados ao estado (função “tanh”). No próximo passo, são combinados estes dois para criar uma atualização para um novo estado celular. Finalmente, é decidido o estado final (função sigmoid “output gate” seguida de uma função “tanh”). Por isso, precisamos normalizar a escala.

Aplicação em R

Para a aplicação em R vamos usar dados da taxa de câmbio brasileira a partir do ano 2000. Por convenção.

brazil <- read_rds("/opt/datasets/brazil.rds")
serie <- brazil %>%
  filter(Subject == "Currency exchange rates, monthly average",
         TIME >= "2000-01-01")

serie %>%
  select(Time = TIME, Value) %>%
  ggplot() + geom_point(aes(x = Time, y = Value), alpha = .5) +
  geom_line(aes(x = Time, y = Value), color = "black", alpha = .5) + 
  theme_minimal()

Antes de iniciar a modelagem, é interessante retirar qualquer tendência da série histórica. Para issso vamos diferenciar a série uma vez e dividi-la em 2 partes. Uma de treino e outra de teste.

# Tamanho de treino
N <- round(.7*nrow(serie))

# Retirando tendencia
serie_diff <- serie %>%
  mutate(Diff = Value - lag(Value)) %>%
  filter(row_number() >= 2) %>%
  mutate(
    Base = if_else(row_number() >= N, "test", "train"),
  )

# Como ficou
serie_diff %>%
  select(TIME, Value, Diff, Base)

## # A tibble: 235 x 4
##    TIME       Value    Diff Base 
##    <date>     <dbl>   <dbl> <chr>
##  1 2000-02-01  1.77 -0.0284 train
##  2 2000-03-01  1.74 -0.0333 train
##  3 2000-04-01  1.77  0.0261 train
##  4 2000-05-01  1.83  0.0598 train
##  5 2000-06-01  1.81 -0.0196 train
##  6 2000-07-01  1.80 -0.0105 train
##  7 2000-08-01  1.81  0.0114 train
##  8 2000-09-01  1.84  0.0300 train
##  9 2000-10-01  1.88  0.0404 train
## 10 2000-11-01  1.95  0.0684 train
## # … with 225 more rows

Além disso, como visto na primeira sessão, precisamos deixar a séria na escala [0,1]. Veja que não precisaremos mudar a escala para toda a série, porém, devemos utilizar os valores de amplitude e de máximo baseados na base de treinamento.

scales <- serie_diff %>%
  filter(Base == "train") %>%
  summarise(min = min(Diff), max = max(Diff), Range = max - min)

serie_norm <- serie_diff %>%
  select(Time = TIME, Value, Diff, Base) %>%
  mutate(DiffPadrao = 2*((Diff - scales$min)/scales$Range - .5))

serie_norm %>% ggplot(aes(x = Time, y = DiffPadrao, color = Base)) +
  geom_point(alpha = .5) +
  geom_line(alpha = .5) + 
  theme_minimal() + 
  scale_color_manual(values = c("black","blue"))

Modelagem

Para a modelagem, vamos utilizar um modelo LSTM AR(1), isto é, que depende somente da observação anterior. Sendo assim, teremos uma variável explitava que será o valor do mês anterior.

# To Model
serie_norm2 <- serie_norm %>%
  transmute(Time, Y = DiffPadrao, X1 = lag(Y), Base) %>% filter(row_number() >= 2)

x_train <- serie_norm2 %>%
  filter(Base == "train") %>%
  select(X1) %>%
  as.matrix()
y_train <- serie_norm2 %>%
  filter(Base == "train") %>%
  select(Y) %>%
  as.matrix()
x_test <- serie_norm2 %>%
  filter(Base == "test") %>%
  select(X1) %>%
  as.matrix()
y_test <- serie_norm2 %>%
  filter(Base == "test") %>%
  select(Y) %>%
  as.matrix()

Agora, vamos preparar os parâmetros para a estimação do modelo.

# Reshape the input to 3-dim
dim(x_train) <- c(length(x_train), 1, 1)

# specify required arguments
X_shape2 = dim(x_train)[2]
X_shape3 = dim(x_train)[3]
batch_size = 1 # must be a common factor of both the train and test samples
units = 5 # can adjust this, in model tuninig phase

Estimando os pesos. Veja que o processo será recursivo.

model <- keras_model_sequential()
model%>%
  layer_lstm(units,
             batch_input_shape = c(batch_size, X_shape2, X_shape3),
             stateful= TRUE) %>%
  layer_dense(units = 1)


model %>% compile(
  loss = 'mean_squared_error',
  optimizer = optimizer_adam(lr= 0.02, decay = 1e-6 ),
  metrics = c('mse')
)

Epochs = 200
for(i in 1:Epochs ){
  model %>% fit(x_train, y_train, epochs = 1,
                batch_size = batch_size,
                verbose = 1,
                shuffle = FALSE)
  cat("Processing",i,"of 100\n")
  model %>% reset_states()
}

Predição

Para a predição, vamos precisar de \(Y_{i-1}\). Também será um processo recursivo.

L = length(x_test)
predictions = numeric(L)

for(i in 1:L){
  if(i==1){
    X = x_test[1]
  }
  else{
    X <- predictions[i-1]
  }
  dim(X) = c(1,1,1)
  yhat = model %>% predict(X, batch_size=batch_size)
  # invert scaling
  yhat = (yhat/2 + .5)*scales$Range + scales$min
  # store
  predictions[i] <- yhat
}

predict_df <- serie_norm %>% filter(row_number() >= N) %>%
  mutate(yhat = predictions)


predict_df %>%
  mutate(Vhat = yhat + lag(Value)) %>%
  select(Time, Value, Vhat) %>%
  pivot_longer(cols = -Time,
               names_to = "Serie",
               values_to = "Valor") %>%
  ggplot(aes(x = Time, y = Valor, color = Serie)) +
  geom_point() +
  geom_line() +
  theme_minimal() + 
  scale_color_manual(values = c("black","blue"))

## Warning: Removed 1 rows containing missing values (geom_point).

## Warning: Removed 1 rows containing missing values (geom_path).

# Erro quadratico medio
erro_lstm <- predict_df %>%
  mutate(Vhat = yhat + lag(Value)) %>%
  summarise(erro = mean((Value - Vhat)^2,na.rm = TRUE))

erro_lstm

## # A tibble: 1 x 1
##     erro
##    <dbl>
## 1 0.0160

O modelo falha em um ponto perto do fim, mas mostra-se bastante eficaz.

Modelos clássicos

Vamos agora testar modelos mais clássicos como suavização exponencial, Arima e Naive.

library(forecast)

## Registered S3 method overwritten by 'xts':
##   method     from
##   as.zoo.xts zoo

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method             from    
##   fitted.fracdiff    fracdiff
##   residuals.fracdiff fracdiff

(mod_arima <- auto.arima(y_train))

## Series: y_train 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1     mean
##       0.3286  -0.1644
## s.e.  0.0738   0.0268
## 
## sigma^2 estimated as 0.05389:  log likelihood=7.71
## AIC=-9.42   AICc=-9.26   BIC=-0.13

prev_arima <- forecast(mod_arima,h = 71)

prev_arima %>% 
  autoplot() + 
  theme_minimal()

(mod_ets <- ets(as.vector(y_train)))

## ETS(A,N,N) 
## 
## Call:
##  ets(y = as.vector(y_train)) 
## 
##   Smoothing parameters:
##     alpha = 0.3123 
## 
##   Initial states:
##     l = -0.1553 
## 
##   sigma:  0.2471
## 
##      AIC     AICc      BIC 
## 378.5056 378.6565 387.7868

prev_ets <- forecast(mod_ets,h = 71)

prev_ets %>% 
  autoplot() + 
  theme_minimal()

predict_df %>% 
  mutate(Vhat = yhat + lag(Value)) %>%
  mutate(yhat_ar = (prev_arima$mean/2 + 0.5)*scales$Range + scales$min,
         yhat_ets = (prev_ets$mean/2 + 0.5)*scales$Range + scales$min) %>% 
  summarise(erro_lstm = mean((Value - Vhat)^2,na.rm = TRUE),
            erro_ar = mean((Value - yhat_ar)^2,na.rm = TRUE),
            erro_ets = mean((Value - yhat_ets)^2,na.rm = TRUE))

## # A tibble: 1 x 3
##   erro_lstm erro_ar erro_ets
##       <dbl>   <dbl>    <dbl>
## 1    0.0160    10.8     10.6

O LSTM ganhou de lavada mas já era de se esperar, pois ele foi projetado para ter alta performance a longo prazo. Vamos comparar agora para as 3 primeiras projeções para que a disputa seja justa.

predict_df %>% 
  mutate(Vhat = yhat + lag(Value)) %>%
  mutate(yhat_ar = (prev_arima$mean/2 + 0.5)*scales$Range + scales$min,
         yhat_ets = (prev_ets$mean/2 + 0.5)*scales$Range + scales$min) %>%
  filter(row_number()<=3) %>% 
  summarise(erro_lstm = mean((Value - Vhat)^2,na.rm = TRUE),
            erro_ar = mean((Value - yhat_ar)^2,na.rm = TRUE),
            erro_ets = mean((Value - yhat_ets)^2,na.rm = TRUE))

## # A tibble: 1 x 3
##   erro_lstm erro_ar erro_ets
##       <dbl>   <dbl>    <dbl>
## 1   0.00503    5.22     5.07

A rede neural continua ganhando, ou seja, é uma grande ferramenta para projeção de séries temporais.