SARIMA no comércio varejista brasileiro (PMC)

Índice de volume de vendas do varejo ampliado por UF, IBGE, base 2022 = 100, sem ajuste sazonal

Autor

Gustavo Gobetti Bertocco

Data de Publicação

25 de outubro de 2026

Sobre o relatório

A base é a PMC (Pesquisa Mensal do Comércio) do IBGE: índice de volume de vendas do comércio varejista ampliado (inclui veículos e materiais de construção), por Unidade da Federação, mensal desde janeiro de 2004, base 2022 = 100, sem ajuste sazonal.

Etapas da análise: olhar os dados, testar raiz unitária, identificar, ajustar o AutoARIMA(season_length=12), diagnosticar os resíduos e prever 24 meses à frente com intervalos de 80% e 95%.

A UF é um parâmetro. Para gerar o relatório de outra UF:

uv run quarto render template-exercicio-sarima.qmd -P UF:"Bahia"
Código 
UF = "Minas Gerais"

Setup

Código 
from datetime import date
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import pandera.polars as pa
import polars as pl
import seaborn as sns
from IPython.display import Markdown, display
from scipy import stats
from statsforecast import StatsForecast
from statsforecast.arima import ndiffs, nsdiffs
from statsforecast.models import ARIMA, AutoARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.stats.diagnostic import acorr_ljungbox
from statsmodels.tsa.seasonal import STL
from statsmodels.tsa.stattools import adfuller, kpss

SEASON = 12  # dados mensais
HORIZON = 24  # meses à frente
LEVELS = [80, 95]
# o Python formata meses em inglês; nomes explícitos evitam depender do locale
MONTHS_PT = [
    "jan", "fev", "mar", "abr", "mai", "jun",
    "jul", "ago", "set", "out", "nov", "dez",
]  # fmt: skip

# Okabe-Ito: paleta categórica segura para daltonismo
OKABE_ITO = [
    "#E69F00",
    "#56B4E9",
    "#009E73",
    "#F0E442",
    "#0072B2",
    "#D55E00",
    "#CC79A7",
]
BLUE, ORANGE, VERMILION, GREY = "#0072B2", "#E69F00", "#D55E00", "#6B6B6B"
sns.set_theme(style="ticks", palette=OKABE_ITO)
plt.rcParams.update({"axes.spines.top": False, "axes.spines.right": False})

# Tabelas completas, sem as linhas de shape e dtype do polars
pl.Config.set_tbl_rows(-1)
pl.Config.set_tbl_hide_dataframe_shape(True)
pl.Config.set_tbl_hide_column_data_types(True)
pl.Config.set_fmt_str_lengths(80)


def month_pt(d) -> str:
    """Data como 'abr/2020'."""
    return f"{MONTHS_PT[d.month - 1]}/{d.year}"


def correlogram(
    values: np.ndarray, ax: plt.Axes, kind: str = "acf", lags: int = 36
) -> None:
    """ACF ou PACF com banda de significância neutra e lags sazonais marcados."""
    plotter = plot_acf if kind == "acf" else plot_pacf
    extra = {} if kind == "acf" else {"method": "ywm"}
    plotter(
        values,
        lags=lags,
        ax=ax,
        title=kind.upper(),
        color=BLUE,
        vlines_kwargs={"colors": BLUE},
        **extra,
    )
    # o statsmodels pinta a banda com a próxima cor da paleta; cinza evita que
    # ela seja lida como dado
    for band in ax.collections[1:]:
        band.set_facecolor(GREY)
        band.set_alpha(0.15)
    for lag in range(SEASON, lags + 1, SEASON):
        ax.axvline(lag, color=VERMILION, lw=0.6, ls=":")

Carregando a base PMC

O CSV é baixado uma vez e guardado em cache como Parquet em data/raw/, então o relatório é gerado da mesma forma mesmo se o release do GitHub estiver fora do ar. O Parquet também preserva o tipo Date.

Código 
DATA_URL = (
    "https://github.com/padsInsper/202634-fa/releases/download/data/pmc_volume_uf.csv"
)
CACHE = Path("data/raw/pmc_volume_uf.parquet")

if CACHE.exists():
    pmc = pl.read_parquet(CACHE)
else:
    pmc = pl.read_csv(DATA_URL).with_columns(pl.col("data").str.to_date("%Y-%m-%d"))
    CACHE.parent.mkdir(parents=True, exist_ok=True)
    pmc.write_parquet(CACHE)

valid_ufs = pmc["uf"].unique().sort().to_list()
if UF not in valid_ufs:
    raise ValueError(f"UF {UF!r} não está na base. Valores válidos: {valid_ufs}")

# formato longo da Nixtla: unique_id, ds, y
df_sf = (
    pmc.filter(pl.col("uf") == UF)
    .select(
        pl.col("uf").alias("unique_id"),
        pl.col("data").alias("ds"),
        pl.col("volume_idx").cast(pl.Float64).alias("y"),
    )
    .sort("ds")
)
df_sf.head()
unique_id ds y
"Minas Gerais" 2004-01-01 47.04283
"Minas Gerais" 2004-02-01 43.47343
"Minas Gerais" 2004-03-01 49.05633
"Minas Gerais" 2004-04-01 48.69024
"Minas Gerais" 2004-05-01 49.3309

