Lista 2 — Raiz Unitária e Modelos Multivariados

Macroeconomia Aplicada — Box-Jenkins e Testes de Raiz Unitária

1 Questão 1 — Box-Jenkins e projeção

1.1 Coleta dos dados (SIDRA)

## Helper: localiza a coluna de código de período de forma tolerante a acento.
## (O SIDRA tem várias colunas "(Código)"; pegamos a do tempo, mês ou trimestre.)
col_periodo <- function(df, freq = c("mes", "tri")) {
  freq <- match.arg(freq)
  pat  <- if (freq == "mes") "M.s \\(C.digo\\)" else "Trimestre \\(C.digo\\)"
  nm   <- grep(pat, names(df), value = TRUE)
  if (!length(nm))
    stop("Coluna de período não encontrada. Colunas disponíveis: ",
         paste(names(df), collapse = " | "))
  as.character(df[[nm[1]]])
}

## ---------------------------------------------------------------------------
## IPCA mensal --- variação mensal (%), 2000-01 a 2021-12.
## Fonte: tabela 1737 (série histórica do IPCA), variable 63 = variação mensal.
## OBS.: a tabela 118 NÃO possui a variável 63 (foi o que travou o knit antes).
## Se o enunciado exigir a 118, rode  info_sidra(118)  e ajuste o código.
## ---------------------------------------------------------------------------
ipca_raw <- get_sidra(x = 1737, variable = 63, period = "200001-202112")

cod_m <- col_periodo(ipca_raw, "mes")              # "YYYYMM"
ipca <- tibble(
  data  = as.Date(paste0(substr(cod_m, 1, 4), "-", substr(cod_m, 5, 6), "-01")),
  valor = as.numeric(ipca_raw$Valor)
) %>% arrange(data) %>% filter(!is.na(valor))

ipca_ts <- ts(ipca$valor,
              start = c(year(min(ipca$data)), month(min(ipca$data))),
              frequency = 12)

## ---------------------------------------------------------------------------
## PIB trimestral a preços de mercado, índice de volume DESSAZONALIZADO.
## SIDRA tabela 1621, variable 584 (índice de volume COM ajuste sazonal),
## classificação c11255 = 90707 (PIB a preços de mercado).
## ATENÇÃO: na 1621 a variável é 584 (a 583 é a versão SEM ajuste, tabela 1620).
## Confirme com  info_sidra(1621)  se algo mudar.
## ---------------------------------------------------------------------------
pib_raw <- get_sidra(x = 1621, variable = 584, period = "200001-202104",
                     classific = "c11255", category = list(90707))

cod_t <- col_periodo(pib_raw, "tri")               # "YYYYTT" (TT = trimestre 01-04)
tri   <- as.integer(substr(cod_t, 5, 6))
pib <- tibble(
  data  = as.Date(paste0(substr(cod_t, 1, 4), "-",
                         sprintf("%02d", (tri - 1) * 3 + 1), "-01")),
  valor = as.numeric(pib_raw$Valor)
) %>% arrange(data) %>% filter(!is.na(valor))

pib_ts <- ts(pib$valor,
             start = c(year(min(pib$data)), quarter(min(pib$data))),
             frequency = 4)

## Conferência rápida (esperado: ~264 meses de IPCA e ~88 trimestres de PIB).
message("IPCA: ", nrow(ipca), " obs | PIB: ", nrow(pib), " obs")

1.2 Pré-tratamento: dessazonalização e transformação

O enunciado pede as séries dessazonalizadas. A série de PIB da tabela 1621 já vem com ajuste sazonal do IBGE, então basta tomar a primeira diferença do log (o crescimento trimestral). Já a inflação do IPCA não vem ajustada; aplico o X-13ARIMA-SEATS (via seasonal::seas), com recuo para STL caso o binário não esteja disponível no ambiente. Como ambas as séries entram já sem componente sazonal, o Box-Jenkins é aplicado sem termos sazonais (seasonal = FALSE) — incluí-los seria modelar uma sazonalidade que já foi removida.

