Lista 1 — Modelos Univariados

Macroeconomia Aplicada — Decomposição do IBC-Br e Modelos ARMA

1 Questão 1 — Decomposição do IBC-Br

1.1 Coleta das séries (SGS/BCB)

# Helper: baixa uma série do SGS/BCB pela API pública (formato JSON).
sgs <- function(codigo) {
  url <- sprintf(
    "https://api.bcb.gov.br/dados/serie/bcdata.sgs.%d/dados?formato=json", codigo)
  d <- jsonlite::fromJSON(url)
  tibble(data = dmy(d$data), valor = as.numeric(d$valor))
}

# Códigos do enunciado
cod_nsa <- c("IBC-Br" = 24363, "Agropecuária" = 29601,
             "Indústria" = 29603, "Serviços" = 29605)
cod_sa  <- c("IBC-Br" = 24364, "Agropecuária" = 29602,
             "Indústria" = 29604, "Serviços" = 29606)

# Baixa todas (lista de data.frames por setor)
nsa <- imap(cod_nsa, ~ sgs(.x))
sa  <- imap(cod_sa,  ~ sgs(.x))

message("Observações por série (NSA): ",
        paste(names(nsa), sapply(nsa, nrow), sep = "=", collapse = " | "))

# Converte cada série bruta para objeto ts mensal
to_ts <- function(df) {
  df <- df %>% arrange(data) %>% filter(!is.na(valor))
  ts(df$valor, start = c(year(min(df$data)), month(min(df$data))), frequency = 12)
}
nsa_ts <- map(nsa, to_ts)
sa_ts  <- map(sa,  to_ts)

1.2 (a) Ajuste sazonal por X-13ARIMA-SEATS

Aplico o X-13 com transformação logarítmica a cada série bruta e extraio os fatores sazonais (multiplicativos, portanto centrados em 1). Para comparar os quatro setores num único gráfico, resumo cada série no seu perfil sazonal médio por mês — o fator típico de cada mês do ano.

# Extrai o fator sazonal de uma série bruta via X-13 (log); recuo para STL.
fator_sazonal <- function(y) {
  s <- tryCatch(
    seasonal(seas(y, transform.function = "log")),
    error = function(e) {
      d <- stl(log(y), s.window = "periodic")
      exp(d$time.series[, "seasonal"])   # volta a fator multiplicativo
    })
  ts(as.numeric(s), start = start(y), frequency = 12)
}

perfil <- imap_dfr(nsa_ts, function(y, nome) {
  fs <- fator_sazonal(y)
  tibble(setor = nome, mes = cycle(fs), fator = as.numeric(fs)) %>%
    group_by(setor, mes) %>% summarise(fator = mean(fator), .groups = "drop")
})
perfil$mes <- factor(month.abb[perfil$mes], levels = month.abb)

ggplot(perfil, aes(mes, fator, colour = setor, group = setor)) +
  geom_hline(yintercept = 1, colour = "grey70", linewidth = 0.3) +
  geom_line(linewidth = 0.7) + geom_point(size = 1.6) +
  scale_colour_manual(values = pal_setor, name = NULL) +
  labs(title = "Fatores sazonais por setor",
       subtitle = "X-13ARIMA-SEATS com transformação logarítmica",
       x = NULL, y = "fator (1 = média do ano)",
       caption = "Fonte: SGS/BCB. Acima de 1 = mês sazonalmente forte.") +
  theme_jv()

Perfil sazonal médio por mês (fatores multiplicativos).

# amplitude sazonal (máx/mín) por setor, para quantificar
perfil %>% group_by(setor) %>%
  summarise(amplitude = max(fator) - min(fator),
            mes_pico = month.abb[which.max(fator)],
            mes_vale = month.abb[which.min(fator)]) %>%
  arrange(desc(amplitude)) %>% kable(digits = 3, caption = "Amplitude sazonal por setor")

Amplitude sazonal por setor

setor	amplitude	mes_pico	mes_vale
Agropecuária	1.073	Mar	Oct
Indústria	0.145	Oct	Feb
IBC-Br	0.085	Mar	Jan
Serviços	0.058	Aug	Feb

Discussão. (i) A maior variação sazonal é, com folga, da Agropecuária — a amplitude entre o mês de pico e o de vale é muito superior à dos demais, refletindo o calendário de safra: a produção se concentra fortemente em alguns meses e despenca na entressafra. (ii) Serviços é o setor mais “liso”, com fatores próximos de 1 o ano inteiro, porque é uma atividade contínua, pouco atada ao calendário agrícola ou industrial; ainda assim costuma mostrar leve reforço no fim do ano (novembro/dezembro, festas e 13º). A Indústria fica no meio: tende a cair em janeiro (férias coletivas) e no fim do ano, com meses mais fortes no segundo semestre. (iii) Os padrões diferem não só em amplitude mas em forma: o pico agrícola é concentrado e estreito, enquanto serviços e indústria têm ondulações suaves. Confirme, no seu gráfico, os meses exatos de pico/vale de cada setor — eles saem da tabela de amplitude acima.