Antes de modelar, uma checagem de qualidade dos dados: um mês faltando ou duplicado deslocaria silenciosamente todos os lags sazonais, então é melhor o relatório falhar aqui.

Código 
class MonthlySeries(pa.DataFrameModel):
    """Uma linha por mês: datas no dia 1, sem lacunas, sem duplicatas, y finito."""

    unique_id: str
    ds: pl.Date = pa.Field(unique=True)
    y: pl.Float64 = pa.Field(nullable=False)

    @pa.check("ds", name="month_start")
    def month_start(cls, data: pa.PolarsData) -> pl.LazyFrame:
        return data.lazyframe.select(pl.col(data.key).dt.day() == 1)

    @pa.check("y", name="finite")
    def finite(cls, data: pa.PolarsData) -> pl.LazyFrame:
        return data.lazyframe.select(pl.col(data.key).is_finite())

    @pa.dataframe_check(name="no_missing_months")
    def no_gaps(cls, data: pa.PolarsData) -> pl.LazyFrame:
        months = pl.col("ds").dt.year() * 12 + pl.col("ds").dt.month()
        return data.lazyframe.select(months.diff().fill_null(1) == 1)


df_sf = MonthlySeries.validate(df_sf)
first, last = df_sf["ds"].min(), df_sf["ds"].max()
Markdown(
    f"**UF: {UF}** · {df_sf.height} meses, de {month_pt(first)} a {month_pt(last)} · "
    "checagens de qualidade aprovadas."
)

UF: Minas Gerais · 266 meses, de jan/2004 a fev/2026 · checagens de qualidade aprovadas.

Análise descritiva

A série

A média móvel centrada de 12 meses remove o ciclo sazonal e deixa a tendência mais visível. Ela serve só para este gráfico: uma janela centrada usa valores futuros, então nunca deve entrar como variável de um modelo.

Código 
plot_df = df_sf.with_columns(
    pl.col("y").rolling_mean(window_size=SEASON, center=True).alias("ma12")
)
fig, ax = plt.subplots()
ax.plot(plot_df["ds"], plot_df["y"], color=BLUE, lw=1, label="Índice de volume")
ax.plot(
    plot_df["ds"], plot_df["ma12"], color=ORANGE, lw=2, label="Média móvel de 12 meses"
)
ax.axhline(100, color=GREY, lw=0.8, ls=":")
ax.set(
    title=f"Volume de vendas do varejo ampliado, {UF}",
    ylabel="Índice (2022 = 100)",
    xlabel=None,
)
ax.legend(frameon=False)
plt.show()

Decomposição (STL)

A STL estima tendência e sazonalidade por regressões locais (Loess). Com robust=True, choques pontuais como abril de 2020 ficam no resíduo em vez de distorcer o padrão sazonal.

Código 
y = df_sf["y"].to_numpy()
dates = df_sf["ds"].to_list()
stl = STL(y, period=SEASON, robust=True).fit()

fig, axes = plt.subplots(4, 1, figsize=(10, 8), sharex=True)
for ax, values, name in zip(
    axes,
    [y, stl.trend, stl.seasonal, stl.resid],
    ["Observado", "Tendência", "Sazonal", "Resíduo"],
    strict=True,
):
    ax.plot(dates, values, color=BLUE, lw=1)
    ax.set_ylabel(name)
axes[3].axhline(0, color=GREY, lw=0.8)
fig.suptitle("Decomposição STL (robusta)")
plt.tight_layout()
plt.show()

# Força da tendência e da sazonalidade (Wang, Smith & Hyndman 2006). O teste
# sazonal do AutoARIMA aplica uma diferença sazonal quando a força sazonal > 0,64.
var_r = np.var(stl.resid)
trend_strength = max(0.0, 1 - var_r / np.var(stl.trend + stl.resid))
seasonal_strength = max(0.0, 1 - var_r / np.var(stl.seasonal + stl.resid))
Markdown(
    f"Força da tendência: **{trend_strength:.2f}** · "
    f"força da sazonalidade: **{seasonal_strength:.2f}** (limiar 0,64)."
)

Força da tendência: 0.98 · força da sazonalidade: 0.82 (limiar 0,64).

Padrão sazonal por mês

O boxplot mostra quais meses ficam sistematicamente acima ou abaixo da tendência, e quão estável é esse efeito ao longo dos anos. Ele usa a série sem tendência (observado menos a tendência da STL): nos níveis brutos, cada caixa mostraria principalmente a tendência (o índice praticamente dobrou no período), escondendo o efeito sazonal.

Código 
by_month = df_sf.with_columns(
    (pl.col("y") - pl.Series(stl.trend)).alias("detrended"),
    pl.col("ds")
    .dt.month()
    .replace_strict(dict(enumerate(MONTHS_PT, start=1)), return_dtype=pl.String)
    .alias("month"),
)
fig, ax = plt.subplots()
sns.boxplot(
    data=by_month.to_pandas(),
    x="month",
    y="detrended",
    order=MONTHS_PT,
    color=OKABE_ITO[1],
    ax=ax,
)
ax.axhline(0, color=GREY, lw=0.8)
ax.set(
    title="Efeito sazonal: índice sem tendência por mês do ano",
    xlabel=None,
    ylabel="Pontos acima / abaixo da tendência",
)
plt.show()