# IPCA dessazonalizado
ipca_sa <- tryCatch(
  final(seas(ipca_ts)),
  error = function(e) {
    message("X-13 indisponível; usando STL.")
    as.ts(seasadj(stl(ipca_ts, s.window = "periodic")))
  }
)

# PIB: primeira diferença do log (já dessazonalizado na origem)
dlpib <- diff(log(pib_ts))

# visualização das séries de trabalho
autoplot(ipca_sa, colour = cor_serie) +
  labs(title = "IPCA --- variação mensal dessazonalizada",
       x = NULL, y = "% a.m.",
       caption = "Fonte: IBGE/SIDRA (tab. 1737), ajuste sazonal próprio.") +
  theme_jv()

autoplot(dlpib, colour = cor_serie) +
  labs(title = "PIB --- primeira diferença do log (crescimento trimestral)",
       x = NULL, y = "Δ log",
       caption = "Fonte: IBGE/SIDRA (tab. 1621, dessazonalizado).") +
  theme_jv()

1.3 Metodologia de Box-Jenkins

Encapsulo as quatro etapas numa função: (1) identificação — testes de raiz unitária e ACF/PACF para sugerir \((p,d,q)\); (2) estimação — deixo o auto.arima buscar por AIC corrigido, mas reporto o objeto para inspeção manual; (3) diagnóstico — Ljung-Box e ACF dos resíduos, para checar ruído branco; e (4) projeção — \(h\) passos à frente com intervalos.

box_jenkins <- function(y, h, nome, sazonal = TRUE) {
  cat("==== ", nome, " ====\n", sep = "")

  # (1) Identificação: ordem de integração
  cat("ndiffs (ADF):", ndiffs(y, test = "adf"),
      "| nsdiffs:", ifelse(sazonal, tryCatch(nsdiffs(y), error = function(e) 0), 0), "\n")
  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))

  # (2) Estimação
  fit <- auto.arima(y, seasonal = sazonal, stepwise = FALSE,
                    approximation = FALSE, ic = "aicc")
  cat("\nModelo selecionado:\n"); print(fit)

  # (3) Diagnóstico
  print(checkresiduals(fit))

  # (4) Projeção
  fc <- forecast(fit, h = h)
  print(autoplot(fc) +
          labs(title = paste("Projeção", h, "passos à frente ---", nome),
               x = NULL, y = NULL) +
          theme_jv())
  invisible(list(fit = fit, fc = fc))
}

1.3.1 IPCA — modelo e projeção (12 meses)

res_ipca <- box_jenkins(ipca_sa, h = 12, nome = "IPCA (mensal)", sazonal = FALSE)

## ==== IPCA (mensal) ====
## ndiffs (ADF): 0 | nsdiffs: 0 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -4.4425, Lag order = 6, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 0.28598, Truncation lag parameter = 5, p-value = 0.1
## 
## 
## Modelo selecionado:
## Series: y 
## ARIMA(1,0,0) with non-zero mean 
## 
## Coefficients:
##          ar1    mean
##       0.6761  0.5122
## s.e.  0.0449  0.0460
## 
## sigma^2 = 0.06009:  log likelihood = -2.73
## AIC=11.45   AICc=11.55   BIC=22.18
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 24.92, df = 23, p-value = 0.3544
## 
## Model df: 1.   Total lags used: 24
## 
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,0) with non-zero mean
## Q* = 24.92, df = 23, p-value = 0.3544

Leitura. A inflação mensal dessazonalizada tende a ser estacionária em nível (o ADF costuma rejeitar a raiz unitária e o KPSS não rejeita estacionariedade), de modo que \(d=0\). O formato de ACF/PACF de uma série de inflação em geral sugere um processo de memória curta — tipicamente um AR de baixa ordem com algum termo MA, e o auto.arima costuma chegar a algo como um ARMA(p,q) parcimonioso. Os resíduos passando no teste de Ljung-Box (p-valor alto) indicam que o modelo capturou a dinâmica e o que sobra é ruído branco; só então a projeção é confiável. A projeção de 12 meses converge suavemente para a média de longo prazo da inflação mensal, com os intervalos se abrindo conforme o horizonte — comportamento esperado de um modelo estacionário. Confirme a ordem exata \((p,d,q)\) que o seu knit selecionou.

