Introdução

Uma das partes mais importantes em data science é poder ajudar na tomada de decisão e na projeção de resultados para empresas. Entender seu faturamento indicadores sociais e financeiros é, geralmente caro, envolve mutias pessoas e leva um tempo considerável.

E se ao invés de tudo isso pudéssemos utilizar o poder da matemática para nos salvar? E se criássemos um algorítimo com poderemos maiores que o da própria Mãe Diná? Parece magia, tem cheiro de magia, se move como se fosse mágico, mas no final do dia é pura ciência!

Neste texto vou mostrar como criar uma rede neural recorrente capaz de olhar para a cotação diária do dolar desde 02/01/2000 até 31/01/2022 e prever os próximos valores!

Para isso eu utilizei a base de dados pública do Banco Central e como sempre disponibilizarei o código fonte completo ao final do texto, que hoje tem um toque especial, vamos comparar a parformance de nossa rede neural com outro modelo de previsão de dados, chamado série temporal.

Os dados podem ser baixados através do link https://www3.bcb.gov.br/sgspub

Preparação dos pacotes

Abaixo está a lista dos pacotes que serão utilizados e o método que faz a ativação ou instalação se necessário.

# ---- Pacotes utilizados ---- 
pacotes <- c("ggplot2", "plotly", "dplyr", "forecast", "reshape2", "keras", "tidyr")

if(sum(as.numeric(!pacotes %in% installed.packages())) != 0){
  instalador = 
    pacotes[!pacotes %in% installed.packages()]
  for(i in 1:length(instalador)) {
    install.packages(instalador, dependencies = T, verbose = F)
    break()
  }
  sapply(pacotes, require, character = T) 
} else {
  sapply(pacotes, require, character = T) 
}
##  ggplot2   plotly    dplyr forecast reshape2    keras    tidyr 
##     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE     TRUE

Análise descritiva dos dados

Como trabalharemos com séries temporais nossa base é “simples” contendo uma coluna com a data e outra com o valor da cotação do dolar naquele dia.

dolar = 
  read.csv("dolar.csv", sep = ";", dec = ",") %>%
  mutate(data = as.Date(data))

head(dolar)
##         data cotacao
## 1 2000-01-03  1.8003
## 2 2000-01-04  1.8329
## 3 2000-01-05  1.8536
## 4 2000-01-06  1.8453
## 5 2000-01-07  1.8273
## 6 2000-01-10  1.8153

Ao utilizarmos o ggplot para visualizarmos a cotação do dolar ao longo do tempo, podemos observar, em azul que há uma tendêndia de crescimento, isso precisará ser tratado mais a frente.

## `geom_smooth()` using formula 'y ~ x'

Preparação dos dados para o modelo

Primeiramente definiremos um ponto de corte para nossos dados, por se tratar de uma série temporai, estabeleceremos este ponto de separação das bases de treino e teste de acordo com uma determinada posição na linha do tempo, desta forma separaremos um período para treinar a rede neural e um período para testarmos se ela é mesmo capaz de prever dados no “futuro”, comparando os dados preditos com a informação real.

Neste caso definiremos 80% da base como o ponto ideal, deixando 20% para testes.

#embora os dados já estejam ordenados por data, vamos forçar a ordenação
#de forma crescente com o comando arrange do pacote dplyr
dolar = 
  dolar %>% 
  arrange(data)

ponto_corte <- round(0.8 * nrow(dolar))

Em seguida precisamos eliminar a tendência de crescimento na cotação, este tipo de fenômeno diminui a acurácia do modelo e a capacidade da rede neural realizar previsões assertivas. Por isso é necessário removê-la e para isso ao invés de trabalharmos com o valor da cotação utilizaremos a diferença entre o valor atual e o valor anterior. Vamos aproveitar também para classificar os dados em treino e teste de acordo com o ponto de corte definido anteriormente

