<jupyter_start><jupyter_text>2 - Análise Exploratória de Séries Temporais - Faturamento TotalProjeto para a disciplina de **Estatística** (Módulo 4) do Data Science Degree (turma de julho de 2020)  Equipe* <NAME>* <NAME>* <NAME>* <NAME>Link para [projeto no Github](https://github.com/flimao/case-previsao-faturamento)  IntroduçãoEste notebook é uma continuação da análise exploratória inicial.Neste notebook, vamos progredir para a análise exploratória de séries temporais.<jupyter_code># importação de bibliotecas
import datetime as dt
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import json

# importação de bibliotecas de análise
from statsmodels.tsa.seasonal import seasonal_decompose
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import acf, pacf
from statsmodels.tsa.arima_process import ArmaProcess
from pmdarima.arima import auto_arima
from pmdarima.arima.arima import ARIMA

# teste para verificar estacionariedade (Dickey-Fuller: https://en.wikipedia.org/wiki/Dickey%E2%80%93Fuller_test)
from statsmodels.tsa.stattools import adfuller

# metricas 
from sklearn.metrics import mean_absolute_percentage_error as smape, mean_squared_error as smse, mean_absolute_error as smae 

# pacote com funções para análise desse projeto
import os
cwd = os.getcwd()
os.chdir("../")
import py_scripts.plots, py_scripts.transform, py_scripts.metrics
os.chdir(cwd)

import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 120
mpl.rcParams['figure.figsize'] = (10, 4)<jupyter_output><empty_output><jupyter_text>Importação dos dados<jupyter_code>ts_raw = pd.read_csv(r'../data/sim_ts_limpo.csv')
tsd, tswide = py_scripts.transform.pipeline(ts_raw)
fat_total = tswide.sum(axis = 'columns').dropna()
fat_total
fat_total.plot(linestyle = '', marker = 'o')
plt.title('Faturamento (R$ bi)')
plt.show()<jupyter_output><empty_output><jupyter_text>Análise ExploratóriaVamos primeiramente analisar o faturamento total contido na série histórica:<jupyter_code>fat_total = tsd['total']
fat_total.describe()
sns.scatterplot(data = fat_total)
plt.title('Série histórica (últimos 4 anos)')
plt.show()<jupyter_output><empty_output><jupyter_text>Parece haver um salto entre 2014 e 2015 no faturamento total. Esse salto é devido ao lançamento de um outro produto, `transporte`. O faturamento deste novo produto é uma ordem de magnitude menor que o faturamento do produto `alimenticio` (como vimos brevemente no gráfico de dados faltantes e veremos com detalhes mais a frente), mas é o suficiente para que seja notado no faturamento total.<jupyter_code>fig = plt.figure(figsize = (6, 8))
sns.boxplot(y = fat_total)
plt.ylabel('Faturamento total (R$ bi)')
plt.title('Boxplot - Série Histórica completa')
plt.show()
sns.histplot(fat_total)
plt.xlabel('Faturamento total (R$ bi)')
plt.title('Histograma - Série Histórica completa')
plt.show()<jupyter_output><empty_output><jupyter_text>No entanto, medidas descritivas de séries temporais devem ser tomadas em relação ao tempo. Vamos separar essas medidas ano a ano:<jupyter_code>n_anos = 4
anos_recentes = fat_total[fat_total.index >= dt.datetime.now() - dt.timedelta(days = n_anos * 365) + pd.tseries.offsets.YearBegin()]
anos_recentes.describe()
sns.scatterplot(data = anos_recentes)
plt.title(f'Série histórica (últimos {n_anos} anos)')
plt.show()
sns.boxplot(y = fat_total, x = fat_total.index.year)
plt.ylabel('Faturamento total (R$ bi)')
plt.title('Boxplot - Série Histórica completa')
plt.show()

sns.boxplot(y = anos_recentes, x = anos_recentes.index.year)
plt.ylabel('Faturamento total (R$ bi)')
plt.title(f'Boxplot - Série Histórica (últimos {n_anos} anos)')
plt.show()<jupyter_output><empty_output><jupyter_text>Parece haver alguns *outliers* em 2021.No entanto, a série de 2021 está incompleta (vai somente até outubro). Historicamente, há um salto no faturamento em agosto, o que pode estar causando essa deturpação das medidas descritivas. Excluindo o ano de 2021...<jupyter_code>anos_recentes_exc2021 = fat_total[(fat_total.index >= '2016') & (fat_total.index < '2021')]
sns.boxplot(y = anos_recentes_exc2021, x = anos_recentes_exc2021.index.year)
plt.ylabel('Faturamento total (R$ bi)')
plt.title(f'Boxplot - Série Histórica (2016-2020)')
plt.show()<jupyter_output><empty_output><jupyter_text>A pandemia se faz notar nos dados apenas com o aumento ligeiro da mediana em relação à distância entre o Q1 e o Q3.<jupyter_code>sns.histplot(x = anos_recentes, hue = anos_recentes.index.year, multiple = 'dodge', shrink = .8, common_norm = False, palette = sns.color_palette()[:4])
plt.xlabel('Faturamento total (R$ bi)')
plt.title(f'Histograma - Série Histórica (últimos {n_anos} anos)')
plt.show()<jupyter_output><empty_output><jupyter_text>Os histogramas ano a ano estão melhor comportados que o histograma da série histórica completa. Fazendo uma análise mês a mês para cada ano...<jupyter_code>hue = fat_total.index.year

palette = []
for i, year in enumerate(hue.unique()):
    if year not in [2014, 2015]:
        palette += ['lightgray']
    else:
        palette += [sns.color_palette()[i]]

ax = sns.lineplot(
    y = fat_total, x = fat_total.index.month, 
    hue = fat_total.index.year,
    palette = palette
)
ax.set_xlabel('Mês')
ax.set_ylabel('Faturamento total (R$ bi)')
ax.set_title(f"Receita por mes do ano")
plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>... nota-se claramente o salto dado de 2014 para 2015 com a entrada do novo produto.  Estacionariedade Para que a série de faturamentos mensais totais possa ser decomposta, é necessário que ela seja estacionária. Não parece ser, mas vamos testar através do teste estatístico de Dickey-Fuller.A hipótese nula do teste de Dickey-Fuller é que a série é um passeio aleatório (random walk):<jupyter_code>testedf = adfuller(fat_total)
pvalor = testedf[1]
alpha = 0.05

print(f'Valor-p: {pvalor:.3%}', end = '')

if pvalor < alpha:
    print(f' < {alpha:.0%}')
    print('    Série de faturamentos mensais é estacionária. Rejeita-se a hipótese de a série ser um passeio aleatório.')
else:
    print(f' > {alpha:.0%}')
    print('    Série de faturamentos mensais é um passeio aleatório. Não podemos rejeitar a hipótese nula.')<jupyter_output>Valor-p: 98.342% > 5%
    Série de faturamentos mensais é um passeio aleatório. Não podemos rejeitar a hipótese nula.<jupyter_text>Isso é evidenciado pela decomposição da série temporal. Decomposição em séries de Fourier<jupyter_code>decomp_total = seasonal_decompose(fat_total)

# plot      
fig, axs = plt.subplots(nrows = 4, figsize = (10, 8), sharex = True)

sns.lineplot(data = fat_total, ax = axs[0])
axs[0].set_title('Faturamento total')

sns.lineplot(data = decomp_total.trend, ax = axs[1])
axs[1].set_ylabel('Tendência')

sns.lineplot(data = decomp_total.seasonal, ax = axs[2])
axs[2].set_ylabel('Sazonalidade')

resid = (decomp_total.resid - decomp_total.resid.mean())/decomp_total.resid.std()
sns.scatterplot(data = resid, ax = axs[3])
axs[3].set_ylabel('Residual')

fig.suptitle(f"Decomposição temporal: faturamento total")
plt.show()<jupyter_output><empty_output><jupyter_text>Como mostrado anteriormente, esta série temporal não é estacionária, o que podemos ver através dos resíduos padronizados do último quadro (onde há um padrão claro oscilatório).  Modelo autorregressivo - Faturamento total Para analisar e prever essa série temporal, é necessário um modelo mais completo. Utilizaremos aqui um modelo autorregressivo integrado de média móvel com sazonalidade - **SARIMA**.OBS.: o modelo completo chama-se SARIMAX; o `X` adicional permite a modelagem de variáveis exógenas. No entanto, não utilizaremos variáveis exógenas neste caso.<jupyter_code># excluindo o período pré-2015

test_begin = '2020-01-01'
fat_modelo = fat_total['2015-01-01':]

total_train = fat_modelo[:test_begin].iloc[:-1]
total_test = fat_modelo[test_begin:]


train_test_split_idx = int(fat_modelo.shape[0] * 0.8 + 1)
total_train = fat_modelo[:train_test_split_idx]
total_test = fat_modelo[train_test_split_idx:]

total_train.plot(label = 'Treino')
total_test.plot(label = 'Teste')
plt.title('Train test split - Faturamento total')
plt.ylabel('Faturamento total (R$ bi)')
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>O modelo SARIMA contém alguns parâmetros, `S(P, D, Q, S)`, `AR(p)`, `I(d)` e `MA(q)`.Para determinarmos o parâmetro `d`, uma boa indicação é o gráfico de autocorrelação:<jupyter_code>fig = plt.figure()
ax = fig.gca()
plot_pacf(fat_modelo, lags = 20, method = 'ywm', ax = ax)
ax.set_xlabel('Lags')
ax.set_title('Autocorrelação parcial - Série de faturamento total')
plt.show()<jupyter_output><empty_output><jupyter_text>Neste caso, uma boa estimativa para o parâmetro `d` é 1 subtraído do número de *lags* em que a correlação é estatisticamente significativa. Neste caso, $d \sim 1$.<jupyter_code>arimas = {}
arimas['total'] = auto_arima(
    y = total_train,
    start_p = 1, max_p = 3,
    d = 2, max_d = 4,
    start_q = 1, max_q = 3,
    start_P = 1, max_P = 3,
    D = None, max_D = 4,
    start_Q = 1, max_Q = 3,
    #max_order = 6,
    m = 12,
    seasonal = True,
    alpha = 0.05,
    stepwise = True,
    trace = True,
    n_fits = 500,
)

modelo_corrente = ARIMA(order = (0, 2, 0), seasonal_order = (1, 1, 1, 12), with_intercept = True).fit(y = total_train)

modelo_funcional = [
    ARIMA(order = (0, 1, 0), seasonal_order = (0, 1, 0, 12), with_intercept = False).fit(y = total_train),
    ARIMA(order = (0, 2, 3), seasonal_order = (2, 1, 1, 12), with_intercept = False).fit(y = total_train),
    ARIMA(order = (0, 2, 0), seasonal_order = (1, 1, 1, 12), with_intercept = False).fit(y = total_train),
    ARIMA(order = (0, 2, 0), seasonal_order = (1, 1, 1, 12), with_intercept = True).fit(y = total_train),
]
arimas['total'] = ARIMA(order = (0, 2, 0), seasonal_order = (1, 1, 1, 12), with_intercept = True).fit(y = total_train)
arimas['total'].summary()<jupyter_output><empty_output><jupyter_text>Métricas para o modelo autorregressivo SARIMAX Primeiramente, podemos avaliar o ajuste visualmente:<jupyter_code>n_test_periods = total_test.shape[0]
arr_preds = arimas['total'].predict(n_test_periods)

idx = pd.date_range(freq = 'MS', start = total_test.index[0], periods = n_test_periods)
preds = pd.Series(arr_preds, index = idx)
preds.name = 'yearly_preds'

preds.plot(label = 'Predição')
total_test.plot(label = 'Conjunto de teste')

plt.legend()
plt.ylabel('Faturamento total (R$ bi)')
plt.title('Predição contra conjunto de teste')
plt.show()<jupyter_output><empty_output><jupyter_text>Vamos aplicar algumas métricas quantitativas ao modelo:<jupyter_code>kwargs_total = dict(
    y_true = total_test,
    y_pred = preds,
    n = total_train.shape[0],
    dof = arimas['total'].df_model()
)

py_scripts.metrics.mostrar_metricas(**kwargs_total)
arimas['total'].arima_res_.data.endog
fat_total['2019-06':'2020-06']<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11889965271431183, 'fraction_numerical': 0.02092853226101261, 'mean_word_length': 5.544856459330144, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 68, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '668448', 'license': 'MIT', 'n_tokens_mistral': 4446, 'n_tokens_neox': 4163, 'n_words': 1049, 'path': 'notebooks_exploration/2-faturamento_total.ipynb', 'repo_name': 'flimao/case-previsao-faturamento', 'size': '11016'}