Artigo 41 — Séries Temporais e Previsão com Prophet e statsmodels

Prof. Ricardo Matos Módulo 7 · Data Science e Machine Learning · Artigo 41 de 52

Introdução

Séries temporais estão em toda parte — preços de ações, temperatura diária, vendas mensais, tráfego de um site, consumo de energia elétrica. Diferente de dados tabulares convencionais, séries temporais têm uma dimensão fundamental: o tempo. A ordem das observações importa, e padrões como tendência, sazonalidade e ciclos precisam ser identificados e modelados. Python tem ferramentas excelentes para isso — desde a análise exploratória até modelos de previsão de produção.

Instalação

pip install statsmodels prophet pandas numpy matplotlib plotly

Criando e Manipulando Séries Temporais com pandas

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime


# Criando série temporal sintética
np.random.seed(42)
datas  = pd.date_range(start="2021-01-01", end="2024-12-31", freq="D")
n      = len(datas)

# Componentes da série
tendencia     = np.linspace(100, 200, n)
sazonalidade  = 20 * np.sin(2 * np.pi * np.arange(n) / 365.25)
sazon_semanal = 10 * np.sin(2 * np.pi * np.arange(n) / 7)
ruido         = np.random.normal(0, 5, n)

valores = tendencia + sazonalidade + sazon_semanal + ruido

# Criando DataFrame com DatetimeIndex
df = pd.DataFrame({"valor": valores}, index=datas)
df.index.name = "data"

print(df.head())
print(f"\nPeríodo: {df.index.min()} → {df.index.max()}")
print(f"Observações: {len(df)}")
print(f"\nEstatísticas:\n{df.describe().round(2)}")

Operações com DatetimeIndex

# Resample — agregação por período
mensal   = df.resample("ME").mean()       # média mensal
semanal  = df.resample("W").sum()         # soma semanal
trimestral = df.resample("QE").agg({
    "valor": ["mean", "min", "max", "std"]
})

print("=== Médias Mensais ===")
print(mensal.head(6).round(2))

# Frequências úteis de resample
# D   — diário
# W   — semanal
# ME  — fim de mês
# MS  — início de mês
# QE  — fim de trimestre
# YE  — fim de ano
# H   — horário
# 15min — 15 minutos

# Seleção por período
jan_2024   = df["2024-01"]
primeiro_t = df["2024-01":"2024-03"]
ultimo_ano = df["2024"]

# Atributos de data no índice
df["ano"]       = df.index.year
df["mes"]       = df.index.month
df["dia"]       = df.index.day
df["dia_semana"] = df.index.dayofweek   # 0=segunda, 6=domingo
df["semana"]    = df.index.isocalendar().week.astype(int)
df["trimestre"] = df.index.quarter

# Estatísticas por componente temporal
media_por_mes    = df.groupby("mes")["valor"].mean()
media_por_dia_sem = df.groupby("dia_semana")["valor"].mean()

print("\n=== Média por Mês ===")
nomes_meses = ["Jan","Fev","Mar","Abr","Mai","Jun",
               "Jul","Ago","Set","Out","Nov","Dez"]
for mes, media in media_por_mes.items():
    print(f"  {nomes_meses[mes-1]:4}: {media:.2f}")

Análise Exploratória de Séries Temporais

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec


def plot_serie_completa(df: pd.DataFrame, coluna: str = "valor"):
    fig = plt.figure(figsize=(16, 10))
    gs  = gridspec.GridSpec(3, 2, figure=fig)

    # Série completa
    ax1 = fig.add_subplot(gs[0, :])
    ax1.plot(df.index, df[coluna], linewidth=0.8, color="#2E75B6", alpha=0.8)
    ax1.set_title("Série Temporal Completa")
    ax1.set_ylabel("Valor")
    ax1.grid(alpha=0.3)

    # Média móvel
    mm_30  = df[coluna].rolling(30).mean()
    mm_90  = df[coluna].rolling(90).mean()
    ax1.plot(df.index, mm_30, color="orange", linewidth=1.5, label="MM 30d")
    ax1.plot(df.index, mm_90, color="red",    linewidth=1.5, label="MM 90d")
    ax1.legend()

    # Sazonalidade anual
    ax2 = fig.add_subplot(gs[1, 0])
    media_mensal = df.groupby("mes")[coluna].mean()
    ax2.bar(nomes_meses, media_mensal.values, color="#2E75B6", alpha=0.8)
    ax2.set_title("Padrão Mensal Médio")
    ax2.set_ylabel("Valor Médio")

    # Sazonalidade semanal
    ax3 = fig.add_subplot(gs[1, 1])
    nomes_dias = ["Seg", "Ter", "Qua", "Qui", "Sex", "Sáb", "Dom"]
    media_diaria = df.groupby("dia_semana")[coluna].mean()
    ax3.bar(nomes_dias, media_diaria.values, color="#375623", alpha=0.8)
    ax3.set_title("Padrão Semanal Médio")

    # Distribuição dos valores
    ax4 = fig.add_subplot(gs[2, 0])
    df[coluna].hist(bins=50, ax=ax4, color="#2E75B6", alpha=0.8)
    ax4.set_title("Distribuição dos Valores")

    # Autocorrelação
    ax5 = fig.add_subplot(gs[2, 1])
    pd.plotting.autocorrelation_plot(df[coluna].dropna(), ax=ax5)
    ax5.set_title("Autocorrelação")
    ax5.set_xlim(0, 60)

    plt.suptitle("Análise Exploratória da Série Temporal",
                 fontsize=14, fontweight="bold")
    plt.tight_layout()
    plt.savefig("eda_serie_temporal.png", dpi=120)
    plt.show()


plot_serie_completa(df)

Decomposição da Série Temporal

from statsmodels.tsa.seasonal import seasonal_decompose, STL


def decompor_serie(serie: pd.Series, modelo: str = "additive", periodo: int = 365):
    """
    Decompõe a série em tendência, sazonalidade e resíduo.
    modelo='additive': Y = T + S + R
    modelo='multiplicative': Y = T × S × R
    """
    decomposicao = seasonal_decompose(
        serie.dropna(),
        model=modelo,
        period=periodo,
        extrapolate_trend="freq"
    )

    fig, axes = plt.subplots(4, 1, figsize=(14, 10), sharex=True)

    componentes = [
        (serie,                      "Original",    "#2E75B6"),
        (decomposicao.trend,         "Tendência",   "#C00000"),
        (decomposicao.seasonal,      "Sazonalidade","#375623"),
        (decomposicao.resid,         "Resíduo",     "#888888"),
    ]

    for ax, (comp, titulo, cor) in zip(axes, componentes):
        ax.plot(comp, color=cor, linewidth=0.8)
        ax.set_ylabel(titulo)
        ax.grid(alpha=0.3)

    axes[0].set_title(f"Decomposição {modelo.capitalize()} — Período={periodo}d")
    plt.tight_layout()
    plt.savefig("decomposicao.png", dpi=120)
    plt.show()

    return decomposicao


# Usando dados mensais para decomposição mais clara
serie_mensal = df["valor"].resample("ME").mean()
decomp = decompor_serie(serie_mensal, modelo="additive", periodo=12)


# STL — mais robusto que seasonal_decompose
stl    = STL(serie_mensal, period=12, robust=True)
res_stl = stl.fit()

fig = res_stl.plot()
fig.set_size_inches(14, 8)
plt.savefig("stl_decomposicao.png", dpi=120)
plt.show()

Teste de Estacionariedade

from statsmodels.tsa.stattools import adfuller, kpss


def testar_estacionariedade(serie: pd.Series, nome: str = "Série"):
    """
    ADF — H0: série tem raiz unitária (não estacionária)
    KPSS — H0: série é estacionária
    """
    print(f"\n{'='*50}")
    print(f"Testes de Estacionariedade — {nome}")
    print('='*50)

    # Augmented Dickey-Fuller
    adf_result = adfuller(serie.dropna(), autolag="AIC")
    print(f"\nAugmented Dickey-Fuller:")
    print(f"  Estatística:  {adf_result[0]:.4f}")
    print(f"  p-valor:      {adf_result[1]:.4f}")
    print(f"  Lag:          {adf_result[2]}")
    print(f"  Estacionária: {'Sim' if adf_result[1] < 0.05 else 'Não'} "
          f"(α=0.05)")

    # KPSS
    kpss_result = kpss(serie.dropna(), regression="c")
    print(f"\nKPSS:")
    print(f"  Estatística:  {kpss_result[0]:.4f}")
    print(f"  p-valor:      {kpss_result[1]:.4f}")
    print(f"  Estacionária: {'Sim' if kpss_result[1] > 0.05 else 'Não'} "
          f"(α=0.05)")

    return adf_result[1] < 0.05