ACF e PACF: série original e diferenças

Primeiro a diferença sazonal, depois a regular. Se a diferença regular vier primeiro, um ciclo sazonal forte pode ser confundido com tendência.

Código 
diffs = df_sf.select(
    pl.col("y"),
    pl.col("y").diff(SEASON).alias("y_d12"),
    pl.col("y").diff(SEASON).diff(1).alias("y_d1_d12"),
)
series = {
    "Original": diffs["y"].to_numpy(),
    "Diferença sazonal (Δ₁₂)": diffs["y_d12"].drop_nulls().to_numpy(),
    "Diferença sazonal + regular (ΔΔ₁₂)": diffs["y_d1_d12"].drop_nulls().to_numpy(),
}

# mathtext nos títulos: a fonte dos gráficos não tem glifos de subscrito Unicode
plot_titles = [
    "Original",
    r"Diferença sazonal ($\Delta_{12}$)",
    r"Diferença sazonal + regular ($\Delta\Delta_{12}$)",
]

fig, axes = plt.subplots(3, 3, figsize=(12, 9))
for row, title, values in zip(axes, plot_titles, series.values(), strict=True):
    row[0].plot(dates[len(dates) - len(values) :], values, color=BLUE, lw=0.9)
    row[0].set_title(title, fontsize=10, loc="left")
    correlogram(values, row[1], "acf")
    correlogram(values, row[2], "pacf")
plt.tight_layout()
plt.show()

As linhas verticais pontilhadas marcam os lags sazonais (12, 24, 36). Leia-os separadamente dos lags curtos 1 a 3: os lags curtos sugerem p e q, e os sazonais sugerem P e Q.

Testes de raiz unitária

Os dois testes têm hipóteses nulas opostas, então são lidos em conjunto:

ADF (H0: raiz unitária) KPSS (H0: estacionária) Conclusão
rejeita não rejeita estacionária
não rejeita rejeita raiz unitária: diferenciar
ambos rejeitam inconclusivo (possivelmente estacionária em tendência)
nenhum rejeita inconclusivo
Código 
def joint_verdict(adf_p: float, kpss_p: float, alpha: float = 0.05) -> str:
    """Leitura conjunta de ADF e KPSS, conforme a tabela acima."""
    adf_rejects, kpss_rejects = adf_p < alpha, kpss_p < alpha
    if adf_rejects and not kpss_rejects:
        return "estacionária"
    if kpss_rejects and not adf_rejects:
        return "raiz unitária: diferenciar"
    return "inconclusivo"


def fmt_kpss_p(p: float) -> str:
    """O statsmodels trunca o p-valor do KPSS no intervalo da tabela [0,01; 0,10]."""
    if p <= 0.01:
        return "≤ 0.01"
    if p >= 0.10:
        return "≥ 0.10"
    return f"{p:.3f}"


rows = []
for name, values in series.items():
    adf_stat, adf_p, *_ = adfuller(values, autolag="AIC")
    kpss_stat, kpss_p, *_ = kpss(values, regression="c", nlags="auto")
    rows.append(
        {
            "série": name,
            "estat. ADF": round(adf_stat, 2),
            "p-valor ADF": round(adf_p, 3),
            "estat. KPSS": round(kpss_stat, 2),
            "p-valor KPSS": fmt_kpss_p(kpss_p),
            "leitura conjunta (5%)": joint_verdict(adf_p, kpss_p),
        }
    )
pl.DataFrame(rows)
série estat. ADF p-valor ADF estat. KPSS p-valor KPSS leitura conjunta (5%)
"Original" -1.9 0.332 1.91 "≤ 0.01" "raiz unitária: diferenciar"
"Diferença sazonal (Δ₁₂)" -2.22 0.2 0.43 "0.064" "inconclusivo"
"Diferença sazonal + regular (ΔΔ₁₂)" -5.9 0.0 0.01 "≥ 0.10" "estacionária"

ndiffs e nsdiffs estimam o número de diferenças automaticamente, com os mesmos testes que o AutoARIMA usa internamente: um teste de força sazonal para D, depois o KPSS na série com diferença sazonal para d.

Código 
D_hat = nsdiffs(y, period=SEASON)
y_seasonal = y[SEASON * D_hat :] - y[: len(y) - SEASON * D_hat] if D_hat else y
d_hat = ndiffs(y_seasonal)
Markdown(
    f"`nsdiffs` → **D = {D_hat}** diferença(s) sazonal(is) · "
    f"`ndiffs` → **d = {d_hat}** diferença(s) regular(es)."
)

nsdiffs → D = 1 diferença(s) sazonal(is) · ndiffs → d = 1 diferença(s) regular(es).

Modelagem: AutoARIMA

O AutoARIMA(season_length=12) escolhe D e d com os testes acima e depois busca p, q, P e Q passo a passo, ficando com o modelo de menor AICc. Ele é uma ferramenta, não uma resposta: o diagnóstico dos resíduos abaixo decide se o modelo é confiável.

Código 
sf = StatsForecast(models=[AutoARIMA(season_length=SEASON)], freq="1mo", n_jobs=1)
sf.fit(df_sf)
model = sf.fitted_[0, 0].model_