dolar_s_tendencia =
  dolar %>%
  #a função lag retorna o valor anterior da cotação do dólar
  mutate(diferenca = cotacao - lag(cotacao)) %>%
  #como a primeira linha não possui um valor anterior ela deve ser removida
  filter(row_number() >= 2) %>%
  mutate(tipo_dado = ifelse(row_number() >= ponto_corte, "teste", "treino"))

O próximo passo é a padronização dos dados em uma escala de -1 a 1, onde todas as diferenças entre as cotações serão recalculados de forma proporcional.

escala_serie = 
  dolar_s_tendencia %>%
  filter(tipo_dado == "treino") %>%
  summarise(menor_diferenca = min(diferenca), 
            maior_diferenca = max(diferenca), 
            aplitude_diferenca = maior_diferenca - menor_diferenca)

dolar_normalizados = 
  dolar_s_tendencia %>%
  select(data, cotacao, diferenca, tipo_dado) %>%
  mutate(diferenca_padronizada = 2*((diferenca - escala_serie$menor_diferenca)/escala_serie$aplitude_diferenca - .5))

Podemos observar que a tendência foi eliminada dos dados e que as bases de treino e teste estão separadas

## `geom_smooth()` using formula 'y ~ x'

Preparação do modelo

Para a preparação do modelo de dados é necessário realizar um segundo lag, desta vez criando a diferença entre os dados, uma diferença entre as diferenças, utilizamos novamente o comando lag e formatamos o nome das colunas para aqueles que serão utilizados no modelo.

Com os dados devidamente normalizados e prontos faremos a separação das bases de treino e teste

#Dados de treino
treino_y = 
  dolar_normalizados_modelo %>%
  filter(tipo_dado == "treino") %>%
  select(var_y) %>%
  as.matrix()

treino_x =
  dolar_normalizados_modelo %>%
  filter(tipo_dado == "treino") %>%
  select(var_x) %>%
  as.matrix()

# Dados para teste de acurácia do modelo
teste_y =
  dolar_normalizados_modelo %>%
  filter(tipo_dado == "teste") %>%
  select(var_y) %>%
  as.matrix()

teste_x = 
  dolar_normalizados_modelo %>%
  filter(tipo_dado == "teste") %>%
  select(var_x) %>%
  as.matrix()

Configuração do modelo

Para que possamos utilizar o pacote keras para a criação de uma série temporal é necessário transformarmos a variável de treino X em uma matriz de 3 dimensões sendo que a posição X tem o tamanho da base e as demais apenas 1 posição

dim(treino_x) <- c(length(treino_x), 1, 1)
batch_size = 1

Na preparação do modelo, informaremos alguns hiperparametros obrigatórios e que serão fundamentais para a criação de nosso modelo de predição:

modelo_lstm = 
  keras_model_sequential() %>%
  layer_lstm(units = 10,
             batch_input_shape = c(batch_size, 1, 1),
             stateful= T,
             dropout = 0.05) %>%
  layer_dense(units = 1)
## Loaded Tensorflow version 2.8.0

O comando compile indica demais configurações do modelo como a taxa de aprendizado (0.05) e decay (1e-6), que indica a taxa na queda do aprendizado da rede neural, isso ajuda a evitar eventuais overffitings, quando o modelo tenta se tornar tão completo que acaba sendo incapaz de prever novos dados. O último parâmetro indica quam métrica será utilizada para avaliar a taxa de melhora a cada iteração do modelo, neste caso o MSE ou erro quadrático médio

modelo_lstm %>% 
  compile(
    loss = 'mean_squared_error',
    optimizer = optimizer_adam(learning_rate = 0.05, decay = 1e-6 ),
    metrics = c('mse'))

Execução do modelo

Neste exemplo a rede será executada 5 vezes, sendo que a cada execução, alguns indicadores eram reiniciados para evitar lixo nas previsões futuras.

