Estudo sobre Gradient Boosting (XGBoost)

 

Gradient Boosting

Bibliotecas necessárias

library(randomForest)
library(tidyr)
library(rpart)
library(glmnet)
library(xgboost)
library(tibble)
library(ggplot2)
library(scales)
library(dplyr)

options(scipen = 999)

Base de dados

source("C:/Users/Thomas/Desktop/FGC/Códigos/Base de dados.R")
df_macro <- df_macro %>% select(-any_of("hiato_inflacao"))

Definindo parâmetros

params <- list(
 objective = "reg:squarederror", # regressão
 eta = 0.05, # learning rate
 max_depth = 4, # profundidade máx de cada árvore
 min_child_weight = 1, # tamanho mínimo de nó
 subsample = 0.8, # fração de dados para cada árvore
 colsample_bytree = 0.8, # fração de variáveis para cada árvore
 gamma = 0, # penalidade por complexidade
 lambda = 1, # regularização L2
 alpha = 0 # regularização L1
)

Treinando modelo

Matriz de desenho

#dummies tratadas automaticamente
df_model <- df_model %>% select(-any_of("hiato_inflacao"))
X <- model.matrix(inflacao ~ . - 1, data = df_macro)
y <- df_macro$inflacao
X_scaled <- scale(X)
#Criar estrutura de dados XGBoost
df_model <- na.omit(df_macro)

df_model$data <- NULL

X <- model.matrix(inflacao ~ . - 1, data = df_model)
y <- df_model$inflacao

X_scaled <- scale(X)

nrow(X_scaled)
length(y)

dtrain <- xgb.DMatrix(data = X_scaled, label = y)

Validação cruzada

xgb_cv <- xgb.cv(
  params = params,
  data = dtrain,
  nrounds = 1000,
  nfold = 5,
  showsd = TRUE,
  metrics = "rmse",
  early_stopping_rounds = 20,
  verbose = 0
)

Gráfico - Evolução do Modelo

log <- xgb_cv$evaluation_log

plot(log$iter, log$test_rmse_mean, type = "l",
     xlab = "Iterações", ylab = "RMSE",
     main = "Curva de aprendizado (XGBoost)")

lines(log$iter, log$train_rmse_mean, lty = 2)

Número ótimo de rodadas

Não se deseja um número muito alto de rodadas em função de nossa limitada capacidade computacional. Contudo, é perceptível a queda do erro quadrático médio com o crescimento das iterações, o que sugere que alto erro para árvore individual.

log <- xgb_cv$evaluation_log

best_round <- which.min(log$test_rmse_mean)

cat("Número ótimo de rounds:", best_round, "\n")

Treinando modelo final

train_idx <- sample(1:nrow(X_scaled), 0.8 * nrow(X_scaled))

X_train <- X_scaled[train_idx, ]
X_test  <- X_scaled[-train_idx, ]

y_train <- y[train_idx]
y_test  <- y[-train_idx]
dtrain <- xgb.DMatrix(data = X_train, label = y_train)
dtest  <- xgb.DMatrix(data = X_test, label = y_test)

xgb_final <- xgb.train(
  params = params,
  data = dtrain,  # agora é o dtrain do split
  best_round <- which.min(xgb_cv$evaluation_log$test_rmse_mean),
  verbose = 0
)

Importância das variáveis

importance_xgb <- xgb.importance(model = xgb_final)

print(head(importance_xgb[, c("Feature", "Gain")], 10))

O modelo Gradient Boosting não pretende revelar os coeficientes associados a cada variável, muito menos o sinal que cada um teria, por outro lado, este modelo reflete apenas o grau de importância que cada variável tem em explicar os resultados de inflação. De fato, o modelo mostra-nos que a variável mais importante neste cenário é a variação do câmbio, seguida por expectativa de inflação e taxa de desemprego.

A taxa de câmbio reflete de maneira quase pleonástica a inflação, pois revela a valorização ou desvalorização da moeda frente a outras divisas estrangeiras. Devido ao fato de ser possível haver inflação mesmo quando a taxa de câmbio permanece estável contra moedas comuns em comparações (o que é usual, inclusive), conclui-se que a flutuação da taxa de câmbio não necessariamente indica inflação. Outro cenário que é interessante é a possibilidade de a taxa de câmbio desvalorizar-se mesmo em cenário deflacionário.