# o statsforecast guarda as ordens como (p, q, P, Q, s, d, D)
p, q, P, Q, s, d, D = model["arma"]
order_label = f"SARIMA({p},{d},{q})({P},{D},{Q})[{s}]"
coefs = pl.DataFrame(
    {
        "termo": list(model["coef"]),
        "estimativa": [round(v, 4) for v in model["coef"].values()],
    }
)
Markdown(
    f"Modelo selecionado: **{order_label}** · AICc = {model['aicc']:.1f} · "
    f"σ² = {model['sigma2']:.2f}"
)

Modelo selecionado: SARIMA(3,1,2)(1,1,1)[12] · AICc = 1357.9 · σ² = 11.40

Código 
coefs
termo estimativa
"ar1" -0.3323
"ar2" 0.1965
"ar3" 0.221
"ma1" -0.1433
"ma2" -0.5179
"sar1" 0.161
"sma1" -0.8747

Diagnóstico dos resíduos

Os resíduos devem se comportar como ruído branco. Em ordem de prioridade: sem autocorrelação, variância constante e (desejável, mas não essencial) normalidade. Os primeiros d + D·s resíduos vêm do período de inicialização das diferenças e são descartados.

Código 
burn_in = d + D * s
resid = np.asarray(model["residuals"])[burn_in:]
resid_dates = dates[burn_in:]

fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes[0, 0].plot(resid_dates, resid, color=BLUE, lw=0.9)
axes[0, 0].axhline(0, color=GREY, lw=0.8)
axes[0, 0].set_title("Resíduos ao longo do tempo", loc="left")

correlogram(resid, axes[0, 1], "acf")
axes[0, 1].set_title("")
axes[0, 1].set_title("ACF dos resíduos", loc="left")

sns.histplot(resid, stat="density", color=OKABE_ITO[1], ax=axes[1, 0])
grid = np.linspace(resid.min(), resid.max(), 200)
axes[1, 0].plot(
    grid,
    stats.norm.pdf(grid, resid.mean(), resid.std()),
    color=VERMILION,
    label="Normal",
)
axes[1, 0].set_title("Distribuição vs normal", loc="left")
axes[1, 0].set_ylabel("Densidade")
axes[1, 0].legend(frameon=False)

stats.probplot(resid, dist="norm", plot=axes[1, 1])
axes[1, 1].get_lines()[0].set(color=BLUE, markersize=3)
axes[1, 1].get_lines()[1].set(color=VERMILION)
axes[1, 1].set_title("")  # remove o título padrão do scipy ("Probability Plot")
axes[1, 1].set_title("Q-Q normal", loc="left")
axes[1, 1].set(xlabel="Quantis teóricos", ylabel="Quantis observados")
plt.tight_layout()
plt.show()

O Ljung-Box testa se as primeiras h autocorrelações são conjuntamente zero. H0 é ruído branco, então um p-valor alto é o resultado bom. model_df = p + q + P + Q desconta os graus de liberdade usados pelos termos ARMA estimados.

Código 
lb = acorr_ljungbox(resid, lags=[12, 24, 36], model_df=p + q + P + Q)
ljung_box = pl.from_pandas(lb.reset_index(names="lag")).with_columns(
    pl.col("lb_stat").round(2),
    pl.col("lb_pvalue").round(4),
    pl.when(pl.col("lb_pvalue") > 0.05)
    .then(pl.lit("compatível com ruído branco"))
    .otherwise(pl.lit("resta autocorrelação"))
    .alias("leitura (5%)"),
)
jb_stat, jb_p = stats.jarque_bera(resid)
ljung_box
lag lb_stat lb_pvalue leitura (5%)
12 11.99 0.035 "resta autocorrelação"
24 28.47 0.0398 "resta autocorrelação"
36 41.43 0.0631 "compatível com ruído branco"

Comparação com o modelo airline

A leitura manual da ACF de ΔΔ₁₂ (picos isolados nos lags 1 e 12) sugere o modelo airline, SARIMA(0,1,1)(0,1,1)₁₂. Como ele tem as mesmas diferenças (d = 1, D = 1) do modelo escolhido, os AICc são comparáveis: diferença menor que 2 é empate, entre 2 e 10 é evidência moderada, acima de 10 é decisiva.

Código 
def fit_summary(data: pl.DataFrame, spec, label: str) -> tuple[dict, dict]:
    """Ajusta um modelo; resume ordens, AICc e Ljung-Box dos resíduos."""
    sf_ = StatsForecast(models=[spec], freq="1mo", n_jobs=1).fit(data)
    fitted = sf_.fitted_[0, 0].model_
    # maiúsculas = ordens sazonais, como na notação SARIMA
    p_, q_, P_, Q_, s_, d_, D_ = fitted["arma"]  # noqa: N806
    res = np.asarray(fitted["residuals"])[d_ + D_ * s_ :]
    pvals = acorr_ljungbox(res, lags=[12, 24, 36], model_df=p_ + q_ + P_ + Q_)
    summary = {
        "modelo": label,
        "ordens": f"({p_},{d_},{q_})({P_},{D_},{Q_})[{s_}]",
        "parâmetros ARMA": p_ + q_ + P_ + Q_,
        "AICc": round(fitted["aicc"], 1),
        **{
            f"LB p lag {h}": round(v, 3)
            for h, v in zip((12, 24, 36), pvals["lb_pvalue"], strict=True)
        },
    }
    return summary, fitted