# Testando a série original e diferenciada
eh_estacionaria = testar_estacionariedade(serie_mensal, "Série Original")

if not eh_estacionaria:
    serie_diff = serie_mensal.diff().dropna()
    testar_estacionariedade(serie_diff, "Série Diferenciada (d=1)")

ARIMA / SARIMA com statsmodels

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from sklearn.metrics import mean_absolute_error, mean_squared_error
import warnings
warnings.filterwarnings("ignore")


def plot_acf_pacf(serie: pd.Series, lags: int = 40):
    """Plota ACF e PACF para identificar parâmetros p e q."""
    fig, axes = plt.subplots(1, 2, figsize=(14, 4))
    plot_acf(serie.dropna(),  lags=lags, ax=axes[0], title="ACF")
    plot_pacf(serie.dropna(), lags=lags, ax=axes[1], title="PACF")
    plt.tight_layout()
    plt.savefig("acf_pacf.png", dpi=120)
    plt.show()


plot_acf_pacf(serie_mensal.diff().dropna())


def treinar_sarima(serie: pd.Series, ordem=(1,1,1), ordem_sazonal=(1,1,1,12)):
    """
    SARIMA(p,d,q)(P,D,Q,s)
    p,d,q: parâmetros AR, diferenciação, MA
    P,D,Q: parâmetros sazonais
    s:     período sazonal
    """
    # Divisão treino/teste
    n_teste  = 12   # últimos 12 meses para teste
    treino   = serie[:-n_teste]
    teste    = serie[-n_teste:]

    # Treinamento
    modelo = SARIMAX(
        treino,
        order=ordem,
        seasonal_order=ordem_sazonal,
        enforce_stationarity=False,
        enforce_invertibility=False
    )
    resultado = modelo.fit(disp=False)

    print(resultado.summary())

    # Previsão no período de teste
    previsoes = resultado.forecast(steps=n_teste)

    # Métricas
    mae  = mean_absolute_error(teste, previsoes)
    rmse = np.sqrt(mean_squared_error(teste, previsoes))
    mape = np.mean(np.abs((teste - previsoes) / teste)) * 100

    print(f"\nMétricas no conjunto de teste:")
    print(f"  MAE:  {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.2f}%")

    # Visualização
    fig, ax = plt.subplots(figsize=(14, 5))
    treino.plot(ax=ax, label="Treino",   color="#2E75B6")
    teste.plot(ax=ax,  label="Real",     color="#375623")
    previsoes.plot(ax=ax, label="Previsão", color="#C00000",
                   linestyle="--", marker="o")

    # Intervalo de confiança para previsão futura
    prev_futuro = resultado.get_forecast(steps=n_teste + 12)
    ic = prev_futuro.conf_int()
    ax.fill_between(ic.index, ic.iloc[:, 0], ic.iloc[:, 1],
                    alpha=0.2, color="#C00000", label="IC 95%")

    ax.set_title(f"SARIMA{ordem}{ordem_sazonal}")
    ax.legend()
    ax.grid(alpha=0.3)
    plt.tight_layout()
    plt.savefig("sarima_previsao.png", dpi=120)
    plt.show()

    return resultado, previsoes


modelo_sarima, prev_sarima = treinar_sarima(
    serie_mensal,
    ordem=(1, 1, 1),
    ordem_sazonal=(1, 1, 0, 12)
)

Prophet: Previsão Moderna do Meta

Prophet é a biblioteca de previsão do Meta — lida automaticamente com sazonalidade, feriados e dados ausentes:

from prophet import Prophet
from prophet.plot import plot_plotly, plot_components_plotly
import pandas as pd
import numpy as np


