markdown: wrap: 72 —

1 Introdução

Em muitas empresas, produtos sazonais ou de baixa rotatividade representam um desafio para a previsão de demanda. Esse tipo de produto frequentemente apresenta comportamento intermitente ou lumpy, com longos períodos de ausência de vendas seguidos por picos inesperados. Esses padrões dificultam o uso de modelos estatísticos tradicionais, que dependem de pressupostos como sazonalidade bem definida ou estacionariedade.

Uma abordagem comum para lidar com esse tipo de série é tentar estabilizar a variância ou aplicar transformações que tornem a série estacionária. Embora essas transformações melhorem o desempenho de modelos como ARIMA ou ETS, elas nem sempre são suficientes para capturar a complexidade das séries. Em muitos casos, mesmo após essas manipulações, os modelos falham em prever corretamente os momentos de maior impacto, como picos de vendas.

Neste estudo, mostramos que modelos de aprendizado de máquina — em especial o XGBoost — podem superar essas limitações quando aplicados de forma global e com previsão recursiva. O objetivo é ir além da compreensão dos padrões históricos e propor uma abordagem preditiva robusta que maximize a acurácia, mesmo em contextos com ruído elevado, baixa frequência e ausência de tendências claras.

2 Sobre os dados

O conjunto de dados utilizado neste estudo refere-se às vendas de itens em quatro lojas distintas ao longo de vários anos. Observa-se um padrão recorrente de interrupção nas vendas entre os meses de janeiro e setembro, caracterizando uma forte sazonalidade operacional. As atividades de vendas são retomadas com intensidade no último trimestre do ano, apresentando picos notáveis durante o período natalino. Essa característica torna o padrão de vendas altamente intermitente e concentrado, o que impõe desafios adicionais à modelagem e previsão da demanda

A estrutura do dataset está como se segue

Exemplo dos dados
Data Vendas Loja
2023-01-01 0 Loja_1
… … …
2024-09-01 100 Loja_1
2024-09-01 86 Loja_2 
vendas_lojas <- read_excel("vendas_lojas.xlsx")
data <- vendas_lojas
data %>%
  plot_time_series(
    .date_var     = date,
    .value        = vendas,
    .facet_vars   = id,
    .facet_ncol   = 2,
    .smooth       = FALSE,
    .legend_show  = TRUE,
    .interactive  = FALSE,
    .title        = "Figura 1 - Vendas por Loja (2022 - 2024)"
  )

3 Avalição dos Dados

3.1 Tipo de Demanda

Antes de prosseguir com a modelagem, é fundamental entender o perfil da demanda com a qual estamos lidando. Para isso, utilizamos duas métricas amplamente reconhecidas na literatura de previsão de demanda: ADI (Average Demand Interval) e CV² (Coeficiente de Variação ao Quadrado)

\(CV^2\): \[(\frac{\sigma}{\mu})^2\]

Onde:

  • \(\sigma\) é o desvio padrão da demanda

  • \(\mu\) é a média da demanda

    \(ADI\): \[(\frac{T}{N})\]

Onde:

  • \(T\) é total de periodos observáveis

  • \(N\) é o total de períodos de demanda positiva diferente de 0 (\(y_t\) > 0)

\[ \text{Demand Profile} = \begin{cases} \text{Smooth Demand} & \text{if } \text{ADI} < 1.32 \text{ and } CV^2 < 0.49 \\ \text{Intermittent Demand} & \text{if } \text{ADI} \geq 1.32 \text{ and } CV^2 < 0.49 \\ \text{Erratic Demand} & \text{if } \text{ADI} < 1.32 \text{ and } CV^2 \geq 0.49 \\ \text{Lumpy Demand} & \text{otherwise} \end{cases} \]

Essa classificação permite uma compreensão antecipada do comportamento da série temporal, orientando a escolha de abordagens adequadas de modelagem, seleção de variáveis explicativas e transformações necessárias. Em especial, perfis intermitentes ou lumpy demandam estratégias específicas que diferem substancialmente daquelas aplicadas à demanda suave.

3.2 Identificação da Estacionariedade da Série

Outro passo essencial na análise de séries temporais é a verificação de sua estacionariedade. Esse diagnóstico é crucial, pois define como podemos esperar que a série se comporte no futuro — em termos de tendência, sazonalidade e variabilidade. Por exemplo, se a série não for estacionária, significa que suas propriedades estatísticas (como média e variância) mudam ao longo do tempo, o que pode indicar a presença de tendências ou mudanças estruturais. Isso afeta diretamente a escolha do modelo e a interpretação dos resultados. De maneira simplificada, uma série estacionária é aquela cujas propriedades estatísticas são constantes ao longo do tempo. Se a média estiver aumentando gradualmente, se os picos e vales tornarem-se mais intensos, ou se a variância apresentar crescimento, é provável que estejamos lidando com uma série não estacionária. Uma forma formal de avaliar a estacionariedade é por meio da detecção de raízes unitárias. A presença de uma raiz unitária indica que a série é não estacionária, ou seja, que 1 é uma solução da equação de processo estocástico autoregressivo que governa a série. Em termos práticos, esse conceito está associado ao fato de que, para uma série ser estacionária, a distribuição conjunta de qualquer subconjunto de observações, como:\[y_t, y_{t+1}, y_{t+2},..., y_{t+k}\] não deve depender do tempo \(t\), para todo lag \(k\).

3.2.1 Teste de Dickey-Fuller Aumentado (ADF)

No R, o teste mais comumente utilizado para verificar a presença de raiz unitária é o Augmented Dickey-Fuller Test (ADF), implementado na função adf.test() do pacote tseries.

adf.test(sua_serie) do pacote tseries

  • Se p-valor < 0,05: rejeita-se a hipótese nula de presença de raiz unitária ⇒ a série é provavelmente estacionária.

  • Se p-valor ≥ 0,05: não há evidência suficiente para rejeitar a hipótese de raiz unitária ⇒ a série possivelmente não é estacionária.

Entretando, a decisão sobre a estacionariedade não deve se basear exclusivamente no resultado do teste estatístico. A inspeção visual da série, da média e variância ao longo do tempo, bem como do histograma das diferenças (primeiras diferenças, por exemplo), é essencial. Se as diferenças estiverem concentradas em torno de zero ou assumirem uma distribuição aproximadamente normal, isso pode indicar estacionariedade após a diferenciação

3.3 Identificação de Correlações

3.3.1

3.3.2 ACF e PACF

Entender a dependência temporal da série é uma etapa essencial tanto na construção de modelos tradicionais (como ARIMA e SARIMA) quanto em modelos de aprendizado de máquina. Através da análise das funções de autocorrelação, conseguimos determinar quais valores passados influenciam diretamente o valor atual, orientando a escolha de defasagens (lags) relevantes para modelagem.