Expectativa de inflação, embora com menor importância, tem papel relevante visto que não apenas revela-se como preditor do mercado ao futuro da moeda como também materializa-se na economia através de ansia por maiores salários, remarcações de preços de outros bens, etc. Assim sendo, se a expectativa de inflação é crescente - em virtude da perda de confiança na moeda - haverá uma pressão constante em favor de se adiantar preços maiores futuros, evitando prejuízos econômicos. Outro aspecto relevante e intuitivo, embora pouco comentado, é o fato de que uma maior expectativa de inflação também prediz maior inflação porque pode se traduzir em redução na demanda por moeda.

Taxa de desemprego pode ser o ponto de maior complexidade na atual interlocação, visto que usualmente a teoria associa maior empregabilidade a maior demanda por bens na economia e, assim, maiores preços e maior inflação. Contudo, deve-se afirmar o básico: inflação trata-se de desvalorização da moeda. Neste sentido, é possível que com maior taxa de desemprego, na realidade, a demanda por moeda caia, o que aumentaria a taxa de inflação.

# checar/gerar importance_xgb
if (!exists("importance_xgb")) {
  if (!exists("xgb_final")) stop("xgb_final não encontrado. Treine o modelo antes.")
  importance_xgb <- xgb.importance(model = xgb_final)
}
if (nrow(importance_xgb) == 0) stop("importance_xgb vazio.")

# localizar colunas (tolerante a case)
col_names <- tolower(names(importance_xgb))
feat_col <- names(importance_xgb)[which(col_names == "feature")[1]]
gain_col <- names(importance_xgb)[which(col_names == "gain")[1]]
if (is.na(feat_col) || is.na(gain_col)) stop("Colunas 'Feature' ou 'Gain' não encontradas em importance_xgb.")


top_n <- 10  # defina quantas variáveis quer mostrar

imp_tbl <- xgb.importance(model = xgb_final)
imp_tbl <- imp_tbl[1:top_n, ]

# gráfico horizontal elegante
ggplot(imp_tbl, aes(x = Feature, y = Gain, fill = Gain)) +
  geom_col(width = 0.7, show.legend = FALSE) +
  geom_text(aes(label = round(Gain, 4)), 
            hjust = -0.1, size = 3.5, color = "black") +
  coord_flip() +
  scale_fill_gradient(low = "#BFD3E6", high = "#08306B") +
  theme_minimal(base_size = 12) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  ) +
  labs(
    title = paste0("Top ", top_n, " - Importância das Variáveis (XGBoost)"),
    y = "Gain"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15)))

Previsões

pred_xgb <- predict(xgb_final, dtest)