auto_summary, _ = fit_summary(df_sf, AutoARIMA(season_length=SEASON), "AutoARIMA")
airline_summary, _ = fit_summary(
    df_sf,
    ARIMA(order=(0, 1, 1), seasonal_order=(0, 1, 1), season_length=SEASON),
    "Airline",
)
pl.DataFrame([auto_summary, airline_summary])
modelo ordens parâmetros ARMA AICc LB p lag 12 LB p lag 24 LB p lag 36
"AutoARIMA" "(3,1,2)(1,1,1)[12]" 7 1357.9 0.035 0.04 0.063
"Airline" "(0,1,1)(0,1,1)[12]" 2 1371.1 0.0 0.0 0.0

Sensibilidade ao choque de 2020

Um único choque extremo entra em todas as autocorrelações calculadas sobre ele. Para separar um problema de estrutura do modelo de um choque não modelado, o AutoARIMA é reajustado com variáveis dummy (intervenção) para março a junho de 2020.

Código 
SHOCK_MONTHS = [date(2020, m, 1) for m in (3, 4, 5, 6)]
df_shock = df_sf.with_columns(
    (pl.col("ds") == d_).cast(pl.Float64).alias(f"covid_{d_:%Y_%m}")
    for d_ in SHOCK_MONTHS
)
shock_summary, shock_model = fit_summary(
    df_shock, AutoARIMA(season_length=SEASON), "AutoARIMA + dummies 2020"
)
# o statsforecast nomeia os coeficientes das exógenas como ex_1, ex_2, ...
shock_effects = ", ".join(
    f"{month_pt(d_)}: {shock_model['coef'][f'ex_{i}']:+.1f}"
    for i, d_ in enumerate(SHOCK_MONTHS, start=1)
)
display(pl.DataFrame([auto_summary, shock_summary]))
Markdown(f"Efeito estimado das dummies (pontos do índice): {shock_effects}.")
modelo ordens parâmetros ARMA AICc LB p lag 12 LB p lag 24 LB p lag 36
"AutoARIMA" "(3,1,2)(1,1,1)[12]" 7 1357.9 0.035 0.04 0.063
"AutoARIMA + dummies 2020" "(2,1,2)(0,1,2)[12]" 6 1302.4 0.305 0.068 0.134

Efeito estimado das dummies (pontos do índice): mar/2020: -7.5, abr/2020: -21.4, mai/2020: -11.9, jun/2020: -6.7.

Previsão: 24 meses à frente

Código 
fc = sf.predict(h=HORIZON, level=LEVELS)
fc_cols = {
    "previsão": "AutoARIMA",
    **{
        f"{side}_{lv}": f"AutoARIMA-{side}-{lv}"
        for lv in LEVELS
        for side in ("lo", "hi")
    },
}

history = df_sf.filter(pl.col("ds") >= pl.date(last.year - 7, 1, 1))
fig, ax = plt.subplots(figsize=(11, 5))
ax.plot(history["ds"], history["y"], color=GREY, lw=1.2, label="Histórico")
ax.fill_between(
    fc["ds"],
    fc[fc_cols["lo_95"]],
    fc[fc_cols["hi_95"]],
    color=BLUE,
    alpha=0.15,
    lw=0,
    label="Intervalo de 95%",
)
ax.fill_between(
    fc["ds"],
    fc[fc_cols["lo_80"]],
    fc[fc_cols["hi_80"]],
    color=BLUE,
    alpha=0.30,
    lw=0,
    label="Intervalo de 80%",
)
ax.plot(fc["ds"], fc[fc_cols["previsão"]], color=BLUE, lw=2, label="Previsão")
ax.axvline(last, color=GREY, lw=0.8, ls=":")
ax.set(ylabel="Índice (2022 = 100)", xlabel=None)
ax.set_title(f"{order_label}: previsão de {HORIZON} meses, {UF}", pad=30)
# legenda em uma linha acima do gráfico, para nunca cobrir os dados
ax.legend(frameon=False, ncols=4, loc="lower left", bbox_to_anchor=(0, 1.0))
plt.show()

A previsão é desenhada como um gráfico em leque: duas faixas aninhadas, com a de 80% (mais escura) dentro da de 95%, para que a incerteza seja lida como parte da previsão e não como uma linha a mais. Só os últimos 7 anos de histórico aparecem, para o período previsto não ficar espremido na borda direita.

