CO\(_2\) Atmosférico em Mauna Loa

Trabalho 1 — Análise de Séries Temporais

# --- Séries temporais e modelagem ---
library(forecast)    # auto.arima, Arima, forecast, ACF/PACF via ggplot
library(tseries)     # adf.test, kpss.test: testes de raiz unitária
library(TSA)         # eacf (ESACF para identificação de ordens)

# --- Visualização ---
library(ggplot2)     # base gráfica
library(patchwork)   # composição de painéis
library(scales)      # formatação de eixos
library(ggfortify)   # autoplot para objetos ts e decompose

# --- Tabelas ---
library(gt)          # tabelas com formatação avançada
library(broom)       # tidy() para extrair coeficientes de modelos

# --- Manipulação ---
library(tidyverse)   # dplyr, tidyr, stringr, purrr

theme_set(theme_minimal(base_size = 12))

# Paleta de cores consistente
cores <- c("#2E86AB", "#E94F37", "#3BB273", "#F4A261", "#8338EC")

# Função auxiliar: ACF/PACF → data.frame com lag em meses (inteiro)
acf_df <- function(x, lag.max, type = "correlation") {
  obj <- acf(x, lag.max = lag.max, type = type, plot = FALSE)
  tibble(
    lag = as.integer(round(obj$lag * frequency(x))),  # converte anos → meses
    acf = as.numeric(obj$acf)
  ) |>
    filter(lag > 0 | type == "correlation")           # remove lag-0 da PACF
}

pacf_df <- function(x, lag.max) {
  obj <- pacf(x, lag.max = lag.max, plot = FALSE)
  tibble(
    lag = as.integer(round(obj$lag * frequency(x))),
    acf = as.numeric(obj$acf)
  )
}

plot_acf <- function(df, ic, title, ylab = "ACF",
                     vlines = NULL, vlab = NULL) {
  p <- ggplot(df, aes(x = lag, y = acf)) +
    geom_hline(yintercept = 0, color = "gray50") +
    geom_segment(aes(xend = lag, yend = 0), color = cores[1], linewidth = 0.7) +
    geom_hline(yintercept = c(-ic, ic),
               linetype = "dashed", color = cores[2], linewidth = 0.7) +
    labs(title = title, x = "Defasagem (meses)", y = ylab) +
    scale_x_continuous(breaks = seq(0, max(df$lag), by = 6))
  if (!is.null(vlines)) {
    p <- p +
      geom_vline(xintercept = vlines,
                 linetype = "dotted", color = cores[4], linewidth = 0.5)
    if (!is.null(vlab))
      p <- p + annotate("text", x = vlines, y = max(df$acf) * 0.9,
                        label = paste0("lag ", vlines),
                        size = 2.8, color = cores[4])
  }
  p
}

Introdução

O dióxido de carbono (\(\text{CO}_2\)) é o principal gás de efeito estufa de origem antrópica, responsável por aproximadamente dois terços do aquecimento global antropogênico. Sua concentração na atmosfera é monitorada continuamente desde março de 1958 pelo Observatório de Mauna Loa, localizado a 3.397 metros de altitude no vulcão Mauna Loa, na ilha de Havaí. A escolha desse local não foi acidental: a altitude elevada, o isolamento geográfico no meio do Oceano Pacífico e a distância de fontes locais significativas de poluição garantem medições representativas da concentração global de \(\text{CO}_2\) na atmosfera, minimizando influências urbanas, industriais ou de variações locais da vegetação.

As medições iniciadas pelo cientista Charles David Keeling geraram a série histórica mais longa, precisa e cientificamente importante de \(\text{CO}_2\) atmosférico disponível — conhecida mundialmente como Curva de Keeling (Keeling et al., 1976). Esta série documenta de forma inequívoca o aumento contínuo das concentrações de \(\text{CO}_2\) desde o início das medições sistemáticas e constitui evidência fundamental nas discussões sobre mudanças climáticas globais.

Este trabalho analisa a série co2 disponível nativamente no R, que contém as concentrações médias mensais de \(\text{CO}_2\) (em partes por milhão, ppm) de janeiro de 1959 a dezembro de 1997 — 468 observações mensais. A série é amplamente estudada na literatura de séries temporais por reunir, em um único conjunto de dados, dois fenômenos de grande interesse estatístico e científico:

1. Tendência crescente de longo prazo: reflexo do acúmulo antropogênico de \(\text{CO}_2\) proveniente principalmente da queima de combustíveis fósseis (petróleo, carvão, gás natural) e, em menor escala, de mudanças no uso da terra e desmatamento.

2. Sazonalidade anual regular: reflexo do ciclo sazonal de fotossíntese e respiração vegetal, predominantemente do hemisfério norte, onde se concentra a maior parte da massa terrestre vegetada do planeta. No inverno boreal (outubro a maio), a respiração das plantas supera a fotossíntese, liberando \(\text{CO}_2\) na atmosfera; no verão boreal (maio a setembro), a fotossíntese intensificada absorve \(\text{CO}_2\), reduzindo sua concentração atmosférica.

Essa combinação de tendência não estacionária e sazonalidade estocástica a torna um caso de referência clássico para a metodologia SARIMA (Seasonal Autoregressive Integrated Moving Average).

Pergunta de interesse: É possível ajustar um modelo SARIMA parcimonioso que capture simultaneamente a tendência e a sazonalidade da série e produza previsões confiáveis para horizontes de até dois anos?

A análise segue rigorosamente a metodologia de Brockwell; Davis (2010) e Morettin; Toloi (2004), e está organizada em quatro etapas: (1) análise descritiva — caracterização da tendência, sazonalidade e componentes da série via decomposição STL; (2) identificação e ajuste — seleção de modelos SARIMA candidatos por critérios de informação (AIC, AICc, BIC) e estimação de parâmetros por máxima verossimilhança; (3) diagnóstico — validação das suposições do modelo via análise de resíduos (autocorrelação, normalidade, homoscedasticidade); e (4) previsão — projeções para 1998-1999 com quantificação rigorosa da incerteza via intervalos de confiança.

---

Descrição dos Dados e Análise Exploratória

A Série Temporal

# Carregar a série nativa do R (1959-1997)
data(co2, package = "datasets")

# Converter para tibble: time(co2) retorna anos decimais (ex: 1959.000, 1959.083...)
co2_df <- tibble(
  tempo = as.numeric(time(co2)),
  ano   = as.integer(floor(tempo)),
  mes   = as.integer(round((tempo - floor(tempo)) * 12)) + 1L,
  ppm   = as.numeric(co2)
) |>
  mutate(mes = ifelse(mes > 12L, 12L, mes))

n      <- length(co2)
inicio <- start(co2)
fim    <- end(co2)
freq   <- frequency(co2)

# Calcular taxa de crescimento (para uso ao longo do documento)
crescimento <- round(
  (mean(co2[time(co2) >= 1997]) - mean(co2[time(co2) < 1960])) / 39, 2
)

cat(sprintf("Período       : %d/%02d a %d/%02d\n", inicio[1], inicio[2], fim[1], fim[2]))

Período       : 1959/01 a 1997/12

cat(sprintf("Observações   : %d (mensais, s = %d)\n", n, freq))

Observações   : 468 (mensais, s = 12)

cat(sprintf("Mínimo / Máx. : %.2f / %.2f ppm\n", min(co2), max(co2)))

Mínimo / Máx. : 313.18 / 366.84 ppm

cat(sprintf("Média global  : %.2f ppm\n", mean(co2)))

Média global  : 337.05 ppm

cat(sprintf("Desvio padrão : %.2f ppm\n", sd(co2)))

Desvio padrão : 14.97 ppm

A série compreende 468 observações mensais de 1959 a 1997, com concentração de \(\text{CO}_2\) variando de 313.18 a 366.84 ppm ao longo das quase quatro décadas de registro. Trata-se de uma série temporal univariada, mensal (escala de tempo: meses), de natureza contínua, coletada pelo Observatório de Mauna Loa (Havaí, EUA) e disponibilizada como conjunto de dados nativo do R (datasets::co2), originalmente publicada por Keeling et al. (1976).

Visualização da Série

ggplot(co2_df, aes(x = tempo, y = ppm)) +
  geom_line(color = cores[1], linewidth = 0.6, alpha = 0.9) +
  geom_smooth(method = "loess", span = 0.15,
              color = cores[2], linewidth = 1.0,
              se = FALSE, linetype = "dashed") +
  scale_x_continuous(breaks = seq(1960, 1995, by = 5)) +
  scale_y_continuous(labels = label_number(suffix = " ppm")) +
  labs(
    title    = expression(CO[2] ~ "Atmosférico -- Mauna Loa Observatory"),
    subtitle = "Linha tracejada: tendência LOESS | Oscilações: sazonalidade anual",
    x        = NULL,
    y        = expression(CO[2] ~ "(ppm)")
  )