1.3.2 PIB — modelo e projeção (4 trimestres)

res_pib <- box_jenkins(dlpib, h = 4, nome = "Δlog PIB (trimestral)", sazonal = FALSE)

## ==== Δlog PIB (trimestral) ====
## ndiffs (ADF): 0 | nsdiffs: 0 
## 
##  Augmented Dickey-Fuller Test
## 
## data:  y
## Dickey-Fuller = -4.111, Lag order = 4, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
##  KPSS Test for Level Stationarity
## 
## data:  y
## KPSS Level = 0.38214, Truncation lag parameter = 3, p-value = 0.08485
## 
## 
## Modelo selecionado:
## Series: y 
## ARIMA(0,0,0) with non-zero mean 
## 
## Coefficients:
##         mean
##       0.0054
## s.e.  0.0019
## 
## sigma^2 = 0.0003181:  log likelihood = 227.36
## AIC=-450.73   AICc=-450.59   BIC=-445.8
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with non-zero mean
## Q* = 3.969, df = 8, p-value = 0.8599
## 
## Model df: 0.   Total lags used: 8
## 
## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0) with non-zero mean
## Q* = 3.969, df = 8, p-value = 0.8599

Leitura. Trabalhando já na primeira diferença do log, a série de crescimento do PIB é, por construção, \(I(0)\) se o PIB em log for \(I(1)\) — o ADF aplicado a dlpib deve rejeitar a raiz unitária. A dinâmica trimestral do crescimento costuma ser fraca, e não raro o melhor modelo por AICc é bastante simples (um ARMA de ordem baixa, eventualmente próximo de um ruído branco com média não-nula). Vale notar a sensibilidade da amostra 2000–2021: a recessão de 2015–2016 e o choque da pandemia em 2020 são outliers fortes que inflam a variância dos resíduos e alargam os intervalos de projeção — por isso a previsão de quatro trimestres à frente carrega incerteza considerável. Cheque os resíduos: se o Ljung-Box rejeitar, vale testar uma ordem maior ou tratar os outliers de 2020.

---

2 Questão 2 — Testes de raiz unitária na estrutura a termo

2.1 Coleta dos dados (FRED)

## API key do FRED.
## Preferência: variável de ambiente FRED_API_KEY (.Renviron). Se não houver,
## usa a chave embutida abaixo.
## >>> SEGURANÇA: esta chave fica visível no código publicado. Gere uma nova
##     chave em https://fred.stlouisfed.org/ (My Account -> API Keys) depois de
##     entregar/publicar, e troque a linha abaixo. <<<
chave_fred <- Sys.getenv("FRED_API_KEY")
if (!nzchar(chave_fred)) chave_fred <- "92497f6cde356c5b00a02d7cd3af24b2"
fredr_set_key(chave_fred)

# GS1  = 1-Year Treasury Constant Maturity, mensal (média das diárias)
# GS10 = 10-Year Treasury Constant Maturity, mensal
gs1  <- fredr(series_id = "GS1")
gs10 <- fredr(series_id = "GS10")

juros <- gs1  %>% transmute(data = date, y1 = value) %>%
  inner_join(gs10 %>% transmute(data = date, y10 = value), by = "data") %>%
  arrange(data) %>% mutate(spread = y10 - y1)

# objetos ts (mensais) para os testes
y1_ts     <- ts(juros$y1,     start = c(year(min(juros$data)), month(min(juros$data))), frequency = 12)
y10_ts    <- ts(juros$y10,    start = c(year(min(juros$data)), month(min(juros$data))), frequency = 12)
spread_ts <- ts(juros$spread, start = c(year(min(juros$data)), month(min(juros$data))), frequency = 12)