# Prophet exige colunas 'ds' (datetime) e 'y' (valor)
df_prophet = serie_mensal.reset_index()
df_prophet.columns = ["ds", "y"]

print(df_prophet.head())
print(f"Shape: {df_prophet.shape}")


def treinar_prophet(
    df: pd.DataFrame,
    periodos: int = 12,
    intervalo_confianca: float = 0.95,
    sazonalidade_anual: bool = True,
    sazonalidade_semanal: bool = False
) -> tuple:
    """Treina modelo Prophet e gera previsão."""

    # Divisão treino/teste
    n_teste = 12
    treino  = df[:-n_teste]
    teste   = df[-n_teste:]

    # Instanciando e configurando
    modelo = Prophet(
        interval_width=intervalo_confianca,
        yearly_seasonality=sazonalidade_anual,
        weekly_seasonality=sazonalidade_semanal,
        daily_seasonality=False,
        seasonality_mode="additive",
        changepoint_prior_scale=0.05,    # flexibilidade da tendência
        seasonality_prior_scale=10.0,    # força da sazonalidade
    )

    # Adicionando feriados brasileiros
    modelo.add_country_holidays(country_name="BR")

    # Sazonalidade mensal customizada
    modelo.add_seasonality(
        name="mensal",
        period=30.5,
        fourier_order=5
    )

    # Treinamento
    modelo.fit(treino)

    # DataFrame futuro
    futuro = modelo.make_future_dataframe(
        periods=periodos + n_teste,
        freq="ME"
    )

    # Previsão
    previsao = modelo.predict(futuro)

    # Métricas no período de teste
    prev_teste = previsao[-n_teste - periodos:-periodos][["ds", "yhat",
                                                           "yhat_lower",
                                                           "yhat_upper"]]
    valores_reais = teste["y"].values

    mae  = mean_absolute_error(valores_reais, prev_teste["yhat"])
    rmse = np.sqrt(mean_squared_error(valores_reais, prev_teste["yhat"]))
    mape = np.mean(np.abs((valores_reais - prev_teste["yhat"]) /
                           valores_reais)) * 100

    print(f"\nMétricas Prophet no teste:")
    print(f"  MAE:  {mae:.2f}")
    print(f"  RMSE: {rmse:.2f}")
    print(f"  MAPE: {mape:.2f}%")

    return modelo, previsao, prev_teste


modelo_prophet, previsao, prev_teste = treinar_prophet(
    df_prophet,
    periodos=24
)

# Visualizações do Prophet
fig1 = modelo_prophet.plot(previsao, figsize=(14, 5))
fig1.suptitle("Previsão Prophet", fontsize=14)
plt.savefig("prophet_previsao.png", dpi=120)
plt.show()

# Componentes da previsão
fig2 = modelo_prophet.plot_components(previsao, figsize=(14, 8))
plt.savefig("prophet_componentes.png", dpi=120)
plt.show()

Comparando Modelos

from sklearn.metrics import mean_absolute_error, mean_squared_error
import pandas as pd


def comparar_modelos(resultados: dict) -> pd.DataFrame:
    """Compara métricas de múltiplos modelos."""
    linhas = []
    for nome, (real, previsto) in resultados.items():
        mae  = mean_absolute_error(real, previsto)
        rmse = np.sqrt(mean_squared_error(real, previsto))
        mape = np.mean(np.abs((real - previsto) / real)) * 100
        linhas.append({
            "Modelo": nome,
            "MAE":    round(mae,  2),
            "RMSE":   round(rmse, 2),
            "MAPE %": round(mape, 2)
        })
    df_comp = pd.DataFrame(linhas).set_index("Modelo")
    return df_comp.sort_values("RMSE")


# Modelo Naive (baseline) — repete o último valor
n_teste   = 12
real_test = serie_mensal[-n_teste:].values
naive_prev = np.repeat(serie_mensal[-n_teste - 1], n_teste)

# Média Móvel simples
mm_prev = serie_mensal.rolling(6).mean()[-n_teste:].values

# Coletando resultados
resultados_modelos = {
    "Naive (baseline)":    (real_test, naive_prev),
    "Média Móvel 6m":      (real_test, mm_prev),
    "SARIMA(1,1,1)(1,1,0,12)": (real_test, prev_sarima.values),
}