rmse_xgb <- sqrt(mean((y_test - pred_xgb)^2))
r2_xgb <- 1 - sum((y_test - pred_xgb)^2) / sum((y_test - mean(y_test))^2)
cat("RMSE XGBoost:", rmse_xgb, "
")
cat("R² XGBoost:", r2_xgb, "
")

Comparação com Random Forest e Decision Tree

Preparação dos modelos

# limpeza e treino compacto para evitar NA em preditores
library(rpart); library(randomForest); library(tibble)
set.seed(123)

# 1) detectar NAs na base usada para model.matrix
if (!exists("df_model")) df_model <- df_macro


# 2) imputação simples: mediana para numéricas, moda para categóricas
num_cols <- sapply(df_model, is.numeric)
cat_cols <- sapply(df_model, function(x) is.factor(x) || is.character(x))

# mediana para numéricas
for (nm in names(df_model)[num_cols]) {
  if (any(is.na(df_model[[nm]]))) {
    df_model[[nm]][is.na(df_model[[nm]])] <- median(df_model[[nm]], na.rm = TRUE)
  }
}

# função moda
mode_val <- function(x) {
  ux <- na.omit(unique(x))
  if (length(ux) == 0) return(NA)
  ux[which.max(tabulate(match(x, ux)))]
}

# moda para categóricas
for (nm in names(df_model)[cat_cols]) {
  if (any(is.na(df_model[[nm]]))) {
    df_model[[nm]][is.na(df_model[[nm]])] <- mode_val(df_model[[nm]])
  }
  # garantir fator
  if (!is.factor(df_model[[nm]])) df_model[[nm]] <- as.factor(df_model[[nm]])
}


# 3) recriar X e y, remover colunas de variância zero antes de escalar
X <- model.matrix(inflacao ~ . - 1, data = df_model)
y <- df_model$inflacao

# remover colunas com sd == 0 (causa NaN ao escalar)
zero_sd <- apply(X, 2, sd, na.rm = TRUE) == 0
if (any(zero_sd)) {
  cat("Removendo colunas com variância zero:", sum(zero_sd), "\n")
  X <- X[, !zero_sd, drop = FALSE]
}

# 4) escalar e split
X_scaled <- scale(X)
train_idx <- sample(seq_len(nrow(X_scaled)), size = floor(0.8 * nrow(X_scaled)))
X_train <- X_scaled[train_idx, , drop = FALSE]
X_test  <- X_scaled[-train_idx, , drop = FALSE]
y_train <- y[train_idx]
y_test  <- y[-train_idx]

Treinando Decision Tree e Random Forest

# 5) treinar Decision Tree e Random Forest
train_df <- data.frame(inflacao = y_train, as.data.frame(X_train))
test_df  <- data.frame(inflacao = y_test,  as.data.frame(X_test))

tree_single <- rpart(inflacao ~ ., data = train_df, method = "anova", control = rpart.control(cp = 0.001))
pred_tree_single <- predict(tree_single, newdata = test_df)
rmse_tree_single <- sqrt(mean((y_test - pred_tree_single)^2))
r2_tree_single <- 1 - sum((y_test - pred_tree_single)^2) / sum((y_test - mean(y_test))^2)

rf_model <- randomForest(x = X_train, y = y_train, ntree = 500, importance = TRUE)
pred_rf <- predict(rf_model, newdata = X_test)
rmse_rf <- sqrt(mean((y_test - pred_rf)^2))
r2_rf <- 1 - sum((y_test - pred_rf)^2) / sum((y_test - mean(y_test))^2)

Calculando Decision Tree e Random Forest

# garantir que X_train/X_test/y_train/y_test existem; caso contrário, recriar split a partir de X_scaled/y
if (!exists("X_train") || !exists("X_test") || !exists("y_train") || !exists("y_test")) {
  train_idx <- sample(1:nrow(X_scaled), 0.8 * nrow(X_scaled))
  X_train <- X_scaled[train_idx, ]
  X_test  <- X_scaled[-train_idx, ]
  y_train <- y[train_idx]
  y_test  <- y[-train_idx]
}

# preparar data.frames para rpart (rpart usa fórmula mais facilmente)
train_df <- data.frame(inflacao = y_train, as.data.frame(X_train))
test_df  <- data.frame(inflacao = y_test,  as.data.frame(X_test))

# 1) Decision Tree (árvore única)
tree_single <- rpart(inflacao ~ ., data = train_df, method = "anova", control = rpart.control(cp = 0.001))
pred_tree_single <- predict(tree_single, newdata = test_df)
rmse_tree_single <- sqrt(mean((y_test - pred_tree_single)^2))
r2_tree_single <- 1 - sum((y_test - pred_tree_single)^2) / sum((y_test - mean(y_test))^2)

# imprimir métricas da árvore
cat("Decision Tree - RMSE:", round(rmse_tree_single, 6), " R2:", round(r2_tree_single, 6), "\n")

# 2) Random Forest
# usar um número moderado de árvores e mtry padrão (sqrt do número de variáveis) para ser parcimonioso
rf_model <- randomForest(x = X_train, y = y_train, ntree = 500, importance = TRUE)
pred_rf <- predict(rf_model, newdata = X_test)
rmse_rf <- sqrt(mean((y_test - pred_rf)^2))
r2_rf <- 1 - sum((y_test - pred_rf)^2) / sum((y_test - mean(y_test))^2)

# imprimir métricas do RF
cat("Random Forest - RMSE:", round(rmse_rf, 6), " R2:", round(r2_rf, 6), "\n")

# (Opcional) importância das variáveis no RF
rf_imp <- importance(rf_model)
print(head(rf_imp[order(-rf_imp[, "IncNodePurity"]), , drop = FALSE], 10))

Resultados

results_comparison <- tribble(
 ~Modelo, ~RMSE, ~R2,
 "Decision Tree", rmse_tree_single, r2_tree_single,
 "Random Forest", rmse_rf, r2_rf,
 "XGBoost", rmse_xgb, r2_xgb
)
print(results_comparison)

Esse resultado reforça a impressão que foi obtida através do gráfico, de que um número incipiente de árvores de decisão era associado a um erro quadrático maior que com uma quantidade crescente de árvores. Ainda, sob a técnica XGBoost, em que cada árvore anterior é usada para melhorar a seguinte - em oposição ao Random Forest que calcula a média do erro de diversas árvores aleatórias - o erro médio é ainda menor.

Observação gráfica

# garantir nomes/tipos corretos e ordenar por RMSE
results_plot <- results_comparison %>%
  mutate(
    Modelo = factor(Modelo, levels = Modelo[order(RMSE)]),
    RMSE = as.numeric(RMSE),
    R2 = as.numeric(R2)
  )

# gráfico horizontal com rótulos e paleta
ggplot(results_plot, aes(x = Modelo, y = RMSE, fill = Modelo)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 4)), 
            hjust = -0.1, size = 3.5, color = "black") +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal(base_size = 12) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  ) +
  labs(
    title = "Comparação de RMSE entre Modelos",
    y = "RMSE"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15)))