3.3.3

3.3.4 Autocorrelação (ACF)

A Função de Autocorrelação (ACF — Autocorrelation Function) mede o grau de correlação linear entre os valores da série e suas defasagens. Em outras palavras, ela indica quanto os valores passados ajudam a explicar o valor presente da série. Para isso, utilizamos o gráfico gerado por:

acf(series)

Esse gráfico apresenta, para cada defasagem \(k\), o valor da correlação entre \(y_t\) e \(y_{t-k}\). Os valores variam de -1 a 1 e os lags considerados estatisticamente significativos são aqueles cujas barras ultrapassam as linhas tracejadas azuis (intervalos de confiança). Esses lags são candidatos fortes para inclusão em modelos autoregressivos.

Exemplo: Se o lag 7 apresenta uma correlação significativa, isso pode indicar uma sazonalidade semanal — ou seja, o valor atual é fortemente influenciado pelo valor de sete dias atrás.

Observações importantes:

  1. O lag 0 sempre terá correlação igual a 1, pois é o valor correlacionado com ele mesmo.

  2. O gráfico de ACF apresenta correlações para lags positivos. Lags negativos são simétricos, mas não são usualmente utilizados em modelos preditivos unidirecionais.

3.3.5 Autocorrelação Parcial (PACF)

A Função de Autocorrelação Parcial (PACF — Partial Autocorrelation Function) mede a correlação entre o valor atual \(y_t\) e um valor defasado \(y_{t-k}\), eliminando o efeito intermediário das defasagens entre eles. Em termos práticos, ela mostra a influência direta de um lag específico, isoladamente.

pacf(serie)

Se, por exemplo, apenas o lag 7 apresenta uma autocorrelação parcial significativa, isso indica que o valor de sete dias atrás possui contribuição direta na explicação do valor atual, enquanto os demais lags não trazem informação adicional relevante.

Essa análise é particularmente útil na identificação da ordem dos termos autoregressivos (AR) ou de médias móveis (MA) em modelos clássicos, mas também pode ser usada na engenharia de variáveis para modelos de machine learning, ao selecionar lags informativos.

#Função para calcular o ADI
calc_ADi <- function(x) {
  demand_pos <- which(x > 0)
  if(length(demand_pos) < 2) return(NA_real_)
  mean(diff(demand_pos))
}

# Função para calcular CV^2 (coeficiente de variação ao quadrado)
calc_CV2 <- function(x) {
  if(length(x) == 0) return(NA_real_)
  m <- mean(x)
  s <- sd(x)
  if(m == 0) return(NA_real_)
  (s / m)^2
}

#Classificação
classify_demand <- function(adi, cv2) {
  if (is.na(adi) || is.na(cv2)) return("Unknown")
  if (adi < 1.32 && cv2 < 0.49) return("Smooth")
  if (adi >= 1.32 && cv2 < 0.49) return("Intermittent")
  if (adi < 1.32 && cv2 >= 0.49) return("Erratic")
  if (adi >= 1.32 && cv2 >= 0.49) return("Lumpy")
  return("Unknown")
}

#Montagem de tabela para cada ID/Grupo
class <- data %>%
  group_by(id) %>%
  summarise(
    ADI = calc_ADi(vendas),
    CV2 = calc_CV2(vendas),
    Class = classify_demand(ADI, CV2)
  )


#Visualizando Qual tipo de série para cada grupo de forma visual
ggplot(class, aes(x = ADI, y = CV2, color = Class)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_text_repel(aes(label = id, colour = "black"), size = 3, max.overlaps = 20) +
  geom_vline(xintercept = 1.32, linetype = "dashed") +
  geom_hline(yintercept = 0.49, linetype = "dashed") +
  scale_color_manual(values = c(
    "Smooth" = "#2ECC71",
    "Intermittent" = "#3498DB",
    "Erratic" = "#E67E22",
    "Lumpy" = "#E74C3C"
  )) +
  labs(
    title = "Figura 2 - Classificação de Séries por ADI e CV²",
    x = "Average Demand Interval (ADI)",
    y = "Coeficiente de Variação ao Quadrado (CV²)",
    color = "Tipo de Demanda"
  ) +
  theme_minimal()

Veja que estamos lidando com demandas lumpy, com características bem parecidas com intermediárias, lembre-se que o aspecto visual é sempre de extrema importância para definir, também, nossa série temporal

#Vamos começar selecionando as colunas da categoria, data e valor e renomear elas
data_tbl <- data %>%
  select(id = id, date = date, value = vendas)

#Gráfico para avaliar as diferenças y - yt-1
##Ele ajuda a identificar a média de mudanças (força) de um dia para outro

#Vamos usar apenas um exemplo já que podemos ver VISUALMENTE que as séries são quase idênticas

data_tbl_loja1 <- data_tbl %>% 
  filter(id == "Loja_1")
diffs <- diff(data_tbl_loja1$value)
qtde_aumentos <- sum(diffs > 0, na.rm = TRUE)
qtde_quedas <- sum(diffs < 0, na.rm = TRUE)
qtde_zeros <- sum(diffs == 0, na.rm = TRUE)
list(
  aumentos = qtde_aumentos,
  quedas = qtde_quedas,
  zero = qtde_zeros
)
## $aumentos
## [1] 175
## 
## $quedas
## [1] 193
## 
## $zero
## [1] 727
diffs <- diff(data_tbl_loja1$value)

dens <- density(diffs)

mean(diffs)
## [1] 0.2547945
plot(dens, 
     main = "Figura 3 - Densidade das Diferenças",
     xlab = "Δ Vendas", 
     ylab = "Densidade",
     col = "tomato", 
     lwd = 2)

#Teste Dicky-Fuller
adf.test(data$vendas)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data$vendas
## Dickey-Fuller = -4.9845, Lag order = 16, p-value = 0.01
## alternative hypothesis: stationary
#Você pode controlar os lags utilizados em definindo k nos args

adf.test(data$vendas, k = 12)
## 
##  Augmented Dickey-Fuller Test
## 
## data:  data$vendas
## Dickey-Fuller = -5.074, Lag order = 12, p-value = 0.01
## alternative hypothesis: stationary
#Isso te ajuda a construir features para o modelo, escolher lags até onde a
##série se mantém estacionária pode ajudar muito o modelo

Veja que o resultado do teste de raiz unitária mostra que não temos evidência de não acreditar que não estamos lidando com dados estacionários, isso é, não esperamos mudanças abruptas em longo prazo, isso em si mesmo já é um insight

Temos também que as diferenças entre vendas dia a dia se concentram em torno de zero com mais variações positivas

library(tseries)
data_tbl_loja1 <- data_tbl %>% 
  filter(id == "Loja_1")

data_tbl_loja2 <- data_tbl %>% 
  filter(id == "Loja_2")

data_tbl_loja3 <- data_tbl %>% 
  filter(id == "Loja_3")

data_tbl_loja4 <- data_tbl %>% 
  filter(id == "Loja_4")

#Exemplo da Loja_1
##as demais seguem a mesma estrutura
p1 <- ggAcf(data_tbl_loja1[["value"]], lag.max = 30) + 
  ggtitle("ACF") + 
  theme_minimal(base_size = 8)

p2 <- ggPacf(data_tbl_loja1[["value"]], lag.max = 30) + 
  ggtitle("PACF") + 
  theme_minimal(base_size = 8)

p1 + p2 + plot_annotation(title = "Figura 4 - ACF e PACF da Série Temporal")

As análises gráficas confirmam que estamos lidando com séries intermitentes, o que limita a eficácia de modelos tradicionais baseados em suposições de continuidade e regularidade. Além disso, observamos estacionariedade nos lags mais elevados, sugerindo que o comportamento da série tende a se manter ao longo do tempo, sem fortes tendências ou mudanças estruturais esperadas.

Nos gráficos de autocorrelação, não há indícios claros de não estacionariedade — apenas uma leve queda nas correlações, possivelmente influenciada por uma frequência maior de vales em comparação a picos.

A função de autocorrelação parcial (PACF) indica que os lags de 1 a 4 possuem correlações significativas com o valor atual, o que será crucial na construção do modelo. Além disso, não foi identificada uma sazonalidade dominante, o que reduz a necessidade de componentes sazonais explícitos.

3.4 Métricas de Avaliação

Para avaliar a eficácia dos modelos propostos, utilizamos diferentes métricas, cada uma com um papel específico:

3.4.1 1. RMSE (Root Mean Squared Error)

O RMSE mede o erro médio das previsões em unidades da variável original. Ele penaliza fortemente grandes erros, sendo adequado para comparar desempenho geral entre modelos.

\[ RMSE = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (y_i - \hat{y}_i)^2} \]