Código 
forecast_table = fc.select(
    pl.col("ds").dt.strftime("%Y-%m").alias("mês"),
    *[pl.col(c).round(1).alias(k) for k, c in fc_cols.items()],
).with_columns((pl.col("hi_95") - pl.col("lo_95")).round(1).alias("amplitude_95"))
forecast_table
mês previsão lo_80 hi_80 lo_95 hi_95 amplitude_95
"2026-03" 104.1 99.8 108.4 97.5 110.7 13.2
"2026-04" 100.7 95.8 105.5 93.2 108.1 14.9
"2026-05" 105.7 100.6 110.8 97.9 113.6 15.7
"2026-06" 103.2 97.6 108.9 94.6 111.9 17.3
"2026-07" 107.6 101.7 113.4 98.7 116.5 17.8
"2026-08" 108.0 101.8 114.1 98.6 117.3 18.7
"2026-09" 106.8 100.5 113.1 97.1 116.5 19.4
"2026-10" 110.6 104.1 117.1 100.6 120.6 20.0
"2026-11" 111.3 104.6 118.1 101.0 121.6 20.6
"2026-12" 121.4 114.5 128.3 110.8 132.0 21.2
"2027-01" 104.1 97.0 111.2 93.2 115.0 21.8
"2027-02" 98.2 90.9 105.5 87.0 109.3 22.3
"2027-03" 106.4 98.6 114.3 94.5 118.4 23.9
"2027-04" 102.7 94.5 110.8 90.2 115.1 24.9
"2027-05" 107.5 99.1 115.9 94.7 120.4 25.7
"2027-06" 105.6 96.9 114.3 92.3 118.9 26.6
"2027-07" 109.7 100.7 118.6 96.0 123.4 27.4
"2027-08" 110.5 101.3 119.7 96.5 124.6 28.1
"2027-09" 108.5 99.1 118.0 94.1 122.9 28.8
"2027-10" 112.3 102.6 121.9 97.5 127.0 29.5
"2027-11" 113.4 103.5 123.3 98.3 128.5 30.2
"2027-12" 123.3 113.3 133.4 107.9 138.7 30.8
"2028-01" 106.0 95.7 116.3 90.3 121.7 31.4
"2028-02" 100.1 89.6 110.6 84.1 116.2 32.1

Discussão

Evidências para apoiar as respostas

Estes números saem diretamente da análise acima. A interpretação é minha.

Código 
width = forecast_table["amplitude_95"]
resid_acf = np.array([np.corrcoef(resid[:-k], resid[k:])[0, 1] for k in (12, 24)])
band = 1.96 / np.sqrt(len(resid))
all_acf = np.array([np.corrcoef(resid[:-k], resid[k:])[0, 1] for k in range(1, 37)])
n_outside = int((np.abs(all_acf) > band).sum())
lb_pvalues = " / ".join(f"{v:.3f}" for v in lb["lb_pvalue"])
worst = np.abs(resid).argmax()
last12_mean = df_sf["y"].tail(12).mean()
fc12_mean = fc[fc_cols["previsão"]].head(12).mean()

evidence = [
    f"- **Sazonalidade / tendência (STL):** força sazonal {seasonal_strength:.2f}, "
    f"força da tendência {trend_strength:.2f}.",
    f"- **Diferenças:** `nsdiffs` D = {D_hat}, `ndiffs` d = {d_hat}; "
    f"o AutoARIMA usou d = {d}, D = {D}.",
    f"- **Modelo selecionado:** {order_label}, com termos AR p = {p}, P = {P} "
    f"e termos MA q = {q}, Q = {Q}.",
    f"- **P-valores do Ljung-Box** (lags 12/24/36): {lb_pvalues}.",
    f"- **Autocorrelação dos resíduos nos lags sazonais:** "
    f"lag 12 = {resid_acf[0]:.2f}, lag 24 = {resid_acf[1]:.2f} "
    f"(banda de 95% ±{band:.2f}); {n_outside} de 36 lags fora da banda "
    f"(por acaso, espera-se cerca de 2).",
    f"- **Normalidade (Jarque-Bera):** p = {jb_p:.3g}; maior resíduo em módulo "
    f"{abs(resid[worst]):.1f} em {month_pt(resid_dates[worst])}.",
    f"- **Nível:** média dos últimos 12 meses observados {last12_mean:.1f}; "
    f"média dos 12 primeiros meses previstos {fc12_mean:.1f}.",
    f"- **Amplitude do intervalo de 95%:** {width[0]:.1f} em h = 1, "
    f"{width[11]:.1f} em h = 12, {width[23]:.1f} em h = 24 pontos.",
]
Markdown("\n".join(evidence))
  • Sazonalidade / tendência (STL): força sazonal 0.82, força da tendência 0.98.
  • Diferenças: nsdiffs D = 1, ndiffs d = 1; o AutoARIMA usou d = 1, D = 1.
  • Modelo selecionado: SARIMA(3,1,2)(1,1,1)[12], com termos AR p = 3, P = 1 e termos MA q = 2, Q = 1.
  • P-valores do Ljung-Box (lags 12/24/36): 0.035 / 0.040 / 0.063.
  • Autocorrelação dos resíduos nos lags sazonais: lag 12 = 0.03, lag 24 = -0.07 (banda de 95% ±0.12); 0 de 36 lags fora da banda (por acaso, espera-se cerca de 2).
  • Normalidade (Jarque-Bera): p = 2.61e-17; maior resíduo em módulo 17.0 em abr/2020.
  • Nível: média dos últimos 12 meses observados 104.1; média dos 12 primeiros meses previstos 106.8.
  • Amplitude do intervalo de 95%: 13.2 em h = 1, 22.3 em h = 12, 32.1 em h = 24 pontos.

Modelo e sazonalidade por UF

O mesmo AutoARIMA(season_length=12) ajustado para todas as UFs, com a força sazonal da STL e o tamanho do pico de dezembro (mediana do índice sem tendência, em pontos).