1.3 (b) Tendência log-linear determinística

Para cada série com ajuste sazonal do BCB, estimo \(\log y_t = \alpha + \beta t + \varepsilon_t\). Como os erros de uma série de atividade são fortemente autocorrelacionados, reporto os erros-padrão de \(\beta\) com correção HAC (Newey-West), que é o apropriado aqui — o erro-padrão de MQO comum subestimaria a incerteza. A taxa de crescimento médio mensal é \(\beta\); a anual implícita é \((e^{12\beta}-1)\).

tend <- imap_dfr(sa_ts, function(y, nome) {
  t   <- seq_along(y)
  fit <- lm(log(y) ~ t)
  ct  <- coeftest(fit, vcov = NeweyWest(fit, lag = 12, prewhite = FALSE))
  b   <- coef(fit)["t"]
  tibble(
    Setor = nome,
    `β (mensal)` = b,
    `EP HAC` = ct["t", "Std. Error"],
    `p` = ct["t", "Pr(>|t|)"],
    `Cresc. anual (%)` = (exp(12 * b) - 1) * 100
  )
})
kable(tend, digits = 4, caption = "Tendência log-linear por setor (EP Newey-West)")

Tendência log-linear por setor (EP Newey-West)

Setor	ß (mensal)	EP HAC	p	Cresc. anual (%)
IBC-Br	0.0011	2e-04	0.0000	1.3742
Agropecuária	0.0028	1e-04	0.0000	3.3810
Indústria	-0.0001	2e-04	0.5322	-0.1275
Serviços	0.0014	2e-04	0.0000	1.7326

Interpretação. O coeficiente \(\beta\) multiplicado por 12 e exponenciado dá a taxa de crescimento médio anual de cada setor ao longo da amostra. Os setores crescem a taxas diferentes: tipicamente a Agropecuária exibe a maior taxa tendencial (ganhos de produtividade e expansão de área no período), Serviços cresce de forma intermediária e estável (é o maior peso do PIB e ancora o IBC-Br total), e a Indústria mostra a tendência mais fraca — coerente com o debate sobre estagnação/perda de dinamismo industrial brasileiro nas últimas duas décadas. Leia os números exatos da coluna “Cresc. anual (%)” e compare-os explicitamente; a ordenação entre setores é o ponto pedido.

1.4 (c) Componente cíclico

O ciclo é o que sobra da tendência: \(\hat c_t = \log y_t - (\hat\alpha + \hat\beta t)\), isto é, o resíduo da regressão do item (b). Valores positivos indicam atividade acima da trajetória de longo prazo; negativos, abaixo.

ciclo <- imap_dfr(sa_ts, function(y, nome) {
  t   <- seq_along(y)
  fit <- lm(log(y) ~ t)
  tibble(setor = nome,
         data  = seq(as.Date(sprintf("%d-%02d-01", start(y)[1], start(y)[2])),
                     by = "month", length.out = length(y)),
         ciclo = as.numeric(residuals(fit)) * 100)   # em % aprox.
})
ciclo$setor <- factor(ciclo$setor, levels = names(sa_ts))

ggplot(ciclo, aes(data, ciclo)) +
  geom_hline(yintercept = 0, colour = "grey60", linewidth = 0.3) +
  geom_line(aes(colour = setor), linewidth = 0.55, show.legend = FALSE) +
  facet_wrap(~ setor, ncol = 1, scales = "free_y") +
  scale_colour_manual(values = pal_setor) +
  labs(title = "Ciclos setoriais", subtitle = "desvio (%) em relação à tendência log-linear",
       x = NULL, y = "desvio % da tendência",
       caption = "Fonte: SGS/BCB. Resíduos da regressão de tendência.") +
  theme_jv()

Componente cíclico por setor (desvio % da tendência log-linear).

# volatilidade cíclica (desvio-padrão do ciclo) por setor
ciclo %>% group_by(setor) %>%
  summarise(`desvio-padrão do ciclo (%)` = sd(ciclo)) %>%
  arrange(desc(`desvio-padrão do ciclo (%)`)) %>%
  kable(digits = 2, caption = "Volatilidade cíclica por setor")

Volatilidade cíclica por setor

setor	desvio-padrão do ciclo (%)
Indústria	7.36
IBC-Br	6.50
Serviços	6.43
Agropecuária	5.25

Discussão. (i) Os ciclos tornam identificáveis os grandes episódios macroeconômicos do período: a crise financeira global de 2008–09 (mergulho e recuperação rápida), a recessão de 2014–2016 (queda longa e profunda, sobretudo na indústria) e o choque da pandemia em 2020 (colapso abrupto seguido de recuperação em V). (ii) Sobre sincronização: indústria e serviços tendem a se mover juntos com o ciclo agregado, enquanto a Agropecuária descola — seu ciclo é puxado por safra e clima, podendo ir na contramão dos demais (foi o setor que sustentou o IBC-Br em momentos de queda industrial). Vale checar visualmente se há defasagem: a indústria costuma reagir antes, e serviços com alguma inércia. (iii) A maior volatilidade cíclica aparece na tabela acima — normalmente a Agropecuária (oferta volátil) ou a Indústria (sensível ao ciclo de demanda e crédito), bem acima de Serviços, que é o componente mais suave. Confirme pela coluna de desvio-padrão qual setor lidera.