# Adicionando Prophet se disponível
prev_prophet_vals = prev_teste["yhat"].values[:n_teste]
resultados_modelos["Prophet"] = (real_test, prev_prophet_vals)

df_comparacao = comparar_modelos(resultados_modelos)
print("=== Comparação de Modelos ===")
print(df_comparacao.to_string())


# Visualização comparativa
fig, ax = plt.subplots(figsize=(14, 6))
datas_teste = serie_mensal[-n_teste:].index

ax.plot(datas_teste, real_test, label="Real", color="black",
        linewidth=2, marker="o")
ax.plot(datas_teste, naive_prev, label="Naive",
        color="gray", linestyle=":")
ax.plot(datas_teste, mm_prev, label="Média Móvel 6m",
        color="orange", linestyle="--")
ax.plot(datas_teste, prev_sarima.values, label="SARIMA",
        color="#2E75B6", linestyle="--", marker="s")
ax.plot(datas_teste, prev_prophet_vals, label="Prophet",
        color="#C00000", linestyle="--", marker="^")

ax.set_title("Comparação de Modelos de Previsão")
ax.legend()
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig("comparacao_modelos.png", dpi=120)
plt.show()

Exemplo Completo: Pipeline de Previsão de Demanda

import pandas as pd
import numpy as np
from prophet import Prophet
from sklearn.metrics import mean_absolute_error
import json
from pathlib import Path
from dataclasses import dataclass
from typing import Optional


@dataclass
class ResultadoPrevisao:
    modelo:        str
    mae:           float
    rmse:          float
    mape:          float
    previsoes:     pd.DataFrame
    modelo_obj:    object


class PipelinePrevisao:
    """Pipeline completo de previsão de séries temporais."""

    def __init__(self, nome: str):
        self.nome    = nome
        self._dados  = None
        self._treino = None
        self._teste  = None

    def carregar(self, df: pd.DataFrame, col_data: str, col_valor: str,
                 freq: str = "ME") -> "PipelinePrevisao":
        serie = df.set_index(col_data)[col_valor]
        serie.index = pd.to_datetime(serie.index)
        self._dados  = serie.resample(freq).mean().dropna()
        print(f"Série carregada: {len(self._dados)} observações "
              f"({self._dados.index.min().date()} → "
              f"{self._dados.index.max().date()})")
        return self

    def dividir(self, n_teste: int = 12) -> "PipelinePrevisao":
        self._treino = self._dados[:-n_teste]
        self._teste  = self._dados[-n_teste:]
        print(f"Treino: {len(self._treino)} obs | Teste: {len(self._teste)} obs")
        return self

    def prever_prophet(self, periodos_futuros: int = 12) -> ResultadoPrevisao:
        df_treino = self._treino.reset_index()
        df_treino.columns = ["ds", "y"]

        modelo = Prophet(
            yearly_seasonality=True,
            weekly_seasonality=False,
            daily_seasonality=False,
            changepoint_prior_scale=0.05,
        )
        modelo.add_country_holidays(country_name="BR")
        modelo.fit(df_treino)

        n_teste = len(self._teste)
        futuro  = modelo.make_future_dataframe(
            periods=n_teste + periodos_futuros, freq="ME"
        )
        previsao = modelo.predict(futuro)

        # Métricas no teste
        prev_teste = previsao.set_index("ds")["yhat"].reindex(
            self._teste.index
        )
        real  = self._teste.values
        pred  = prev_teste.values
        mae   = mean_absolute_error(real, pred)
        rmse  = np.sqrt(np.mean((real - pred) ** 2))
        mape  = np.mean(np.abs((real - pred) / real)) * 100

        print(f"\n[Prophet] MAE={mae:.2f} | RMSE={rmse:.2f} | MAPE={mape:.2f}%")

        return ResultadoPrevisao(
            modelo="Prophet",
            mae=round(mae, 2),
            rmse=round(rmse, 2),
            mape=round(mape, 2),
            previsoes=previsao[["ds", "yhat", "yhat_lower", "yhat_upper"]],
            modelo_obj=modelo
        )

    def plotar_resultado(self, resultado: ResultadoPrevisao):
        fig, ax = plt.subplots(figsize=(14, 6))

        # Dados históricos
        self._treino.plot(ax=ax, label="Treino", color="#2E75B6", linewidth=1)
        self._teste.plot(ax=ax,  label="Real",   color="#375623",
                         linewidth=2, marker="o")

        # Previsão
        prev = resultado.previsoes.set_index("ds")
        prev_periodo = prev.loc[self._teste.index[0]:]

        ax.plot(prev_periodo.index, prev_periodo["yhat"],
                color="#C00000", linestyle="--",
                linewidth=2, marker="^", label=resultado.modelo)

        ax.fill_between(
            prev_periodo.index,
            prev_periodo["yhat_lower"],
            prev_periodo["yhat_upper"],
            alpha=0.2, color="#C00000", label="IC 95%"
        )

        ax.set_title(
            f"Pipeline de Previsão — {self.nome}\n"
            f"MAE={resultado.mae:.2f} | MAPE={resultado.mape:.2f}%"
        )
        ax.legend()
        ax.grid(alpha=0.3)
        plt.tight_layout()
        plt.savefig(f"pipeline_{self.nome.lower()}.png", dpi=120)
        plt.show()

    def exportar(self, resultado: ResultadoPrevisao, pasta: str = "previsoes"):
        Path(pasta).mkdir(exist_ok=True)
        caminho = Path(pasta) / f"{self.nome}_previsao.csv"
        resultado.previsoes.to_csv(caminho, index=False)

        meta = {
            "serie":    self.nome,
            "modelo":   resultado.modelo,
            "mae":      resultado.mae,
            "rmse":     resultado.rmse,
            "mape":     resultado.mape,
        }
        with open(Path(pasta) / f"{self.nome}_meta.json", "w") as f:
            json.dump(meta, f, indent=2)

        print(f"Exportado: {caminho}")