Código 
all_ufs = pmc.select(
    pl.col("uf").alias("unique_id"),
    pl.col("data").alias("ds"),
    pl.col("volume_idx").cast(pl.Float64).alias("y"),
).sort("unique_id", "ds")
sf_all = StatsForecast(models=[AutoARIMA(season_length=SEASON)], freq="1mo", n_jobs=1)
sf_all.fit(all_ufs)

rows = []
for uid, fitted in zip(sf_all.uids, sf_all.fitted_[:, 0], strict=True):
    p_, q_, P_, Q_, s_, d_, D_ = fitted.model_["arma"]
    uf_df = all_ufs.filter(pl.col("unique_id") == uid)
    y_uf = uf_df["y"].to_numpy()
    stl_uf = STL(y_uf, period=SEASON, robust=True).fit()
    var_ru = np.var(stl_uf.resid)
    seasonal_strength_uf = 1 - var_ru / np.var(stl_uf.seasonal + stl_uf.resid)
    december = uf_df["ds"].dt.month().to_numpy() == 12
    rows.append(
        {
            "UF": uid,
            "modelo": f"({p_},{d_},{q_})({P_},{D_},{Q_})[{s_}]",
            "d": d_,
            "força sazonal": round(seasonal_strength_uf, 2),
            "pico de dez.": round(float(np.median((y_uf - stl_uf.trend)[december])), 1),
            "crescimento": round(y_uf[-12:].mean() / y_uf[:12].mean(), 1),
            "σ resíduos": round(float(np.sqrt(fitted.model_["sigma2"])), 1),
        }
    )
pl.DataFrame(rows).sort("pico de dez.", descending=True).with_columns(
    (pl.col("UF") == UF).alias("minha UF")
)
UF modelo d força sazonal pico de dez. crescimento σ resíduos minha UF
"Rio de Janeiro" "(1,1,1)(1,1,1)[12]" 1 0.84 24.2 1.5 4.2 false
"Bahia" "(2,1,1)(0,1,2)[12]" 1 0.78 22.3 1.8 5.0 false
"Alagoas" "(1,1,1)(0,1,2)[12]" 1 0.81 21.3 3.0 4.4 false
"Rio Grande do Norte" "(0,1,2)(1,1,2)[12]" 1 0.83 20.9 2.4 4.1 false
"Pará" "(1,1,1)(0,1,2)[12]" 1 0.83 20.6 3.0 4.1 false
"Rio Grande do Sul" "(2,0,1)(2,1,1)[12]" 0 0.84 20.2 1.9 4.0 false
"Pernambuco" "(1,1,1)(0,1,2)[12]" 1 0.75 20.1 2.3 4.7 false
"Amapá" "(2,0,0)(2,1,0)[12]" 0 0.59 19.6 3.2 6.2 false
"Sergipe" "(0,1,1)(0,1,2)[12]" 1 0.79 19.6 2.3 4.9 false
"Paraíba" "(1,1,1)(0,1,1)[12]" 1 0.73 18.4 3.0 4.7 false
"Piauí" "(1,1,1)(0,1,2)[12]" 1 0.76 17.8 2.8 4.8 false
"Ceará" "(0,1,3)(1,1,2)[12]" 1 0.68 17.4 2.9 4.8 false
"Acre" "(1,1,1)(0,1,1)[12]" 1 0.71 17.2 5.0 5.3 false
"São Paulo" "(2,1,1)(0,1,1)[12]" 1 0.77 17.2 1.7 4.0 false
"Amazonas" "(2,0,2)(1,1,1)[12]" 0 0.7 17.0 2.7 4.8 false
"Paraná" "(1,1,1)(1,1,1)[12]" 1 0.79 16.4 1.9 4.1 false
"Santa Catarina" "(2,0,1)(1,1,1)[12]" 0 0.8 16.1 2.9 3.4 false
"Maranhão" "(1,1,1)(0,1,2)[12]" 1 0.73 15.9 3.7 4.6 false
"Distrito Federal" "(1,1,1)(0,1,2)[12]" 1 0.71 15.4 1.4 5.5 false
"Rondônia" "(1,1,1)(0,1,2)[12]" 1 0.71 14.9 3.9 5.4 false
"Minas Gerais" "(3,1,2)(1,1,1)[12]" 1 0.82 14.0 2.0 3.4 true
"Mato Grosso do Sul" "(1,1,2)(0,1,2)[12]" 1 0.81 12.5 2.3 3.6 false
"Espírito Santo" "(1,1,2)(1,1,2)[12]" 1 0.55 11.6 3.1 5.3 false
"Goiás" "(0,1,1)(1,1,1)[12]" 1 0.65 11.6 2.0 4.9 false
"Mato Grosso" "(3,0,0)(0,1,2)[12]" 0 0.77 10.3 2.5 3.5 false
"Roraima" "(2,0,1)(0,1,1)[12]" 0 0.68 9.0 4.5 3.1 false
"Tocantins" "(2,0,2)(1,1,1)[12]" 0 0.61 8.0 5.8 5.2 false

Crescimento: média dos últimos 12 meses dividida pela média dos 12 primeiros (2004).