Motivo de uso: permite comparar modelos globais (ex: XGBoost) e tradicionais (ex: ARIMA) de forma padronizada, sensível a grandes erros.

3.4.2 2. MAAPE (Mean Arctangent Absolute Percentage Error)

O MAAPE é uma métrica robusta para dados com valores pequenos ou intermitentes, pois suaviza os efeitos dos erros relativos extremos que afetam o MAPE tradicional.

\[ MAAPE = \frac{1}{n} \sum_{i=1}^{n} \arctan\left(\left|\frac{y_i - \hat{y}_i}{y_i}\right|\right) \]

Motivo de uso: ideal para séries com comportamento lumpy ou intermitente, pois evita distorções comuns no MAPE.

3.4.3 3. \(R^2\) (Coeficiente de Determinação)

O \(R^2\) indica a proporção da variância da variável dependente que é explicada pelo modelo.

\[ R^2 = 1 - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]

Motivo de uso: mede a qualidade do ajuste no conjunto de treino e teste, útil principalmente em modelos de machine learning.

3.4.4 4. Trade-off de Viés e Variância

Modelos simples tendem a ter alto viés e baixa variância, enquanto modelos complexos apresentam baixo viés e alta variância. O ideal é encontrar o equilíbrio que minimize o erro total:

\[ Erro\ total = Viés^2 + Variância + \text{Ruído irreducível} \]

Motivo de uso: ao avaliar modelos com aprendizado de máquina, é essencial entender como a complexidade afeta a generalização do modelo.

3.5 Treinando Modelos Tradicionais

Como ponto de partida para a comparação entre abordagens, iniciamos a análise com o ajuste de modelos locais com métodos tradicionais de previsão, que tradicionalmente oferecem maior acurácia ao serem personalizados para cada série individual. Neste estudo, utilizamos três modelos amplamente consolidados na literatura de séries temporais:

  • ARIMA (AutoRegressive Integrated Moving Average)

  • ETS (Error, Trend, Seasonal Models)

  • Prophet (modelo aditivo com tendência e sazonalidade, desenvolvido pelo Facebook)

O objetivo aqui é avaliar o desempenho desses modelos em séries intermitentes e com forte seasson, como observado nos dados em questão.

A implementação foi realizada com o pacote modeltime do R, que oferece uma interface unificada para treinamento, calibração e comparação de modelos de séries temporais. Este pacote facilita a integração entre diferentes frameworks e permite realizar avaliações com métricas padronizadas, como RMSE, MAE, MAPE e R².

Agora com os dados transformados, vamos criar uma tabela extendida, passo para o workflow do modeltime para previsão de modelos com funções locais

#Agora, seguindo o workflow do Modeltime para construção de local models
##vamos criar um nested_dataset

nested_data_tbl <- data_tbl %>%
  
  # 1. Extending: Vamos predizer 52 dias no futuro
  extend_timeseries(
    .id_var        = id, #nosso id (categoria -> nome da coluna)
    .date_var      = date, #nosso date (data -> nome da coluna)
    .length_future = 365 #quantos passos no futuro vamos prever 
  ) %>%
  
  # 2. Nesting: Vamos agrupar por id e criar dataset futuro
  #    selecionamos quantidade de passos futuros vamos prever
  #    selecionamos o tamanho do dataset original
  nest_timeseries(
    .id_var        = id,
    .length_future = 365,
    .length_actual = length_actual
  ) %>%
  
  # 3. Splitting: Por fim criar um split de dados de treino e teste
  #    isso precisa ser feito para medir acurácia pelo teste
  #    e o resto para treino
  split_nested_timeseries(
    .length_test = 365
  )

Criado os dados extendidos, vamos definir todos os modelos, para esse vamos usar a função workflow() para cada modelo, cada workflow vai receber um add_model() onde vamos definir quais modelos:

  • prophet_reg("regression") para prophet
  • arima_reg() %>% set_engine("auto_arima") para AUTOARIMA (define parâmetros sozinho)
  • exp_smoothing() %>% set_engine("ets") para o ETS

E vamos adicionar em cada um, uma receita criada pela função recipe() onde vamos detalhar o fit()value ~ date que diz que vamos prevêr receita por data apenas e colocar uma função extract_nested_train_split() para dividir nosso dataset extendido em dados de treino e teste

No fim de cada workflow vamos adicionar o add_recipe() e passar o nome da recipe

#Modelo Prophet do Modeltime
rec <- recipe(value ~ date, extract_nested_train_split(nested_data_tbl)) 

wflw_prophet <- workflow() %>%
  add_model(
    prophet_reg("regression", seasonality_yearly = TRUE) %>% 
      set_engine("prophet")
  ) %>%
  add_recipe(rec)