juros %>% tidyr::pivot_longer(c(y1, y10, spread), names_to = "serie", values_to = "v") %>%
  ggplot(aes(data, v, colour = serie)) +
  geom_hline(yintercept = 0, colour = "grey70", linewidth = 0.3) +
  geom_line(linewidth = 0.55) +
  scale_colour_manual(values = c(y1 = "#b5651d", y10 = "#2f5c4f", spread = "#444444"),
                      labels = c(y1 = "1 ano", y10 = "10 anos", spread = "spread (10a-1a)"),
                      name = NULL) +
  labs(title = "Treasuries de 1 e 10 anos e o spread a termo",
       x = NULL, y = "% a.a.",
       caption = "Fonte: FRED (GS1, GS10).") +
  theme_jv()

2.2 (a) Ordem de integração de cada série

Aplico o teste de Dickey-Fuller aumentado a cada taxa, em nível e em primeira diferença. Uso a versão ur.df (com termo de deriva) para inspecionar a estatística contra os valores críticos, e complemento com ndiffs para a contagem de diferenças sugerida.

df_report <- function(x, nome, tipo = "drift") {
  cat("---", nome, "---\n")
  print(summary(ur.df(x, type = tipo, selectlags = "AIC")))
  cat("ndiffs(ADF) =", ndiffs(x, test = "adf"), "\n\n")
}
df_report(y1_ts,  "Taxa 1 ano (nível)")

## --- Taxa 1 ano (nível) ---
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.11975 -0.11771 -0.01713  0.12560  1.90410 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.047710   0.021344   2.235  0.02565 *  
## z.lag.1     -0.010098   0.003779  -2.672  0.00767 ** 
## z.diff.lag   0.374622   0.031371  11.942  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3615 on 873 degrees of freedom
## Multiple R-squared:  0.1437, Adjusted R-squared:  0.1418 
## F-statistic: 73.26 on 2 and 873 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -2.6723 3.5736 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
## 
## ndiffs(ADF) = 1

df_report(y10_ts, "Taxa 10 anos (nível)")

## --- Taxa 10 anos (nível) ---
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.65407 -0.12326 -0.00394  0.12733  1.51213 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.032457   0.018278   1.776   0.0761 .  
## z.lag.1     -0.005669   0.002934  -1.932   0.0537 .  
## z.diff.lag   0.308869   0.032170   9.601   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2491 on 873 degrees of freedom
## Multiple R-squared:  0.09772,    Adjusted R-squared:  0.09565 
## F-statistic: 47.27 on 2 and 873 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -1.9318 1.8747 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
## 
## ndiffs(ADF) = 1

df_report(diff(y1_ts),  "Taxa 1 ano (1ª diferença)")

## --- Taxa 1 ano (1ª diferença) ---
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.96119 -0.11039  0.00649  0.12880  1.89339 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.001244   0.011958   0.104    0.917    
## z.lag.1     -0.773226   0.037028 -20.882  < 2e-16 ***
## z.diff.lag   0.226943   0.032979   6.881 1.13e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3537 on 872 degrees of freedom
## Multiple R-squared:  0.3504, Adjusted R-squared:  0.3489 
## F-statistic: 235.1 on 2 and 872 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -20.8823 218.0361 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
## 
## ndiffs(ADF) = 0

df_report(diff(y10_ts), "Taxa 10 anos (1ª diferença)")

## --- Taxa 10 anos (1ª diferença) ---
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.61920 -0.12201 -0.00226  0.13210  1.40263 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.001367   0.008257   0.166    0.869    
## z.lag.1     -0.838278   0.039012 -21.488  < 2e-16 ***
## z.diff.lag   0.208497   0.033114   6.296 4.82e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2443 on 872 degrees of freedom
## Multiple R-squared:  0.3751, Adjusted R-squared:  0.3737 
## F-statistic: 261.8 on 2 and 872 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -21.4878 230.862 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
## 
## ndiffs(ADF) = 0