Concentração mensal de \(\text{CO}_2\) em Mauna Loa (1959–1997). Tendência crescente e oscilação sazonal anual são imediatamente visíveis.

A figura evidencia as duas componentes dominantes da série:

• Tendência: crescimento monotônico e aproximadamente linear ao longo das quase quatro décadas, a uma taxa média de 1.23 ppm/ano (de 313 ppm em 1959 para 367 ppm em 1997).
• Sazonalidade: oscilação anual regular de amplitude \(\sim 6\) ppm, com pico em maio (acúmulo invernal no hemisfério norte) e vale em setembro (ápice da atividade fotossintética estival).

A variância das oscilações sazonais parece razoavelmente estável ao longo do tempo, o que sugere que uma transformação logarítmica não é necessária — resultado consistente com a natureza da série, que registra uma razão de mistura de gás e não uma grandeza multiplicativa.

Decomposição Clássica

stl_fit <- stl(co2, s.window = "periodic", robust = TRUE)

autoplot(stl_fit) +
  labs(
    title    = expression("Decomposição STL --" ~ CO[2] ~ "Mauna Loa"),
    subtitle = "Tendência estimada por LOESS; sazonalidade periódica"
  ) +
  theme(strip.text = element_text(face = "bold"))

Decomposição STL da série co2: tendência, sazonalidade e resíduo. A barra cinza à direita de cada painel indica a escala relativa de cada componente.

A decomposição STL (Cleveland et al., 1990) confirma a estrutura visual:

• Tendência: perfeitamente monotônica, sem inflexões relevantes — sugestiva de diferenciação simples de ordem 1 para induzir estacionariedade.
• Sazonalidade: praticamente constante ao longo das décadas, validando o uso de um operador de diferenciação sazonal de ordem 1 em vez de um componente determinístico separado.
• Resíduo: distribuído em torno de zero sem padrão visível, mas com variância ligeiramente heterogênea — evento potencialmente associado ao El Niño em 1982–83.

Padrão Sazonal

meses_abr <- c("Jan","Fev","Mar","Abr","Mai","Jun",
                "Jul","Ago","Set","Out","Nov","Dez")

p1 <- co2_df |>
  mutate(mes_f = factor(mes, levels = 1:12, labels = meses_abr)) |>
  ggplot(aes(x = mes_f, y = ppm, fill = mes_f)) +
  geom_boxplot(alpha = 0.75, outlier.size = 0.8,
               outlier.alpha = 0.4, show.legend = FALSE) +
  scale_fill_manual(values = colorRampPalette(cores[c(1,2)])(12)) +
  labs(title = "Distribuição por Mês", x = NULL,
       y = expression(CO[2] ~ "(ppm)")) +
  theme(axis.text.x = element_text(size = 8))

p2 <- co2_df |>
  mutate(mes_f = factor(mes, labels = meses_abr)) |>
  ggplot(aes(x = mes, y = ppm, group = ano, color = ano)) +
  geom_line(alpha = 0.55, linewidth = 0.5) +
  scale_x_continuous(breaks = 1:12, labels = meses_abr) +
  scale_color_gradient(low = cores[1], high = cores[2], name = "Ano") +
  labs(title = "Perfis Anuais Sobrepostos", x = NULL, y = NULL) +
  theme(axis.text.x = element_text(size = 8),
        legend.key.height = unit(0.4, "cm"))

p1 + p2

Padrão sazonal de \(\text{CO}_2\): (esq.) boxplot por mês; (dir.) perfis anuais sobrepostos. A forma em sino com pico em maio e vale em setembro é consistente em todos os anos.

O padrão sazonal é extremamente regular: o \(\text{CO}_2\) aumenta de outubro a maio (menor atividade fotossintética no hemisfério norte) e diminui de maio a setembro. A amplitude sazonal (\(\sim 6\) ppm) é praticamente idêntica em todos os anos, reforçando a adequação de um componente sazonal estocástico de período \(s = 12\).

Taxa de Crescimento Anual

co2_df |>
  group_by(ano) |>
  summarise(media_anual = mean(ppm), .groups = "drop") |>
  mutate(delta = c(NA, diff(media_anual))) |>
  filter(!is.na(delta)) |>
  ggplot(aes(x = ano, y = delta)) +
  geom_col(fill = cores[3], alpha = 0.85, width = 0.7) +
  geom_smooth(method = "loess", se = FALSE,
              color = cores[2], linewidth = 1.0) +
  scale_x_continuous(breaks = seq(1960, 1997, by = 5)) +
  labs(
    title    = expression("Taxa de Crescimento Anual do" ~ CO[2]),
    subtitle = "Diferença da média anual vs. ano anterior",
    x        = "Ano",
    y        = "ppm / ano"
  )

Taxa de crescimento anual do \(\text{CO}_2\) (ppm/ano): média por ano após remoção da sazonalidade mensal. O crescimento é consistente, com aceleração visível após 1980.

---

Identificação do Modelo

Verificação de Estacionariedade

A série original claramente não é estacionária (fraca): apresenta média crescente no tempo e, portanto, \(E(Y_t)\) varia com \(t\). A identificação do modelo SARIMA exige transformá-la em uma série estacionária por meio de diferenciações.

Função de Autocorrelação da Série Original

ic_orig <- 1.96 / sqrt(n)

p_acf <- plot_acf(
  acf_df(co2, lag.max = 48) |> filter(lag > 0),
  ic    = ic_orig,
  title = "ACF -- Série Original"
)

p_pacf <- plot_acf(
  pacf_df(co2, lag.max = 48),
  ic    = ic_orig,
  title = "PACF -- Série Original",
  ylab  = "PACF"
)

p_acf + p_pacf

ACF e PACF da série original (sem diferenciação). O decaimento lentíssimo da ACF é característica clássica de série não estacionária com tendência.

A ACF decai para zero de forma extremamente lenta (ainda significativa na defasagem 48), confirmando a não estacionariedade. O padrão senoidal com período 12 evidencia a sazonalidade.

Diferenciação Simples (\(d = 1\))

co2_d1   <- diff(co2, differences = 1)
n_d1     <- length(co2_d1)
ic_d1    <- 1.96 / sqrt(n_d1)
tempo_d1 <- as.numeric(time(co2_d1))