O tempo de duração depende de uma série de fatores: quantidade de dados, complexidade do modelo proposto, capacidade computacional, utilização ou não de GPU (placas de vídeo podem ser utilizadas para realização de cálculos complexos, você pode saber mais em https://www.cienciaedados.com/gpu-e-deep-learning/)

O processo geração é recursivo, ou seja, ele utiliza dados dele mesmo para aprimorar-se.

Predição (hora de ver se esse negócio funciona mesmo)

Assim como no processo de criação do modelo o processo de predição também é recursivo e deve ser realizado uma vez para cada linha da base de treino que reservamos, neste caso 1117 linhas

Ao final é possível criar um novo data frame combinando os dados já normalizados com os valores preditos

df_predicoes = 
  dolar_normalizados %>% 
  #Observe que somente as informações separadas para teste foram utilizadas
  filter(tipo_dado == "teste") %>%
  mutate(yhat = predicoes, # variável y com o valor previsto (diferença entre as cotações)
         valor_predito = yhat + lag(cotacao)) # cotação real somada a diferença predita

E possível verificar que os valores previstos e a cotação real estão muito próximos indicando que o modelo parece ser preciso em suas predições

df_predicoes %>%
  select(data, cotacao, valor_predito) %>%
  pivot_longer(cols = -data,
               names_to = "dolar",
               values_to = "Valor") %>%
  ggplot(aes(x = data, y = Valor, color = dolar)) +
  geom_line() +
  theme_minimal()
## Warning: Removed 1 row(s) containing missing values (geom_path).

Para verificarmos de forma mais objetiva a qualidade do modelo criado utilizaremos o Erro Quadrático médio (EQM) ou MSE (em inglês). Ele faz a comparação da média dos valores previstos e os valores reais e quando mais próximo de zero melhor

df_predicoes %>%
  group_by() %>%
  summarise(eqm = mean((cotacao - valor_predito)^2,na.rm = T))
## # A tibble: 1 x 1
##       eqm
##     <dbl>
## 1 0.00850

Comparação com outro modelo de série temporal

Vamos verificar agora a qualidade da rede neural recorrente na previsão de valores e modelos baseados em séries temporais, outra excelente forma de se prever valores financeiros e muito utilizada no mercado.

Assim como para a rede neural os modelos abaixo dependem de testes de verificação de tendência nos dados, algo que já fizemos e não precisaremos repetir aqui, só precisamos garantir que utilizaremos os mesmos dados em todos os modelos.

Para isso utilizaremos dois modelos, o Arima, muito comum e repleto de opções se parametrização. Para isso é necessário combinar dois comandos, o primeiro é o auto.arima, que avalia a base de dados e propõe o melhor modelo possível.

Em seguida utilizaremos o comando forecastpara realizar predição de dados, observe o parâmetro h = 1117, ele indica que o modelo deve prever 1117 dias com base nos dados existentes

predicoes_arima = 
  forecast(auto.arima(treino_y), h = 1117)

Ao verificarmos como os valores preditos pela série temporal se encaixam nos valores reais temos um resultado aparentemente muito preciso

df_predicoes %>% 
  mutate(valor_predito_arima = valor_predito + (predicoes_arima$mean/2 + 0.5) * escala_serie$aplitude_diferenca + escala_serie$menor_diferenca) %>%
  select(data, cotacao, valor_predito_arima) %>%
  ggplot() +
  geom_line(aes(y = cotacao, x=data, color="cotacao")) + 
  geom_line(aes(y = valor_predito_arima, x= data, color="valor_predito_arima")) + 
  theme_minimal()
## Warning: Removed 1 row(s) containing missing values (geom_path).

Vamos verificar agora o erro quadrático médio dos três modelos e avaliar qual performa melhor (chega mais próximo de zero por ter pouca diferença entre os valores preditos e valores reais)

df_predicoes %>% 
  mutate(valor_predito_arima = (predicoes_arima$mean/2 + 0.5)*escala_serie$aplitude_diferenca + escala_serie$menor_diferenca) %>%
  summarise(lstm = mean((cotacao - valor_predito) ^ 2,na.rm = T),
            arima = mean((cotacao - valor_predito_arima)^2,na.rm = T)) %>%
  melt()
## No id variables; using all as measure variables
##   variable        value
## 1     lstm  0.008501429
## 2    arima 21.067978785

Incrível! A rede neural recorrente praticamente zerou a diferença entre os valores e foi capaz de realizar predições muito confiáveis. Olhando para estes valores de forma visual, a diferença entre os modelos fica ainda mais evidente

df_predicoes %>% 
  mutate(valor_predito_arima = (predicoes_arima$mean/2 + 0.5)*escala_serie$aplitude_diferenca + escala_serie$menor_diferenca) %>%
  summarise(lstm = mean((cotacao - valor_predito) ^ 2,na.rm = TRUE),
            arima = mean((cotacao - valor_predito_arima)^2,na.rm = TRUE)) %>%
  melt() %>%
  ggplot(aes(y = value, x=variable, fill = factor(variable))) +
  geom_bar(stat = "identity") + 
  geom_label(aes(label = (round(value,3))), hjust = 1.2, color = "white", size = 5) +
  coord_flip() +
  labs(title="Erro Quadrático Médio (EQM) ou MSE", 
       x = "bases",
       y = "") +
  theme_minimal()
## No id variables; using all as measure variables

# Outra métrica que pode ser utilizada para a comparação entre os modelos é a raiz quadrática média ela possui algumas particularidades em comparação com o EQM (ou MSE) mas por não ser o foco deste texto vamos utilizá-la como um indicador complementar e o resultado é confirma o que já vimos até aqui a rede neural é superior a técnica de série temporal, para este conjunto de dados.

df_predicoes %>% 
  mutate(valor_predito_arima = (predicoes_arima$mean/2 + 0.5) * escala_serie$aplitude_diferenca + escala_serie$menor_diferenca) %>%
  na.omit() %>%
  summarise(lstm = sqrt(mean((cotacao - valor_predito)^2)),
            arima = sqrt(mean((cotacao - valor_predito_arima)^2))) %>%
  melt()
## No id variables; using all as measure variables
##   variable      value
## 1     lstm 0.09220319
## 2    arima 4.59106598
df_predicoes %>% 
  mutate(valor_predito_arima = (predicoes_arima$mean/2 + 0.5) * escala_serie$aplitude_diferenca + escala_serie$menor_diferenca) %>%
  na.omit() %>%
  summarise(lstm = sqrt(mean((cotacao - valor_predito)^2)),
            arima = sqrt(mean((cotacao - valor_predito_arima)^2))) %>%
  melt() %>%
  ggplot(aes(y = value, x=variable, fill = factor(variable))) +
  geom_bar(stat = "identity") + 
  geom_label(aes(label = (round(value,3))), hjust = 1.2, color = "white", size = 5) +
  coord_flip() +
  labs(title="RMSE", 
       x = "modelos",
       y = "") +
  theme_minimal()
## No id variables; using all as measure variables

Estes modelos foram criados utilizando versões mais resumidas e simples de seus respectivos pacote e podem ser otimizados para entregar predições ainda mais precisas. É importante ter em mente que uma série temporal pode apresentar excelentes resultados se otimizada e configurada corretamente. Este exemplo tem como objetivo apresentar a técnica da rede neural recorrente e como ela pode ser uma excelente opção, mesmo sem otimização.

Materiais para consulta

Link para o código fonte e base de dados: https://1drv.ms/u/s!AtsEI2mTf0TYkp8E2f-YA6QEPuBa6Q?e=hiXVgQ

Outros sites para consulta

  1. https://www.linkedin.com/pulse/rmse-ou-mae-como-avaliar-meu-modelo-de-machine-learning-rezende/?originalSubdomain=pt

  2. https://robjhyndman.com/expsmooth/

  3. https://www.pedronl.com/post/previsao-de-series-temporais-com-o-r/

  4. https://forecasters.org/wp-content/uploads/gravity_forms/7-621289a708af3e7af65a7cd487aee6eb/2015/07/Smyl_Slawek_ISF2015.pdf