Regressões Lineares

Modelo OLS

# OLS (lm) usando os dados já escalados
if (!exists("X_train") || !exists("X_test") || !exists("y_train") || !exists("y_test")) {
  stop("X_train/X_test/y_train/y_test não encontrados. Execute o pré-processamento e split antes.")
}

# transformar X_train/X_test em data.frames para usar com lm
train_df_ols <- data.frame(inflacao = y_train, as.data.frame(X_train))
test_df_ols  <- data.frame(inflacao = y_test,  as.data.frame(X_test))

set.seed(123)
ols_model <- lm(inflacao ~ ., data = train_df_ols)

pred_ols <- predict(ols_model, newdata = test_df_ols)

rmse_ols <- sqrt(mean((y_test - pred_ols)^2))
r2_ols   <- 1 - sum((y_test - pred_ols)^2) / sum((y_test - mean(y_test))^2)

cat("OLS - RMSE:", round(rmse_ols, 6), " R2:", round(r2_ols, 6), "\n")

Modelo Ridge

# garantir glmnet carregado
library(glmnet)

# glmnet espera matrizes; X_train/X_test já são matrizes escaladas
if (!is.matrix(X_train)) X_train_mat <- as.matrix(X_train) else X_train_mat <- X_train
if (!is.matrix(X_test))  X_test_mat  <- as.matrix(X_test)  else X_test_mat  <- X_test

set.seed(123)
cv_ridge <- cv.glmnet(x = X_train_mat, y = y_train, alpha = 0, nfolds = 5, type.measure = "mse")

lambda_ridge <- cv_ridge$lambda.min

ridge_model <- glmnet(x = X_train_mat, y = y_train, alpha = 0, lambda = lambda_ridge)

pred_ridge <- predict(ridge_model, newx = X_test_mat, s = lambda_ridge)
pred_ridge <- as.numeric(pred_ridge)

rmse_ridge <- sqrt(mean((y_test - pred_ridge)^2))
r2_ridge   <- 1 - sum((y_test - pred_ridge)^2) / sum((y_test - mean(y_test))^2)

cat("Ridge - lambda:", signif(lambda_ridge,4), " RMSE:", round(rmse_ridge, 6), " R2:", round(r2_ridge, 6), "\n")

Modelo Lasso

set.seed(123)
cv_lasso <- cv.glmnet(x = X_train_mat, y = y_train, alpha = 1, nfolds = 5, type.measure = "mse")