#Modelo Arima do Modeltime
wflw_arima <- workflow() %>%
  add_model(
    arima_reg() %>%
      set_engine("auto_arima")
  ) %>%
  add_recipe(rec)

#Modelo ETS do Modeltime
wflw_ets <- workflow() %>%
  add_model(
    exp_smoothing() %>%
      set_engine("ets")
  ) %>%
  add_recipe(rec)

Feito isso, o próximo passo do workflow é montar uma tabela de modelos onde vamos chamar a função modeltime_nested_fit() e nela vamos passar o arg nested_data onde vamos colocar nossa tabela, e adicionar nossos modelos

Após isso vamos treinar o modelo e capturar suas acurácias com extract_nested_test_accuracy() e plotar em uma tabela usando table_modeltime_accuracy()

#Agora vamos juntar esses modelos dentro de uma tabela para o workflow 
nested_modeltime_tbl <- modeltime_nested_fit(
  # Nested data 
  nested_data = nested_data_tbl,
  
  # Add workflows models
  wflw_prophet,
  wflw_arima,
  wflw_ets
)
Resultados da Acurácia dos Modelos Tradicionais
ID Modelo RMSE RSQ
Loja_1 Prophet 100.35 0.70
Loja_1 Arima 264.87 0.00
Loja_1 ETS 196.40 -
Loja_2 Prophet 102.70 0.67
Loja_2 Arima 285.70 0.00
Loja_2 ETS 205.66 -
Loja_3 Prophet 98.62 0.71
Loja_3 Arima 191.97 0.00
Loja_3 ETS 184.19 -
Loja_4 Prophet 95.95 0.72
Loja_4 Arima 258.95 0.00
Loja_4 ETS 179.11 - 
#Plotando resultado
nested_modeltime_tbl %>% 
  extract_nested_test_forecast() %>%
  group_by(id) %>%         # garante que só séries com pelo menos 200 pontos fiquem
  plot_modeltime_forecast(
    .facet_ncol  = 2,
    .interactive = FALSE,
    .title = "Figura 5 - Previsão teste dos Modelos Tradicionais"
  )

3.5.1 Resultados dos Modelos Tradicionais

A avaliação dos modelos locais tradicionais — ARIMA, ETS e Prophet — revelou limitações importantes ao lidar com séries de demanda intermitente e altamente ruidosas. Mesmo com um destaque com Propeht é correto afirmar que esses modelos, por sua natureza, tendem a suavizar os dados e se ajustar fortemente à média, o que compromete a capacidade de capturar variações esparsas e picos isolados de demanda.

Diante dessas limitações, optamos por explorar abordagens baseadas em aprendizado de máquina, especialmente métodos de boosting, que possuem maior flexibilidade para capturar padrões não lineares e aprender a partir dos erros residuais dos modelos anteriores.

4 Modelos Globais

Como discutido por Montero-Manso e Hyndman (2021), modelos globais se diferenciam dos modelos locais por aplicarem uma única função de previsão a todas as subséries do conjunto de dados. Essa abordagem contrasta com os modelos locais, que ajustam uma função específica para cada série individualmente. Modelos globais têm demonstrado, em muitos contextos, desempenho comparável — e até superior — ao dos modelos locais, especialmente quando há padrões compartilhados entre as séries.

Essa vantagem ocorre porque modelos globais são capazes de capturar dinâmicas comuns e estruturas complexas presentes em múltiplas séries, permitindo maior generalização e estabilidade nos resultados, mesmo em situações de dados esparsos ou intermitentes.

4.1 Seleção e Pré-processamento das Séries para o Modelo Global

Antes de ajustar o modelo global, é fundamental realizar uma curadoria cuidadosa das séries que serão incluídas, garantindo que possuam informações suficientes e características compatíveis com o aprendizado conjunto. As seguintes estratégias foram adotadas para essa filtragem:

  • Remoção de séries descontinuadas: produtos que deixaram de vender após um curto período ou que apresentaram apenas vendas pontuais foram descartados.

  • Exclusão de séries com pouco histórico ou datas de início e fim muito distintas em relação às demais.

As regras específicas adotadas foram:

  1. Threshold de dias com demanda positiva: Selecionamos apenas séries com mais de 30 observações com \(y_t > 0\), garantindo uma base mínima de dados informativos para o modelo aprender.

  2. Padronização dos intervalos de tempo: Todas as séries foram ajustadas para possuir o mesmo intervalo de datas, com preenchimento dos períodos faltantes com zero, assegurando consistência temporal.

  3. Filtro por pico mínimo de demanda: \[max_t y_t >= 10\] Selecionamos apenas séries com pelo menos um dia com demanda significativa:Filtro por média de atividade: \[\frac{1}{T} \Sigma y_t > 1\], onde \(T\) é total de pontos de datas isso calcula a média de movimentação por período o que garante que vamos selecionar apenas séries ativas. Esse critério assegura que a série apresenta atividade consistente ao longo do tempo.

4.2 Modelagem com XGBoost

Dentre os métodos baseados em boosting, o XGBoost (Extreme Gradient Boosting) destaca-se como uma das técnicas mais robustas e amplamente utilizadas em tarefas de previsão. O princípio central do boosting é o ajuste iterativo de modelos fracos, de forma que cada novo modelo aprenda a corrigir os erros dos anteriores. Isso permite capturar não apenas o sinal da série, mas também seus resíduos, tornando o modelo especialmente útil em contextos com ruído elevado. Além do XGBoost, o framework inclui variantes como AdaBoost e Gradient Boosting Machines (GBM). No entanto, o XGBoost é amplamente reconhecido por seu desempenho e eficiência computacional, sendo uma escolha popular em competições e aplicações do mundo real. Apesar de seu poder preditivo, o XGBoost também apresenta alta complexidade, o que exige atenção especial à sua configuração para evitar overfitting. Entre os hiperparâmetros mais relevantes, destacam-se:

  • tree (nrounds) esse parâmetro controla quantas iterações o modelo vai executar para reduzir o erro. Mais árvores podem melhorar a performance no treino, mas aumentam o risco de overfitting se não houver controle adequado por eta, max_depth e técnicas de regularização.

  • max_depth termina o nível máximo de splits que cada árvore pode ter. Profundidades maiores permitem que o modelo capture interações complexas entre variáveis, mas também aumentam o risco de overfitting.

  • eta (learning rate): atua como um fator de encolhimento aplicado aos pesos $\alpha_m$ de cada modelo base. Esse parâmetro ajuda a reduzir o risco de overfitting ao frear o impacto de cada árvore adicionada, tornando o processo de aprendizagem mais gradual e controlado.

Esses e outros hiperparâmetros (como subsample, colsample_bytree, entre outros) podem ser ajustados com rigor, preferencialmente por meio de validação cruzada ou estratégias de tuning automatizado, para garantir um equilíbrio adequado entre viés e variância.

