Pacotes utilizados

library(tidyverse)
library(readr)
library(forecast)
library(tidyquant)
library(timetk)
library(sweep)
library(reshape2)
library(tseries)
library(earth)
library(keras)
library(gridExtra)

Banco de dados

O banco de dados utilizado para o projeto foi de um desafio de previsão de demanda de itens de lojas. Esse desafio consta no site do kaggle, onde foi extraído o arquivo em formato csv

Introdução e Objetivos

O objetivo desse trabalho é de avaliar o desempenho de cada modelo para as previsões de cada item de cada loja estudada. Como existe no banco de dados, 10 lojas com cada loja apresentando 50 items em estoque, então teremos cercas de 500 séries a serem modeladas. Por questões de tempo computacional, vai ser feito apenas uma previsão (um passo a frente) para cada série. No final, vai ser obtida uma tabela com cada Erro quadrátivo médio para cada um dos 5 modelos trabalhados.

Análise descritiva

Antes de aplicarmos essa sequência de modelos, é útil fazer um estudo descritivo antes, afim de evideciarmos o comportamento dos dados trabalhados. Olhando o histograma abaixo, podemos perceber que existe uma certa assimetria a direita, com maiores concetrações para as vendas abaixo de 50 unidades.

g1<-train %>% ggplot(aes(x=sales)) + 
  geom_histogram(binwidth = 12, color="black", fill="white") +
  labs(x="Vendas", y="Frequência")
g1

Podemos fazer uma vizualização da série, selecionando apenas o item 1 da loja 1. Foi feito também uma média móvel para vizualizar não só os estoques dia a dia, mas semanalmente e mensalmente. Podemos dizer que existe um comportamento similar entre cada ano.

#Vizualization of the store 1
st_store1<-train %>% 
  filter(store==1, item==1) %>% 
  mutate(sales_ma_weekly = ma(sales,order=7), sales_ma_month = ma(sales,order=30)) %>% 
  ggplot() +
  geom_line(aes(x=date, y=sales, col="Sales"))+
  geom_line(aes(x=date, y=sales_ma_weekly, col= "Weekly"))+
  geom_line(aes(x=date, y=sales_ma_month, col="Monthly"))+
  labs(y="Time Series of Sales",x="Time")    

st_store1

Modelos Arima , ETS e SNAIVE.

Nesta seção será mostrada como foi aplicado os modelos acima. Primeramente foi feito as bases de treino e teste.

train2<- train %>% 
  filter(date <= "2017-12-30") %>%  mutate(Flag = "Train")
test2<-  train %>% 
  filter(date > "2017-12-30") %>%  mutate(Flag = "Test")

data_all<- train2 %>%
  bind_rows(test2) %>%
  transmute(
    store,
    item,
    date,
    log_sales = log(sales + 1e-6),
    Flag)
 
train <- data_all %>% filter(Flag == "Train") %>% dplyr::select(-Flag)
test <- data_all %>% filter(Flag == "Test") %>% dplyr::select(-Flag)

Utilizando a função nest() do pacote tidyverse, podemos reorganizar o dataset anterior, no qual foi feito um agrupamento pela loja e item, no qual, em cada linha desse banco dataset apresenta a série temporal da loja e item em análise. Podemos observar isso abaixo.

dataset <- train %>%
  group_by(store, item) %>%
  nest() %>% 
  mutate(data.ts = map(.x       = data, 
                       .f       = tk_ts, 
                       select   = -date, 
                       start    = 2013,
                       freq     = 365))
dataset

## # A tibble: 500 x 4
## # Groups:   store, item [500]
##    store  item           data data.ts
##    <dbl> <dbl> <list<df[,2]>> <list> 
##  1     1     1    [1,825 × 2] <ts>   
##  2     2     1    [1,825 × 2] <ts>   
##  3     3     1    [1,825 × 2] <ts>   
##  4     4     1    [1,825 × 2] <ts>   
##  5     5     1    [1,825 × 2] <ts>   
##  6     6     1    [1,825 × 2] <ts>   
##  7     7     1    [1,825 × 2] <ts>   
##  8     8     1    [1,825 × 2] <ts>   
##  9     9     1    [1,825 × 2] <ts>   
## 10    10     1    [1,825 × 2] <ts>   
## # … with 490 more rows