lambda_lasso <- cv_lasso$lambda.min

lasso_model <- glmnet(x = X_train_mat, y = y_train, alpha = 1, lambda = lambda_lasso)

pred_lasso <- predict(lasso_model, newx = X_test_mat, s = lambda_lasso)
pred_lasso <- as.numeric(pred_lasso)

rmse_lasso <- sqrt(mean((y_test - pred_lasso)^2))
r2_lasso   <- 1 - sum((y_test - pred_lasso)^2) / sum((y_test - mean(y_test))^2)

cat("Lasso - lambda:", signif(lambda_lasso,4), " RMSE:", round(rmse_lasso, 6), " R2:", round(r2_lasso, 6), "\n")

# coeficientes não nulos (variáveis selecionadas)
coef_lasso <- as.matrix(coef(lasso_model, s = lambda_lasso))
nonzero <- coef_lasso[coef_lasso[,1] != 0, , drop = FALSE]
nonzero <- data.frame(Feature = rownames(nonzero), Coef = nonzero[,1], row.names = NULL)
nonzero <- nonzero[order(-abs(nonzero$Coef)), ]
print(head(nonzero, 20))

A estimação via Lasso, que subjulga o poder de predição de diversas variáveis para contextos mais realistas, representa síntese próxima do resultado obtido através do gradient boostind. Novamente variação do câmbio aparece como variável mais importante em termos de explicação, mas com a diferença de que outras variáveis de maior dificuldade interpretativa também apareceram, como output gap e credibilidade.

Agregação de Resultados

# garantir que métricas dos outros modelos existem
required_vars <- c("rmse_tree_single","r2_tree_single","rmse_rf","r2_rf","rmse_xgb","r2_xgb")
missing_vars <- required_vars[!sapply(required_vars, exists)]
if (length(missing_vars) > 0) {
  stop("Faltam métricas de alguns modelos: ", paste(missing_vars, collapse = ", "))
}

# criar/atualizar tabela de resultados
results_comparison <- tribble(
 ~Modelo, ~RMSE, ~R2,
 "Decision Tree", rmse_tree_single, r2_tree_single,
 "Random Forest", rmse_rf, r2_rf,
 "XGBoost", rmse_xgb, r2_xgb,
 "OLS", rmse_ols, r2_ols,
 "Ridge", rmse_ridge, r2_ridge,
 "Lasso", rmse_lasso, r2_lasso
)

print(results_comparison)

# gráfico ordenado por RMSE
results_plot <- results_comparison %>%
  mutate(
    Modelo = factor(Modelo, levels = Modelo[order(RMSE)]),
    RMSE = as.numeric(RMSE),
    R2 = as.numeric(R2)
  )

library(ggplot2)
ggplot(results_plot, aes(x = Modelo, y = RMSE, fill = Modelo)) +
  geom_col(width = 0.6, show.legend = FALSE) +
  geom_text(aes(label = round(RMSE, 4)), 
            hjust = -0.1, size = 3.5, color = "black") +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  theme_minimal(base_size = 12) +
  theme(
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold"),
    panel.grid.major.y = element_blank()
  ) +
  labs(
    title = "Comparação de RMSE entre Modelos (incluindo OLS/Ridge/Lasso)",
    y = "RMSE"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15)))

Conclusão

Para esta base de dados - marcada por relações lineares - os modelos de maior complexidade associados a árvores de decisão aparecem como mais ineficientes. Tal resultado podia ser esperado, porquanto podia-se ter em vista, logo de imediato, que esses modelos proporcionam maior confiabilidade em contextos não lineares, portanto diferentes do nosso.

Ainda assim, ter a lista de variáveis ordenadas pelo grau de importância é um dado importante e pode fornecer boas sugestões para interpretações econômicas ao economista que analisa tal questão. De todo modo, o ponto de maior estranheza trata-se do erro menor associado ao modelo mais simples de todos, o OLS. Contudo, isto não nos parece enfraquecer a relevância econométrica dos modelos Ridge e Lasso, visto que permitem, no primeiro caso a seleção das variáveis mais importantes, e no segundo caso, a parcimônia dos coeficientes dessas variáveis.