# Usando o pipeline
df_input = df[["valor"]].reset_index()
df_input.columns = ["data", "valor"]

pipeline = PipelinePrevisao("vendas_mensais")
pipeline.carregar(df_input, "data", "valor")
pipeline.dividir(n_teste=12)

resultado = pipeline.prever_prophet(periodos_futuros=12)
pipeline.plotar_resultado(resultado)
pipeline.exportar(resultado)

Resumo

• Séries temporais têm componentes de tendência, sazonalidade, ciclo e resíduo
• resample() agrega séries por período; rolling() e expanding() calculam estatísticas móveis
• Decomposição com seasonal_decompose e STL separa os componentes para análise individual
• ADF e KPSS testam estacionariedade — requisito do ARIMA; use diferenciação para tornar estacionária
• ACF e PACF identificam os parâmetros p e q do modelo ARIMA/SARIMA
• SARIMA incorpora sazonalidade explícita com os parâmetros (P, D, Q, s)
• Prophet lida automaticamente com sazonalidade múltipla, feriados e dados ausentes
• Compare modelos com MAE, RMSE e MAPE — sempre avalie em conjunto de teste real
• Um pipeline robusto separa claramente carga, divisão, treinamento, avaliação e exportação

Referências e Leituras Complementares

• statsmodels — séries temporais — https://www.statsmodels.org/stable/tsa.html
• Prophet — documentação oficial — https://facebook.github.io/prophet/docs/quick_start.html
• pandas — séries temporais — https://pandas.pydata.org/docs/user_guide/timeseries.html
• statsmodels — SARIMAX — https://www.statsmodels.org/stable/statespace.html
• HYNDMAN, Rob J.; ATHANASOPOULOS, George. Forecasting: Principles and Practice. 3. ed. OTexts, 2021. — referência completa e gratuita em https://otexts.com/fpp3/
• NIELSEN, Aileen. Practical Time Series Analysis. O'Reilly Media, 2019. — abordagem prática com Python cobrindo ARIMA, Prophet e deep learning para séries temporais.

Prof. Ricardo Matos — Dominando o Python em 1 Ano · Artigo 41 de 52 Próximo: Artigo 42 — MLOps: versionamento de modelos, MLflow e deploy de ML

you asked

Sim

claude response

Dominando o Python em 1 Ano