p_d1 <- tibble(t = tempo_d1, y = as.numeric(co2_d1)) |>
  ggplot(aes(x = t, y = y)) +
  geom_line(color = cores[1], linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
  scale_x_continuous(breaks = seq(1960, 1997, by = 5)) +
  labs(title    = "Série Diferenciada (d = 1)",
       x = NULL, y = expression(nabla ~ CO[2] ~ "(ppm)"))

p_acf_d1 <- plot_acf(
  acf_df(co2_d1, lag.max = 48) |> filter(lag > 0),
  ic    = ic_d1,
  title = "ACF -- Após d = 1"
)

p_d1 / p_acf_d1

Série após diferenciação simples (\(d = 1\)). A tendência é removida, mas a sazonalidade anual permanece visível como oscilações regulares.

Após \(\nabla Y_t = Y_t - Y_{t-1}\), a tendência é eliminada, mas a ACF ainda apresenta picos sazonais nas defasagens 12, 24, 36 — indicando que uma diferenciação sazonal adicional é necessária.

Diferenciação Sazonal (\(d = 1\), \(D = 1\), \(s = 12\))

co2_d1s1   <- diff(diff(co2, lag = 12), differences = 1)
n_d1s1     <- length(co2_d1s1)
ic_d1s1    <- 1.96 / sqrt(n_d1s1)
tempo_d1s1 <- as.numeric(time(co2_d1s1))

p_d1s1 <- tibble(t = tempo_d1s1, y = as.numeric(co2_d1s1)) |>
  ggplot(aes(x = t, y = y)) +
  geom_line(color = cores[1], linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray60") +
  scale_x_continuous(breaks = seq(1960, 1997, by = 5)) +
  labs(title = "Série Diferenciada (d = 1, D = 1, s = 12)",
       x = NULL,
       y = expression(nabla ~ nabla[12] ~ CO[2] ~ "(ppm)"))

p_acf_d1s1 <- plot_acf(
  acf_df(co2_d1s1, lag.max = 48) |> filter(lag > 0),
  ic       = ic_d1s1,
  title    = "ACF -- Após d = 1, D = 1",
  vlines   = c(12, 24, 36),
  vlab     = TRUE
)

p_pacf_d1s1 <- plot_acf(
  pacf_df(co2_d1s1, lag.max = 48),
  ic       = ic_d1s1,
  title    = "PACF -- Após d = 1, D = 1",
  ylab     = "PACF",
  vlines   = c(12, 24, 36),
  vlab     = TRUE
)

p_d1s1 / (p_acf_d1s1 + p_pacf_d1s1)

Série após diferenciações simples e sazonal (\(d = 1\), \(D = 1\), \(s = 12\)). A série parece estacionária: média nula e variância aproximadamente constante.

Após \(\nabla \nabla_{12} Y_t = (1-B)(1-B^{12})Y_t\), a série é visualmente estacionária: flutuação em torno de zero sem tendência aparente e sem oscilação sazonal sistemática.

Testes Formais de Raiz Unitária

# suppressWarnings: os warnings indicam apenas que o p-valor é menor/maior
# que o valor impresso na tabela de referência do teste (sem impacto no resultado)
testar_estacionariedade <- function(serie, nome) {
  adf  <- suppressWarnings(adf.test(serie, alternative = "stationary"))
  kpss <- suppressWarnings(kpss.test(serie, null = "Level"))
  tibble(
    Serie               = nome,
    ADF_p               = round(adf$p.value,  4),
    KPSS_p              = round(kpss$p.value, 4),
    ADF_H0              = ifelse(adf$p.value  < 0.05, "Sim (estacionária)", "Não"),
    KPSS_H0             = ifelse(kpss$p.value < 0.05, "Sim (não estac.)",   "Não")
  )
}

testes <- bind_rows(
  testar_estacionariedade(co2,      "Original"),
  testar_estacionariedade(co2_d1,   "d = 1"),
  testar_estacionariedade(co2_d1s1, "d = 1, D = 1")
)

testes |>
  gt() |>
  tab_header(
    title    = md("**Tabela 1.** Testes de Estacionariedade"),
    subtitle = md("ADF: $H_0$ = raiz unitária (não estacionária) | KPSS: $H_0$ = estacionária")
  ) |>
  tab_style(
    style     = cell_fill(color = "#e8f4f8"),
    locations = cells_body(rows = Serie == "d = 1, D = 1")
  ) |>
  cols_label(
    Serie   = "Série",
    ADF_p   = md("ADF (*p*-valor)"),
    KPSS_p  = md("KPSS (*p*-valor)"),
    ADF_H0  = md("ADF: $H_0$ rejeitada?"),
    KPSS_H0 = md("KPSS: $H_0$ rejeitada?")
  )

Tabela 1. Testes de Estacionariedade
ADF: \(H_0\) = raiz unitária (não estacionária) | KPSS: \(H_0\) = estacionária
Série	ADF (p-valor)	KPSS (p-valor)	ADF: \(H_0\) rejeitada?	KPSS: \(H_0\) rejeitada?
Original	0.2269	0.01	Não	Sim (não estac.)
d = 1	0.0100	0.10	Sim (estacionária)	Não
d = 1, D = 1	0.0100	0.10	Sim (estacionária)	Não

Os testes formais confirmam o diagnóstico gráfico:

• Série original: ADF não rejeita \(H_0\) (raiz unitária presente); KPSS rejeita \(H_0\) (evidência contra estacionariedade). Conclusão: não estacionária.
• Após \(d=1\): situação intermediária — tendência removida, mas sazonalidade ainda viola a estacionariedade.
• Após \(d=1\), \(D=1\): ADF rejeita \(H_0\) (\(p < 0{,}05\)) e KPSS não rejeita \(H_0\). Conclusão: série estacionária, pronta para a identificação de \(p, q, P, Q\).

Leitura da ACF e PACF para Identificação de Ordens

p_acf_f <- plot_acf(
  acf_df(co2_d1s1, lag.max = 48) |> filter(lag > 0),
  ic       = ic_d1s1,
  title    = expression("ACF --" ~ nabla ~ nabla[12] ~ Y[t]),
  vlines   = c(12, 24, 36),
  vlab     = TRUE
)

p_pacf_f <- plot_acf(
  pacf_df(co2_d1s1, lag.max = 48),
  ic       = ic_d1s1,
  title    = expression("PACF --" ~ nabla ~ nabla[12] ~ Y[t]),
  ylab     = "PACF",
  vlines   = c(12, 24, 36),
  vlab     = TRUE
)

p_acf_f / p_pacf_f

ACF e PACF da série estacionária \(\nabla\nabla_{12}Y_t\) com destaque nos lags sazonais (múltiplos de 12). A leitura dessas funções guia a escolha das ordens \(p\), \(q\), \(P\), \(Q\).

A leitura das funções de autocorrelação da série estacionária \(\nabla\nabla_{12}Y_t\) segue o princípio de Brockwell; Davis (2010): a ACF e PACF das defasagens não-sazonais (\(h = 1, 2, \ldots, 11\)) guiam a escolha de \(p\) e \(q\); as defasagens sazonais (\(h = 12, 24, \ldots\)) guiam a escolha de \(P\) e \(Q\).

Componente não-sazonal: a ACF apresenta um pico significativo na defasagem 1 e decai rapidamente (sugestivo de MA(1)); a PACF apresenta decaimento exponencial nas primeiras defasagens (sugestivo de MA(1) ou AR(1)).

Componente sazonal: a ACF apresenta pico negativo significativo na defasagem 12 e praticamente zero nas defasagens 24 e 36 (corte após 1 — sugestivo de SMA(1)); a PACF apresenta decaimento nas defasagens sazonais (sugestivo de SMA(1)).

Candidatos iniciais:

• \(\text{SARIMA}(1, 1, 1)(0, 1, 1)_{12}\)
• \(\text{SARIMA}(0, 1, 1)(0, 1, 1)_{12}\) — modelo “Airline” de Box et al. (2015)
• \(\text{SARIMA}(1, 1, 1)(1, 1, 1)_{12}\)

---

Ajuste de Modelos SARIMA

Definição do Modelo SARIMA

Um processo \(\{Y_t\}\) é \(\text{SARIMA}(p, d, q)(P, D, Q)_s\) se a série diferenciada \[W_t = (1-B)^d(1-B^s)^D Y_t\] é um processo ARMA causador satisfazendo: \[\Phi(B)\tilde{\Phi}(B^s) W_t = \Theta(B)\tilde{\Theta}(B^s)\varepsilon_t, \quad \varepsilon_t \sim \mathcal{RB}(0, \sigma^2),\] onde \(\Phi(z) = 1 - \varphi_1 z - \cdots - \varphi_p z^p\), \(\tilde{\Phi}(z) = 1 - \tilde{\varphi}_1 z - \cdots - \tilde{\varphi}_P z^P\), \(\Theta(z) = 1 + \theta_1 z + \cdots + \theta_q z^q\) e \(\tilde{\Theta}(z) = 1 + \tilde{\theta}_1 z + \cdots + \tilde{\theta}_Q z^Q\) (Brockwell; Davis, 2010, cap. 9).

Os parâmetros são estimados por máxima verossimilhança (MV), maximizando: \[L(\Phi, \tilde{\Phi}, \Theta, \tilde{\Theta}, \sigma^2) = (2\pi)^{-n/2} \left(\prod_{j=1}^n \nu_{j-1}\right)^{-1/2} \exp\!\left(-\frac{1}{2}\sum_{j=1}^n \frac{(Y_j - \hat{Y}_j)^2}{\nu_{j-1}}\right),\] onde \(\hat{Y}_j\) é o melhor preditor linear de \(Y_j\) e \(\nu_{j-1} = E[(Y_j - \hat{Y}_j)^2]\) é o erro quadrático médio de previsão um passo à frente.

Seleção por AIC, AICc e BIC

# CSS-ML: inicializa via mínimos quadrados condicionais, depois refina por MV.
# Mais robusto numericamente que ML puro para certos candidatos.
ajustar_sarima <- function(ordens, serie = co2) {
  tryCatch(
    Arima(serie,
          order    = ordens[1:3],
          seasonal = list(order = ordens[4:6], period = 12),
          method   = "CSS-ML"),
    error = function(e) {
      # fallback para CSS se CSS-ML também falhar
      tryCatch(
        Arima(serie,
              order    = ordens[1:3],
              seasonal = list(order = ordens[4:6], period = 12),
              method   = "CSS"),
        error = function(e2) NULL
      )
    }
  )
}

candidatos <- list(
  "SARIMA(0,1,1)(0,1,1)[12]" = c(0,1,1,0,1,1),
  "SARIMA(1,1,1)(0,1,1)[12]" = c(1,1,1,0,1,1),
  "SARIMA(0,1,1)(1,1,1)[12]" = c(0,1,1,1,1,1),
  "SARIMA(1,1,1)(1,1,1)[12]" = c(1,1,1,1,1,1),
  "SARIMA(2,1,1)(0,1,1)[12]" = c(2,1,1,0,1,1),
  "SARIMA(0,1,2)(0,1,1)[12]" = c(0,1,2,0,1,1),
  "SARIMA(1,1,0)(0,1,1)[12]" = c(1,1,0,0,1,1)
)

fits <- lapply(candidatos, ajustar_sarima)

# Remover modelos que falharam
fits <- Filter(Negate(is.null), fits)

tab_comp <- map_dfr(names(fits), function(nm) {
  f <- fits[[nm]]
  tibble(
    Modelo = nm,
    k      = length(f$coef),
    LogLik = round(f$loglik, 2),
    AIC    = round(f$aic,    2),
    AICc   = round(f$aicc,   2),
    BIC    = round(f$bic,    2)
  )
}) |>
  arrange(AICc) |>
  mutate(
    dAIC  = round(AIC  - min(AIC),  2),
    dAICc = round(AICc - min(AICc), 2),
    dBIC  = round(BIC  - min(BIC),  2)
  )

tab_comp |>
  gt() |>
  tab_header(
    title    = md("**Tabela 2.** Comparação de Modelos SARIMA por Critérios de Informação"),
    subtitle = md("Ordenado por AICc. $\\Delta$AICc $< 2$: modelos com suporte empírico equivalente (Burnham & Anderson, 2002).")
  ) |>
  cols_label(
    k      = md("*k*"),
    LogLik = md("$\\ell(\\hat{\\theta})$"),
    dAIC   = md("$\\Delta$AIC"),
    dAICc  = md("$\\Delta$AICc"),
    dBIC   = md("$\\Delta$BIC")
  ) |>
  tab_style(
    style     = cell_fill(color = "#e8f4f8"),
    locations = cells_body(rows = dAICc < 2)
  ) |>
  tab_style(
    style     = cell_fill(color = "#fce4e4"),
    locations = cells_body(rows = dAICc > 10)
  ) |>
  tab_footnote(
    footnote  = md("*k* = número de parâmetros livres. Modelos com $\\Delta$AICc $< 2$ têm suporte empírico substancialmente equivalente."),
    locations = cells_column_labels(columns = k)
  )

Tabela 2. Comparação de Modelos SARIMA por Critérios de Informação
Ordenado por AICc. \(\Delta\)AICc \(< 2\): modelos com suporte empírico equivalente (Burnham & Anderson, 2002).
Modelo	k1	\(\ell(\hat{\theta})\)	AIC	AICc	BIC	\(\Delta\)AIC	\(\Delta\)AICc	\(\Delta\)BIC
SARIMA(2,1,1)(0,1,1)[12]	4	-83.91	177.83	177.96	198.43	0.00	0.00	7.91
SARIMA(1,1,1)(0,1,1)[12]	3	-85.03	178.07	178.16	194.55	0.24	0.20	4.03
SARIMA(0,1,1)(0,1,1)[12]	2	-86.08	178.16	178.21	190.52	0.33	0.25	0.00
SARIMA(0,1,2)(0,1,1)[12]	3	-85.55	179.09	179.18	195.57	1.26	1.22	5.05
SARIMA(1,1,1)(1,1,1)[12]	4	-84.88	179.76	179.89	200.36	1.93	1.93	9.84
SARIMA(0,1,1)(1,1,1)[12]	3	-85.96	179.92	180.01	196.40	2.09	2.05	5.88
SARIMA(1,1,0)(0,1,1)[12]	2	-89.09	184.18	184.23	196.54	6.35	6.27	6.02
1 k = número de parâmetros livres. Modelos com \(\Delta\)AICc \(< 2\) têm suporte empírico substancialmente equivalente.

aa <- auto.arima(co2, stepwise = FALSE, approximation = FALSE,
                 max.p = 3, max.q = 3, max.P = 2, max.Q = 2,
                 d = 1, D = 1, ic = "aicc")

cat("Modelo selecionado por auto.arima:\n")

Modelo selecionado por auto.arima:

print(aa)

Series: co2 
ARIMA(0,1,3)(0,1,1)[12] 

Coefficients:
          ma1      ma2      ma3     sma1
      -0.3394  -0.0180  -0.0973  -0.8538
s.e.   0.0475   0.0497   0.0467   0.0256

sigma^2 = 0.08518:  log likelihood = -83.43
AIC=176.86   AICc=177   BIC=197.47

cat(sprintf("\nAIC  = %.2f | AICc = %.2f | BIC = %.2f\n",
            aa$aic, aa$aicc, aa$bic))

AIC  = 176.86 | AICc = 177.00 | BIC = 197.47

A Tabela 2 mostra que os três primeiros modelos apresentam ΔAICc inferior a 0.3, configurando um empate técnico segundo o critério de Burnham; Anderson (2002): modelos com diferença de AICc menor que 2 têm suporte empírico substancialmente equivalente. Nesse cenário, o SARIMA(2,1,1)(0,1,1)₁₂ apresenta o menor AICc absoluto (177.96), seguido pelo (1,1,1)(0,1,1)₁₂ (AICc = 178.16, ΔAICc = 0.20) e pelo (0,1,1)(0,1,1)₁₂ (AICc = 178.21, ΔAICc = 0.25).

Diante do empate técnico (ΔAICc < 0.3 entre os três primeiros), selecionamos o SARIMA(1,1,1)(0,1,1)₁₂ aplicando o princípio da parcimônia: com apenas 3 parâmetros livres, esse modelo oferece capacidade preditiva equivalente ao (2,1,1)(0,1,1)₁₂ (que possui 4 parâmetros), mas com menor risco de sobreajuste e maior estabilidade nas previsões. O auto.arima com busca exaustiva sugere o SARIMA(0,1,3)(0,1,1)₁₂ (AICc = 177.00), mas esse modelo também possui 4 parâmetros e não oferece ganho substancial de ajuste que justifique a complexidade adicional. Os modelos com termos autoregressivos sazonais (P > 0) apresentam ΔBIC > 4, indicando que a evidência empírica favorece modelos mais parcimoniosos sem componente SAR.

Modelo Final: \(\text{SARIMA}(1,1,1)(0,1,1)_{12}\)

fit_final <- fits[["SARIMA(1,1,1)(0,1,1)[12]"]]

coef_df <- tidy(fit_final) |>
  rename(Parametro  = term,
         Estimativa = estimate,
         EP         = std.error) |>
  mutate(
    IC_inf    = round(Estimativa - 1.96 * EP, 4),
    IC_sup    = round(Estimativa + 1.96 * EP, 4),
    z         = round(Estimativa / EP, 3),
    p_valor   = round(2 * pnorm(-abs(z)), 5),
    Parametro = recode(Parametro,
      "ar1"  = "phi[1] (AR nao-sazonal)",
      "ma1"  = "theta[1] (MA nao-sazonal)",
      "sma1" = "theta_tilde[1] (SMA sazonal, lag 12)"
    )
  ) |>
  select(Parametro, Estimativa, EP, IC_inf, IC_sup, z, p_valor)

coef_df |>
  gt() |>
  fmt_number(columns = c(Estimativa, EP, IC_inf, IC_sup, z), decimals = 4) |>
  fmt_scientific(columns = p_valor, decimals = 2) |>
  tab_header(
    title    = md("**Tabela 3.** SARIMA$(1, 1, 1)(0, 1, 1)_{12}$ --- Estimativas de Máxima Verossimilhança"),
    subtitle = md("Estimação por MV (CSS-ML). ICs assintóticos: $\\hat{\\beta} \\pm 1{,}96 \\cdot \\text{EP}$.")
  ) |>
  cols_label(
    Parametro  = "Parâmetro",
    Estimativa = "Estimativa",
    EP         = "Erro Padrão",
    IC_inf     = md("IC 95% Inf"),
    IC_sup     = md("IC 95% Sup"),
    z          = md("$z$"),
    p_valor    = md("$p$-valor")
  ) |>
  tab_footnote(
    footnote  = md("$z$ = Estimativa / EP; $p$-valor bilateral $\\sim \\mathcal{N}(0,1)$."),
    locations = cells_column_labels(columns = z)
  ) |>
  tab_style(
    style     = cell_text(weight = "bold"),
    locations = cells_body(columns = Parametro)
  )

Tabela 3. SARIMA\((1, 1, 1)(0, 1, 1)_{12}\) — Estimativas de Máxima Verossimilhança
Estimação por MV (CSS-ML). ICs assintóticos: \(\hat{\beta} \pm 1{,}96 \cdot \text{EP}\).
Parâmetro	Estimativa	Erro Padrão	IC 95% Inf	IC 95% Sup	\(z\)1	\(p\)-valor
phi[1] (AR nao-sazonal)	0.2399	0.1430	−0.0405	0.5202	1.6770	9.35 × 10−2
theta[1] (MA nao-sazonal)	−0.5710	0.1237	−0.8134	−0.3286	−4.6170	0.00
theta_tilde[1] (SMA sazonal, lag 12)	−0.8516	0.0256	−0.9017	−0.8014	−33.2870	0.00
1 \(z\) = Estimativa / EP; \(p\)-valor bilateral \(\sim \mathcal{N}(0,1)\).

O modelo ajustado pode ser escrito explicitamente como: \[\underbrace{(1 - \hat{\varphi}_1 B)}_{\text{AR}(1)} \underbrace{(1 - B)(1 - B^{12})}_{\text{diferenciações}} Y_t = \underbrace{(1 + \hat{\theta}_1 B)}_{\text{MA}(1)} \underbrace{(1 + \hat{\tilde{\theta}}_1 B^{12})}_{\text{SMA}(1)} \varepsilon_t, \quad \varepsilon_t \sim \mathcal{RB}(0, \hat{\sigma}^2).\]

Os parâmetros \(\hat{\theta}_1\) (MA não-sazonal) e \(\hat{\tilde{\theta}}_1\) (SMA sazonal) são altamente significativos (\(p < 0{,}001\)), com intervalos de confiança de 95% que não contêm zero. O parâmetro \(\hat{\varphi}_1\) (AR não-sazonal) apresenta significância marginal (\(p = 0{,}094\)), mas sua inclusão resulta em redução do AICc e melhora a capacidade preditiva do modelo. O parâmetro SMA (\(\hat{\tilde{\theta}}_1 \approx -0{,}85\)) captura a dependência negativa entre as inovações sazonais sucessivas, efeito característico de séries com sazonalidade estocástica estável (Box et al., 2015).

Verificação de Invertibilidade e Causalidade

phi   <- c(1, -fit_final$coef["ar1"])
theta <- c(1,  fit_final$coef["ma1"])

raiz_ar <- polyroot(phi)
raiz_ma <- polyroot(theta)

tibble(
  Componente  = c("AR nao-sazonal", "MA nao-sazonal"),
  Raiz_Re     = round(c(Re(raiz_ar), Re(raiz_ma)), 4),
  Raiz_Im     = round(c(Im(raiz_ar), Im(raiz_ma)), 4),
  Mod_raiz    = round(c(Mod(raiz_ar), Mod(raiz_ma)), 4),
  Fora_UC     = c(Mod(raiz_ar) > 1, Mod(raiz_ma) > 1),
  Propriedade = c("Causalidade", "Invertibilidade")
) |>
  gt() |>
  tab_header(
    title    = md("**Tabela 4.** Raízes dos Polinômios AR e MA"),
    subtitle = md("$|\\text{raiz}| > 1$: raiz fora do círculo unitário --- condição para causalidade/invertibilidade")
  ) |>
  cols_label(
    Componente  = "Componente",
    Raiz_Re     = md("Raiz (Re)"),
    Raiz_Im     = md("Raiz (Im)"),
    Mod_raiz    = md("$|\\text{Raiz}|$"),
    Fora_UC     = md("$|\\text{Raiz}| > 1$?"),
    Propriedade = "Propriedade"
  ) |>
  tab_style(
    style     = cell_fill(color = "#e8f4f8"),
    locations = cells_body(rows = Fora_UC == TRUE)
  )

Tabela 4. Raízes dos Polinômios AR e MA
\(|\text{raiz}| > 1\): raiz fora do círculo unitário — condição para causalidade/invertibilidade
Componente	Raiz (Re)	Raiz (Im)	\(|\text{Raiz}|\)	\(|\text{Raiz}| > 1\)?	Propriedade
AR nao-sazonal	4.1686	0	4.1686	TRUE	Causalidade
MA nao-sazonal	1.7514	0	1.7514	TRUE	Invertibilidade

Ambas as raízes estão fora do círculo unitário (\(|\text{raiz}| > 1\)), garantindo que o modelo é causal e invertível (Brockwell; Davis, 2010, cap. 3).

---

Diagnóstico do Modelo Ajustado

Análise dos Resíduos

Os resíduos do modelo ajustado devem se comportar como um ruído branco \(\mathcal{RB}(0, \sigma^2)\): média zero, variância constante e sem autocorrelação.

res       <- residuals(fit_final)
res_clean <- na.omit(as.numeric(res))
n_res     <- length(res_clean)

tibble(t = as.numeric(time(res)), r = as.numeric(res)) |>
  filter(!is.na(r)) |>
  ggplot(aes(x = t, y = r)) +
  geom_line(color = cores[1], linewidth = 0.5) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  geom_hline(yintercept = c(-2, 2) * sd(res_clean),
             linetype = "dotted", color = cores[2], linewidth = 0.6) +
  scale_x_continuous(breaks = seq(1960, 1997, by = 5)) +
  labs(title    = "Resíduos do Modelo Final",
       subtitle = "Linhas pontilhadas: ±2 desvios padrão",
       x = NULL, y = "Resíduo (ppm)")

Resíduos do modelo SARIMA\((1,1,1)(0,1,1)_{12}\). Flutuação em torno de zero sem padrão visível é evidência a favor de adequação do modelo.

ic_res <- 1.96 / sqrt(n_res)

p_acf_res <- plot_acf(
  acf_df(res_clean, lag.max = 48) |> filter(lag > 0),
  ic    = ic_res,
  title = "ACF dos Resíduos"
)

p_pacf_res <- plot_acf(
  pacf_df(res_clean, lag.max = 48),
  ic    = ic_res,
  title = "PACF dos Resíduos",
  ylab  = "PACF"
)

p_acf_res + p_pacf_res

ACF e PACF dos resíduos. Ausência de autocorrelação significativa indica que o modelo capturou toda a estrutura de dependência linear da série.

A ACF e a PACF dos resíduos não apresentam autocorrelações significativas em nenhuma defasagem (dentro das bandas de confiança de 95%: \(\pm 1{,}96/\sqrt{n}\)). O modelo capturou a estrutura de dependência linear da série. Visualmente, os resíduos flutuam em torno de zero com variância aproximadamente constante ao longo do tempo, com exceção de um valor atípico isolado em torno de 1983-84 (período associado ao El Niño), sem comprometer a adequação geral do modelo.

Normalidade dos Resíduos

p_hist <- tibble(r = res_clean) |>
  ggplot(aes(x = r)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 30, fill = cores[1],
                 color = "white", alpha = 0.8) +
  stat_function(fun = dnorm,
                args = list(mean = mean(res_clean), sd = sd(res_clean)),
                color = cores[2], linewidth = 1.0) +
  labs(title = "Histograma dos Resíduos",
       x = "Resíduo (ppm)", y = "Densidade")

p_qq <- tibble(r = res_clean) |>
  ggplot(aes(sample = r)) +
  stat_qq(color = cores[1], alpha = 0.6, size = 0.9) +
  stat_qq_line(color = cores[2], linewidth = 0.9) +
  labs(title = "Q-Q Normal dos Resíduos",
       x = "Quantis teóricos", y = "Quantis amostrais")

p_hist + p_qq

Diagnóstico de normalidade dos resíduos: histograma com curva normal sobreposta (esq.) e gráfico Q-Q (dir.).

sw <- shapiro.test(res_clean)

mu3     <- mean((res_clean - mean(res_clean))^3)
mu4     <- mean((res_clean - mean(res_clean))^4)
sk      <- mu3 / sd(res_clean)^3
ku      <- mu4 / sd(res_clean)^4 - 3
jb_stat <- (n_res / 6) * (sk^2 + ku^2 / 4)
jb_p    <- pchisq(jb_stat, df = 2, lower.tail = FALSE)

tibble(
  Teste        = c("Shapiro-Wilk", "Jarque-Bera"),
  Estatistica  = round(c(sw$statistic, jb_stat), 4),
  p_valor      = round(c(sw$p.value, jb_p), 5),
  H0_rejeitada = c(sw$p.value < 0.05, jb_p < 0.05)
) |>
  gt() |>
  tab_header(
    title    = md("**Tabela 5.** Testes de Normalidade dos Resíduos"),
    subtitle = md("$H_0$: resíduos seguem distribuição normal")
  ) |>
  fmt_scientific(columns = p_valor, decimals = 2) |>
  cols_label(
    Estatistica  = "Estatística",
    p_valor      = md("$p$-valor"),
    H0_rejeitada = md("$H_0$ rejeitada (5%)?")
  ) |>
  tab_style(
    style     = cell_fill(color = "#fce4e4"),
    locations = cells_body(rows = H0_rejeitada == TRUE)
  )

Tabela 5. Testes de Normalidade dos Resíduos
\(H_0\): resíduos seguem distribuição normal
Teste	Estatística	\(p\)-valor	\(H_0\) rejeitada (5%)?
Shapiro-Wilk	0.9969	5.32 × 10−1	FALSE
Jarque-Bera	1.9176	3.83 × 10−1	FALSE

Os testes formais não rejeitam a hipótese de normalidade dos resíduos (Shapiro-Wilk: \(p = 0{,}532\); Jarque-Bera: \(p = 0{,}383\)), indicando que a distribuição dos resíduos é consistente com a normal. Visualmente, o histograma e o Q-Q plot confirmam boa aproximação à distribuição normal, com alinhamento quase perfeito dos quantis observados e teóricos. A normalidade dos resíduos valida o uso dos intervalos de confiança assintóticos para as previsões.

Teste de Ljung-Box (Ausência de Autocorrelação Residual)

k_params <- length(fit_final$coef)

lb_res <- map_dfr(c(6, 12, 18, 24, 36, 48), function(h) {
  lb <- Box.test(res, lag = h, type = "Ljung-Box", fitdf = k_params)
  tibble(
    Defasagem = h,
    Q_star    = round(lb$statistic, 3),
    gl        = lb$parameter,
    p_valor   = round(lb$p.value, 5)
  )
})

lb_res |>
  gt() |>
  tab_header(
    title    = md("**Tabela 6.** Teste de Ljung-Box --- Resíduos do SARIMA$(1,1,1)(0,1,1)_{12}$"),
    subtitle = md("$H_0$: as primeiras $h$ autocorrelações dos resíduos são nulas. g.l. $= h - k$.")
  ) |>
  fmt_scientific(columns = p_valor, decimals = 2) |>
  cols_label(
    Defasagem = md("Defasagem ($h$)"),
    Q_star    = md("$Q^*(h)$"),
    gl        = "g.l.",
    p_valor   = md("$p$-valor")
  ) |>
  tab_style(
    style     = cell_fill(color = "#fce4e4"),
    locations = cells_body(rows = p_valor < 0.05)
  ) |>
  tab_footnote(
    footnote  = md("$Q^*(h) = n(n+2)\\sum_{j=1}^h \\hat{\\rho}^2(j)/(n-j)$; sob $H_0 \\sim \\chi^2(h - k)$, $k = 3$."),
    locations = cells_column_labels(columns = Q_star)
  )

Tabela 6. Teste de Ljung-Box — Resíduos do SARIMA\((1,1,1)(0,1,1)_{12}\)
\(H_0\): as primeiras \(h\) autocorrelações dos resíduos são nulas. g.l. \(= h - k\).
Defasagem (\(h\))	\(Q^*(h)\)1	g.l.	\(p\)-valor
6	5.288	3	1.52 × 10−1
12	11.621	9	2.36 × 10−1
18	12.568	15	6.36 × 10−1
24	20.832	21	4.69 × 10−1
36	34.423	33	4.00 × 10−1
48	46.592	45	4.07 × 10−1
1 \(Q^*(h) = n(n+2)\sum_{j=1}^h \hat{\rho}^2(j)/(n-j)\); sob \(H_0 \sim \chi^2(h - k)\), \(k = 3\).

O teste de Ljung-Box não rejeita \(H_0\) em nenhuma das defasagens testadas (\(p > 0{,}10\) em todos os casos), confirmando que os resíduos não apresentam autocorrelação serial detectável — validando o ajuste do modelo.

Sumário dos Diagnósticos

p_lb12 <- round(Box.test(res, lag = 12, type = "Ljung-Box",
                          fitdf = k_params)$p.value, 4)

tibble(
  Diagnostico = c("Média dos resíduos",
                  "Desvio padrão estimado",
                  "Assimetria amostral",
                  "Excesso de curtose",
                  "Ljung-Box (h = 12)",
                  "Shapiro-Wilk"),
  Valor       = c(format(round(mean(res_clean), 5), nsmall = 5),
                  format(round(sd(res_clean),   4), nsmall = 4),
                  as.character(round(sk, 3)),
                  as.character(round(ku, 3)),
                  sprintf("p = %.4f", p_lb12),
                  sprintf("p = %.4f", sw$p.value)),
  Conclusao   = c("Média aprox. zero",
                  "---",
                  if (abs(sk) < 0.5) "Assimetria baixa" else "Assimetria moderada",
                  if (abs(ku) < 1)   "Curtose próx. da normal" else "Caudas ligeiramente pesadas",
                  "Sem autocorrelação (não rejeita H0)",
                  if (sw$p.value > 0.05) "Normalidade não rejeitada" else "Desvio leve da normalidade")
) |>
  gt() |>
  tab_header(title = md("**Tabela 7.** Sumário dos Diagnósticos do Modelo Final")) |>
  cols_label(Diagnostico = "Diagnóstico", Conclusao = "Conclusão") |>
  tab_style(
    style     = cell_text(weight = "bold"),
    locations = cells_body(columns = Diagnostico)
  )

Tabela 7. Sumário dos Diagnósticos do Modelo Final
Diagnóstico	Valor	Conclusão
Média dos resíduos	0.01721	Média aprox. zero
Desvio padrão estimado	0.2873	---
Assimetria amostral	-0.151	Assimetria baixa
Excesso de curtose	0.084	Curtose próx. da normal
Ljung-Box (h = 12)	p = 0.2355	Sem autocorrelação (não rejeita H0)
Shapiro-Wilk	p = 0.5316	Normalidade não rejeitada

---

Previsão

Previsões para os Próximos 24 Meses

O melhor previsor linear de \(Y_{n+h}\) dado \(\{Y_1, \ldots, Y_n\}\), denotado \(P_n Y_{n+h}\), é obtido do algoritmo das inovações e minimiza o erro quadrático médio de previsão (Brockwell; Davis, 2010, cap. 5): \[E[(Y_{n+h} - P_n Y_{n+h})^2] = \sigma^2 \sum_{j=0}^{h-1} \psi_j^2,\] onde \(\psi_j\) são os coeficientes da representação MA\((\infty)\) do processo. Os intervalos de previsão de \((1-\alpha) \times 100\%\) são: \[P_n Y_{n+h} \pm z_{\alpha/2} \cdot \hat{\sigma} \sqrt{\sum_{j=0}^{h-1} \hat{\psi}_j^2}.\]

h    <- 24
prev <- forecast(fit_final, h = h, level = c(80, 95))

# Horizonte temporal: sequência de anos decimais após o fim da amostra
t_prev <- seq(fim[1] + fim[2] / 12, by = 1/12, length.out = h)

prev_df <- tibble(
  tempo   = t_prev,
  central = as.numeric(prev$mean),
  lo80    = as.numeric(prev$lower[, 1]),
  hi80    = as.numeric(prev$upper[, 1]),
  lo95    = as.numeric(prev$lower[, 2]),
  hi95    = as.numeric(prev$upper[, 2])
)

obs_recente <- co2_df |> filter(ano >= 1993)

ggplot() +
  geom_line(data = obs_recente,
            aes(x = tempo, y = ppm),
            color = cores[1], linewidth = 0.8) +
  geom_ribbon(data = prev_df,
              aes(x = tempo, ymin = lo95, ymax = hi95),
              fill = cores[2], alpha = 0.15) +
  geom_ribbon(data = prev_df,
              aes(x = tempo, ymin = lo80, ymax = hi80),
              fill = cores[2], alpha = 0.25) +
  geom_line(data = prev_df,
            aes(x = tempo, y = central),
            color = cores[2], linewidth = 1.0) +
  geom_vline(xintercept = 1998,
             linetype = "dashed", color = "gray40", linewidth = 0.7) +
  annotate("text", x = 1997.65, y = 358,
           label = "Fim obs.", angle = 90, size = 3, color = "gray40") +
  scale_x_continuous(breaks = seq(1993, 2000, by = 1)) +
  scale_y_continuous(labels = label_number(suffix = " ppm")) +
  labs(
    title    = expression("Previsão SARIMA" * (1 * "," * 1 * "," * 1) *
                          (0 * "," * 1 * "," * 1)[12] * " -- 24 meses"),
    subtitle = "Faixa escura: IC 80% | Faixa clara: IC 95% | Linha vermelha: previsão pontual",
    x        = NULL,
    y        = expression(CO[2] ~ "(ppm)")
  )

Previsões SARIMA para 24 meses (1998–1999) com intervalos de 80% e 95%. A série observada de 1993 a 1997 é mostrada à esquerda da linha de corte.

meses_pt <- c("Janeiro","Fevereiro","Março","Abril","Maio","Junho",
              "Julho","Agosto","Setembro","Outubro","Novembro","Dezembro")

prev_df |>
  mutate(
    Ano = as.integer(floor(tempo)),
    Mes = meses_pt[as.integer(round((tempo - floor(tempo)) * 12)) + 1L],
    .before = central
  ) |>
  select(-tempo) |>
  gt(groupname_col = "Ano") |>
  tab_header(
    title    = md("**Tabela 8.** Previsões Mensais de $\\text{CO}_2$ (ppm) --- 1998 e 1999"),
    subtitle = md("Modelo SARIMA$(1,1,1)(0,1,1)_{12}$")
  ) |>
  fmt_number(columns = c(central, lo80, hi80, lo95, hi95), decimals = 2) |>
  cols_label(
    Mes     = "Mês",
    central = "Prev. Pontual",
    lo80    = "IC 80% Inf",
    hi80    = "IC 80% Sup",
    lo95    = "IC 95% Inf",
    hi95    = "IC 95% Sup"
  ) |>
  tab_style(
    style     = cell_fill(color = "#f0f7ff"),
    locations = cells_row_groups()
  ) |>
  tab_footnote(
    footnote  = md("ICs: $\\hat{Y}_{n+h} \\pm z_{\\alpha/2} \\cdot \\hat{\\sigma}\\sqrt{\\sum_{j=0}^{h-1}\\hat{\\psi}_j^2}$."),
    locations = cells_column_labels(columns = lo95)
  )

Tabela 8. Previsões Mensais de \(\text{CO}_2\) (ppm) — 1998 e 1999
Modelo SARIMA\((1,1,1)(0,1,1)_{12}\)
Mês	Prev. Pontual	IC 80% Inf	IC 80% Sup	IC 95% Inf1	IC 95% Sup
1998
Janeiro	365.18	364.81	365.56	364.61	365.75
Fevereiro	365.97	365.52	366.42	365.28	366.66
Março	366.82	366.31	367.32	366.05	367.58
Abril	368.16	367.62	368.71	367.33	369.00
Maio	368.73	368.14	369.32	367.83	369.63
Junho	368.04	367.42	368.67	367.09	369.00
Julho	366.54	365.88	367.20	365.53	367.55
Agosto	364.49	363.79	365.18	363.43	365.54
Setembro	362.63	361.90	363.35	361.52	363.73
Outubro	362.75	362.00	363.51	361.60	363.91
Novembro	364.19	363.40	364.97	362.99	365.38
Dezembro	365.60	364.79	366.41	364.36	366.84
1999
Janeiro	366.65	365.80	367.51	365.35	367.96
Fevereiro	367.49	366.60	368.38	366.13	368.85
Março	368.35	367.43	369.27	366.94	369.76
Abril	369.70	368.75	370.66	368.24	371.16
Maio	370.27	369.29	371.25	368.77	371.78
Junho	369.58	368.57	370.60	368.03	371.13
Julho	368.08	367.04	369.12	366.49	369.68
Agosto	366.03	364.96	367.10	364.39	367.66
Setembro	364.17	363.07	365.27	362.49	365.85
Outubro	364.30	363.17	365.42	362.58	366.01
Novembro	365.73	364.58	366.88	363.97	367.49
Dezembro	367.14	365.96	368.32	365.34	368.94
1 ICs: \(\hat{Y}_{n+h} \pm z_{\alpha/2} \cdot \hat{\sigma}\sqrt{\sum_{j=0}^{h-1}\hat{\psi}_j^2}\).

Análise da Incerteza de Previsão

tibble(
  h       = 1:h,
  largura = prev_df$hi95 - prev_df$lo95
) |>
  ggplot(aes(x = h, y = largura)) +
  geom_line(color = cores[1], linewidth = 0.9) +
  geom_point(color = cores[1], size = 2) +
  scale_x_continuous(breaks = seq(2, 24, by = 2)) +
  labs(
    title    = "Largura do IC 95% de Previsão por Horizonte",
    subtitle = "Incerteza acumulada: de aprox. 0,5 ppm (h = 1) a aprox. 2,0 ppm (h = 24)",
    x        = "Horizonte h (meses)",
    y        = "Largura do IC 95% (ppm)"
  )

Largura dos intervalos de previsão (IC 95%) em função do horizonte \(h\). O crescimento reflete o acúmulo de incerteza com a distância no tempo.

A largura do intervalo de previsão cresce com o horizonte \(h\), reflexo do acúmulo de incerteza inerente à previsão de longo prazo. Para \(h=1\) mês, o IC 95% tem largura de aproximadamente \(\pm 0{,}5\) ppm; para \(h=24\) meses, a incerteza acumula-se a cerca de \(\pm 2{,}0\) ppm. Em contexto de valores de \(\text{CO}_2\) na faixa de 360–370 ppm, esse nível de incerteza representa um erro relativo inferior a 0,6%, demonstrando a excelente capacidade preditiva do modelo.

Modelo de Regressão com Erros ARMA (Sinal + Ruído)

Como alternativa ao SARIMA, o modelo de regressão com erros ARMA (Brockwell; Davis, 2010, cap. 6) decompõe explicitamente a série em um componente determinístico (sinal) e um componente estocástico (ruído): \[Y_t = \beta_0 + \beta_1 t + \sum_{k=1}^{K} \left[\alpha_k \cos\!\left(\frac{2\pi k t}{12}\right) + \delta_k \sin\!\left(\frac{2\pi k t}{12}\right)\right] + W_t, \quad W_t \sim \text{ARMA}(p, q),\] onde os termos trigonométricos capturam a sazonalidade determinística (usamos \(K=2\) harmônicos fundamentais para parcimônia) e \(W_t\) é um processo ARMA estacionário com média zero.

t_idx <- 1:length(co2)

# Matriz de regressores: tendência linear + harmônicos fundamentais (k=1,2)
# Usamos apenas os 2 primeiros harmônicos para evitar sobreparametrização
fourier_terms <- cbind(
  cos1 = cos(2 * pi * 1 * t_idx / 12),
  sin1 = sin(2 * pi * 1 * t_idx / 12),
  cos2 = cos(2 * pi * 2 * t_idx / 12),
  sin2 = sin(2 * pi * 2 * t_idx / 12)
)
xreg <- cbind(trend = t_idx, fourier_terms)

# Ajuste direto: ARMA(1,1) nos resíduos da regressão
fit_reg <- Arima(co2, order = c(1, 0, 1), xreg = xreg, method = "CSS-ML")

cat("Modelo de regressão com erros ARMA(1,1):\n")

Modelo de regressão com erros ARMA(1,1):

print(fit_reg)

Series: co2 
Regression with ARIMA(1,0,1) errors 

Coefficients:
         ar1      ma1  intercept   trend     cos1    sin1    cos2     sin2
      0.9946  -0.3378   313.3255  0.1068  -1.7198  2.1957  0.7707  -0.0042
s.e.  0.0044   0.0514     2.0433  0.0059   0.0294  0.0295  0.0185   0.0185

sigma^2 = 0.1039:  log likelihood = -132.12
AIC=282.23   AICc=282.63   BIC=319.57

cat(sprintf("\nAIC  = %.2f | AICc = %.2f | BIC = %.2f\n",
            fit_reg$aic, fit_reg$aicc, fit_reg$bic))

AIC  = 282.23 | AICc = 282.63 | BIC = 319.57

tibble(
  Modelo     = c("SARIMA$(1,1,1)(0,1,1)_{12}$",
                  "Regressão + ARMA(1,1) (sinal+ruído)"),
  Parametros = c(length(fit_final$coef), length(fit_reg$coef)),
  AIC        = round(c(fit_final$aic,    fit_reg$aic),    2),
  AICc       = round(c(fit_final$aicc,   fit_reg$aicc),   2),
  BIC        = round(c(fit_final$bic,    fit_reg$bic),    2),
  sigma2     = round(c(fit_final$sigma2, fit_reg$sigma2), 5)
) |>
  gt() |>
  tab_header(
    title    = md("**Tabela 9.** Comparação: SARIMA vs. Regressão com Erros ARMA"),
    subtitle = "Ambas as abordagens modelam tendência e sazonalidade com suposições distintas."
  ) |>
  fmt_markdown(columns = Modelo) |>
  cols_label(
    Parametros = "Parâmetros",
    sigma2     = md("$\\hat{\\sigma}^2$")
  ) |>
  tab_style(
    style     = cell_fill(color = "#e8f4f8"),
    locations = cells_body(rows = AICc == min(AICc))
  )

Tabela 9. Comparação: SARIMA vs. Regressão com Erros ARMA
Ambas as abordagens modelam tendência e sazonalidade com suposições distintas.
Modelo	Parâmetros	AIC	AICc	BIC	\(\hat{\sigma}^2\)
SARIMA\((1,1,1)(0,1,1)_{12}\)	3	178.07	178.16	194.55	0.08560
Regressão + ARMA(1,1) (sinal+ruído)	8	282.23	282.63	319.57	0.10391

O SARIMA\((1,1,1)(0,1,1)_{12}\) apresenta AICc inferior ao modelo de regressão com erros ARMA(1,1), sendo preferido pelo critério de seleção. Além disso, o SARIMA é mais parcimonioso e trata a sazonalidade como estocástica, o que é mais realista para séries longas nas quais o padrão sazonal pode variar ligeiramente entre anos. A abordagem de regressão tem valor quando o interesse recai sobre a interpretação direta de parâmetros como a taxa de crescimento anual (via \(\hat{\beta}_1\)) ou a amplitude e fase dos harmônicos sazonais.

---

Resultados

Síntese da Análise

Análise descritiva: A série cresce a uma taxa média de \(\sim\) 1.23 ppm/ano, com sazonalidade anual de amplitude \(\sim 6\) ppm, pico em maio e vale em setembro. A decomposição STL confirmou que tendência e sazonalidade são os componentes dominantes.

Identificação e seleção de modelos: A combinação de diferenciações \(d=1\) (remove tendência) e \(D=1\) (remove sazonalidade estocástica) com período \(s=12\) foi necessária e suficiente para induzir estacionariedade, confirmado pelos testes ADF e KPSS. A comparação por AICc revelou três modelos com suporte empírico equivalente (\(\Delta\)AICc \(< 0{,}3\)). Selecionamos o \(\text{SARIMA}(1,1,1)(0,1,1)_{12}\) por equilibrar parcimônia e capacidade preditiva, com 3 parâmetros: \(\hat{\varphi}_1 = 0{,}24\) (\(p = 0{,}094\)), \(\hat{\theta}_1 = -0{,}57\) (\(p < 0{,}001\)) e \(\hat{\tilde{\theta}}_1 = -0{,}85\) (\(p < 0{,}001\)).

Diagnóstico: Os resíduos não apresentam autocorrelação detectável (Ljung-Box, \(p > 0{,}10\) para \(h \leq 48\)), têm média próxima de zero e distribuição aproximadamente normal. O modelo é causal e invertível.

Previsão: Para dezembro de 1998, o modelo prevê 365.6 ppm; para dezembro de 1999, 367.14 ppm (IC 95%: 365.34–368.94 ppm). O crescimento monotônico e a sazonalidade anual são fielmente reproduzidos nas previsões, com incerteza crescente mas contida (largura máxima de \(\sim \pm 2\) ppm em \(h=24\)).

---

Conclusão

Este trabalho demonstrou que a série de \(\text{CO}_2\) de Mauna Loa pode ser descrita de forma precisa e parcimoniosa pelo modelo \(\text{SARIMA}(1,1,1)(0,1,1)_{12}\). Entre vários candidatos com desempenho equivalente (\(\Delta\)AICc \(< 0{,}3\)), esse modelo foi selecionado por equilibrar parcimônia (3 parâmetros) e capacidade preditiva. O modelo captura a tendência de longo prazo via diferenciação simples, a sazonalidade anual via diferenciação sazonal com período 12, e a dependência residual de curto prazo via os componentes AR(1) não-sazonal (\(\hat{\varphi}_1 = 0{,}24\)), MA(1) não-sazonal (\(\hat{\theta}_1 = -0{,}57\)) e SMA(1) sazonal (\(\hat{\tilde{\theta}}_1 = -0{,}85\)).

O diagnóstico confirma a adequação do ajuste: resíduos sem autocorrelação detectável, média nula e distribuição próxima da normal. As previsões para 1998 e 1999 preservam o padrão esperado de crescimento monotônico e sazonalidade regular, com intervalos de incerteza estreitos que refletem a alta regularidade da série.

A comparação com o modelo de regressão com erros ARMA mostrou que ambas as abordagens são válidas, mas o SARIMA é preferido por critérios de informação e parcimônia. Em contextos onde a interpretação explícita da taxa de crescimento anual ou da amplitude sazonal for de interesse, o modelo de regressão oferece uma leitura mais direta desses parâmetros.

Do ponto de vista substantivo, a análise quantifica a tendência de aumento do \(\text{CO}_2\) atmosférico ao longo de quase quatro décadas: um crescimento de \(\sim\) 1.23 ppm/ano que reforça a urgência das discussões sobre mudanças climáticas (Keeling et al., 1976).

---

Referências

BOX, George E. P. et al. Time Series Analysis: Forecasting and Control. 5th. ed. Hoboken, NJ: Wiley, 2015.

BROCKWELL, Peter J.; DAVIS, Richard A. Introduction to Time Series and Forecasting. 2nd. ed. New York: Springer, 2010.

BURNHAM, Kenneth P.; ANDERSON, David R. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach. 2nd. ed. New York: Springer, 2002.

CLEVELAND, Robert B. et al. STL: A seasonal-trend decomposition procedure based on LOESS. Journal of Official Statistics, v. 6, n. 1, p. 3–73, 1990.

KEELING, Charles D. et al. Atmospheric carbon dioxide variations at Mauna Loa Observatory, Hawaii. Tellus, v. 28, n. 6, p. 538–551, 1976.

MORETTIN, Pedro A.; TOLOI, Clélia M. C. Análise de Séries Temporais. São Paulo: Editora Edgar Blücher, 2004.

---

Apêndice: Códigos

Todos os códigos utilizados neste trabalho estão disponíveis ao longo do documento com opção de expansão (code-fold: true). Para visualizá-los, clique em “Code” ao lado de cada bloco, ou use o botão “</> Code” no topo da página para expandir todos simultaneamente.

sessionInfo()

R version 4.5.0 (2025-04-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United Kingdom.utf8 
[2] LC_CTYPE=English_United Kingdom.utf8   
[3] LC_MONETARY=English_United Kingdom.utf8
[4] LC_NUMERIC=C                           
[5] LC_TIME=English_United Kingdom.utf8    

time zone: America/Sao_Paulo
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] lubridate_1.9.5  forcats_1.0.1    stringr_1.6.0    dplyr_1.2.1     
 [5] purrr_1.2.2      readr_2.2.0      tidyr_1.3.2      tibble_3.3.1    
 [9] tidyverse_2.0.0  broom_1.0.13     gt_1.3.0         ggfortify_0.4.19
[13] scales_1.4.0     patchwork_1.3.2  ggplot2_4.0.3    TSA_1.3.1       
[17] tseries_0.10-61  forecast_9.0.2  

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.58          htmlwidgets_1.6.4  lattice_0.22-6    
 [5] tzdb_0.5.0         quadprog_1.5-8     vctrs_0.7.3        tools_4.5.0       
 [9] generics_0.1.4     curl_7.1.0         parallel_4.5.0     xts_0.14.2        
[13] pkgconfig_2.0.3    Matrix_1.7-3       RColorBrewer_1.1-3 S7_0.2.2          
[17] lifecycle_1.0.5    compiler_4.5.0     farver_2.1.2       leaps_3.2         
[21] litedown_0.9       sass_0.4.10        htmltools_0.5.9    yaml_2.3.12       
[25] pillar_1.11.1      nlme_3.1-168       fracdiff_1.5-4     commonmark_2.0.0  
[29] tidyselect_1.2.1   locfit_1.5-9.12    digest_0.6.39      stringi_1.8.7     
[33] labeling_0.4.3     splines_4.5.0      fastmap_1.2.0      grid_4.5.0        
[37] colorspace_2.1-2   cli_3.6.6          magrittr_2.0.5     base64enc_0.1-6   
[41] dichromat_2.0-0.1  withr_3.0.2        backports_1.5.1    timechange_0.4.0  
[45] TTR_0.24.4         rmarkdown_2.31     quantmod_0.4.28    otel_0.2.0        
[49] timeDate_4052.112  gridExtra_2.3      hms_1.1.4          zoo_1.8-15        
[53] urca_1.3-4         evaluate_1.0.5     knitr_1.51         markdown_2.0      
[57] mgcv_1.9-4         rlang_1.2.0        Rcpp_1.1.1-1.1     glue_1.8.1        
[61] xml2_1.5.2         jsonlite_2.0.0     R6_2.6.1           fs_2.1.0