# ============================================================
# CDI mensal | Pré-teste (log & sazonalidade) | SARIMA x ARIMA | Previsão
# ============================================================
# Repositório CRAN e timeout
options(repos = c(CRAN = "https://cloud.r-project.org"))
options(timeout = max(getOption("timeout"), 120))
suppressPackageStartupMessages({
library(rbcb); library(httr); library(jsonlite); library(dplyr)
library(lubridate); library(xts); library(zoo)
library(forecast); library(tseries); library(knitr)
})
# ---- 1) Coleta robusta do CDI diário (% a.a.) ----
get_cdi_daily_resilient <- function(start_date = "2015-01-01", end_date = Sys.Date()) {
# 1) Tenta série 4389 (CDI diário % a.a.)
df <- NULL
try({ df <- rbcb::get_series(4389, start_date = start_date, end_date = end_date) }, silent = TRUE)
if (!is.null(df) && nrow(df)) {
names(df) <- c("date","value"); df$date <- as.Date(df$date)
return(dplyr::arrange(df, date))
}
# 2) Tenta série 12 (alternativa % a.a.)
df2 <- NULL
try({ df2 <- rbcb::get_series(12, start_date = start_date, end_date = end_date) }, silent = TRUE)
if (!is.null(df2) && nrow(df2)) {
names(df2) <- c("date","value"); df2$date <- as.Date(df2$date)
return(dplyr::arrange(df2, date))
}
# 3) Fallback HTTP com RETRY
sgs_url <- function(id, di, df) {
paste0(
"https://api.bcb.gov.br/dados/serie/bcdata.sgs.", id,
"/dados?formato=json&dataInicial=", format(as.Date(di), "%d/%m/%Y"),
"&dataFinal=", format(as.Date(df), "%d/%m/%Y")
)
}
fetch_chunk <- function(id, di, df) {
resp <- httr::RETRY(
verb = "GET", url = sgs_url(id, di, df),
times = 4, pause_base = 0.7, pause_cap = 3,
httr::user_agent("Mozilla/5.0 (Macintosh) R/httr fallback"),
httr::accept_json(), httr::timeout(60)
)
if (httr::http_error(resp)) return(tibble::tibble(date = as.Date(character()), value = numeric()))
txt <- httr::content(resp, as = "text", encoding = "UTF-8")
jj <- tryCatch(jsonlite::fromJSON(txt), error = function(e) NULL)
if (is.null(jj) || !NROW(jj)) return(tibble::tibble(date = as.Date(character()), value = numeric()))
tibble::tibble(
date = as.Date(jj$data, "%d/%m/%Y"),
value = suppressWarnings(as.numeric(gsub(",", ".", jj$valor)))
)
}
fetch_range <- function(id, start_date, end_date, by = "1 month") {
pts <- seq(as.Date(start_date), as.Date(end_date), by = by)
if (tail(pts,1) < as.Date(end_date)) pts <- c(pts, as.Date(end_date))
out <- vector("list", length(pts) - 1)
for (i in seq_along(out)) {
di <- pts[i]; df <- pts[i+1] - 1
out[[i]] <- fetch_chunk(id, di, df)
Sys.sleep(0.15)
}
dplyr::bind_rows(out) |>
dplyr::filter(!is.na(value)) |>
dplyr::distinct(date, .keep_all = TRUE) |>
dplyr::arrange(date)
}
df3 <- fetch_range(4389, start_date, end_date); if (nrow(df3)) return(df3)
df4 <- fetch_range(12, start_date, end_date); if (nrow(df4)) return(df4)
stop("Não foi possível obter o CDI (4389/12) mesmo com fallbacks.")
}
cdi_daily <- get_cdi_daily_resilient("2015-01-01", Sys.Date())
cdi_daily <- dplyr::filter(cdi_daily, !is.na(value))
# ---- 2) CDI diário (% a.a.) -> CDI mensal (% a.m., média do mês do fator equivalente) ----
cdi_daily_xts <- xts::xts(cdi_daily$value, order.by = cdi_daily$date) # % a.a.
cdi_m_equiv <- (1 + cdi_daily_xts/100)^(1/12) - 1 # fração a.m.
cdi_m_frac <- xts::apply.monthly(cdi_m_equiv, mean, na.rm = TRUE) # média no mês
cdi_m <- 100 * cdi_m_frac # % a.m.
cdi_m_ts <- ts(
as.numeric(cdi_m),
start = c(lubridate::year(start(zoo::as.zoo(cdi_m))),
lubridate::month(start(zoo::as.zoo(cdi_m)))),
frequency = 12
)
# Histórico para gráficos finais (definido ANTES das figuras finais)
hist_df <- data.frame(
date = as.Date(zoo::index(cdi_m)),
cdi = as.numeric(cdi_m)
)
cat("Janela da amostra mensal: ",
format(min(hist_df$date), "%Y-%m"), " → ", format(max(hist_df$date), "%Y-%m"),
"\n", sep = "")
## Janela da amostra mensal: 2015-01 → 2025-11
# Exploratória rápida
plot.ts(cdi_m_ts, main = "CDI mensal (% a.m.)", ylab = "% a.m.", xlab = "Tempo")
# ---- 3) Pré-teste: log e sazonalidade ----
y <- log1p(cdi_m_ts/100) # estabiliza a variância
d_sug <- forecast::ndiffs(y, alpha = 0.05, test = "adf")
D_sug <- forecast::nsdiffs(y, m = 12, test = "ocsb")
cat("Diferenças sugeridas: d =", d_sug, " | D =", D_sug, "\n")
## Diferenças sugeridas: d = 2 | D = 0
adf_lvl <- tseries::adf.test(y)
cat("ADF no nível (log CDI): p-valor =", signif(adf_lvl$p.value, 3), "\n")
## ADF no nível (log CDI): p-valor = 0.494
if (D_sug >= 1) {
yD <- stats::na.omit(diff(y, lag = 12))
adf_D <- tseries::adf.test(yD)
cat("ADF após D=1: p-valor =", signif(adf_D$p.value, 3), "\n")
}
stl_fit <- stats::stl(ts(as.numeric(y), frequency = 12, start = start(y)),
s.window = "periodic", robust = TRUE)
plot(stl_fit, main = "STL em y (log do CDI % a.m.)")
s_comp <- stl_fit$time.series[, "seasonal"]
r_comp <- stl_fit$time.series[, "remainder"]
forca <- 1 - var(r_comp, na.rm = TRUE) / var(r_comp + s_comp, na.rm = TRUE)
cat("Força de sazonalidade (0~1):", round(forca, 3), "\n")
## Força de sazonalidade (0~1): 0.046
# ---- ADF sazonal para justificar D = 0 ----
y_seas_diff <- stats::na.omit(diff(y, lag = 12, differences = 1))
adf_seas <- tseries::adf.test(y_seas_diff)
cat("ADF na diferença sazonal (lag=12): p-valor =", signif(adf_seas$p.value, 3), "\n")
## ADF na diferença sazonal (lag=12): p-valor = 0.0151
print(adf_seas)
##
## Augmented Dickey-Fuller Test
##
## data: y_seas_diff
## Dickey-Fuller = -3.9284, Lag order = 4, p-value = 0.01513
## alternative hypothesis: stationary
# ---- 4) SARIMA (com sazonalidade) ----
sarima_fit <- forecast::auto.arima(
y, seasonal = TRUE, D = D_sug, stepwise = FALSE, approximation = FALSE
)
print(sarima_fit)
## Series: y
## ARIMA(2,2,2)(0,0,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 ma2 sma1
## -0.6421 -0.6424 0.0077 0.6003 -0.2028
## s.e. 0.2146 0.1262 0.2271 0.1290 0.1158
##
## sigma^2 = 2.205e-08: log likelihood = 955.81
## AIC=-1899.62 AICc=-1898.93 BIC=-1882.46
sarima_ord <- forecast::arimaorder(sarima_fit)
sarima_aic <- AIC(sarima_fit); sarima_bic <- BIC(sarima_fit)
# ---- 5) ARIMA (sem sazonalidade) ----
arima_fit <- forecast::auto.arima(
y, seasonal = FALSE, stepwise = FALSE, approximation = FALSE
)
print(arima_fit)
## Series: y
## ARIMA(2,2,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.5640 -0.6219 -0.0508 0.6710
## s.e. 0.1386 0.1026 0.1330 0.0908
##
## sigma^2 = 2.247e-08: log likelihood = 954.28
## AIC=-1898.56 AICc=-1898.07 BIC=-1884.26
arima_ord <- forecast::arimaorder(arima_fit)
arima_aic <- AIC(arima_fit); arima_bic <- BIC(arima_fit)
# ---- 6A) Teste de Razão de Verossimilhança (ARIMA vs SARIMA) ----
ll_ar <- as.numeric(logLik(arima_fit))
ll_sa <- as.numeric(logLik(sarima_fit))
# df = número de parâmetros adicionais do SARIMA em relação ao ARIMA
k_df <- length(coef(sarima_fit)) - length(coef(arima_fit))
k_df <- ifelse(k_df < 1, 1, k_df) # proteção de segurança
LR <- -2 * (ll_ar - ll_sa)
p_lr <- pchisq(LR, df = k_df, lower.tail = FALSE)
cat(sprintf("LR test ARIMA vs SARIMA: LR = %.3f | df = %d | p-valor = %.4f\n",
LR, k_df, p_lr))
## LR test ARIMA vs SARIMA: LR = 3.054 | df = 1 | p-valor = 0.0805
# Tabela opcional para incluir no relatório
kable(
data.frame(Modelo = c("ARIMA", "SARIMA"), LogLik = c(ll_ar, ll_sa)),
caption = "Log-likelihood dos modelos"
)
Log-likelihood dos modelos
| ARIMA |
954.2809 |
| SARIMA |
955.8079 |
# ---- 6) Comparação AIC/BIC (evita & no LaTeX) ----
comp <- data.frame(
Modelo = c(
sprintf("SARIMA(%d,%d,%d)(%d,%d,%d)[12]",
sarima_ord["p"],sarima_ord["d"],sarima_ord["q"],
sarima_ord["P"],sarima_ord["D"],sarima_ord["Q"]),
sprintf("ARIMA(%d,%d,%d)",
arima_ord["p"],arima_ord["d"],arima_ord["q"])
),
AIC = c(round(sarima_aic,2), round(arima_aic,2)),
BIC = c(round(sarima_bic,2), round(arima_bic,2))
)
kable(comp, caption = "Comparação de AIC e BIC entre SARIMA e ARIMA")
Comparação de AIC e BIC entre SARIMA e ARIMA
| SARIMA(2,2,2)(0,0,1)[12] |
-1899.62 |
-1882.46 |
| ARIMA(2,2,2) |
-1898.56 |
-1884.26 |
# Como AIC e BIC são muito próximos e a força de sazonalidade é baixa (~0.05),
# e o LR test não rejeita ARIMA como modelo suficiente, escolhemos ARIMA.
melhor_modelo <- "ARIMA"
cat("\nDecisão: AIC e BIC praticamente iguais; sazonalidade fraca; LR test não significativo.\n")
##
## Decisão: AIC e BIC praticamente iguais; sazonalidade fraca; LR test não significativo.
cat("=> Modelo final escolhido por critério estatístico e parcimônia: ARIMA.\n\n")
## => Modelo final escolhido por critério estatístico e parcimônia: ARIMA.
cat("Modelo final selecionado:", melhor_modelo, "\n")
## Modelo final selecionado: ARIMA
fit_final <- if (melhor_modelo == "ARIMA") arima_fit else sarima_fit
ord_final <- forecast::arimaorder(fit_final)
modelo_str <- if (melhor_modelo == "ARIMA") {
sprintf("ARIMA(%d,%d,%d)", ord_final["p"],ord_final["d"],ord_final["q"])
} else {
sprintf("SARIMA(%d,%d,%d)(%d,%d,%d)[12]",
ord_final["p"],ord_final["d"],ord_final["q"],
ord_final["P"],ord_final["D"],ord_final["Q"])
}
# Definir h **antes** de imprimir
h <- 12
cat("Modelo final:", modelo_str, "| h =", h, "\n")
## Modelo final: ARIMA(2,2,2) | h = 12
# ---- 7) FAC/FACP do MODELO ESCOLHIDO (série estacionária usada pelo modelo) ----
# --- Função correta para estacionarizar conforme ord_final (usa base::diff) ---
build_stationary_series <- function(y, ord_final) {
# pega d e D com fallback seguro
d <- as.integer(ifelse(is.na(ord_final["d"]), 0, ord_final["d"]))
D <- as.integer(ifelse(is.na(ord_final["D"]), 0, ord_final["D"]))
m <- 12L
y_stat <- y
# diferença sazonal primeiro (se houver)
if (D > 0) {
y_stat <- stats::na.omit(base::diff(y_stat, lag = m, differences = D))
}
# diferença não sazonal depois (se houver)
if (d > 0) {
y_stat <- stats::na.omit(base::diff(y_stat, differences = d))
}
# garante classe ts
y_stat <- ts(as.numeric(y_stat), frequency = frequency(y))
return(y_stat)
}
y_stat <- build_stationary_series(y, ord_final)
par(mfrow=c(1,2), mar=c(4.5,4,3.5,1))
acf(y_stat, lag.max = 36,
main = paste("FAC -", modelo_str),
ylab = "Autocorrelação",
xlab = "Defasagem",
cex.main = 0.9)
pacf(y_stat, lag.max = 36,
main = paste("FACP -", modelo_str),
ylab = "Autocorrelação Parcial",
xlab = "Defasagem",
cex.main = 0.9)
par(mfrow=c(1,1))
# ---- 8) Diagnósticos do modelo final (resíduos) ----
res_final <- residuals(fit_final)
par(mfrow=c(1,2), mar=c(4.5,4,3.5,1))
acf(res_final, main="FAC dos resíduos", ylab="Autocorrelação", xlab="Defasagem", cex.main=0.9)
acf(res_final^2,main="FAC dos resíduos²", ylab="Autocorrelação", xlab="Defasagem", cex.main=0.9)
par(mfrow=c(1,1))
lb <- Box.test(res_final, lag=24, type="Ljung-Box")
lb2 <- Box.test(res_final^2, lag=24, type="Ljung-Box")
jb <- tseries::jarque.bera.test(res_final)
cat("\n--- Teste Ljung-Box (resíduos) ---\n")
##
## --- Teste Ljung-Box (resíduos) ---
print(lb)
##
## Box-Ljung test
##
## data: res_final
## X-squared = 18.696, df = 24, p-value = 0.768
cat("\n--- Ljung-Box (resíduos²) ---\n")
##
## --- Ljung-Box (resíduos²) ---
print(lb2)
##
## Box-Ljung test
##
## data: res_final^2
## X-squared = 16.79, df = 24, p-value = 0.8575
cat("\n--- Jarque-Bera (normalidade) ---\n")
##
## --- Jarque-Bera (normalidade) ---
print(jb)
##
## Jarque Bera Test
##
## data: res_final
## X-squared = 68.418, df = 2, p-value = 1.443e-15
diag_tab <- data.frame(
Teste = c("Ljung-Box(res)", "Ljung-Box(res²)", "Jarque-Bera"),
`p-valor` = signif(c(lb$p.value, lb2$p.value, jb$p.value), 3)
)
kable(diag_tab, caption = "Testes de diagnóstico do modelo final")
Testes de diagnóstico do modelo final
| Ljung-Box(res) |
0.768 |
| Ljung-Box(res²) |
0.857 |
| Jarque-Bera |
0.000 |
# ---- 9) Previsão determinística (12 meses) ----
fc_final <- forecast::forecast(fit_final, h = 12)
last_ym <- zoo::as.yearmon(tail(zoo::index(cdi_m), 1))
pred_final <- data.frame(
date = as.Date(seq(last_ym + 1/12, by = 1/12, length.out = 12)),
mean = (exp(as.numeric(fc_final$mean)) - 1) * 100,
lo95 = (exp(as.numeric(fc_final$lower[,"95%"])) - 1) * 100,
hi95 = (exp(as.numeric(fc_final$upper[,"95%"])) - 1) * 100
)
# Gráfico: histórico + previsão determinística
plot(hist_df$date, hist_df$cdi, type="l", xlab="Tempo", ylab="% a.m.",
main=paste0("CDI mensal + previsão (", modelo_str, ")"))
lines(pred_final$date, pred_final$mean, lwd=2)
polygon(c(pred_final$date, rev(pred_final$date)),
c(pred_final$lo95, rev(pred_final$hi95)),
col=rgb(0,0,0,0.15), border=NA)
cat("\nPrevisão concluída com o modelo: ", modelo_str, "\n", sep="")
##
## Previsão concluída com o modelo: ARIMA(2,2,2)