A estratégia agora é:

Criar 3 novas colunas usando a função mutate().
Depois, utilizar a função map(), usando cada linha do dataset que é no formato <ts>, ou seja, série temporal.
Aplicar os modelos Arima(), ets() e snaive().

dataset2 <- dataset %>% mutate(
  model_arima = map(data.ts, Arima),
  model_ets = map(data.ts, ets), 
  model_snaive = map(data.ts, snaive)) %>%
  select(store, item, model_arima, model_ets, model_snaive)

dataset2

## # A tibble: 500 x 5
## # Groups:   store, item [500]
##    store  item model_arima model_ets model_snaive
##    <dbl> <dbl> <list>      <list>    <list>      
##  1     1     1 <fr_ARIMA>  <ets>     <forecast>  
##  2     2     1 <fr_ARIMA>  <ets>     <forecast>  
##  3     3     1 <fr_ARIMA>  <ets>     <forecast>  
##  4     4     1 <fr_ARIMA>  <ets>     <forecast>  
##  5     5     1 <fr_ARIMA>  <ets>     <forecast>  
##  6     6     1 <fr_ARIMA>  <ets>     <forecast>  
##  7     7     1 <fr_ARIMA>  <ets>     <forecast>  
##  8     8     1 <fr_ARIMA>  <ets>     <forecast>  
##  9     9     1 <fr_ARIMA>  <ets>     <forecast>  
## 10    10     1 <fr_ARIMA>  <ets>     <forecast>  
## # … with 490 more rows

Nos códigos abaixo, podemos fazer também as previsões para cada linha. Posteriormente fazer 3 bancos pra cada modelo, denominados por actual.pred_arima, actual.pred_ets, actual.pred_snaive.

#data with the forecast by row
tst_nest <- test %>%
  group_by(store, item) %>%
  nest() %>% 
  inner_join(dataset2, by = c("store", "item")) %>%
  mutate(fcast.arima = map( model_arima, forecast, h=1) , 
         fcast.ets =   map( model_ets, forecast, h=1),
         fcast.snaive = map(model_snaive, forecast, h=1))

actual.pred_arima <- tst_nest %>%
  mutate(pred_arima = map(fcast.arima, sw_sweep, fitted = F, timetk_idx = TRUE)) %>%
  unnest(pred_arima) %>% filter(key=="forecast") %>% 
  select(store, item, index , log_sales) %>% 
  mutate(real_pred_arima = round(exp(log_sales),0)) %>% 
  bind_cols(test2) %>%
  select(store, item, index, sales, real_pred_arima)

actual.pred_ets <- tst_nest %>%
  mutate(pred_ets = map(fcast.ets, sw_sweep, fitted = F, timetk_idx = TRUE)) %>%
  unnest(pred_ets) %>% filter(key=="forecast") %>% 
  select(store, item, index, log_sales) %>% 
  mutate(real_pred_ets = round(exp(log_sales),0)) %>% 
  bind_cols(test2) %>%
  select(store, item, index, sales, real_pred_ets) 

actual.pred_snaive <- tst_nest %>%
  mutate(pred_snaive = map(fcast.snaive, sw_sweep, fitted = F, timetk_idx = TRUE)) %>%
  unnest(pred_snaive) %>% filter(key=="forecast") %>% 
  select(store, item, index , log_sales) %>% 
  mutate(real_pred_snaive = round(exp(log_sales),0)) %>% 
  bind_cols(test2) %>%
  select(store, item, index, sales, real_pred_snaive)

No final, calcula-se os EQM´s para cada modelo, utlizando os data.frame´s recém formados.

EQM_ARIMA <-mean((actual.pred_arima$sales - actual.pred_arima$real_pred_arima )^2)
EQM_ETS <-mean((actual.pred_ets$sales - actual.pred_ets$real_pred_ets)^2)
EQM_SNAIVE<-mean((actual.pred_snaive$sales - actual.pred_snaive$real_pred_snaive)^2)

