Análise de séries temporais para previsão de demanda
Click here to read this article in English.
A análise de séries temporais tem por objetivo a identificação de padrões, previsão de tendências e tomada de decisões, que podem ser aplicadas nas mais diversas áreas uma vez que auxilia o entendimento do comportamento de algo ao longo do tempo.
Nesse artigo mostro o caminho que eu fiz para encontrar o melhor modelo preditivo na demanda por vinhos entre os modelos com: abordagem Naive, Média Móvel, Holt’s Linear Trend Model, ARIMA, Prophet e PyCaret.
* Observação
Este é um relato completo do estudo, incluindo códigos e metodologias utilizadas. Caso queira, publiquei um resumo mais direto, onde trago apenas os resultados dessa pesquisa:
Para acessar o artigo resumido acesse aqui.
Sumário
- Sobre o Projeto
1.1. Previsão de Demanda
1.2. Séries Temporais
1.3. Business Understanding - Objetivo Geral
2.1. Objetivos Específicos - Obtenção dos Dados
- Dicionário de Variáveis
- Importação dos Dados e das Bibliotecas
- Análise Exploratória dos Dados
6.1. Por que utilizar Pandas Profiling?
6.2. Gerar Relatório
6.2.1. Relatório de Produtos
6.2.2. Relatório de Vendas - Tratamento dos Dados
7.1. Dataset de Produtos
7.2. Dataset de Vendas
7.3. Junção dos Datasets - Análise de Dados
8.1. Remover datas de 2016 e 2017
8.2. Vinho mais caro e mais barato - Feature Engineering
- Análise do Negócio
10.1. Vinhos: mais vendidos x maior receita
10.2. Produtores: mediana de preço x quantidade de produtos
10.3. Produtores: mais vendidos x maior receita
10.4. Região: quantidade de produto x receita - Previsão de Demanda
11.1. Preparação dos dados
11.2. Série Estacionária
11.2.1. Teste Augmented Dickey Fuller (ADF)
11.2.2. Transformar Série
11.3. Dividir o conjunto em dados de treino e de validação
11.4. Fazer previsões
11.4.1. Abordagem Naive
11.4.2. Média Móvel
11.4.3. Holt’s Linear Trend Model
11.4.4. ARIMA
11.4.5. Prophet
11.4.6. PyCaret - Avaliação dos Modelos
12.1. Reverter transformações para série estacionária
12.2. Gráficos
12.3. Métricas de Avaliação - Conclusão
1. Sobre o Projeto
1.1. Previsão de Demanda
Ser capaz de prever a demanda por um produto ou serviço é uma técnica estratégica que pode auxiliar a tomada de decisões, o planejamento, a gestão de estoque, como a otimização de recursos, por exemplo, e até mesmo a satisfação do cliente. São muitas as aplicações em que esse tipo de ferramenta pode ser aplicada e em empresas dos mais diversos setores.
Nesse estudo, vamos nos basear nos requisitos exigidos em uma oferta real de trabalho que foi postada no LinkedIn. O contrato previa apresentar uma solução de previsão das vendas para ajudar na gestão de estoque da empresa. É apontado ainda que seja considerado os desafios peculiares à indústria do vinho.
Todo projeto para realização de uma previsão de demanda exige a análise de dados históricos, modelagem estatística e o uso de algoritmos de aprendizado de máquina. Por isso, eles informam que irão fornecer dados das vendas passadas.
1.2. Séries Temporais
Esse tipo de problema é típico de ser solucionado com a análise de séries temporais, que são aqueles dados distribuídos em intervalos sequenciais e regulares ao longo do tempo. Isso significa que os dados são coletados em intervalos de tempo como: a cada hora, diariamente, mensalmente, etc. Além disso, cada novo dado é dependente dos dados anteriores. Se não fosse a conexão entre os dados e a variável temporal, o problema poderia ser solucionado com regressão linear (Analytics Vidhya, 2018).
1.3. Business Understanding
Todo projeto exige um entendimento do caso em estudo. Compreender as peculiaridades auxilia o cientista de dados na escolha dos melhores parâmetros e de algoritmos de aprendizado de máquina para que se possa construir um modelo preditivo que se encaixe no caso específico.
Nesse caso estamos lidando com um produto principal: vinhos. Eles são considerados bebidas vivas, ou seja, estão em constante evolução. Por esse motivo, eles precisam de condições ideais para que mantenham suas características e não estraguem.
Entre os cuidados necessários estão o armazenamento e o transporte, para que o estoque não seja perdido, e isso faz com que o objetivo desse projeto seja ainda mais importante. Esses cuidados podem fazer a diferença entre a empresa ter lucro ou prejuízo no longo prazo.
A luz solar precisa ser evitada, e os vinhos precisam ser mantidos em um ambiente com umidade e temperatura controladas e constantes. É mais importante que a temperatura seja constante do que qual é a temperatura. Normalmente o alvo é entre 10 e 15 graus, variando de acordo com o tipo de vinho.
Também é importante evitar movimentos como tremores, vibrações ou qualquer tipo de perturbação ao líquido, o que nos indica que o gerenciamento de estoque precisa ser feito de forma otimizada, a fim de não mexer nos vinhos mais do que o necessário. Se possível, é recomendado ainda que sejam mantidos na posição horizontal, o que também pode precisar ser levado em consideração ao planejar o depósito onde os vinhos ficarão.
Por esses motivos, é preciso pensar em otimizar a rotatividade, para que o vinho passe o mínimo de tempo possível em uma posição que possa o deteriorar, garantindo que seja entregue aos clientes ainda em um estado excelente de consumo.
Apesar de alguns vinhos ficarem melhores com o passar dos anos, a grande maioria é feita para consumo “imediato”, desde o momento em que chegam ao mercado, até no máximo 1 a 3 anos. Por isso, é vital não ficar com vinhos em estoque por um período maior do que isso.
2. Objetivo Geral
Desenvolver um algoritmo com base em aprendizado de máquina para prever a demanda por vinhos em uma loja especializada nesse produto.
2.1. Objetivos Específicos
- Realizar uma análise exploratória nos dados para conhecer o data set e verificar presença de anormalidades que necessitem de tratamentos adequados.
- Criar modelos de machine learning com diversos tipos de aprendizado de máquina.
- Avaliar modelos para encontrar aquele que tem o melhor desempenho.
3. Obtenção dos Dados
Os dados utilizados nesse projeto foram disponibilizados pelo Rafael Duarte, cientista de dados e enófilo, que criou esse conjunto dada a dificuldade de se encontrarem dados de vendas desse setor. Ele pontua que:
- A carta de vinhos, ou seja, os produtos, são reais, baseados na oferta real de um e-commerce de vinhos aqui no Brasil. Nomes, safras, e valores são 100% reais, convertidos pra dólar para ter um apelo mais internacional.
- O conjunto de dados de venda é baseado em uma competição do Kaggle. Nele haviam 5 anos de vendas diárias, distribuídas por 10 lojas, com catálogo de 50 produtos. Esse conjunto foi alterado para 3 anos de vendas diárias, distribuídas em 3 lojas e com 219 produtos diferentes em estoque.
Os dados estão separados em dois arquivos: um com os dados históricos das vendas e o outro com informações sobre os vinhos. Ambos foram salvos em nuvem para o caso de indisponibilidade ou alterações futuras e, serão esses arquivos salvos por mim que eu utilizo no código, para que seja mantida a integridade do estudo.
Vale destacar que, por se tratar de um conjunto com dados sintéticos, ainda que com base em dados reais, é possível nos depararmos com problemas que não seriam encontrados com dados reais. Isso deverá ser um ponto de atenção durante toda a análise.
4. Dicionário de Variáveis
A compreensão do conjunto de dados passa pela checagem das variáveis disponíveis nele para que se possa realizar uma boa análise. Segue-se um resumo dos atributos encontrados e seus respectivos significados. Como temos dois conjuntos produtos.csv e vendas.csv, suas variáveis encontram-se em duas tabelas distintas.
produtos.csvem ordem alfabética
country: país de origem do vinhoitem_id: número de identificação do itemkind: tipo do vinho:
-sparkling: espumante
-rose sparkling: espumante rosé
-white: branco
-rosé: rosé
-red: tintoname: nome do vinhoprice_brl: preço em reaisprice_usd: preço em dólar americanoproducer: nome do produtor do vinhoregion: região de produção do vinhovintage: ano da safra
O conjunto de vendas possui também o atributo item_id pelo qual pode-se retornar ao conjunto de produtos e obter mais detalhes sobre o item comprado pelo cliente.
vendas.csvem ordem alfabética
date: data da vendaitem: número de identificação do itemstore: lojasales: quantidade vendida
5. Importação dos Dados e das Bibliotecas
Ao iniciar um projeto é necessário instalar pacotes, importar as bibliotecas que possuem funções específicas a serem utilizadas nas linhas de código seguintes e realizar as configurações necessárias para a saída do código. Também, se prossegue com a importação do dataset, salvando-o em uma variável específica para que seja usada posteriormente.
# instalar pacote adicional
!pip install pycaret -q # auto machine-learning
!pip install pandas_profiling -q # produção de relatório de análise exploratória
!pip install pmdarima -q # auto arima
!pip install pycaret -q # pycaret# importar as bibliotecas necessárias
import pandas as pd # manipulação de dados
import numpy as np # manipulação de arrays
import matplotlib.pyplot as plt # visualização de dados
import seaborn as sns # visualização estatística dos dados
import datetime as dt # manipulação de datas e horários
from pandas_profiling import ProfileReport # produção de relatório de análise exploratória
from statsmodels.tsa.stattools import adfuller # teste adf (série estacionária)
from sklearn.metrics import mean_absolute_error as mae # métrica de avaliação
from sklearn.metrics import mean_absolute_percentage_error as mape # métrica de avaliação
from statsmodels.tsa.seasonal import seasonal_decompose # decompor série temporal
from statsmodels.tsa.holtwinters import Holt # método de previsão Holt's Linear Trend Model
from pmdarima.arima import auto_arima # autoarima
from prophet import Prophet # séries temporais
from pycaret.time_series import *
import warnings # notificações
warnings.filterwarnings('ignore') # configurar notificações para serem ignoradas
# configurações adicionais
plt.style.use('ggplot')
sns.set_style('dark')# importar os conjuntos de dados e atribuí-los em variáveis
df_products_raw = pd.read_csv('https://www.dropbox.com/scl/fi/t8wr73843nxvcw949twi3/products.csv?rlkey=qrms57xtev2rzqzpxuqoc9yif&dl=1')
df_sales_raw = pd.read_csv('https://www.dropbox.com/scl/fi/hlrspaogkay44qe45bbf4/sales.csv?rlkey=j9rlkn2r97j3swfpg9irckvgm&dl=1')6. Análise Exploratória dos Dados
Essa é uma etapa essencial em projetos de ciência de dados onde se busca compreender melhor os dados seja por identificar padrões, outliers, possíveis relações entre as variáveis, etc. Nesse estudo, se exploraram informações que fossem relevantes para orientar as respostas dos objetivos indicados anteriormente (ver Objetivo Geral e Objetivos Específicos).
Para isso vou gerar um relatório que resume os dados, com o uso da biblioteca ProfileReport. A partir dele, se houver necessidade, será feito uma análise mais profunda. Contudo, ele já nos fornecerá informações suficientes para identificar anomalias como outliers e desbalanceamento nos dados.
6.1. Por que utilizar Pandas Profiling?
A biblioteca Pandas Profiling oferece a função ProfileReport que cria automaticamente um relatório completo de análise exploratória de dados. Ele inclui as seguintes informações:
- Visão Geral
Com a quantidade de registros e atributos do dataset, bem como valores ausentes e duplicados. - Estatística Descritiva
Para cada variável irá calcular: média, mediana, valor máximo e mínimo, desvio padrão, quartis, entre outros dados estatísticos. - Distribuição
Para as variáveis numéricas, ele irá criar gráficos como histogramas, boxplot, densidade, para que seja visualmente fácil compreender a distribuição dos dados. - Frequência
Para os atributos numéricos será feita a identificação dos principais valores únicos, e suas respectivas contagens. - Correlação
É gerada uma matriz de correlação para verificação do relacionamento entre as variáveis.
Por isso, o uso do ProfileReport agiliza a rotina da análise exploratória dos dados, contudo, não exime o cientista de dados de realizar o estudo aprofundado do relatório para compreender e interpretar os dados, além de identificar eventuais problemas que necessitam de tratamento.
6.2. Gerar Relatório
Nessa etapa vou gerar os relatórios e, em seguida, imprimir a visualização deles para que se possa analisar as estatísticas do conjunto.
# criar relatórios
report_products = ProfileReport(df_products_raw)
report_sales = ProfileReport(df_sales_raw)6.2.1. Relatório de Produtos
# visualizar relatório do dataset 'products'
report_products.to_notebook_iframe()No relatório gerado com o código acima percebe-se que:
- o conjunto é formado por 219 registros e 9 atributos
- pelas primeiras e últimas entradas do dataset ele parece estar bem preenchido e sem anormalidades
- não há dados ausentes
- não há dados duplicados
- existe 1 coluna numérica e as demais são categóricas, o que merece um checagem uma vez que por ter variáveis de preço elas podem não estar no tipo de variável correta
- foram detectados 4 atributos com alta cardinalidade, ou seja, com muitos valores distintos, o que é condizente com os atributos em questão (
name,producer,price_brleprice_usd) - foram detectas correlações fortes entre
item_ideproducer,countryeitem_id,regioneitem_id,vintageekind,producerekind - as variáveis
item_id,vintage,price_brleprice_usdestão como tipostringmas não é o tipo mais adequado para elas - na variável
countrytemos ao todo 6 países diferentes, onde a França representa 70% do conjunto, seguido pela Itália com 10% e Espanha com 8% - na variável
regiontemos Burgundy representando 32% do conjunto, seguido por Bordeaux 21% e Rhone 6% - na variável
vintagetemos ao todo 17 safras diferentes, sendo a safra de 2018 representando 30% do conjunto, seguido por 2017 com 22%, e tanto 2015 como 2014, cada um com 8.7% de representação - ainda sobre o atributo
vintagenota-se a presença do acrônimo NV que significa Non-Vintage, será necessário alterar esse valor para 0 para poder alterar todo o atributo em tipoint. - na variável
kindtemos o vinho tinto representando 60% do conjunto, seguido pelo branco com 31.5% e o espumante com 4.6% - na variável
producertemos ao todo 58 produtores diferentes, sendo o mais presente nesse conjunto o Domaine Ponsot com 5.5%, seguido por La Chablisienne com 4.6% e Domaine Matrot com 4.1%
Dessa forma, os tratamentos necessários envolvem a alteração de tipos das variáveis item_id para str, price_brl e price_usd para tipo float. Na variável vintage antes de converter para int, é preciso converter os valores NV para 0.
A alteração de tipo em item_id servirá para evitar algum tipo de associação hierárquica pelo modelo de machine-learning.
6.2.2. Relatório de Vendas
# visualizar relatório do dataset 'sales'
report_sales.to_notebook_iframe()No relatório gerado com o código acima percebe-se que:
- o conjunto é formado por 720071 registros e 4 atributos
- pelas primeiras e últimas entradas do dataset ele parece estar bem preenchido e sem anormalidades
- não há dados ausentes
- não há dados duplicados
- existem 2 colunas numéricas e 2 categóricas, o que merece um checagem uma vez que por ter variável de data ela pode não estar no tipo correto
- as datas estão em registros diários, mas não seguem um padrão e não estão em formato
datetime - olhando pelos registros iniciais e finais, aparentemente, o período de registro nesse conjunto é de 1º de janeiro de 2018 até 31 de dezembro de 2020
- além disso, as datas se repetem para cada item vendido
- a variável
storeestá em formatostringe está balanceada - a variável
itemestá em formatofloate temos 219 valores distintos, o que corresponde à quantidade de itens do conjunto de produtos - a variável
salesestá em formatofloat
Precisamos converter o atributo item para str e, sales para int. Também será necessário converter a variável date para o tipo datetime.
Como mencionado no relatório de produtos, a alteração de tipo em item servirá para evitar algum tipo de associação hierárquica pelo modelo de machine-learning.
7. Tratamento dos Dados
Nessa etapa serão realizados os procedimentos necessários para prosseguirmos com a modelagem do algoritmo e, que foram identificados na seção anterior.
7.1. Dataset de Produtos
No dataset de produtos serão feitos os seguintes tratamentos:
- alterar
item_idparastr - alterar
price_brleprice_usdparafloat - na variável
vintageconverter valoresNVpara 0 - alterar
vintageparaint
Começamos realizando uma cópia do dataset para identificar que temos um conjunto que foi alterado e, ainda, mantemos o conjunto original intacto.
# fazer cópia do dataset 'products'
df_products_clean = df_products_raw.copy()# alterar 'item_id' para tipo 'str'
df_products_clean['item_id'] = df_products_clean['item_id'].astype('str')
# checar alteração
type(df_products_clean['item_id'][0])
'''
str
'''# alterar 'price_brl' e 'price_usd' para 'float'
#df_products_clean['price_brl'] = df_products_clean['price_brl'].astype('float')
#df_products_clean['price_usd'] = df_products_clean['price_usd'].astype('float')
# checar alteração
type(df_products_clean['price_brl'][0])
'''
str
'''No código acima temos um erro na conversão de tipo pois foi encontrado valor com vírgula, dessa forma, primeiro vou retirar a vírgula e depois, converter o tipo da variável para float.
# excluir vírgulas
df_products_clean.replace(',', '', regex=True, inplace=True)
# alterar 'price_brl' e 'price_usd' para 'float
df_products_clean['price_brl'] = pd.to_numeric(df_products_clean['price_brl'])
df_products_clean['price_usd'] = pd.to_numeric(df_products_clean['price_usd'])
# checar alteração
type(df_products_clean['price_brl'][0])
'''
numpy.float64
'''# converter valores 'NV' do atributo 'vintage' para 0
df_products_clean['vintage'].replace('NV', 0, inplace=True)
df_products_clean.head()
# alterar 'vintage' para tipo 'int'
df_products_clean['vintage'] = df_products_clean['vintage'].astype('int64')
# checar alteração
type(df_products_clean['vintage'][0])
'''
numpy.int64
'''7.2. Dataset de Vendas
No dataset de vendas serão feitos os seguintes tratamentos:
- alterar
itemparastr - alterar
salesparaint - alterar
dateparadatetime.
Começamos realizando uma cópia do dataset para identificar que temos um conjunto que foi alterado e, ainda, mantemos o conjunto original intacto.
# fazer cópia do dataset 'sales'
df_sales_clean = df_sales_raw.copy()# alterar 'item_id' para tipo 'str'
df_sales_clean['item'] = df_sales_clean['item'].astype('str')
# checar alteração
type(df_sales_clean['item'][0])
'''
str
'''# alterar 'sales' para 'int'
df_sales_clean['sales'] = df_sales_clean['sales'].astype('int64')
# checar alteração
type(df_sales_clean['sales'][0])
'''
numpy.int64
'''# alterar 'date' para 'datetime'
df_sales_clean['date'] = pd.to_datetime(df_sales_clean['date'])
# checar alteração
type(df_sales_clean['date'][0])
'''
pandas._libs.tslibs.timestamps.Timestamp
'''7.3. Junção dos Datasets
Para facilitar uma análise mais minuciosa desses dados, vou unir os dois datasets products e vendas.
# renomar atributo 'item' para 'item_id' e ficar igual ao dataset 'products'
df_sales_clean.rename(columns= {'item': 'item_id'}, inplace=True)# unir os datasets 'sales' e 'products' de acordo com o atributo 'item_id'
df_merged = df_products_clean.merge(df_sales_clean, on='item_id', how='right')
# imprimir tamanho do dataset
print(df_merged.shape)
# verificar as primeiras entradas do dataset
df_merged.head()
Podemos confirmar que a junção foi bem sucedida pois temos os mesmos 720071 registros contidos no dataset sales e, quanto aos atributos, tínhamos 9 em products e 4 em sales, como ambos continham a variável de relação item_id, ela só é contabilizada uma vez e, por isso, esperava-se encontrar 12 atributos, que é a quantidade exata mostrada acima.
Vamos também fazer a checagem dos registros finais do dataset, bem como se o tipo das variáveis foi mantido.
# verificar últimas entradas em 'df_merged'
df_merged.tail()
# verificar tipo das variáveis
df_merged.info()
8. Análise de Dados
Uma vez que foi realizado o tratamento dos dados, podemos gerar o relatório pelo Pandas Profiling novamente, já que temos atributos relevantes como os que se referem aos valores monetários e vendas, que antes, por estarem no tipo errado, não foi possível realizar nenhuma análise.
# criar relatório
report = ProfileReport(df_merged)
# visualizar relatório
report.to_notebook_iframe()No relatório gerado com o código acima percebe-se que:
- o conjunto é formado por 720071 registros e 12 atributos
- pelas primeiras e últimas entradas do dataset ele parece estar bem preenchido e sem anormalidades
- não há dados ausentes
- não há dados duplicados
- existe 1 coluna
datetime, 5 numéricas e 6 categóricas - assim como antes foi detectada a cardinalidade de alguns atributos e a alta correlação entre alguns deles
- o atributo
vintageagora possui 3.7% dos campos com valor 0, devido ao tratamento realizado nessa variável - em
datetemos alguns períodos de 2016 e 2017 que surgiram e não possuem registros
Cabe destacar algumas estatísticas importantes:
As demais informações referentes à safra, país, região e produtores, permanecem as mesmas.
8.1. Remover datas de 2016 e 2017
Como o dataset refere-se à dados de 3 anos, provavelmente esses dados são remanescentes da construção sintética. Vou proceder com a exclusão desses dados referentes à 2016 e 2017.
# remover dados referentes à 2016 e 2017
df_clean = df_merged.loc[(df_merged['date'].dt.year != 2016) & (df_merged['date'].dt.year != 2016)]# verificar exclusão, ordenando os dados do mais antigo para o mais recente
df_clean.sort_values(by='date')
df_clean.head()
8.2. Vinho mais caro e mais barato
Vamos identificar o vinho mais caro e o mais barato.
# identificar vinho mais caro
df_clean[df_clean['price_usd'] == 1901.73].iloc[0]
# identificar vinho mais barato
df_clean[df_clean['price_usd'] == 9.13].iloc[0]
9. Feature Engineering
O processo de feature engineering consiste na criação de novas variáveis a partir dos dados que já existem com o objetivo de melhorar o desempenho de um modelo de machine-learning.
Nesse caso, como temos a variável de tempo e faremos uma análise de Série Temporal, é interessante separá-la em quantos atributos individuais conseguirmos, ou seja, em dia do mês, dia da semana, dia do ano, semana do ano, mês, trimestre e ano.
# criar atributos a partir da data
df_clean['Year'] = df_clean.date.dt.year # ano
df_clean['Quarter'] = df_clean.date.dt.quarter # trimestre
df_clean['Month'] = df_clean.date.dt.month # mês
df_clean['Week'] = df_clean.date.dt.week # semana do ano
df_clean['Weekday'] = df_clean.date.dt.weekday # dia da semana
df_clean['Day'] = df_clean.date.dt.day # dia do mês
df_clean['Dayofyear']= df_clean.date.dt.dayofyear # dia do ano
# verificar alterações
df_clean.head()
Agora que temos informações do dia da semana (0 corresponde à segunda-feira e, assim por diante até terminar em 6, que se refere ao domingo) vou construir um novo atributo para distinguir dia da semana e final de semana.
# criar atributo 'Weekend' com valores igual a 0
df_clean['Weekend'] = 0
# converter valores correspondentes ao final de semana para 1
## sábado = 5
## domingo = 6
df_clean.loc[(df_clean.Weekday == 5) | (df_clean.Weekday == 6), 'Weekend'] = 1
# verificar alteração
df_clean.head()
# verificar registros finais
df_clean.tail()
Vamos verificar visualmente a quantidade média de vendas durante a semana e aos finais de semana.
# plotar gráfico com valor médio de vendas durante a semana x final de semana
df_clean.groupby('Weekend').sales.mean().plot.bar();
Nota-se que são valores médios muito próximos.
Outra informação que podemos extrair é o valor total por dia que cada produto gerou de receita. Para isso podemos multiplicar o atributo de preço do vinho que é dado por price_usd e por price_brl, pela quantidade de unidades vendidas, dado por sales.
# criar novo atributo 'total_amount_usd' e total_amount_brl' para informar a receita
df_clean['total_amount_usd'] = df_clean['price_usd'] * df_clean['sales']
df_clean['total_amount_brl'] = df_clean['price_brl'] * df_clean['sales']
# verificar alterações
df_clean.head()
Vale checar as principais estatísticas desses novos atributos.
# imprimir relatório estatístico
df_clean.describe().round(2)
Em total_amount_usd e total_amount_brl observa-se que o valor médio gasto é de 12138 usd ou 70158, contudo dado o alto valor de desvio padrão, temos uma mediana com valor mais baixo, de 5051 usd ou 29198 brl. Isso pode ser visto ainda pela diferença entre o 3º quartil e o valor máximo, indicando a presença de outliers. Contudo, vou manter esses valores no dataset uma vez que eles eventualmente vão ocorrer e interferir igualmente futuramente, dessa forma, é importante que façam parte dos dados para a análise de série temporal que será realizada.
10. Análise do Negócio
Até esse momento, fizemos a análise exploratória dos dados com o objetivo de conhecer o conjunto de dados disponibilizado e checar a presença de anormalidades que pudessem prejudicar tanto o entendimento do dataset quanto à construção do modelo preditivo. Fizemos o devido tratamento e adicionaram-se variáveis relevantes para uma análise de Série Temporal, como o desmembramento do atributo referente às datas.
Contudo, ainda não foi possível ter uma compreensão mais profunda dos dados, com algo que seja relevante para a empresa conhecer melhor sua demanda e como usar isso para melhorar suas vendas, ofertas, estoque e consequentemente aumentar seus lucros.
Nessa seção vou analisar especificamente a oferta, para compreender as regiões que mais vendem em relação à volume e valor dos produtos. Esse estudo se faz relevante devido ao fato de que podemos ter um vinho que venda 10 garrafas diariamente e que custe 10 dólares, o que corresponderia a 700 dólares por semana. Por outro lado, podemos ter uma garrafa a 700 dólares e que tenha uma unidade sendo vendida por semana. Portanto, em termos de receita, são equivalentes. Logo, qual deles seria melhor ter em estoque? Seria possível vender mais do vinho de 10 dólares? Ou do de 700? São essas perguntas que tentaremos responder com uma análise mais minuciosa.
10.1. Vinhos: mais vendidos x maior receita
Primeiro, vou construir dois conjuntos diferentes:
- os 20 vinhos que geraram mais receita (e, então, coletar a informação de unidades vendidas)
- os 20 vinhos mais vendidos por unidades (e, então, coletar a informação de receita obtida)
Com isso, posso plotar um gráfico de dispersão que mostre a relação entre a quantidade de receita gerada versus a quantidade de unidades vendidas dos 20 vinhos. São dois conjuntos separados uma vez que, ranqueando por receita nem sempre o vinho estará entre o top 20 com unidades mais vendidas e vice e versa.
# top 20 vinhos que mais geraram receita
top_amount = df_clean.groupby('name')['total_amount_usd'].sum().sort_values(ascending=False)
top_amount_sales = df_clean.groupby('name')['sales'].sum().sort_values(ascending=False)
top_amount_merged = pd.merge(top_amount, top_amount_sales, on='name', how='left')
top_amount_merged['rank'] = range(1, (len(top_amount)+1))
top_amount_merged[:20]
# top 20 vinhos mais vendidos por unidade
top_sales = df_clean.groupby('name')['sales'].sum().sort_values(ascending=False)
top_sales_amount= df_clean.groupby('name')['total_amount_usd'].sum().sort_values(ascending=False)
top_sales_merged= pd.merge(top_sales, top_sales_amount, on='name', how='left')
top_sales_merged['rank'] = range(1, (len(top_amount)+1))
top_sales_merged[:20]
# plotar gráfico de dispersão
# os 20 vinhos mais vendidos por unidade x maior receita
# configurar gráfico
top20_amount = top_amount_merged.iloc[0:20]
top20_sales = top_sales_merged.iloc[0:20]
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(x=top20_amount['total_amount_usd'].values, y=top20_amount['sales'].values, c='red', s=70, alpha=0.4)
ax.scatter(x=top20_sales['total_amount_usd'].values, y=top20_sales['sales'].values, c='blue', s=20, alpha=0.4)
ax.set_title('Top 20 vinhos mais vendidos')
ax.set_xlabel('Receita')
ax.set_ylabel('Unidades vendidas')
ax.legend(['Receita', 'Unidades Vendidas'], loc='upper center')
plt.show()
Os pontos dos grupos foram plotados em cores e tamanhos distintos para facilitar a leitura e, com isso, pode-se observar no gráfico acima que temos dois pontos isolados, no lado direito e no alto, que corresponde ao top 1 de unidades vendidas e top 1 de receitas. Logo, é um vinho que tem muitas unidades vendidas e que representa a maior receita da empresa.
Ao voltar nas listas acima, podemos confirmar essa informação e verificar que se trata do vinho ‘Domaine Ponsot Chapelle-Chambertin Grand Cru’, ou, pode-se confirmar essa informação com o código abaixo.
# imprimir top 1 vinho com mais unidades vendidas e o top 1 que gerou maior receita
print('O top 1 vinho com mais unidades vendidas: {}\n'
'O top 1 vinho que gerou mais receita: {}'.format(top_sales_merged.index[0], top_amount_merged.index[0]))
'''
O top 1 vinho com mais unidades vendidas: Domaine Ponsot Chapelle-Chambertin Grand Cru
O top 1 vinho que gerou mais receita: Domaine Ponsot Chapelle-Chambertin Grand Cru
'''Vou checar qual o preço unitário desse vinho:
# identificar preço do vinho mais vendido e que mais gera receita
top1_price = df_products_clean.loc[df_products_clean.name == 'Domaine Ponsot Chapelle-Chambertin Grand Cru']
print('O vinho Domaine Ponsot Chapelle-Chambertin Grand Cru custa US$ {} por unidade.'.format(top1_price['price_usd'].values[0]))
'''
O vinho Domaine Ponsot Chapelle-Chambertin Grand Cru custa US$ 656.06 por unidade.
'''Conforme visto na análise exploratória dos produtos temos que esse é um dos vinhos mais caros que a empresa vende, uma vez que seu preço encontra-se já no percentil 95, o que significa que é um valor maior do que 95% dos outros preços de vinhos do nosso conjunto de dados.
Outro exemplo que merece ser destacado no gráfico de dispersão é o segundo ponto mais alto no gráfico, azul e do lado esquerdo. Isso nos mostra que, apesar dele ser o top 2 vinho em unidades vendidas, ele gera uma receita muito mais baixa se comparada aos demais.
Vamos checar isso.
# imprimir top 2 vinho com mais unidades vendidas
print('O top 2 vinho com mais unidades vendidas: {}'.format(top_sales_merged.index[1]))
'''
O top 2 vinho com mais unidades vendidas: Potensac
'''# checar a colocação dele em relação ao top 20 de receita
# buscar ranking do vinho top2 'Potensac'
top2_in_amount = top_amount_merged.loc[top_amount_merged.index == 'Potensac']
print(top2_in_amount['rank'])
'''
name
Potensac 33
Name: rank, dtype: int64
'''# identificar preço do vinho top2 em unidades vendidas
top2_in_amount_price = df_products_clean.loc[df_products_clean.name == 'Potensac']
print('O vinho Potensac custa US$ {} por unidade.'.format(top2_in_amount_price['price_usd'].values[0]))
'''
O vinho Potensac custa US$ 112.32 por unidade.
'''Temos que o vinho ‘Potensac’ que é o 2º em quantidade de garrafas vendidas está em 33º lugar quando se trata de receita e custa 112 dólares. Esse valor, de acordo com a nossa análise exploratória dos produtos está entre a mediana, de 88 usd e o 3º quartil, de 164 usd, o que nos indica que é um preço mediano, mais para caro. Ou seja, não é um dos vinhos mais em contas para ser comprado entre as opções disponíveis pela empresa.
10.2. Produtores: mediana de preço x quantidade de produtos
Outra relação interessante de se olhar é a quantidade de produtos diferentes que cada produtor oferece. A depender do caso, pode ser relevante traçar uma estratégia para se conseguir comprar em maior quantidade e obter um desconto maior, aumentando o lucro em alguns produtos, por exemplo.
Lembrando que a empresa atualmente tem negócio com 58 produtores diferentes.
Vou construir dois conjuntos diferentes:
- os 20 produtores com maior mediana de preço por garrafa de vinho
- os 20 produtores com maior quantidade de produtos comercializados pela empresa, por produtor
Utilizou-se a mediana por ser uma estatística menos sensível aos outliers.
A partir dessas informações vou plotar um gráfico de dispersão que mostre a relação entre as maiores medianas de preço por produtor versus a quantidade de produtos comercializados por produtor.
# produtores: quantidade de produto x preço mediano
producer_median_price = df_clean.groupby('producer')['price_usd'].median().sort_values(ascending=False)
producer_items = df_clean.groupby('producer')['item_id'].nunique().sort_values(ascending=False)
# juntar informações
producer_mp = pd.merge(producer_median_price, producer_items, on='producer', how='left')
producer_it = pd.merge(producer_items, producer_median_price, on='producer', how='left')# imprimir lista com top 20 produtores por mediana de preço
producer_mp_top20 = producer_mp.iloc[:20]
print(producer_mp_top20)
# imprimir lista com top 20 produtores por quantidade de produtos
producer_it_top20 = producer_it.iloc[:20]
print(producer_it_top20)
# plotar gráfico de dispersão
# configurar gráfico
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(x=producer_mp_top20['price_usd'].values, y=producer_mp_top20['item_id'].values, c='red', s=70, alpha=0.4)
ax.scatter(x=producer_it_top20['price_usd'].values, y=producer_it_top20['item_id'].values, c='blue', s=20, alpha=0.4)
ax.set_title('Top 20 produtores')
ax.set_xlabel('Mediana de preço')
ax.set_ylabel('Quantidade de produtos')
ax.legend(['por Mediana de Preço', 'Quantidade de produto'], loc='upper center')
plt.show()
plt.show()
Podemos observar dois pontos sobrepostos, na extrema direita e alto do gráfico. Isso quer dizer que provavelmente temos um produtor cuja mediana de preço é a mais alta e que ele também é o que possui maior quantidade de produtos. Abaixo, faremos a identificação desses dois pontos para confirmar essa hipótese.
# identificar produtor top1 em mediana de preço e o top1 em quantidade de produtos
print('O top 1 produtor com maior mediana de preço: {}\n'
'O top 1 produtor com maior quantidade de produtos: {}'.format(producer_median_price.index[0], producer_items.index[0]))
'''
O top 1 produtor com maior mediana de preço: Domaine Ponsot
O top 1 produtor com maior quantidade de produtos: Domaine Ponsot
'''Confirmamos que o produtor com a maior mediana de preço também é o que possui a maior quantidade de produtos. Ademais, conforme visto anteriormente, é esse produtor que possui o vinho com maior receita gerada para a empresa e que é o vinho também mais vendido (por garrafa).
Com isso, pode-se pensar em diversas estratégias para a empresa, como, por exemplo:
- identificar outros produtores na mesma região do Domaine Ponsot
- buscar por vinhos semelhantes ao Domaine Ponsot Chapelle-Chambertin Grand Cru
- oferecer outros produtos do Domaine Ponsot que ainda não estão sendo oferecidos pela empresa
- criar experiências desse produtor como kits de degustação, por exemplo
10.3. Produtores: mais vendidos x maior receita
Assim como verificamos os vinhos mais vendidos e que mais geram receitas, podemos fazer a mesma checagem com relação aos produtores. Isso é interessante para, com outras informações como as preferências de consumo dos clientes, se possam criar estratégias que irão consequentemente aumentar as vendas.
Para essa análise vou construir dois conjuntos:
- os 20 produtores que geraram mais receita
- os 20 produtores com maior quantidade de itens vendidos
Em seguida, vou criar um gráfico de dispersão que mostre a relação entre a quantidade de receita gerada versus os produtores com maior quantidade de unidades vendidas. São dois conjuntos separados uma vez que, ranqueando por receita nem sempre o produtor estará entre o top 20 com unidades mais vendidas e vice e versa.
# top 20 produtores que mais geraram receita
top_p_amount= df_clean.groupby('producer')['total_amount_usd'].sum().sort_values(ascending=False)
top_p_sales = df_clean.groupby('producer')['sales'].sum().sort_values(ascending=False)
top_p_amount_merged= pd.merge(top_p_amount, top_p_sales, on='producer', how='left')
top_p_amount_merged['rank'] = range(1, (len(top_p_amount)+1))
top_p_amount_merged[:20]
# top 20 produtores mais vendidos
top_p_sales = df_clean.groupby('producer')['sales'].sum().sort_values(ascending=False)
top_p_amount= df_clean.groupby('producer')['total_amount_usd'].sum().sort_values(ascending=False)
top_p_sales_merged = pd.merge(top_p_sales, top_p_amount, on='producer', how='left')
top_p_sales_merged['rank'] = range(1, (len(top_p_amount)+1))
top_p_sales_merged[:20]
# plotar gráfico de dispersão
# os 20 produtores com maiores vendas x maiores receita
# configurar gráfico
top20_p_amount = top_p_amount_merged.iloc[0:20]
top20_p_sales = top_p_sales_merged.iloc[0:20]
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(x=top20_p_amount['total_amount_usd'].values, y=top20_p_amount['sales'].values, c='red', s=70, alpha=0.4)
ax.scatter(x=top20_p_sales['total_amount_usd'].values, y=top20_p_sales['sales'].values, c='blue', s=20, alpha=0.4)
ax.set_title('Top 20 produtores mais vendidos')
ax.set_xlabel('Receita')
ax.set_ylabel('Unidades vendidas')
ax.legend(['Receita', 'Unidades Vendidas'], loc='upper center')
plt.show()
Mais uma vez, podemos observar dois pontos destoantes dos demais, na extremidade direita e no alto do gráfico. Vamos prosseguir com a confirmação.
# imprimir top 1 produtor com mais unidades vendidas e o top 1 que gerou maior receita
print('O top 1 produtor com mais unidades vendidas: {}\n'
'O top 1 produtor que gerou mais receita: {}'.format(top_p_sales_merged.index[0], top_p_amount_merged.index[0]))
'''
O top 1 produtor com mais unidades vendidas: Domaine Ponsot
O top 1 produtor que gerou mais receita: Domaine Ponsot
'''E, novamente, dá ele: o Domaine Ponsot. Se observamos novamente as listas de ranking acima, de fato, a distância dele para o segundo colocado são muito grandes.
Além disso, destaca-se que o Chateau Latour, que está na segunda posição como produtor que mais gerou receita no período, entretanto, não aparece na lista dos top 20 com maior quantidade de vinhos vendidos. E, se voltarmos na análise da seção anterior, ele está em 5ª colocação entre os produtores com maior mediana de preço (500 usd/garrafa), mas também não figura entre os 20 com maior quantidade de produtos oferecidos. Logo, podemos inferir que, apesar de poucos produtos oferecidos e poucas unidades vendidas, sua segunda posição com maior receita pode se dar devido ao fato do alto valor de seus vinhos.
10.4. Região: quantidade de produto x receita
Podemos comparar ainda a quantidade de produto oferecido e a receita gerada, por região de produção dos vinhos.
Novamente, vou gerar 2 conjuntos de dados:
- as 20 regiões que geraram mais receita
- as 20 regiões maior quantidade de itens oferecidos
Por fim, vou construir o gráfico de dispersão que mostre a relação entre elas.
# top 20 regiões que mais geraram receita
top_r_amount= df_clean.groupby('region')['total_amount_usd'].sum().sort_values(ascending=False)
top_r_sales = df_clean.groupby('region')['item_id'].nunique().sort_values(ascending=False)
top_r_amount_merged= pd.merge(top_r_amount, top_r_sales, on='region', how='left')
top_r_amount_merged['rank'] = range(1, (len(top_r_amount)+1))
top_r_amount_merged[:20]
# top 20 produtores mais vendidos
top_r_sales = df_clean.groupby('region')['item_id'].nunique().sort_values(ascending=False)
top_r_amount= df_clean.groupby('region')['total_amount_usd'].sum().sort_values(ascending=False)
top_r_sales_merged = pd.merge(top_r_sales, top_r_amount, on='region', how='left')
top_r_sales_merged['rank'] = range(1, (len(top_r_amount)+1))
top_r_sales_merged[:20]
# plotar gráfico de dispersão
# as 20 regiões com maiores vendas x maiores receita
# configurar gráfico
fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(x=top_r_amount_merged['total_amount_usd'].values, y=top_r_amount_merged['item_id'].values, c='red', s=70, alpha=0.4)
ax.scatter(x=top_r_sales_merged['total_amount_usd'].values, y=top_r_sales_merged['item_id'].values, c='blue', s=20, alpha=0.4)
ax.set_title('Top 20 regiões mais vendidas')
ax.set_xlabel('Receita')
ax.set_ylabel('Unidades vendidas')
ax.legend(['Receita', 'Unidades Vendidas'], loc='upper center')
plt.show()
Agora, vemos 4 pontos que se destacam, 2 a 2 sobrepostos e, ao checar os rankings acima, confirmamos que os 2 primeiros lugares das regiões que mais geram receitas, também são as 2 regiões que mais possuem produtos diferentes disponíveis. Vale lembrar que a correlação de itens com região já havia sido apontada durante a análise de dados.
Outro ponto interessante é que Ribera del Duero, com apenas 5 produtos, fica em 3º lugar em montante de receita.
11. Previsão de Demanda
Para essa seção serão necessárias as seguintes etapas:
- preparação dos dados para atender os requisitos das bibliotecas
- checagem da série (estacionária ou não-estacionária)
- se não-estacionária, aplicar as transformações necessárias - dividir o conjunto em treino e validação
- fazer previsões com diferentes algoritmos
11.1. Preparação dos dados
Na maioria das bibliotecas, como o Prophet e o PyCaret, o conjunto de dados pode conter apenas as informações de data e os valores, que no nosso caso são as vendas. Ainda, para o Prophet é necessário que os atributos sejam nomeados de forma específica, então, já vou deixar nesse formato também.
# criar variável 'df_ts' para receber os dados formatados
# inserir os dados de venda, por dia (por isso, soma)
df_ts = df_clean.groupby('date',as_index=False)['sales'].sum()
# verificar
df_ts.head()
# renomear atributos
df_ts.columns = ['ds', 'y']
# definir data como index
df_ts.index = pd.to_datetime(df_ts['ds'], format="%Y-%m-%d")
# remover o atributo data que ficou repetido
df_ts.drop('ds', axis=1, inplace=True)
# verificar alterações
df_ts.head()
11.2. Série Estacionária
Para que possamos usar nossos dados para realizar previsões é requisito de alguns algoritmos que estejamos trabalhando com uma série estacionária, o que significa que o conjunto de dados deve possuir características estatísticas constantes em relação ao tempo. Ou, em outras palavras, deve ter a média, variância e covariância permanecendo constante no intervalo de tempo (Imagem 1). Isso se deve ao fato de que a maioria dos métodos estatísticos assume a premissa de se tratar de uma série estacionária para a realização dos cálculos.
Séries Temporais são dados distribuídos em intervalos sequenciais e regulares ao longo do tempo.
Por ter a dependência da variável tempo as técnicas tradicionais não são adequadas para esse tipo de problema. E, não considerar a dependência em relação ao tempo pode mudar o significado dos dados e levar a conclusões equivocadas!
11.2.1. Teste Augmented Dickey Fuller (ADF)
Para checar se um determinado conjunto é uma série estacionária existem diversos métodos que podemos aplicar. Contudo, os testes estatísticos são os mais confiáveis. Entre eles, o mais conhecido é o teste Augmented Dickey Fuller (ADF), que para o teste de hipóteses vamos considerar aqui:
- Hipótese Nula (H0) ▶ p-value> 0.05
A série não é estacionária, logo, os dados possuem alguma dependência em relação ao tempo. - Hipótese Alternativa (H1) ▶ p-value<= 0.05
Rejeita a hipótese nula, ou seja, a série é estacionária.
Nota
O nível de significância (p-value) é um valor determinado de acordo com cada caso para determinar um limite para aceitar ou rejeitar a hipótese nula. Assim, quanto menor for o valor dado por p-value, mais forte são as evidências contra a hipótese nula. Por exemplo, se o resultado obtido em p-value fosse de 0.04, logo, conseguiríamos rejeitar a hipótese nula e defender que a série é estacionária com um nível de confiança de 96%.
Portanto, se o valor ficasse acima do threshold, não conseguimos rejeitar a hipótese nula, mas isso não é uma verdade absoluta. Significa apenas que não encontramos evidências suficientes para rejeitá-la.
# extrair apenas os valores
X = df_ts.y
# aplicar ADF
result = adfuller(X)
# imprimir resultados
print('Augumented Dickey-Fuller')
print('Statistical Test: {:.4f}'.format(result[0]))
print('P-Value: {:.10f}'.format(result[1]))
print('Critical Values: ')
for key, value in result[4].items():
print('\t{}: {:.4f}'.format(key, value))
'''
Augumented Dickey-Fuller
Statistical Test: -2.3601
P-Value: 0.1533303754
Critical Values:
1%: -3.4365
5%: -2.8642
10%: -2.5682
'''O valor calculado para o p-value desse conjunto é de 0.1533, logo, acima threshold de 0.05 que foi determinado. Portanto, precisamos transformar essa série em estacionária.
11.2.2. Transformar Série
O processo de transformação da série em estacionária requer a remoção da tendência e da sazonalidade nos dados. Há diversas formas de se fazer isso, nesse caso, primeiro vou aplicar log para reduzir a magnitude dos valores da série e, então, vou subtrair a média móvel de um determinado período em relação ao log da série.
Vou checar a série, novamente aplicando o teste ADF e, em seguida aplicar também o método de diferenciação. Com uma nova verificação do teste ADF.
Para iniciar, primeiro vou plotar os dados brutos para termos uma noção visual do seu estado original. Como estamos lidando com dados diários, vou usar uma média móvel de 30, que corresponderia a um mês.
# plotar série temporal original
## definir média móvel de 30 períodos
ma = df_ts.rolling(30).mean()
## configurar gráfico
fig, ax = plt.subplots(figsize=(10,5))
df_ts.plot(ax=ax, c='#004a8f', legend=False)
ma.plot(ax=ax, c='#E3F700', legend=False)
plt.tight_layout();
Agora, vou aplicar o log, a fim de reduzir a magnitude, ou seja, a assimetria dos nossos dados. Na sequência, vou subtrair a média móvel do período de 30 dias em relação ao log da série temporal.
# aplicar log
df_log = np.log(df_ts)
# calcular log da média móvel de 30
ma_log = df_log.rolling(30).mean()
# plotar gráfico
fig, ax = plt.subplots(figsize=(10,5))
df_log.plot(ax=ax, c='#004a8f', legend=False)
ma_log.plot(ax=ax, c='#E3F700', legend=False)
plt.tight_layout();
Note que o gráfico anterior, sem a aplicação do log, o eixo y ia até 52 mil. Agora, podemos ver acima que o log reduziu o eixo y consideravelmente: para 10.85.
# subtrair a média do log dos dados
df_sub = (df_log - ma_log).dropna()
# obs: é necessário o dropna porque os 30 primeiros valores estarão em branco
# e esses dados em branco podem prejudicar a transformação
ma_sub = df_sub.rolling(30).mean()
# desvio padrão
std_sub = df_sub.rolling(30).std()
# plotar gráfico
fig, ax = plt.subplots(figsize=(10,5))
df_sub.plot(ax=ax, c='#004a8f', legend=False)
ma_sub.plot(ax=ax, c='#E3F700', legend=False)
std_sub.plot(ax=ax, c='#F57169', legend=False)
plt.tight_layout();
Repetiremos o teste ADF para verificar se agora temos uma série estacionária.
# extrair apenas os valores
X_sub = df_sub.y
# aplicar ADF
result_sub = adfuller(X_sub)
# imprimir resultados
print('Augumented Dickey-Fuller')
print('Statistical Test: {:.4f}'.format(result_sub[0]))
print('P-Value: {:.10f}'.format(result_sub[1]))
print('Critical Values: ')
for key, value in result_sub[4].items():
print('\t{}: {:.4f}'.format(key, value));
'''
Augumented Dickey-Fuller
Statistical Test: -4.8494
P-Value: 0.0000437098
Critical Values:
1%: -3.4366
5%: -2.8643
10%: -2.5682
'''Com o resultado de p-value de 0.0000437, agora sim, já conseguiríamos descartar a hipótese nula e considerar a série como estacionária. Lembrando que o valor anterior obtido de p-value foi de 0.1533, que ficava acima do threshold determinado de 0.05.
Agora, já poderíamos passar para a previsão de demanda, contudo, vou aplicar também a diferenciação para ter um nível de confiança ainda maior nesses dados.
A diferenciação é outra técnica que ajuda a tornar uma série temporal em estacionária. Nela, calculamos a diferença entre duas observações possíveis, logo, esse cálculo é dado por:
Ou seja, dada uma observação, pega-se a observação imediatamente anterior a ela e subtrai.
# aplicar diferenciação
# diff = 1
# porque é de 1 período que queremos calcular a diferenciação
df_diff = df_sub.diff(1)
ma_diff = df_diff.rolling(30).mean()
std_diff = df_diff.rolling(30).std()
fig, ax = plt.subplots(figsize=(10,5))
df_diff.plot(ax=ax, c='#004a8f', legend=False)
ma_diff.plot(ax=ax, c='#E3F700', legend=False)
std_diff.plot(ax=ax, c='#F57169', legend=False)
plt.tight_layout();
# extrair apenas os valores
X_diff = df_diff.y.dropna().values
# aplicar ADF
result_diff = adfuller(X_diff)
# imprimir resultados
print('Augumented Dickey-Fuller')
print('Statistical Test: {:.4f}'.format(result_diff[0]))
print('P-Value: {:.20f}'.format(result_diff[1]))
print('Critical Values: ')
for key, value in result_diff[4].items():
print('\t{}: {:.4f}'.format(key, value));
'''
Augumented Dickey-Fuller
Statistical Test: -10.9003
P-Value: 0.00000000000000000012
Critical Values:
1%: -3.4366
5%: -2.8643
10%: -2.5682
'''Note como conseguimos reduzir o valor de p-value a ponto de termos que aumentar de 10 para 20 casas decimais para termos dimensão do seu valor.
Abaixo, vou plotar os gráficos construídos acima, para checar visualmente a transformação dos dados em série temporal ao longo de cada etapa. Bem como, o valor obtido no teste ADF em cada uma delas.
# configurar gráfico
fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(20,5))
# dados originais
ax0.set_title('Dados Originais\n ADF {:.4}'.format(result[1]))
ax0.plot(df_ts, c='#004a8f')
ax0.plot(ma, c='#E3F700')
# aplicação de log
ax1.set_title('Técnica Log\n ADF {:.4}'.format(result_sub[1]))
ax1.plot(df_sub, c='#004a8f')
ax1.plot(ma_sub, c='#E3F700')
ax1.plot(std_sub, c='#F57169')
# aplicação de diferenciação
ax2.set_title('Técnica Diferenciação\n ADF {:.4}'.format(result_diff[1]))
ax2.plot(df_diff, c='#004a8f')
ax2.plot(ma_diff, c='#E3F700')
ax2.plot(std_diff, c='#F57169')
plt.tight_layout();
11.3. Dividir o conjunto em dados de treino e de validação
Buscamos criar um modelo genérico, isto é, que seja capaz de atender dados do mundo real da melhor maneira possível (sem under-fitting e nem over-fitting). Por isso realizamos a divisão do nosso conjunto de dos em 2 subconjuntos, um de treino e outro de validação. Dessa forma, os dados de treino serão usados para o modelo aprender e, os dados de validação serão utilizados para avaliar o desempenho desse modelo. Além disso, ela deve ocorrer de forma aleatória para não termos uma amostra enviesada.
Não existe regra para o tamanho da divisão, contudo, temos que nos atentar para o fato de que, quanto maior o tempo que estamos tentando prever, maiores as chances do modelo errar. Além disso, o tempo a ser previsto irá depender da utilização desses dados dentro da empresa. Por isso, considerei uma previsão de 120 dias para ser útil em uma projeção trimestral da empresa.
Primeiro, vamos verificar nossa série temporal após as transformações.
# checar série temporal
df_diff.info()
'''
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 1067 entries, 2018-01-30 to 2020-12-31
Data columns (total 1 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 y 1066 non-null float64
dtypes: float64(1)
memory usage: 16.7 KB
'''Nota-se que temos 1067 entradas e 1066 delas estão preenchidas. Logo, temos um dado nulo. Vou removê-lo para não interferir nas previsões.
# remover dado nulo
df_diff.dropna(inplace=True)
# verificar série temporal
df_diff.info()
'''
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 1066 entries, 2018-01-31 to 2020-12-31
Data columns (total 1 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 y 1066 non-null float64
dtypes: float64(1)
memory usage: 16.7 KB
'''Se temos 1066 pontos de dados, logo, nosso conjunto de treino conterão os dados de 0 a 945, enquanto o conjunto de validação ficará com os últimos 30 dias, isto é, dos registros 946 ao 1066.
# dividir conjunto em treino='train' e validação='valid'
train = df_diff[:946]
valid = df_diff[946:]# verificação dos dados de treino
print('Dados de treino\n'
'\t início: {}\n'
'\t fim: {}'.format((train.index[0]), (train.index[-1])))
# verificação dos dados de validação
print('\nDados de Validação\n'
'\t início: {}\n'
'\t fim: {}'.format((valid.index[0]), (valid.index[-1])))
'''
Dados de treino
início: 2018-01-31 00:00:00
fim: 2020-09-02 00:00:00
Dados de Validação
início: 2020-09-03 00:00:00
fim: 2020-12-31 00:00:00
'''11.4. Fazer previsões
A título de comparação serão realizadas diferentes previsões com as técnicas listadas abaixo para que possamos verificar o desempenho das mesmas.
- abordagem Naive¹
- Média Móvel¹
- Holt’s Linear Trend Model¹
- ARIMA²
- Prophet²
- PyCaret¹
[¹] nessas técnicas serão utilizados os dados limpos e processados. Ou seja, sem que a série tenha sido transformada em estacionária: ds_ts, train_ts e valid_ts.
[²] nessas técnicas serão utilizados os dados tratados, ou seja, com a série convertida em estacionária: df_diff, train e valid.
Vamos relembrar a formatação do conjunto de dados limpo e processo, denominado df_ts, e prepará-lo, ou seja, dividí-lo em dados de treino e dados de validação.
# imprimir primeiras 5 entradas do conjunto 'df_ts'
df_ts.head()
# checar série temporal
df_ts.info()
'''
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 1096 entries, 2018-01-01 to 2020-12-31
Data columns (total 1 columns):
# Column Non-Null Count Dtype
--- ------ -------------- -----
0 y 1096 non-null int64
dtypes: int64(1)
memory usage: 17.1 KB
'''# dividir conjunto em treino='train_ts' e validação='valid_ts'
train_ts = df_ts[:976]
valid_ts = df_ts[976:]
# copiar dados de validação
y_hat_ts = valid_ts.copy()# verificação dos dados de treino
print('Dados de treino\n'
'\t início: {}\n'
'\t fim: {}'.format((train_ts.index[0]), (train_ts.index[-1])))
# verificação dos dados de validação
print('\nDados de Validação\n'
'\t início: {}\n'
'\t fim: {}'.format((valid_ts.index[0]), (valid_ts.index[-1])))
'''
Dados de treino
início: 2018-01-01 00:00:00
fim: 2020-09-02 00:00:00
Dados de Validação
início: 2020-09-03 00:00:00
fim: 2020-12-31 00:00:00
'''11.4.1. Abordagem Naive
Com esse método a previsão é feita com base no último dado disponível. Ele serve como uma linha de base inicial, ou seja, um valor de referência para melhorias futuras. É uma inferência simples para resolver o problema e que sabemos que a partir dele podemos melhorar.
# copiar o último valor no conjunto 'train' e atribuir à 'y_hat' em um novo atributo 'naive'
y_hat_ts['naive'] = train_ts.iloc[-1].values[0]
# checar alteração
y_hat_ts.head()
# plotar gráfico
fig, ax = plt.subplots(figsize=(8,4))
train_ts.plot(ax=ax, c='#004a8f')
valid_ts.plot(ax=ax, c='#91aac9')
# criar linha de previsão com y_hat['naive']
y_hat_ts['naive'].plot(ax=ax, c='#F422D1', linewidth=1.8)
# legendas
ax.legend(['Train', 'Valid', 'Naive'])
plt.show()
11.4.2. Média Móvel
Se na abordagem Naive fazemos a previsão com base apenas no último valor, na média móvel utilizamos um intervalo de tempo personalizado. Com isso, as médias móveis suavizam curvas e reduzem o ruído (dispersão).
Vamos fazer a previsão para a nossa série temporal, com base em uma média móvel de 3 dias, dessa forma, ao invés de analisarmos 1 único ponto como ocorre com o Naive, teremos mais 2 informações para agregar.
# calcular a média móvel dos últimos 3 valores disponível no conjunto
# salvar em um novo atributo denominado 'mm3'
y_hat_ts['mm3'] = train_ts.y.rolling(3).mean().iloc[-1]
# checar alteração
y_hat_ts.head()
# plotar gráfico
fig, ax = plt.subplots(figsize=(8,4))
train_ts.plot(ax=ax, c='#004a8f')
valid_ts.plot(ax=ax, c='#91aac9')
# criar linha de previsão com y_hat['mm3']
y_hat_ts['mm3'].plot(ax=ax, c='#F422D1', linewidth=1.8)
# legendas
ax.legend(['Train', 'Valid', 'Média Móvel 3'])
plt.show()
11.4.3. Holt’s Linear Trend Model
Observe que tanto a abordagem Naive quanto a Média Móvel são linhas retas. Isso acontece porque perde-se o senso de tendência. Nesse sentido, o método de previsão utilizando o Holt’s consegue trazer esse parâmetro de volta porque ele trabalha não só com o nível, como as técnicas anteriores, mas com a tendência das séries, o que o leva a ter resultados melhores em relação às anteriores.
# salvar os componentes da série temporal em uma variável
result = seasonal_decompose(train_ts)
# plotar componentes
fig, (ax0, ax1, ax2, ax3) = plt.subplots(4, 1, figsize=(8,10))
result.observed.plot(ax=ax0, c='#91aac9')
result.trend.plot(ax=ax1, c='#004a8f')
result.seasonal.plot(ax=ax2, c='#91aac9')
result.resid.plot(ax=ax3, c='#91aac9')
plt.tight_layout()
Podemos observar acima, em azul mais escuro, o gráfico que mostra a tendência dessa série temporal.
# gerar os valores Holt para essa série (dados de treino = 'train')
# salvar eles em uma nova variável 'holt'
# smoothing_level= ajuste de altura
# smoothing_slope= ajuste de inclinação
# fazer previsão de acordo com o tamanho do conjunto de validação
y_hat_ts['holt'] = Holt(train_ts.y).fit(smoothing_level=0.005,
smoothing_slope=0.01).forecast(len(valid_ts))
# checar alteração
y_hat_ts.head()
# plotar gráfico
fig, ax = plt.subplots(figsize=(8,4))
train_ts.plot(ax=ax, c='#004a8f')
valid_ts.plot(ax=ax, c='#91aac9')
# criar linha de previsão com y_hat['holt']
y_hat_ts['holt'].plot(ax=ax, c='#F422D1', linewidth=1.8)
# legendas
ax.legend(['Train', 'Valid', 'Holt'])
plt.show()
11.4.4. ARIMA
Esse é um dos modelos mais usados para fazer previsões de séries temporais, devido ao fato de que ele consegue capturar conjuntos diferentes de estruturas temporais da série.
Seu acrônimo ARIMA é o conjunto de três valores necessários para criar o modelo:
- Auto Regression (p)
- Integrated (d)
- Moving Average (q)
Vale lembrar que para o ARIMA é necessário que a série seja estacionária, para isso usaremos aqui os dados de treino train e os dados de validação valid, provenientes da divisão do conjunto df_diff.
# instanciar modelo
model = auto_arima(train, # conjunto de dados de treino
suppress_warnings=True) # suprimir alertas
# treinar o modelo
model.fit(train)
# realizar previsões
forecast = model.predict(n_periods=len(valid))
forecast = pd.DataFrame(forecast,index = valid.index,columns=['Prediction'])# plotar gráfico
fig, ax = plt.subplots(figsize=(8,4))
train.plot(ax=ax, c='#004a8f')
valid.plot(ax=ax, c='#91aac9')
# criar linha de previsão com 'forecast'
forecast.plot(ax=ax, c='#F422D1', linewidth=1.8)
# legendas
ax.legend(['Train', 'Valid', 'ARIMA'])
plt.show()
Vou plotar um gráfico apenas com os dados que foram previstos pelo modelo versus os dados reais.
# plotar gráfico
# concatenar dados de previsão e de validação
ax = pd.concat([valid, forecast], axis=1).plot()
# definir cores das linhas
ax.lines[0].set_color('#91aac9') # 'valid'
ax.lines[1].set_color('#004a8f') # 'forecast'
# legendas
ax.legend(['Valid', 'ARIMA'])
plt.show();
11.4.5. Prophet
A biblioteca Prophet, desenvolvida pela empresa Meta (antes denominada Facebook), é simples, fácil de usar e, ainda assim, consegue ser uma ferramenta muito robusta.
Vale lembrar que para o Prophet é necessário que a série seja estacionária, para isso usaremos aqui os dados de treino train e os dados de validação valid (resultado da divisão do conjunto df_diff).
# transformar dados no formato do Prophet
train_p = train.copy()
train_p.reset_index(inplace=True)
valid_p = valid.copy()
valid_p.reset_index(inplace=True)# instanciar modelo
b = Prophet()
# treinar o modelo
b.fit(train_p)
# criar o data frame
future = b.make_future_dataframe(periods=len(valid_p))
# fazer previsões com o data frame 'future' (criado no item anterior)
forecast_p = b.predict(future)
# imprimir
forecast_p.head()Na saída do código acima, o Prophet gera várias informações como resultado. A previsão encontra-se no atributo yhat. Em yhat_lower e yhat_upper temos os dados dos limites inferiores e superiores, respectivamente.
Abaixo, vou gerar a visualização das previsões.
# plotar previsões
b.plot(forecast_p).savefig('forecast.png')
Com o Prophet, podemos facilmente plotar os componentes da nossa previsão, separados em 3 gráficos:
- Tendência Geral
- Tendência Semanal
- Tendência Anual
# plotar componentes da previsão
# geral, semanal, anual
b.plot_components(forecast_p).savefig('components.png')
Assim como vimos em ARIMA, vou plotar aqui os dados de validação e os dados previstos, para que possamos compará-los.
# plotar gráfico
fig, ax = plt.subplots()
ax.plot(valid_p.ds, valid_p.y, c='#91aac9')
ax.plot(forecast_p.ds[946:], forecast_p.yhat[946:], c='#004a8f')
# legendas
ax.legend(['Valid', 'Prophet'], loc='lower left')
# definir tamanho da fonte do eixo x
ax.tick_params(axis='x', labelsize=6)
plt.tight_layout();
Apesar da visualização já nos dar algumas informações, vale a pena plotar uma outra visualização em que possamos checar também os limites previstos pelo modelo, dado por yhat_lower e yhat_upper.
# definir função de comparação dos dados
def make_comparison_dataframe(historical, forecast):
return forecast_p.set_index('ds')[['yhat', 'yhat_lower', 'yhat_upper']].join(historical)
# comparar dados reais com o previstos
cmp_df = make_comparison_dataframe(df_diff, forecast_p)# visualizar resultados
fig, ax = plt.subplots(figsize=(10,5))
# dados para plotar
plt.plot(cmp_df['yhat'])
plt.plot(cmp_df['yhat_lower'])
plt.plot(cmp_df['yhat_upper'])
plt.plot(cmp_df['y'])
# configuração
ax.plot(cmp_df['yhat_lower'], c='#A826FF', ls='--', linewidth=0.5, label='yhat_lower')
ax.plot(cmp_df['yhat_upper'], c='#F422D1', ls='--', linewidth=0.5, label='yhat_upper')
ax.plot(cmp_df['y'], c='#bfbfbf', label='y')
ax.plot(cmp_df['yhat'], c='#004a8f', linewidth=0.8, label='yhat')
# legendas
ax.legend()
plt.tight_layout()
No gráfico acima pode-se notar que o modelo está em um ponto médio dos dados, com o limite inferior mais ajustado do que o superior. Se quisermos realizar previsões mais otimistas, deve-se utilizar o yhat_upper como referência. Já o contrário, ou seja, para previsões mais pessimistas, pode-se usar o yhat_lower.
11.4.6. PyCaret
PyCaret é uma biblioteca de uso simplificado que auxilia o processo de construção de modelos de machine-learning por meio de funções automatizadas de etapas comuns ao desenvolvimento desses modelos. Isso faz dela uma ferramenta extremamente útil para rapidamente construir modelos de aprendizado de máquina.
Vale ressaltar que ferramentas de auto machine-learning como essa não substitui o trabalho do cientista de dados. A utilização dela deve ser feita com o intuito de dar suporte e agilidade ao processo da construção de modelos.
Vou utilizar aqui os dados originais limpos e processados. Ou seja, sem que a série tenha sido transformada em estacionária. Dessa forma, evitamos um processo a mais que alguns algoritmos automaticamente já conseguem lidar.
O processo envolve 7 etapas:
- configuração
- criação e comparação de modelos
- criação do melhor modelo encontrado na etapa anterior
- ajuste de hiperparâmetros, de acordo com a métrica de avaliação MAPE
- visualização do modelo
- previsão no conjunto de testes
- finalização do modelo
A partir disso, podemos fazer a previsão final e gerar um gráfico para visualizar o valor previsto pelo modelo e os dados reais.
A primeira etapa consiste em passar os parâmetros necessários para configuração dos modelos.
# 1. configurar o time series
s = setup(train_ts, # conjunto de dados
fh = 120, # quantidade de pontos para serem previstos
session_id = 123) # seed para gerar resultados reproduzíveisAgora, serão criados e treinados todos os modelos disponíveis na biblioteca do PyCaret, e então, avaliados. Os resultados são ranqueados e o melhor modelo de cada métrica é destacado.
# 2. criar e comparar modelos
best = compare_models()
Acima, foram criados e testados 27 modelos diferentes. Sendo que, os que apresentaram os melhores resultados foram o Huber w/ Cond. Deseasonalize & Detrending e o K Neighbors w/ Cond. Deseasonalize & Detrending.
Podemos confirmar o melhor algoritmo criado com o código abaixo:
# 3. criar modelo com o melhor algoritmo detectado
# identificar o melhor modelo
print(best)
'''
BaseCdsDtForecaster(fe_target_rr=[WindowSummarizer(lag_feature={'lag': [7, 6, 5,
4, 3, 2,
1]},
n_jobs=1)],
regressor=HuberRegressor(), sp=7, window_length=7)
'''Portanto, vamos dar continuidade utilizando o Huber w/ Cond. Deseasonalize & Detrending.
# instanciar e criar o modelo com o melhor algoritmo detectado
huber = create_model('huber_cds_dt')
Note que as médias das métricas mostradas acima são as mesmas para o modelo criado na lista de todos os algoritmos. Isso porque o modelo é criado e avaliado da mesma forma que foi feito anteriormente.
# 4. ajustar hiperparâmetros de acordo como valor do MAPE
tuned_huber = tune_model(huber, optimize='MAPE')
A média obtida com o MAPE foi um pouco maior aqui do que o modelo sem ajustes, contudo, observe que em 2 das 3 avaliações o modelo teve um desempenho melhor. Portanto, vou prosseguir com o modelo ajustado.
Vamos visualizar os dados previstos versus os dados reais no gráfico abaixo.
# 5. visualizar modelo
# plotar gráfico com as previsões
plot_model(tuned_huber, plot='forecast')
Antes de finalizar o modelo, fazemos uma verificação na qual realizamos uma previsão, com o conjunto de testes separado pela função setup na primeira etapa de criação do modelo. Dessa forma, checamos se não há divergências nos resultados apresentados.
# 6. fazer previsões no conjunto de testes
holdout_pred = predict_model(tuned_huber)
O modelo performou um pouco melhor do que o resultado obtido na etapa 4, de 0.323 para 0.0318 e, dentro do desvio padrão de 0.0074 .
Na última etapa do processo de criação do modelo, finalizamos ele com a função finalize_model que irá treinar o modelo com o conjunto completo que foi informado na função setup. Ou seja, ele irá incluir o conjunto que havia sido separado para os testes.
# 7. finalizar modelo
final_huber = finalize_model(tuned_huber)Podemos imprimir o modelo final, para checar os parâmetros utilizados no modelo.
# verificar parâmetros
print(final_huber)
'''
ForecastingPipeline(steps=[('forecaster',
TransformedTargetForecaster(steps=[('model',
BaseCdsDtForecaster(fe_target_rr=[WindowSummarizer(lag_feature={'lag': [7,
6,
5,
4,
3,
2,
1]},
n_jobs=1)],
regressor=HuberRegressor(),
sp=7,
window_length=7))]))])
'''Agora sim, podemos proceder à previsão.
# fazer previsão
forecast_pycaret = predict_model(final_huber, fh=120)# verificar as últimas entradas do conjunto para confirmar as datas
forecast_pycaret.tail()
# criar atributo 'ds' com os valores das datas
forecast_pycaret['ds'] = forecast_pycaret.index.to_timestamp()
# renomear atributos
forecast_pycaret.columns = ['y', 'ds']
# definir 'ds' como index
forecast_pycaret.index = forecast_pycaret.ds
# excluir atributo 'ds' que ficou duplicado
forecast_pycaret.drop('ds', axis=1, inplace=True)
# verificar alterações
forecast_pycaret.head()
Com os dados da previsão, vou plotar em um gráfico junto aos dados reais para que possamos compará-los.
# plotar gráfico
# concatenar dados de previsão e de validação
ax = pd.concat([valid_ts, forecast_pycaret], axis=1).plot()
# definir cores das linhas
ax.lines[0].set_color('#91aac9') # 'valid'
ax.lines[1].set_color('#004a8f') # 'forecast'
# legendas
ax.legend(['Valid', 'PyCaret'])
plt.show();
12. Avaliação dos Modelos
A última etapa da construção de modelos é a avaliação de desempenho, ou seja, o quanto um modelo acertou e errou nas suas previsões.
Para isso, vou trazer 2 subseções: a primeira com gráficos, para uma análise rápida dos modelos e, a segunda, com duas métricas de avaliação que medem os erros dos modelos e, assim, termos um valor estatístico e, portanto, mais preciso de qual modelo obteve o melhor desempenho.
12.1. Reverter transformações para série estacionária
Antes de prosseguir com a avaliação, como nos modelos ARIMA e Prophet utilizamos dados com tratamentos de log, soma da média móvel e diferenciação para torná-la uma série estacionária, para comparar os valores de erro com os outros modelos devemos reverter esse tratamento. Do contrário estaremos comparando dimensões diferentes.
Primeiro, vou preparar os dados do Prophet, separando apenas os dados previstos e corrigindo o índice para as datas nos conjuntos de previsões e no de validação.
# criar variável com dados de previsão do Prophet
forecast_prophet = forecast_p[['ds', 'yhat']][946:]
# configurar datas como índice
forecast_prophet.set_index('ds', inplace=True)
valid_p.set_index('ds', inplace=True)Passamos então à reversão da diferenciação, da soma da média móvel e do log. Para que possamos ter as saídas necessárias (colunas nomeadas em ds e y), foi preciso renomear os atributos e, também remover os dados nulos para evitar erros.
# reverter a diferenciação
valid_revert = valid.cumsum()
valid_p_revert = valid_p.cumsum()
forecast_revert = forecast.cumsum()
forecast_p_revert = forecast_prophet.cumsum()
# renomear atributos
forecast_revert.columns = ['y']
forecast_p_revert.columns = ['y']
# reverter subtração da média
valid_revert = valid_revert + ma_log
valid_p_revert = valid_p_revert + ma_log
forecast_revert = forecast_revert + ma_log
forecast_p_revert = forecast_p_revert + ma_log
# reverter log
valid_revert = np.exp(valid_revert)
valid_p_revert = np.exp(valid_p_revert)
forecast_revert = np.exp(forecast_revert)
forecast_p_revert = np.exp(forecast_p_revert)
# remover dados nulos
valid_revert.dropna(inplace=True)
valid_p_revert.dropna(inplace=True)
forecast_revert.dropna(inplace=True)
forecast_p_revert.dropna(inplace=True)
# verificar reversões
valid_revert.head()
12.2. Gráficos
Uma forma simples de avaliar os modelos criados são por meio visual, ou seja, gráficos. Abaixo, plotei os modelos Naive, Média Móvel e Holt’s juntos e com os dados reais para que possamos compará-los.
# configurar gráfico
fig, ax = plt.subplots(figsize=(10,5))
ax.set_title('Previsões x Validação')
ax.plot(train_ts, c='#bfbfbf')
ax.plot(valid_ts, c='#91aac9')
# abordagem Naive
ax.plot(y_hat_ts['naive'], c='#F422D1', linewidth=1.8)
# média móvel 3
ax.plot(y_hat_ts['mm3'], c='#004a8f', linewidth=1.8)
# holt's
ax.plot(y_hat_ts['holt'], c='#A826FF', linewidth=1.8)
# legendas
ax.legend(['Train', 'Valid', 'Naive', 'Média Móvel 3', 'Holt'])
plt.tight_layout();
Podemos ver acima como a média móvel encontra-se um pouco mais centralizada em relação aos dados e como o Holt, por possuir uma inclinação, está corretamente posicionado na direção da tendência que as vendas se seguiram.
Em um segundo gráfico, trago as previsões dos modelos ARIMA, Prophet e PyCaret (bem como dos dados reais), uma vez que como tentam prever ponto a ponto, eles são mais complexos e plota-los junto aos modelos acima traria uma poluição visual muito grande, dificultando o entendimento. Dessa forma, temos um gráfico mais fácil de ser compreendido e mais fácil de comparar modelos semelhantes.
# configurar gráfico
fig, ax = plt.subplots(figsize=(10,5))
ax.set_title('Previsões x Validação')
ax.plot(valid_revert, c='#bfbfbf')
# ARIMA
ax.plot(forecast_revert.y, c='#F422D1', linewidth=1.5)
# Prophet
ax.plot(forecast_p_revert.y, c='#A826FF', linewidth=1.5)
# PyCaret
ax.plot(forecast_pycaret.y, c='#FB8B24', linewidth=1.5)
# legendas
ax.legend(['Validação', 'ARIMA', 'Prophet', 'PyCaret'], loc='lower left')
# definir tamanho da fonte do eixo x
ax.tick_params(axis='x', labelsize=6)
plt.tight_layout();
Olhando o gráfico acima, vemos que inicialmente os modelos são muito semelhantes, até o Prophet divergir próximo do ponto de 01–11–2020. O PyCaret parece não acompanhar tanto a tendência quanto o ARIMA e o Prophet. Por isso, parece que visualmente o modelo do Prophet se ajusta melhor aos dados reais. Mas vamos tirar a prova a seguir.
12.3. Métricas de Avaliação
Para avaliar os modelos criados anteriormente, ou seja, identificar qual deles possui o melhor desempenho, serão utilizadas duas métricas: o Mean Absolute Error (MAE) e o Mean Absolute Percentage Error (MAPE).
O MAE é o valor absoluto do erro na previsão em relação à série real, calculado a partir da média das somas do valor absoluto da magnitude dos erros. Por isso, quanto menor o seu valor, melhor é o modelo. Seu cálculo é dado por:
Já o MAPE mostra o quanto as previsões diferem do valor real em porcentagem, ou seja, é a porcentagem equivalente ao MAE.
# avaliar abordagem Naive
naive_mae = mae(y_hat_ts.y, y_hat_ts.naive)
naive_mape = mape(y_hat_ts.y, y_hat_ts.naive)
# avaliar média móvel 3
mm3_mae = mae(y_hat_ts.y, y_hat_ts.mm3)
mm3_mape = mape(y_hat_ts.y, y_hat_ts.mm3)
# avaliar holt's
holt_mae = mae(y_hat_ts.y, y_hat_ts.holt)
holt_mape = mape(y_hat_ts.y, y_hat_ts.holt)
# avaliar arima
arima_mae = mae(valid_revert, forecast_revert)
arima_mape = mape(valid_revert, forecast_revert)
# avaliar prophet
prophet_mae = mae(valid_p_revert, forecast_p_revert)
prophet_mape = mape(valid_p_revert, forecast_p_revert)
# avaliar pycaret
pycaret_mae = mae(valid_ts, forecast_pycaret)
pycaret_mape = mape(valid_ts, forecast_pycaret)
print(f'{"MÉTRICAS DE AVALIAÇÃO":^40}')
print('-'*40)
print('\t\t MAE \t\t MAPE')
print('Abordagem Naive: {:.2f} \t {:.4f}\n'
'Média Móvel 3: {:.2f} \t {:.4f}\n'
'Holt: {:.2f} \t {:.4f}\n'
'ARIMA: {:.2f} \t {:.4f}\n'
'Prophet: {:.2f} \t {:.4f}\n'
'PyCaret: {:.2f} \t {:.4f}\n'.format(naive_mae, naive_mape,
mm3_mae, mm3_mape,
holt_mae, holt_mape,
arima_mae, arima_mape,
prophet_mae, prophet_mape,
pycaret_mae, pycaret_mape))
'''
MÉTRICAS DE AVALIAÇÃO
----------------------------------------
MAE MAPE
Abordagem Naive: 1631.18 0.0357
Média Móvel 3: 1428.87 0.0309
Holt: 1370.60 0.0290
ARIMA: 1089.11 0.0233
Prophet: 1346.82 0.0287
PyCaret: 1756.81 0.0383
'''Conseguimos verificar que o melhor modelo gerado foi pelo algoritmo ARIMA, com uma taxa de erro de 2.33% apenas. Embora o Prophet não tenha ficado muito longe, com 2.87% de taxa de erro. Há que se pontuar que o tempo de execução: o Prophet é muito mais ágil comparado ao ARIMA e isso pode ser um ponto importante a ser considerado, a depender do uso do algoritmo.
Cabe destacar ainda, o resultado obtido pelo Holt, de 2.90%, que, por ser um algoritmo mais simples e ter ficado com resultado muito próximo ao Prophet, de acordo com a aplicação prática, pode ser uma excelente alternativa.
13. Conclusão
Esse estudo teve por objetivo central o desenvolvimento de um algoritmo capaz de prever a demanda de vinhos para uma loja especialidade nesse produto.
Com o uso da biblioteca Pandas Profiling foi possível obter em poucos segundos um panorama do conjunto de dados, com a geração de um relatório completo, que auxiliou na análise exploratória.
Após o devido tratamento dos dados, bem como a aplicação de feature engineering que desmembrou as informações temporais em outras 7 variáveis, bem como a adição de 2 variáveis relativas às receitas da empresa, foram realizadas diversas análises para compreensão do negócio. Então, prosseguiu-se para a elaboração de algoritmos preditivos da demanda.
Por estarmos lidando com uma série não-estacionária, de acordo com o teste Augumented Dickey Fuller (ADF), para atingir melhores resultados em alguns algoritmos, foi realizada a transformação da série em estacionária com o uso das técnicas de aplicação de log para redução da magnitude dos valores, subtração da média móvel de 30 períodos e a diferenciação.
A previsão teve como parâmetro 120 dias e foram usados os seguintes métodos: abordagem Naive, Média Móvel, Holt’s Linear Trend Model, ARIMA, Prophet e PyCaret. Sendo que nesse último foram gerados 27 modelos preliminares, do qual se selecionou o melhor para prosseguir com os ajustes de hiperparâmetros e previsão dos dados.
Como resultado, de acordo com o MAPE, o modelo gerado pelo ARIMA obteve o melhor desempenho, com uma taxa de erro de apenas 2.33%, seguido pelo Prophet com 2.87% e do Holt com 2.90%.
Para aplicação do modelo há que se considerar o objetivo ao qual ele servirá, uma vez que dessas opções há algoritmos que rodam com mais agilidade do que outros.
Saiba Mais
Esse estudo encontra-se disponível nas plataformas Google Colab e GitHub. Basta clicar nos links das imagens abaixo para ser redirecionado.