Veja mais parâmetros em: https://xgboost.readthedocs.io/en/latest/parameter.html

4.3 Criação de Features

Para modelos de machine learning alguns passos anteriores da construção do modelo são importantes

Precisamos construír variáveis que vão fornecer insights ao modelo de ML para a previsão de passos futuros. Vamos explicar os mais tradicionais métodos

4.3.1 1. Lags

Podemos entregar para cada tempo \(y_t\) o seu valor anterior no lag \(k\) como \(y_{t-k}\). Isso mostra para o modelo que podemos prever hoje como uma função de lags anteriores. Aqui os gráficos PACF pode ajudar a identificar quais lags são mais significativos para o trabalho

Lags curtos, 1 a 6 passos passados, ajudam a identificar correlações imediatas, lags semanais \(y_{t-7}\), \(y_{t-14}\) ajudam a identificar padrões semanais e lags longos como acima de \(y_{t-21}\) ajudam a capturar padrões longos. Não é recomendado usar todos porque inflar o modelo pode acabar enviesando ou overfittando o modelo, então use o ACF e PACF para definir quais lags usar

4.3.2 2. Janelas Móveis

Ajudam a suavizar a série e explicar a tendência da série. Temos algumas opções como:

  1. Log Trend \[L_t = ln(t)\]. Ajuda a identificar padrões completos da série

  2. Médias Móveis curtas \[MA_{t,n} = \frac{1}{n} \Sigma y_i\] onde n = 2 a 7, evidencia padrões de mudança semanais

  3. Médias Móveis de termos médios, mesma formula de médias móveis curtas, ajuda a identificar padrões de tempos médios como 14 a 21 dias

4.3.3 3. Diferenciações de Primeira e Segunda Ordem

Diferenciações ajudam a identificar direção de mudanção e força das mudanças da série

Diferenciação de primeira ordem \[\Delta y_t = y_t - y_{t-1}\], ajuda ao modelo entender mudanças de direções de um dia para outro

Diferenciação de segunda ordem \[\Delta^2 y_t = \Delta y_t - \Delta y_{t-1} = (y_t - y_{t-1}) - (y_{t-1} - y_{t-2})\], mostram a força dessas mudanças

Em séries estacionárias ou com poucas mudanças de um dia pra outro, podem inflar modelo com viés equivocado de mudanças e gerar falhas estruturais

4.3.4 4. Sazonalidade (Senos e Cossenos)

Aplicaçãos desses cálculos provê insights poderosos ao modelos para capturar padrões ciclicos que se repetem em tempos regulares

\[\sin\frac{2\pi ft}{360}\] e \[\cos\frac{2\pi ft}{360}\]

4.3.5 5. Features de Calendário

Identificar como qual mês a data é, se é feriado (0 para não) (1 para sim), identificar semana, ciclos de férias, etc, podem dar informações para o modelo desenhar melhor as estrtuturas de relacionamento

Para isso criamos variáveis binárias ou multiclasse

4.4 Modelo Recursivo Global

Na modelagem preditiva de séries temporais com machine learning, a construção de features baseadas em defasagens (lags), janelas móveis (rolling windows) e variáveis sazonais é essencial. No entanto, se aplicarmos essas features em cenários padrões, previsão multistep, o modelo não terá as informações para trabalhar janelas mais a frente, por exemplo, \(y_{t+2}\) não terá as informações das features de \(y_{t+1}\), para isso vamos abordar a estratégia de modelo recursivo.

4.4.1 Forecast Recursivo

Modelos recursivos realizam previsões passo a passo, onde a saída prevista em \(y_{t+1}\) é usada como entrada (via lags e janelas) para prever \(y_{t+2}\), e assim por diante. Isso é especialmente importante quando as features utilizadas no treinamento dependem de valores históricos da própria variável alvo.

Se o modelo não for executado de forma recursiva, ele tentará utilizar valores futuros reais (inexistentes) para gerar as features — o que resultará em valores ausentes (NaNs) nas previsões, inviabilizando o processo. Por isso, no nosso fluxo de trabalho com o pacote modeltime, é essencial garantir que o modelo seja ajustado e avaliado de forma recursiva, replicando as condições reais de previsão.

4.4.2 Estratégia de Transformações (Features)

Com base nas análises prévias da série (ACF, PACF e diagnósticos de estacionariedade), definimos as seguintes transformações para o modelo global recursivo:

Lags curtos: Apesar da ausência de padrões persistentes em longas defasagens, observamos correlação significativa até o lag 14, identificado no PACF. Portanto, incluiremos esse lag no conjunto de preditores.

Janelas móveis (rolling features): Utilizaremos a média móvel de curto prazo (ex.: 7 dias) para capturar pequenas tendências sazonais de curta duração — como oscilações semanais. Essa escolha é sustentada pela leve queda observada na função de autocorrelação, mesmo com a série sendo predominantemente estacionária.

Variáveis sazonais de calendário: Incorporaremos variáveis de calendário como mês, semana e distância até o Natal, para ajudar o modelo a identificar os efeitos sazonais pronunciados no período de fim de ano — que, conforme identificado, representam os principais picos de demanda no histórico.

Essas transformações visam oferecer ao modelo informações relevantes e realistas para capturar tanto padrões locais quanto efeitos compartilhados entre as séries, respeitando a lógica de previsão recursiva.

4.4.3 Treinamento e Validação

Antes de avançarmos para a previsão de múltiplos passos à frente, é fundamental avaliar a capacidade de generalização do modelo, garantindo que ele não apenas se ajuste bem aos dados históricos, mas também tenha bom desempenho em dados futuros. Essa etapa é crucial para evitar overfitting, especialmente ao utilizar algoritmos complexos como os métodos de boosting.

4.4.3.1 Overfitting em Séries Temporais

O overfitting ocorre quando o modelo aprende padrões específicos e até ruídos dos dados de treino, resultando em excelente desempenho histórico, mas baixa capacidade preditiva. Isso é particularmente comum em séries estacionárias ou altamente ruidosas, onde pequenas flutuações podem ser interpretadas como padrões reais. Modelos potentes, como o XGBoost, são especialmente suscetíveis a esse problema devido à sua alta flexibilidade.

Para mitigar esse risco, adotamos uma estratégia de validação temporal rigorosa, comparando o desempenho do modelo em dados de treino e teste.

4.4.3.2 Estratégia de Validação

A validação foi realizada com o auxílio do pacote modeltime, que oferece suporte direto para técnicas robustas de amostragem temporal.

Utilizamos a função time_series_split() para particionar os dados em conjunto de treino e teste, preservando a ordenação temporal — aspecto essencial na modelagem de séries temporais.