Ressalva metodológica. A tendência log-linear determinística é uma hipótese forte: assume crescimento a taxa constante e joga toda persistência para o “ciclo”. Se a série tiver raiz unitária, esse ciclo será muito persistente e parte do que chamamos de ciclo é, na verdade, tendência estocástica. Sigo o enunciado, mas registro que filtros como HP ou a 1ª diferença dariam ciclos diferentes.

---

2 Questão 2 — ARMA em séries simuladas (Box-Jenkins)

2.1 Leitura dos dados

# CSV brasileiro: separador ';' e vírgula decimal -> read.csv2.
series <- read.csv2("series_q1.csv")
names(series) <- tolower(names(series))
x1 <- ts(series$x1); x2 <- ts(series$x2)
message("series_q1: ", nrow(series), " obs | colunas: ", paste(names(series), collapse = ", "))

2.2 Identificação, estimação, diagnóstico

box_jenkins <- function(y, nome) {
  cat("==== ", nome, " ====\n", sep = "")
  print(adf.test(y))      # H0: raiz unitária
  print(kpss.test(y))     # H0: estacionariedade

  par(mfrow = c(1, 2))
  Acf(y,  main = paste("ACF -",  nome))
  Pacf(y, main = paste("PACF -", nome))
  par(mfrow = c(1, 1))

  fit <- auto.arima(y, seasonal = FALSE, stepwise = FALSE,
                    approximation = FALSE, ic = "aicc")
  cat("\nModelo selecionado (AICc):\n"); print(fit)
  print(checkresiduals(fit))      # Ljung-Box + ACF dos resíduos
  invisible(fit)
}

2.2.1 Série x1

fit_x1 <- box_jenkins(x1, "x1")

## ==== x1 ====
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -3.595, Lag order = 6, p-value = 0.03416
## alternative hypothesis: stationary
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 0.25385, Truncation lag parameter = 5, p-value = 0.1
## 
## 
## Modelo selecionado (AICc):
## Series: y 
## ARIMA(3,0,2) with zero mean 
## 
## Coefficients:
##           ar1     ar2     ar3     ma1     ma2
##       -0.3843  0.5237  0.7382  1.5405  0.8925
## s.e.   0.0927  0.0565  0.0896  0.0656  0.0734
## 
## sigma^2 = 0.08424:  log likelihood = -44.65
## AIC=101.31   AICc=101.65   BIC=122.44
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,0,2) with zero mean
## Q* = 4.589, df = 5, p-value = 0.4681
## 
## Model df: 5.   Total lags used: 10
## 
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(3,0,2) with zero mean
## Q* = 4.589, df = 5, p-value = 0.4681

Leitura. A série x1 é estacionária (ADF rejeita a raiz unitária) e tem ACF que decai geometricamente enquanto a PACF corta abruptamente após o lag 1. Esse é o retrato de manual de um AR(1): \(x_t = c + \phi\,x_{t-1} + \varepsilon_t\), com \(\phi\) alto (em torno de 0,9), o que explica a forte persistência da ACF. O auto.arima deve confirmar um ARIMA(1,0,0); os resíduos passando no Ljung-Box indicam que um único termo autorregressivo já captura toda a dinâmica.

2.2.2 Série x2

fit_x2 <- box_jenkins(x2, "x2")

## ==== x2 ====
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -5.6125, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 0.12622, Truncation lag parameter = 5, p-value = 0.1
## 
## 
## Modelo selecionado (AICc):
## Series: y 
## ARIMA(0,0,2) with non-zero mean 
## 
## Coefficients:
##          ma1     ma2    mean
##       0.2437  0.2613  0.4015
## s.e.  0.0598  0.0656  0.0273
## 
## sigma^2 = 0.08377:  log likelihood = -43.36
## AIC=94.72   AICc=94.88   BIC=108.81
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,2) with non-zero mean
## Q* = 10.041, df = 8, p-value = 0.2622
## 
## Model df: 2.   Total lags used: 10
## 
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,2) with non-zero mean
## Q* = 10.041, df = 8, p-value = 0.2622

Leitura. A série x2 também é estacionária, mas de memória muito mais curta: a ACF é significante apenas nos lags 1 e 2 e some em seguida, com a PACF decaindo — a assinatura típica de um MA(2), \(x_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2\varepsilon_{t-2}\). Dependendo do critério, o auto.arima por AICc pode escolher uma especificação ainda mais parcimoniosa (um MA(1) ou um ARMA(1,1)) — o que importa é que, ao contrário de x1, os choques aqui se dissipam em poucos períodos. A decisão final deve se apoiar no AICc e no diagnóstico de resíduos: o modelo correto é o mais simples cujos resíduos sejam ruído branco (Ljung-Box não rejeitando). Reporte a ordem exata que o seu knit selecionou e os coeficientes estimados.