Modelos MARS (Multivariate Adaptative Regression Splines)

O modelo Multivariate Adaptive Regression Splines (MARS), é um modelo de regressão não paramétrica que utiliza-se de vários modelos de regressão linear em todo o domínio dos valores das variáveis preditoras. Isso é feito particionando os dados e executando um modelo de regressão linear em cada partição diferente.

O algoritmos do modelo MARS é uma extensão de modelos lineares que não faz suposições sobre o relacionamento entre a variável resposta e as variáveis preditoras. Enquanto os Modelos Lineares Generalizados e os Modelos Aditivos Generalizados assumem que os coeficientes das variáveis preditivas são de forma invariante, o algoritmo MARS especificamente não leva em consideração esse fator.

Vantagens:

Funciona bem com um grande número de variáveis preditoras.
Detecta automaticamente interações entre variáveis.
É um algoritmo eficiente e rápido, apesar de sua complexidade.
Robusto para outliers.

Desvantagens:

Susceptível ao overfitting.
Mais difícil de entender e interpretar do que outros métodos.
Não é eficiente para dados faltantes.

Como estamos lidando agora com um objeto de class <mts>, então para o caso desse modelo, não foi possível fazer a aplicação junto com a dos 3 modelos anteriores. Por isso foi feito separadamente. Seguindo um processo parecido do que foi feito anteriormente, definindo um novo banco de dados.

Olhando para o data_all_MARS, podemos notar que foi adicionado 3 novas variáveis, chamadas de year, month e wday. Essas novas variáveis, vão ser utilizadas como variáveis explicativas para o modelo.

data_all_mars <- train2 %>%
  bind_rows(test2) %>%
  transmute(
    Flag,
    store,
    item,
    log_sales = log(sales + 1e-6),
    date,
    year = scale(as.integer(year(date))),
    month = scale(month(date)),
    wday = scale(wday(date, week_start = 1))
  )

train_mars <- data_all_mars %>% filter(Flag == "Train") %>% select(-Flag)
test_mars <- data_all_mars %>% filter(Flag == "Test") %>% select(-Flag)

Fazendo o mesmo procedimento com a função nest(), podemos reorganizar o data_all_MARS para a mesma estrutura conhecida e já utilizada anteriormente. Podemos notar agora, que houve uma mudança na estrutura de cada linha desse banco, no qual, cada linha é constituída de um data.frame, com estrutrua <mts>, ou seja, uma série temporal multivariada.

dataset_mars <- train_mars %>%
  group_by(store, item) %>%
  nest() %>% 
  mutate(data.ts = map(.x       = data, 
                       .f       = tk_ts, 
                       select   = -date, 
                       start    = 2013,
                       freq     = 365))
dataset_mars

## # A tibble: 500 x 4
## # Groups:   store, item [500]
##    store  item           data data.ts
##    <dbl> <dbl> <list<df[,5]>> <list> 
##  1     1     1    [1,825 × 5] <mts>  
##  2     2     1    [1,825 × 5] <mts>  
##  3     3     1    [1,825 × 5] <mts>  
##  4     4     1    [1,825 × 5] <mts>  
##  5     5     1    [1,825 × 5] <mts>  
##  6     6     1    [1,825 × 5] <mts>  
##  7     7     1    [1,825 × 5] <mts>  
##  8     8     1    [1,825 × 5] <mts>  
##  9     9     1    [1,825 × 5] <mts>  
## 10    10     1    [1,825 × 5] <mts>  
## # … with 490 more rows

Aplicando agora o modelo MARS, utilizando o a função earth().

dataset2_mars <- dataset_mars %>% mutate(model = 
  map(data.ts, ~ earth(log_sales ~ year + month + wday + year:wday , data = .x))) %>%
  select(store, item, model)

Realizando a predição de cada linha do dataset2_mars e criando um data.frame com os valores reais da fase de test e as suas respectivas predições, denominado por actual.pred_snaive.