#antes precisamos mudar o id para factor
data_tbl <- data_tbl %>%
  dplyr::mutate(id = as.factor(id))

splits <- time_series_split(
  data_tbl, 
  assess     = "1 year", 
  cumulative = TRUE
)

Agora vamos criar um recipe criando as features que esperamos colocar em nosso modelo

#No recipe, indicamos os dados como training(splits) para usar nos dados de treino
rec_obj <- recipe(value ~ ., training(splits)) %>%
  #step_mutate_at(id) para tirarmos o id (lembra que criaremos um modelo global)
  step_mutate_at(id, fn = droplevels)%>%
  #step_lag() cria os lags que desejamos
  step_lag(value, lag = c(1,2,3,4,14)) %>%
  #step_mutate_rollmean criamos a media de janela móvel
  step_mutate(
    rollmean = slider::slide_dbl(lag(value, 14), mean, .before = 11, .complete = FALSE)
  ) %>%
  #step_rm() nos excluimos a coluna date
  step_rm(date) %>%
  #step_zv retiramos todas as features sem variação, útil para reduzir a dimensionalidade
  step_zv(all_predictors()) %>%
  #por fim usamos o step_summy() para, caso exista, transformamos features string em dummys
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

Agora vamos treinar nosso modelo com machine learning, utilizando o mesmo fluxo aplicado aos modelos tradicionais. Como já temos todas as features (como lags, janelas móveis e variáveis de calendário) criadas tanto para os dados de treino quanto de teste, não há necessidade de aplicar forecast recursivo neste ponto.

4.4.3.2.1 Fluxo de Modelagem
  1. Criação do Workflow
  2. Adicionar à Tabela do Modeltime modeltime_table()
  3. Calibrar o Modelo

4.4.4 Viés x Variância

O parâmetro eta (também conhecido como learning rate) controla o ritmo de aprendizado do modelo a cada iteração. Valores muito altos de eta podem fazer o modelo aprender rápido demais, aumentando o risco de overfitting — isto é, o modelo se ajusta excessivamente aos dados de treino e perde capacidade de generalização. Já valores muito baixos exigem mais árvores para alcançar uma boa performance, podendo causar underfitting se o número de árvores for insuficiente.

Viés elevado ocorre quando o modelo é simples demais e não consegue capturar a complexidade dos dados, resultando em underfitting. Por outro lado, variância elevada acontece quando o modelo é muito complexo, ajustando-se excessivamente aos dados de treino, o que causa overfitting e prejudica o desempenho em dados novos.

Em vez de analisar o trade-off entre viés e variância apenas em função do eta, estamos agora testando diversas combinações dos parâmetros eta, tree_depth (profundidade das árvores) e trees (quantidade de árvores). Para avaliar o impacto dessas combinações no desempenho do modelo, utilizamos a métrica RMSE nos conjuntos de treino e teste.

Dessa forma, buscamos identificar a combinação de hiperparâmetros que melhor equilibra a complexidade do modelo, reduzindo o erro e garantindo boa capacidade de generalização.

# Grid expandido de hiperparâmetros
param_grid <- expand.grid(
  eta        = c(0.1, 0.25, 0.5),
  tree_depth = c(3, 6, 9),
  trees      = c(500, 1000, 1500)
)

# Lista para armazenar modelos
modelos_treinados <- list()

for (i in seq_len(nrow(param_grid))) {
  
  params <- param_grid[i, ]
  
  modelo <- boost_tree(
    mode = "regression",
    learn_rate = params$eta,
    trees = params$trees,
    tree_depth = params$tree_depth
  ) %>%
    set_engine("xgboost")
  
  wf <- workflow() %>%
    add_model(modelo) %>%
    add_recipe(rec_obj)
  
  fit_modelo <- wf %>% fit(training(splits))
  
  modelos_treinados[[i]] <- list(
    modelo = fit_modelo,
    params = params
  )
}

metricas_list <- list()

for (i in seq_along(modelos_treinados)) {
  modelo_info <- modelos_treinados[[i]]
  modelo_fit <- modelo_info$modelo
  params <- modelo_info$params
  
  model_tbl <- modeltime_table(modelo_fit)
  
  calib_tbl <- model_tbl %>%
    modeltime_calibrate(new_data = testing(splits), id = "id")
  
  metricas <- calib_tbl %>%
    modeltime_accuracy(metric_set = extended_forecast_accuracy_metric_set()) %>%
    mutate(
      eta        = params$eta,
      tree_depth = params$tree_depth,
      trees      = params$trees
    )
  
  metricas_list[[i]] <- metricas
}

# Combinar resultados
metricas_df <- bind_rows(metricas_list)


plot_hiperparametro <- function(df, param) {
  ggplot(df, aes_string(x = param, y = "rmse")) +
    geom_point(color = "blue") +
    stat_summary(fun = mean, geom = "line", aes(group = 1), color = "red") +
    labs(title = paste(param, "vs RMSE"), x = param, y = "RMSE") +
    theme_minimal()
}

p1 <- plot_hiperparametro(metricas_df, "eta")
p2 <- plot_hiperparametro(metricas_df, "tree_depth")
p3 <- plot_hiperparametro(metricas_df, "trees")

# Combine com patchwork
(p1 | p2) / (p3) +
  plot_annotation(
    title = "Figura 6 - Tunning de Hiperparâmetros XGBOOST vs RMSE",
    theme = theme(plot.title = element_text(face = "bold", hjust = 0.5))
  )

Combinações de Hiperparâmetros XGBOOST
eta tree_depth trees RMSE
0.10 9 500 81.76
0.10 9 1000 81.76
0.10 9 1500 82.27

Durante os testes, foi observado que a combinação learning_rate (ou eta) a de 0.1, tree_depth de 9 e 500 trees

O RMSE alcançado foi de ~81, mas podemos notar também que o hiperparâmetro tree não influência de forma demasiada o resultado do modelo

Com isso, os hiperparâmetros finais escolhidos foram:

wflw_xgb <- workflow() %>%
  add_model(
    boost_tree(
      mode = "regression",
      learn_rate = 0.1,    # eta
      tree_depth = 9,      # profundidade da árvore
      trees = 500         # número de árvores
    ) %>%
    set_engine("xgboost")
  ) %>%
  add_recipe(rec_obj) %>%
  fit(training(splits))


model_tbl <- modeltime_table(
  wflw_xgb
)

calib_tbl <- model_tbl %>%
  modeltime_calibrate(
    new_data = testing(splits), 
    id       = "id",
    quiet = FALSE
  )

calib_tbl_treino <- model_tbl %>%
  modeltime_calibrate(
    new_data = training(splits),
    id =    "id"
  )

