Modelos_Previsao
Importando os dados:
library(tidyverse)
library(forecast)
library(tidyquant)
library(timetk)
library(modeltime)
library(tidymodels)
library(workflows)
df <- read_csv2("dados/DadosProduto.csv") %>%
mutate(DATA = as.Date(DATA, format = "%d/%m/%Y %H:%M")) %>%
mutate(QUANTIDADE = ifelse(QUANTIDADE < 0, 0, QUANTIDADE)) %>%
mutate(l_y = log(QUANTIDADE + 1)) # um teste com variavel em log
df %>% head()## # A tibble: 6 x 3
## DATA QUANTIDADE l_y
## <date> <dbl> <dbl>
## 1 2017-01-04 228 5.43
## 2 2017-01-05 228 5.43
## 3 2017-01-06 108 4.69
## 4 2017-01-07 276 5.62
## 5 2017-01-08 0 0
## 6 2017-01-09 0 0tk_get_frequency realiza uma análise espectral da série e retorna o período de frequência dominante da série. O resultado parece indicar uma série semanal. Alternativametne pode-se utilizar a função getfrequency().
df %>%
tk_index() %>%
tk_get_frequency(period = "auto")## [1] 7Converter os dados de diário para semanal usando tq_transmute:
# converter para semana
df_semanal <- df %>% tq_transmute(select = l_y,
mutate_fun = apply.weekly,
FUN = sum)
df_semanal %>% head()## # A tibble: 6 x 2
## DATA l_y
## <date> <dbl>
## 1 2017-01-08 21.2
## 2 2017-01-15 20.3
## 3 2017-01-22 23.0
## 4 2017-01-29 31.2
## 5 2017-02-05 23.0
## 6 2017-02-12 28.0Gráfico da série temporal
A série semanal parece ter alguns padrões de tendência e talvez de sazonalidade.
df_semanal %>%
plot_time_series(DATA, l_y)# codigo alternativo
#df_semanal %>%
# ggplot(aes(x = DATA, y = QUANTIDADE)) +
# geom_line()Para investigar vamos usar plot_stl_diagnostics. Alternativamente pode-se utilizar a função decompose.
df_semanal %>%
plot_stl_diagnostics(DATA, l_y, .frequency = "auto", .trend = "auto")Modelos de Previsão
Particionando os dados
Particionando os dados entre testing set e training set. Vamos utilizar os últimos 6 meses como o testing set (assess = "6 months") e todos os períodos anteriores como training set (cumulative = TRUE). O ideal é utilizar um procedimento de cross-validation.
ts_split <- df_semanal %>%
time_series_split(assess = "6 months", # últimos 6 meses como testing set
cumulative = TRUE) # todos os períodos anteriores como training setO resultado do split do banco de dados:
ts_split %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(DATA, l_y)Modelos
Temos dois tipos de modelos: (1) os modelos clássicos de séries temporais, como o ARIMA e o Prophet, que usarão os valores passados da série para produzir ajuste do modelo e (2) modelos de machine learning que incorporam outras variáveis explicativas, sobretudo derivadas da variável data.
Para facilizar o processo de rodar vários modelos vamos criar um workflow com três passos:
- Definir o modelo
- Pre-processar os dados, criando as variáveis
- Ajustar o modelo sob os dados
Os modelos de TS convencionais sempre usarão a mesma receita, de quantidade (y) em função dos seus valores passados (DATA).
recipe_auto <- recipe(l_y ~ DATA, data = training(ts_split))Para os demais modelos, usaremos um conjunto de variáveis derivadas de DATA. Com step_timeseries_signature convertemos a variável DATA em várias variáveis derivadas. Com step_date(DATA, featurs = "month") criamos uma variável de mês, que captura sazonalidade mensal. Com step_mutate(date_num = as.numeric(DATA)) criamos uma variável numérica e crescente para cada período no tempo. Esta última variável captura tendências na série.
recipe_spec <- recipe(l_y ~ DATA, data = training(ts_split)) %>%
step_timeseries_signature(DATA) %>%
step_fourier(DATA, period = c(12, 24, 36, 48), K=2) %>%
step_rm(contains("iso"), contains("minute"), contains("hour"),
contains("am.pm"), contains("xts")) %>%
step_date(DATA, features = "month", ordinal = FALSE) %>%
step_mutate(date_num = as.numeric(DATA)) %>%
step_normalize(date_num) %>%
step_normalize(contains("index.num")) %>%
step_dummy(contains("lbl"), one_hot = TRUE) %>%
step_dummy(all_nominal())As variáveis criadas por curiosidade:
bake(prep(recipe_spec), new_data = training(ts_split))## # A tibble: 184 x 67
## DATA l_y DATA_index.num DATA_year DATA_half DATA_quarter
## <date> <dbl> <dbl> <int> <int> <int>
## 1 2017-01-08 21.2 -1.72 2017 1 1
## 2 2017-01-15 20.3 -1.70 2017 1 1
## 3 2017-01-22 23.0 -1.68 2017 1 1
## 4 2017-01-29 31.2 -1.66 2017 1 1
## 5 2017-02-05 23.0 -1.64 2017 1 1
## 6 2017-02-12 28.0 -1.62 2017 1 1
## 7 2017-02-19 27.9 -1.61 2017 1 1
## 8 2017-02-26 28.5 -1.59 2017 1 1
## 9 2017-03-05 18.0 -1.57 2017 1 1
## 10 2017-03-12 21.5 -1.55 2017 1 1
## # ... with 174 more rows, and 61 more variables: DATA_month...7 <int>,
## # DATA_day <int>, DATA_second <int>, DATA_wday <int>, DATA_mday <int>,
## # DATA_qday <int>, DATA_yday <int>, DATA_mweek <int>, DATA_week <int>,
## # DATA_week2 <int>, DATA_week3 <int>, DATA_week4 <int>, DATA_mday7 <int>,
## # DATA_sin12_K1 <dbl>, DATA_cos12_K1 <dbl>, DATA_sin12_K2 <dbl>,
## # DATA_cos12_K2 <dbl>, DATA_sin24_K1 <dbl>, DATA_cos24_K1 <dbl>,
## # DATA_sin24_K2 <dbl>, DATA_cos24_K2 <dbl>, DATA_sin36_K1 <dbl>,
## # DATA_cos36_K1 <dbl>, DATA_sin36_K2 <dbl>, DATA_cos36_K2 <dbl>,
## # DATA_sin48_K1 <dbl>, DATA_cos48_K1 <dbl>, DATA_sin48_K2 <dbl>,
## # DATA_cos48_K2 <dbl>, DATA_month...36 <fct>, date_num <dbl>,
## # DATA_month.lbl_01 <dbl>, DATA_month.lbl_02 <dbl>, DATA_month.lbl_03 <dbl>,
## # DATA_month.lbl_04 <dbl>, DATA_month.lbl_05 <dbl>, DATA_month.lbl_06 <dbl>,
## # DATA_month.lbl_07 <dbl>, DATA_month.lbl_08 <dbl>, DATA_month.lbl_09 <dbl>,
## # DATA_month.lbl_10 <dbl>, DATA_month.lbl_11 <dbl>, DATA_month.lbl_12 <dbl>,
## # DATA_wday.lbl_1 <dbl>, DATA_wday.lbl_2 <dbl>, DATA_wday.lbl_3 <dbl>,
## # DATA_wday.lbl_4 <dbl>, DATA_wday.lbl_5 <dbl>, DATA_wday.lbl_6 <dbl>,
## # DATA_wday.lbl_7 <dbl>, DATA_month...37_fev <dbl>,
## # DATA_month...37_mar <dbl>, DATA_month...37_abr <dbl>,
## # DATA_month...37_mai <dbl>, DATA_month...37_jun <dbl>,
## # DATA_month...37_jul <dbl>, DATA_month...37_ago <dbl>,
## # DATA_month...37_set <dbl>, DATA_month...37_out <dbl>,
## # DATA_month...37_nov <dbl>, DATA_month...37_dez <dbl>Modelo 1: ARIMA
O primeiro modelo é um ARIMA simples, que usará o auto.arima() para definir os parâmetros ótimos do modelo. Assim, o workflow para o arima adiciona a receita recipe_auto, depois definimos o modelo e por fim, ajustamos o modelo aos dados.
modelo_fit_arima <- arima_reg() %>%
set_engine(engine = "auto_arima")
workflow_fit_arima <- workflow() %>%
add_recipe(recipe_auto) %>%
add_model(modelo_fit_arima) %>%
fit(training(ts_split))Modelo 2: ARIMA BOOSTED (XGBOOST)
É uma forma de utilizar o modelo XGBoost para modelar os resíduos de um modelo ARIMA. Aqui vamos adicionar adicionar novas variáveis derivadas da coluna utilizando a receita recipe_spec.
modelo_fit_arima_boosted <- arima_boost(
min_n = 2,
learn_rate = 0.015) %>%
set_engine(engine = "auto_arima_xgboost")
workflow_fit_arima_xgboost <- workflow() %>%
add_model(modelo_fit_arima_boosted) %>%
add_recipe(recipe_spec) %>%
fit(training(ts_split))Modelo 3: EXPONENTIAL SMOOTHING
Criar um modelo de Estado espaço Exponential Smoothing do tipo Error-Trend-Season.
modelo_fit_ets <- exp_smoothing() %>%
set_engine(engine = "ets")
workflow_fit_ets <- workflow() %>%
add_model(modelo_fit_ets) %>%
add_recipe(recipe_auto) %>%
fit(training(ts_split))Modelo 4: PROPHET
Modelo Prophet desenvolvido pelo facebook.
modelo_fit_prophet <- prophet_reg(seasonality_yearly = TRUE) %>%
set_engine(engine = "prophet")
workflow_fit_prophet <- workflow() %>%
add_model(modelo_fit_prophet) %>%
add_recipe(recipe_auto) %>%
fit(training(ts_split))Modelo 5: PROPHET BOOSTED
modelo_fit_prophet_boost <- prophet_boost(learn_rate = 0.1) %>%
set_engine(engine = "prophet_xgboost")
workflow_fit_prophet_boost <- workflow() %>%
add_model(modelo_fit_prophet_boost) %>%
add_recipe(recipe_spec) %>%
fit(training(ts_split))Modelo 6: Modelo Linear TSLM
workflow_fit_lm <- linear_reg() %>%
set_engine("lm") %>%
fit(l_y ~ as.numeric(DATA) + factor(month(DATA, label = TRUE), ordered = FALSE),
data = training(ts_split))Modelo 7: Modelo de Random Forest
modelo_fit_rf <- rand_forest(trees = 500, min_n = 50) %>%
set_engine(engine = "randomForest")
worflow_fit_rf <- workflow() %>%
add_model(modelo_fit_rf) %>%
add_recipe(recipe_spec) %>%
fit(training(ts_split))Adicionar Modelos Ajustados
tabela_modelos <- modeltime_table(
workflow_fit_arima,
workflow_fit_arima_xgboost,
workflow_fit_ets,
workflow_fit_prophet,
workflow_fit_prophet_boost,
worflow_fit_rf
)
tabela_modelos## # Modeltime Table
## # A tibble: 6 x 3
## .model_id .model .model_desc
## <int> <list> <chr>
## 1 1 <workflow> ARIMA(2,1,3)(0,0,1)[13]
## 2 2 <workflow> ARIMA(4,1,1)(1,0,0)[13] W/ XGBOOST ERRORS
## 3 3 <workflow> ETS(A,N,A)
## 4 4 <workflow> PROPHET
## 5 5 <workflow> PROPHET W/ XGBOOST ERRORS
## 6 6 <workflow> RANDOMFORESTCalibrando o Modelo
Aqui são calculados as métricas de performance do modelo a partir das previsões feitas para fora da amostra da base de treinamento. Assim, são quantificado os erros e estimados os intervalos de confiança.
tabela_calibracao <- tabela_modelos %>%
modeltime_calibrate(new_data = testing(ts_split))
tabela_calibracao## # Modeltime Table
## # A tibble: 6 x 5
## .model_id .model .model_desc .type .calibration_da~
## <int> <list> <chr> <chr> <list>
## 1 1 <workflow> ARIMA(2,1,3)(0,0,1)[13] Test <tibble [25 x 4~
## 2 2 <workflow> ARIMA(4,1,1)(1,0,0)[13] W/ XGBOOS~ Test <tibble [25 x 4~
## 3 3 <workflow> ETS(A,N,A) Test <tibble [25 x 4~
## 4 4 <workflow> PROPHET Test <tibble [25 x 4~
## 5 5 <workflow> PROPHET W/ XGBOOST ERRORS Test <tibble [25 x 4~
## 6 6 <workflow> RANDOMFOREST Test <tibble [25 x 4~Visualizar a Previsão
Com os dados calibrados, podemos visualziar as previsões sobre o período de teste.
tabela_calibracao %>%
modeltime_forecast(
new_data = testing(ts_split),
actual_data = df_semanal) %>%
plot_modeltime_forecast()Medidas de Performance
O Modelo Prophet parece produzir o menor valor de Mean-Absolute-Error (MAE). Mas ainda muito alto.
tabela_calibracao %>%
modeltime_accuracy() %>%
table_modeltime_accuracy()Reajustar o Modelo
Reajusta ro modelo para o banco de dados completo e realizar previsões para as datas futuras.
tabela_refit <- tabela_calibracao %>%
modeltime_refit(data = df_semanal)
tabela_refit %>%
modeltime_forecast(h = "3 months", actual_data = df_semanal) %>%
plot_modeltime_forecast()Resample
Nosso resultado pode ser sensível ao uso de uma janela específica de treinamento e teste. Podemos avaliar se os resultados se modificam quando utilizamos janelas diferentes. Usamos a função time_serse_cv para produzir 4 splits test/training diferentes, com o treinamento sempre possuindo 12 meses de duração, e o período de teste, 6 meses.
resamples_tscv <- df_semanal %>%
time_series_cv(initial = "12 months", # training window
assess = "6 months", # assessment window
skip = "6 months", # shift between resample sets
slice_limit = 4) # no. of resamplesAs quatro janelas podem ser observadas abaixo:
resamples_tscv %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(DATA, l_y)Usando o tabela_modelos e o resample_tscv cada modelo é internamente refitado para cada novo training set.
library(modeltime.resample)
resamples_fitted <- tabela_modelos %>%
modeltime_fit_resamples(
resamples = resamples_tscv,
control = control_resamples(verbose = FALSE)
)Avaliando os resultados
O modelo linear parece muito sensível ao recorte temporal, produzindo um MAE muito alto no slice 3.
resamples_fitted %>%
plot_modeltime_resamples(
.point_size = 3,
.point_alpha = 0.8,
.interactive = FALSE
)
Tabela de Precisão
A tabela abaixo mostra os resultados médios para os 4 recortes. De maneira geral, o modelo ETS tem o menor MAE médio e o menor RMSE médio.
resamples_fitted %>%
modeltime_resample_accuracy(summary_fns = mean) %>%
table_modeltime_accuracy()