tst_nest_mars <- test_mars %>%
  group_by(store, item) %>%
  nest() %>% 
  inner_join(dataset2_mars, by = c("store", "item")) %>%
  mutate(pred = map2(data, model, ~ exp(predict(.y, newdata = .x))))

actual.pred_mars <- tst_nest_mars %>%
  unnest(pred) %>%
  pull(pred) %>%
  data.frame(real_pred_mars = .) %>% round(0) %>% 
  bind_cols(test2) %>% 
  select(store, item, date, sales, real_pred_mars) %>% 
  as_tibble()

Calculando o EQM.

EQM_MARS <- mean((actual.pred_mars$sales - actual.pred_mars$real_pred_mars)^2)

Modelo LSTM

A LSTM é bem adequada para classificar, processar e prever séries temporais com intervalos de tempo de duração desconhecida. A insensibilidade relativa ao comprimento do “gap” dá uma vantagem à LSTM em relação a RNNs (Redes neurais recursais) tradicionais. A principal vantagem do modelo LSTM, é a capacidade de lembrar momentos distantes no passado, conseguindo perceber o seu ponto de partida após um número grande de previsões (Exemplo: 1000 passos a frente).

data_all_LSTM <- train2 %>%
  bind_rows(test2) %>%
  transmute(
    Flag,
    store,
    item,
    log_sales = log(sales + 1e-6),
    date
    ) %>% dplyr:: select(store, item, date, log_sales, Flag)


dataset <- data_all_LSTM %>%
  group_by(store, item) %>%
  nest() %>% 
  mutate(data.ts = map(.x       = data, 
                       .f       = tk_ts, 
                       select   = -date, 
                       start    = 2013,
                       freq     = 365))

A ideia principal que utilizei foi de aplicar um laço de “for” para cada série do dataset, e aplicar os procedimentos do modelo LSTM. Cada série vai ser aplicada uma diferenciação e padronizar. Essa operação pode ser muito custosa, já que estamos trabalhando com 500 séries diferentes.

N <- 1825
lista_pred<-1:500 %>% map(~rep(0,1))
n<-nrow(dataset)
for (i in 1:n){
  lg<-dataset$data[[i]] 
  lg1<- lg %>% as.data.frame() %>% dplyr::select(date, log_sales, Flag)
  
  #Defining the diff_log_sales
  serie_diff <- lg1 %>% 
    mutate(Diff = log_sales - lag(log_sales)) %>%
    dplyr::filter(row_number() >= 2)

  scales <- serie_diff %>%
    filter(Flag == "Train") %>%
    summarise(min = min(Diff), max = max(Diff), Range = max - min)
  
  serie_norm <- serie_diff %>%
    mutate(DiffPadrao = 2*((Diff - scales$min)/scales$Range - .5)) 
  
  serie_norm2 <- serie_norm %>%
    transmute(date, Y = DiffPadrao, X1 = lag(Y), Flag, log_sales) %>% 
    filter(row_number() >= 2)
  
  x_train <- serie_norm2 %>%
    filter(Flag == "Train") %>%
    dplyr::select(X1) %>%
    as.matrix()
  
  y_train <- serie_norm2 %>%
    filter(Flag == "Train") %>%
    dplyr::select(Y) %>%
    as.matrix()
  
  x_test <- serie_norm2 %>%
    filter(Flag == "Test") %>%
    dplyr::select(X1) %>%
    as.matrix()
  
  y_test <- serie_norm2 %>%
    filter(Flag == "Test") %>%
    dplyr::select(Y) %>%
    as.matrix()
  
  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 

  model <- keras_model_sequential()
  model %>%
    layer_lstm(units=30, batch_input_shape = c(batch_size, X_shape2, X_shape3), stateful= TRUE) %>%
    layer_dropout(rate= 0.2) %>% 
    layer_dense(units=20) %>% 
    layer_dense(units=10) %>% 
    layer_dense(units=5) %>% 
    layer_dense(units = 1)
  
  
  model %>% compile(
    loss = 'mean_squared_error',
    optimizer = optimizer_adam(lr= 0.02, decay = 1e-6 ),
    metrics = c('mse')
  )
  #In real training Epochs must be equal 30
  Epochs = 20
  #cat("Serie temporal de numéro ", i, "\n")
  for(j in 1:Epochs ){
    model %>% fit(x_train, y_train, epochs = 1, batch_size = batch_size,
                  verbose = 0, shuffle = FALSE)
    model %>% reset_states()
  }
  
  L <- length(x_test)
  predictions = numeric(L)
  X_real = lg$log_sales[nrow(lg)]
  X = x_train[length(x_train)]
  for(k in 1:L){
    #cat(" Valor X inicial:",X)
    dim(X) = c(1,1,1)
    y_sdt = model %>% predict(X, batch_size=batch_size)
    # invert scaling
    yhat = (y_sdt/2 + .5)*scales$Range + scales$min
    X<-y_sdt
    # store
    predictions[k] <- yhat + X_real
    X_real <- predictions[k]
  }
  
  predict_df <- serie_norm %>%
    filter(row_number() >= N) %>%
    mutate(yhat = predictions) %>% 
    mutate(Vhat = yhat + log_sales)

lista_pred[[i]]<- as.vector(predict_df$yhat)
}