Com o modelo treinado, vamos ver nossos resultados com o modeltime_accuracy() e table_modeltime_accuracy()

Primeiro vamos ver o modelo global

Resultado do modelo Global com Tunning
.model RMSE RSQ
XGBOOST 81.33 0.80

Dado no modelo de forma global, tivemos uma queda no RMSE comparado ao modelos tradicionais

metrics_train <- calib_tbl_treino %>%
  modeltime_accuracy(
    acc_by_id = TRUE,
    metric_set = extended_forecast_accuracy_metric_set()
  ) %>%
  mutate(dataset = "Treino")

metrics_test <- calib_tbl %>%
  modeltime_accuracy(
    acc_by_id = TRUE,
    metric_set = extended_forecast_accuracy_metric_set()
  ) %>%
  mutate(dataset = "Teste")

metrics_compare <- bind_rows(metrics_train, metrics_test)


ggplot(metrics_compare, aes(x = id, y = maape, fill = dataset)) +
  geom_col(position = "dodge") +
  labs(title = "Figura 7 -Comparação MAAPE entre Treino e Teste", y = "MAAPE", x = "ID") +
  theme_minimal()

ggplot(metrics_compare, aes(x = id, y = rsq, fill = dataset)) +
  geom_col(position = "dodge") +
  labs(title = "Figura 8 - Comparação R^2 entre Treino e Teste", y = "RSQ", x = "ID") +
  theme_minimal()

Resultado por Loja do XGBOOST parametrizado
ID MAAPE RMSE RSQ DATASET
Loja_1 104.8654 10.09058536 0.9968078 Treino
Loja_2 104.6524 0.02688104 1.0000000 Treino
Loja_3 104.6518 0.01930264 1.0000000 Treino
Loja_4 104.6515 0.01914229 1.0000000 Treino
Loja_1 120.4828 94.01560831 0.7249804 Teste
Loja_2 117.7007 75.03682879 0.8208426 Teste
Loja_3 118.9113 76.93570501 0.8260017 Teste
Loja_4 117.8049 77.90996864 0.8115844 Teste

Podemos ver que o MAAPE mostrou um leve aumento, o que é esperado, mas mantém uma boa generalização. O \(R^2\) é comum cair, é uma métrica que costuma sempre ser alta, como 100%, no treino e cair no teste, mas mesmo com a queda temos um resultado bem promissor de uma média de 80%.

Mas comparando aos modelos tradicionais, vemos que há uma grande melhora, comparando o RMSE entre ambos, o XGBOOST teve uma média de erro de ~80 no teste, já os modelos tradicionais ficaram com uma média de ~100 no Prophet, estamos falando de uma queda de 20% no erro.

Em suma, apesar de métricas como MAAPE parecerem elevadas, em séries lumpy é mais importante observar tendência de melhora entre treino e teste e capacidade de adaptação ao padrão irregular do dado.

Vamos agora visualizar as previsões usando as funções do modeltime:

p <- calib_tbl %>%
  modeltime_forecast(
    new_data    = testing(splits),
    actual_data = data_tbl,
    conf_by_id  = TRUE
  ) %>%
  group_by(id) %>%
  plot_modeltime_forecast(
    .facet_ncol  = 2,
    .interactive = FALSE,
    .conf_interval_show = TRUE
  )
  
p <- p + 
  labs(
    title    = "Figura 9 - Previsões de teste com XGBOOST - Etapa de Treino",
    subtitle = "ETA: 0.1, trees: 500, max_depth:9",
    x        = "Data",
    y        = "Vendas"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 10),
    legend.title  = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1,size = 6)
  )

suppressWarnings(print(p))

5 Modelo Final

Agora vamos, enfim, construir o modelo final com XGBoost.

Diferente das etapas anteriores, não utilizaremos mais o recipe(). Em vez disso, vamos:

  • Criar uma função de engenharia de features,

  • Gerar um dataset para datas futuras, e

  • Ajustar o modelo para que possa lidar com previsões recursivas.

Por isso, o workflow também será diferente nesta etapa.

Vamos começar criando a função responsável pela geração das features.

#Vamos construir nossas funções para a engenharia de features
lag_roll_transformer_grouped <- function(data){
  data %>%
    group_by(id) %>%
    tk_augment_lags(value, .lags = c(1,2,3,4, 14)) %>% #criando os lags
    tk_augment_slidify(
      .value   = contains("lag14"),   # janela móveis de média de termo médio
      .f       = ~mean(.x, na.rm = TRUE),
      .period  = 12,
      .partial = TRUE
    ) %>%
    fastDummies::dummy_cols(remove_selected_columns = FALSE) %>%
    ungroup()
}

#Criando features de calendário
add_calendar_vars <- function(data) {
  data %>%
    mutate(
      month_num      = month(date), #indicando número do mês
      week_iso       = isoweek(date), #indicando semana ISO
      # Semana do mês: número da semana relativa ao 1º dia do mês
      week_of_month  = floor((day(date) - 1) / 7) + 1
    )
}

Agora vamos criar nosso modelo utilizando o XGBOOST. Mas antes vamos desenvolver nossas features aplicando as funções criadas

#Vamos prever 120 dias
FORECAST_HORIZON <- 180


#Antes vamos expandir nosso dataset
#Vamos criar dados para até agosto de 2026, isso porque sabemos que são datas 
##com vendas zeradas


novas_datas <- seq.Date(from = as.Date("2025-01-01"), to = as.Date("2025-08-31"), by = "day")
ids <- unique(data_tbl$id)
dados_novos <- expand.grid(id = ids, date = novas_datas)
dados_novos <- dados_novos %>%
  mutate(value = 0)

#Unindo datasets
data_tbl_exp <- bind_rows(data_tbl, dados_novos) %>%
  arrange(id, date)


#Vamos criar o extended dataset por id
m4_extended <- data_tbl_exp %>%
  group_by(id) %>%
  future_frame(
    .length_out = FORECAST_HORIZON,
    .bind_data  = TRUE
  ) %>%
  ungroup()
## .date_var is missing. Using: date
#Vamos aplicar as regras nos dados
m4_lags <- m4_extended %>%
  lag_roll_transformer_grouped()

#Vamos splitar em dados de treino e aplicar variaveis de calendario
train_data <- m4_lags %>%
  drop_na() %>%
  add_calendar_vars()

#Vamos splitar em teste e aplicar variaveis de calendario
future_data <- m4_lags %>%
  filter(is.na(value)) %>%
  add_calendar_vars()

Depois de transformar os dados precisamos selecionar o modelo que queremos. No modeltime para usarmos o XGBOOST vamos chamar a função boost_tree() dentro dela vamos definir alguns args para os hiperparâmetros do modelo e o mode, que para forecast será “regression”