Respostas

  1. O ACF/PACF da série original sugeriu sazonalidade? E tendência? Bate com o que você viu no gráfico de linha?

    Sim, as duas coisas. A ACF decai muito lentamente (0,89 no lag 1 e ainda cerca de 0,45 no lag 36), a assinatura de uma série não estacionária, com tendência; a PACF, com um pico dominante no lag 1 (0,89), reforça a raiz unitária. Sobre esse decaimento aparecem picos nos lags 12 e 24 (0,83 e 0,64, acima dos vizinhos), e a PACF tem picos fortes nos lags 12 e 13: sazonalidade anual. Isso bate com o gráfico de linha, em que o índice quase dobra entre 2004 e 2013 (de cerca de 52 para 96), recua na recessão de 2015–2016 e retoma depois, com pico todo dezembro e vale em fevereiro (+14 e −10 pontos em relação à tendência). A STL confirma: força da tendência 0,98 e sazonal 0,82.

  2. Quantas diferenças (regular e sazonal) o AutoARIMA aplicou? E os componentes AR/MA? Conseguiu observar isso no ACF/PACF?

    Uma de cada (d = 1, D = 1), coerente com os testes: na série original o ADF não rejeita e o KPSS rejeita; após Δ₁₂ o resultado fica inconclusivo e a ACF ainda decai devagar até perto do lag 23; só em ΔΔ₁₂ os dois testes concordam em estacionariedade, como indicam também nsdiffs e ndiffs. Na ΔΔ₁₂, os picos negativos nos lags 1 (−0,39) e 12 (−0,37) sugerem, pela leitura manual, o modelo airline SARIMA(0,1,1)(0,1,1)₁₂. O AutoARIMA escolheu um modelo mais rico, SARIMA(3,1,2)(1,1,1)₁₂, que capta também a estrutura dos lags 2 a 4 (PACF significativa em 1 e 2, ACF em 3 e 4), menos óbvia no gráfico. A comparação justifica a complexidade: o airline tem AICc 13 pontos pior (evidência decisiva) e deixa resíduos autocorrelacionados (Ljung-Box p < 0,001). O MA sazonal de −0,88, próximo de −1, sugere uma sazonalidade quase determinística.

  3. Os resíduos parecem white noise? Algum padrão remanescente sazonal?

    Quase, mas não totalmente. Nenhum dos 36 lags da ACF dos resíduos sai da banda de ±0,12, e os lags sazonais 12 e 24 são pequenos (0,03 e −0,07): não resta padrão sazonal. Ainda assim o Ljung-Box rejeita ruído branco nos lags 12 e 24 (p = 0,035 e 0,040), porque testa a soma de muitas autocorrelações pequenas. A análise de sensibilidade mostra a causa principal: com dummies para março a junho de 2020, os p-valores sobem para 0,31, 0,07 e 0,13, e o efeito estimado de abril é de cerca de −21 pontos. Ou seja, o problema é o choque da pandemia, não a estrutura do modelo. A normalidade é rejeitada (caudas pesadas), o que afeta os intervalos de previsão, não a previsão pontual.

  4. A previsão pareceu razoável visualmente (captura a sazonalidade, fica bem ancorada no histórico, intervalos abrem muito com o horizonte)?

    Sim, com ressalvas. A previsão reproduz a sazonalidade (picos de dezembro em cerca de 121 e 123, vales em fevereiro) e começa ancorada no histórico: a média dos 12 primeiros meses previstos (106,8) fica perto da média dos últimos 12 observados (104,1). Mas ela projeta crescimento contínuo só por inércia: o modelo não conhece juros, crédito ou renda. Os intervalos de 95% abrem de 13 para 32 pontos em 24 meses, o esperado com d = 1 e D = 1, em que a incerteza se acumula; no segundo ano a amplitude já supera a própria oscilação sazonal (cerca de 24 pontos entre fevereiro e dezembro). Com resíduos de caudas pesadas, esses intervalos tendem a ser otimistas; o próximo passo seria validar fora da amostra (MASE contra o SeasonalNaive).

  5. Que outras UFs você imagina que dariam um SARIMA bem diferente do seu (e por quê)? Pense em estados de regiões diferentes, com perfil econômico distinto.

    A tabela por UF mostra que o mês do pico não muda: as 27 UFs têm pico em dezembro e vale em fevereiro. O que muda é a intensidade e a dinâmica. Rio de Janeiro (+24), Bahia (+22), Alagoas e Rio Grande do Norte (+21) têm o pico de Natal mais forte; Tocantins (+8), Roraima (+9), Mato Grosso (+10) e Goiás (+12) têm o mais fraco. Minha hipótese é que, no Centro-Norte, a renda do agronegócio e do setor público se distribui melhor ao longo do ano, enquanto no litoral o 13º salário e o turismo de verão concentram o consumo em dezembro. Na tendência, as 7 UFs que ficaram com d = 0 (entre elas Tocantins, Roraima e Mato Grosso) cresceram mais em média (3,4× contra 2,6× desde 2004), já que a diferença sazonal absorve uma tendência regular; mas não é regra: o Acre cresceu 5× e manteve d = 1. Já economias pequenas e concentradas, como o Amapá (o maior σ dos resíduos, 6,2), geram modelos mais instáveis. Minas, diversificada, fica no meio: pico moderado (+14) e um dos menores ruídos (σ = 3,4).