Utilizando a função unlist() para transformar a lista das predições dos 500 modelos em um vetor coluna. Posteriormente foi transformado esse vetor em um data.frame().

data_pred_lstm<- unlist(lista_pred) %>% data.frame(pred_lstm = .)

actual.pred_LSTM <- test2 %>%  
  bind_cols(data_pred_lstm) %>%
  mutate(pred_real_lstm = round(exp(data_pred_lstm$pred_lstm),0) ) %>% 
  dplyr::select(-pred_lstm)

Calculando o EQM do modelo final.

EQM_LSTM<- mean((actual.pred_LSTM$sales- actual.pred_LSTM$pred_real_lstm)^2 )

Comparativo

Olhando de forma geral, os modelos Arima, MARS e LSTM foram os que tiveram maior destaque para previsões com 1 passo a frente. O dois primeiros modelos citados são de fácil implementação e com um custo computacional menor que o modelo LSTM.

table_results<-data.frame(c(EQM_ARIMA, EQM_ETS, EQM_SNAIVE, EQM_MARS, EQM_LSTM))
dimnames(table_results)<- list(paste(c("ARIMA", "ETS", "SNAIVE", "MARS", "LSTM")), 
                               paste(c("Erro quadrático médio")))

table_results

##        Erro quadrático médio
## ARIMA                 70.636
## ETS                  179.422
## SNAIVE               145.952
## MARS                  63.410
## LSTM                  49.014

Escolhendo apenas os 3 melhores modelos e utilizando gráficos de dipersão dos Valores Preditos \(\times\) Valor Observado. Os modelos selecionados foram o ARIMA, MARS, LSTM. Podemos observar que o modelo que teve melhor desempenho, com valores preditos mais próximos dos reais, foi o modelo LSTM, sendo mais eficaz para previsões com 1 passo à frente.

g1 <- actual.pred_arima %>% ggplot(aes(x = sales, y = real_pred_arima)) + geom_point(size = 1) +
  labs(y = "Previsão - Arima", x = "Valor Obervado") +
  geom_abline(intercept = 0, col="red") +
  ggtitle("Análise da Previsão Utilizando o Modelo Arima")

g2<- actual.pred_mars %>% ggplot(aes(x = sales, y = real_pred_mars)) + geom_point(size = 1) +
  labs(y = "Previsão - MARS", x = "Valor Obervado") +
  geom_abline(intercept = 0, col="red") +
  ggtitle("Análise da Previsão Utilizando o Modelo MARS")

g3 <- actual.pred_LSTM %>% ggplot(aes(x = sales, y = pred_real_lstm)) + geom_point(size = 1) +
  labs(y = "Previsão - LSTM", x = "Valor Obervado") +
  geom_abline(intercept = 0, col="red") +
  ggtitle("Análise da Previsão Utilizando o Modelo LSTM")

gridExtra::grid.arrange(g1,g2,g3, ncol = 1, nrow = 3)

Projeto Final - Aplicação LSTM

Arthur Machado

10/12/2019