depois vamos aplicar um engine para o modelo (importante passo no workflow do modeltime), para isso vamos aplicar a função set_engine() e escolher o “xgboost”

Após isso vamos aplica um fit() e escolher os campos, aqui definimos value ~ . que significa que vamos prever a coluna value para todas as variáveis criadas criando a coluna month, e excluíndo a coluna date que não é necessária

E por fim adicionamos mais um bloco com a função recursive onde vamos indicar qual o “id” caso tenhamos grupos, o dataset futuro, quantos passos futuros também, e a função de transformação

Depois disso criamos uma tabela de modelos (outro passo do workflow) com modeltime_table() e treinamos nosso modelo com modeltime_calibrate() especificando o id para dar garantias de que vamos usar a função global para cada id e indicando o new_data como os dados futuros

Vamos também definir eta como 0.1, tree_depth = 3, trees = 500

set.seed(123)
#Construindo XGBOOST -> MOdelo de ML para predição
model_fit_xgb_recursive_1 <- boost_tree(
  mode = "regression",
  learn_rate = 0.1,       # nosso eta
  tree_depth = 9,          # substitui max_depth
  trees = 500             # ajuste do número de árvores
  
) %>%
  set_engine("xgboost") %>%
  fit(
    value ~ .
    + month(date, label = TRUE) 
    + as.numeric(date) 
    - date, 
    data = train_data
  ) %>%
  recursive(
    id         = "id", # Precisamos usar a função Recursiva aqui e identificar o id
    transform  = lag_roll_transformer_grouped,
    # Vamos colocar os dados de treino como dataset para a previsão recursiva
    train_tail = panel_tail(train_data, id, FORECAST_HORIZON)
  ) # aqui o *train_data* são dados futuros, quando definimos isso prevemos dados para frente do ja2 visto
#poderiamos, claro, no split, pegar dados ja existentes e apenas testar o modelo

#Criar tabela de modelos (no caso só temos 1)
model_tbl <- modeltime_table(
  model_fit_xgb_recursive_1
)

#Treinar e calibrar modelo
calibration_tbl <- model_tbl %>%
  modeltime_calibrate(
    new_data = train_data,
    id = "id"
  )

Agora precisamos plotar o modelo chamando, para a tabela calibrado a função modeltime_forecast() e definir, principalmente, os args new_data que seria nossa tabela future_data e o arg actual_data que seria nossa tabela data_tbl_exp

E no fim chamamos a função plot_modeltime_forecast() para ver o gráfico do resulto

Para prever qualquer resultado, bastaria trazer um dataframe como o future_data e adicionar dentro da função modeltime_forecast(), a partir dessa função podemos extrair as previsões, gráficos, IC, etc

#Avalidando resultados
p <- calibration_tbl %>%
  modeltime_forecast(
    new_data    = future_data,
    actual_data = data_tbl_exp,
    conf_interval = 0.95,
    conf_method   = "conformal_split", #gera IC baseado em quantis
    conf_by_id    = TRUE,
    keep_data   = TRUE
  ) %>%
  group_by(id) %>%
  plot_modeltime_forecast(
    .interactive        = FALSE,
    .conf_interval_show = TRUE,
    .facet_ncol         = 2
  )

p <- p + 
  labs(
    title    = "Figura 10 - Previsões com XGBOOST - Forecast Recursivo",
    subtitle = "ETA: 0.1, trees: 500, max_depth: 9",
    x        = "Data",
    y        = "Vendas"
  ) +
  theme_minimal(base_size = 10) +
  theme(
    plot.title    = element_text(face = "bold", size = 15),
    plot.subtitle = element_text(size = 10),
    legend.title  = element_blank(),
    axis.text.x = element_text(angle = 45, hjust = 1, size = 6)
  )

suppressWarnings(print(p))

6 Transfer Learning (Hard)

Uma das estratégias utilizadas para permitir previsões em diferentes níveis da hierarquia (como lojas, produtos ou regiões) é o Transfer Learning — uma técnica onde um modelo treinado em uma tarefa é reutilizado para inferência em outra, mesmo que essa nova tarefa envolva dados não observados no treinamento.

Entre as abordagens existentes, o Hard Transfer Learning é uma das mais diretas: utiliza-se o modelo pré-treinado sem qualquer readequação ou fine-tuning. Ou seja, os parâmetros aprendidos $\theta$ são mantidos fixos e aplicados diretamente a novos dados.

  • Modelo pré-treinado, onde \(x_i\) representa as features que conhecemos, \(L\) a função de perda e \(\theta\) os parâmetros criados. Ou seja, definimos \(\theta\) que minimize a função custo \(L\) para as features \(x_i\)

\[ f_\theta = arg min_\theta \Sigma_{(x_i,y_i)} L(f(x_i,\theta),y_i) \]

  • Posteriormente, ao se deparar com novos dados, usamos \(\theta\) treinado para prevêr novos pontos com novos dados \(x_j\) não vistos

    \[ \hat{y}_j = f(x_j, \theta) \]

7 Resultados

Os experimentos demonstram que o uso do pacote modeltime com abordagens de recursividade e aprendizado de máquina é eficaz para capturar padrões complexos em séries temporais ruidosas e de difícil modelagem, como as do tipo lumpy ou intermitente.

Os resultados evidenciam ganhos significativos de desempenho quando comparamos modelos baseados em aprendizado de máquina com os modelos tradicionais. Por exemplo, na Loja_1, o RMSE caiu de aproximadamente 100 (melhor modelo tradicional) para cerca de 80 com o modelo ajustado de machine learning — uma melhora de 20%.

As outras lojas tiveram um resultado semelhante com quedas mais acentuadas.

Embora modelos tradicionais como ARIMA, ETS e Prophet mantenham sua importância — especialmente em séries estacionárias e com estruturas mais simples — eles tendem a ter desempenho inferior quando o objetivo é prever com maior precisão em horizontes curtos ou em contextos de alta variabilidade.

Modelos tradicionais geralmente assumem que a média da série se mantém no longo prazo, o que pode ser útil para detectar mudanças de sinal ou tendência. No entanto, para tarefas com foco em previsão granular e acurada, os modelos baseados em aprendizado de máquina demonstram uma vantagem clara, tanto em robustez quanto em adaptabilidade.

Inclusive vale destacar que modelos globais podem superar modelos locais, a custo de uma redução minima de acurácia.

8 Referências

Chaves, M. C. S. (2024). Time Series Forecasting via Machine Learning in Large-Scale E-Commerce Using Transfer Learning and Nonparametric Bootstrap. April 24, 2024.

Nielsen, A. (1900). Análise Prática de Séries Temporais.

Business Science. (2024). Modeltime: Time Series Forecasting for Business Science. Acessado em 24 de julho de 2025, de https://business-science.github.io/modeltime/index.html