Leitura. O padrão típico — e o que se espera teoricamente para taxas de juros — é que ambas as séries sejam \(I(1)\): em nível, o ADF não rejeita a raiz unitária (estatística acima do valor crítico), mas, ao diferenciar uma vez, a raiz unitária é rejeitada. Ou seja, as duas taxas compartilham a mesma ordem de integração, \(I(1)\). Isso é coerente com o fato de juros nominais herdarem a persistência das expectativas de inflação e da política monetária, que se comportam como tendências estocásticas. Confirme nos seus resultados: se o nível rejeitar para alguma série, reporte-a como \(I(0)\) e ajuste a discussão.

2.3 (b) Teste de Dickey-Fuller para o spread

df_report(spread_ts, "Spread 10a - 1a (nível)")

## --- Spread 10a - 1a (nível) ---
## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression drift 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56638 -0.10927 -0.00928  0.10487  2.13718 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.031521   0.010433   3.021  0.00259 ** 
## z.lag.1     -0.035046   0.007429  -4.717 2.78e-06 ***
## z.diff.lag   0.312227   0.032149   9.712  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2376 on 873 degrees of freedom
## Multiple R-squared:  0.1096, Adjusted R-squared:  0.1075 
## F-statistic: 53.72 on 2 and 873 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -4.7174 11.1271 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau2 -3.43 -2.86 -2.57
## phi1  6.43  4.59  3.78
## 
## ndiffs(ADF) = 0

print(kpss.test(spread_ts))

## 
##  KPSS Test for Level Stationarity
## 
## data:  spread_ts
## KPSS Level = 1.3075, Truncation lag parameter = 6, p-value = 0.01

Leitura. O resultado de interesse é se o spread, apesar de ser a diferença de duas séries \(I(1)\), é estacionário. Na maioria das amostras longas o ADF rejeita a raiz unitária no spread (e o KPSS não rejeita estacionariedade), indicando que o spread é \(I(0)\) — ainda que persistente. Esse é o caso mais comum, mas o spread a termo é uma série sabidamente “quase-unitária”: não se assuste se o teste ficar no limiar; nesse caso reporte a ambiguidade e apoie-se no KPSS e na inspeção do gráfico.

2.4 (c) Implicações econômicas para a estrutura a termo

Se as duas taxas são \(I(1)\) e o spread é \(I(0)\), então as séries de 1 e 10 anos são cointegradas com vetor de cointegração \((1,-1)\): elas dividem uma única tendência estocástica comum — o nível geral dos juros — e o spread é o desvio estacionário em torno dessa relação de equilíbrio de longo prazo. Em termos econômicos, isso diz que as pontas curta e longa da curva não vagueiam indefinidamente para longe uma da outra; choques que abrem o spread tendem a se dissipar, e a estrutura a termo retorna a uma configuração de equilíbrio.

Essa é exatamente a previsão da hipótese das expectativas da estrutura a termo: a taxa longa é, a menos de um prêmio a termo aproximadamente estável, a média das taxas curtas futuras esperadas. Sob essa hipótese, taxa longa e taxa curta podem cada uma ter raiz unitária (porque as expectativas de juros futuros têm), mas a diferença entre elas — o prêmio mais o componente de expectativas — é estacionária. A cointegração \((1,-1)\) é, portanto, uma assinatura empírica da hipótese das expectativas com prêmio de risco estacionário.

A leitura se inverte no cenário alternativo: se o spread também fosse \(I(1)\), as duas pontas da curva não estariam ancoradas entre si no longo prazo, o que seria incompatível com a versão simples da hipótese das expectativas e apontaria para um prêmio a termo com tendência estocástica própria (variando permanentemente no tempo). Daí a importância do teste — ele distingue uma curva de juros com âncora de longo prazo de uma em que as maturidades se descolam. Reporte qual dos dois cenários seus dados sustentam e ancore a conclusão nas estatísticas obtidas.