Estimaci´on de la Matriz de Varianzas y Covarianzas con Bootstrap para regresión Cuantilica

knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

# Cargar librerías necesarias
library(quantreg)  # Para regresión cuantil
data(engel)
# Definir la cantidad de réplicas bootstrap
B <- 100000  # Número de muestras bootstrap
n <- nrow(engel)  # Tamaño de la muestra original

# Inicializar matriz para almacenar los coeficientes estimados
beta_boot <- matrix(NA, nrow = B, ncol = 2)  # Dos columnas para Intercepto y Beta1

# Estimación original del modelo cuantil para tau = 0.5
modelo_original <- rq(foodexp ~ income, data = engel, tau = 0.5)
beta_hat <- coef(modelo_original)  # Coeficientes originales
summary(modelo_original,se = "boot", bsmethod= "xy", R=100000)

## 
## Call: rq(formula = foodexp ~ income, tau = 0.5, data = engel)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 81.48225 27.22146    2.99331  0.00306
## income       0.56018  0.03485   16.07430  0.00000

# Bootstrap
# set.seed(123)  # Fijar semilla para reproducibilidad
for (b in 1:B) {
  # Seleccionar muestra bootstrap con reemplazo
  sample_indices <- sample(1:n, size = n, replace = TRUE)
  engel_boot <- engel[sample_indices, ]  # Nueva muestra bootstrap

  # Ajustar modelo de regresión cuantil en la muestra bootstrap
  modelo_boot <- rq(foodexp ~ income, data = engel_boot, tau = 0.5)
  
  # Almacenar los coeficientes estimados
  beta_boot[b, ] <- coef(modelo_boot)
}

# Cálculo de la varianza bootstrap
var_bootstrap <- apply(beta_boot, 2, var)
desv_est_bootstrap <- sqrt(var_bootstrap)  # Desviación estándar bootstrap

# Resultados
cat("Estimación original del modelo:\n")

## Estimación original del modelo:

print(beta_hat)

## (Intercept)      income 
##  81.4822474   0.5601806

cat("\nVarianza bootstrap de los coeficientes:\n")

## 
## Varianza bootstrap de los coeficientes:

print(var_bootstrap)

## [1] 7.410242e+02 1.216356e-03

cat("\nDesviación estándar bootstrap de los coeficientes:\n")

## 
## Desviación estándar bootstrap de los coeficientes:

print(desv_est_bootstrap)

## [1] 27.22175996  0.03487629

# Histograma de las distribuciones bootstrap de los coeficientes
par(mfrow = c(1,2))  # Para visualizar dos gráficos en la misma ventana
hist(beta_boot[, 1], col = "lightblue", main = "Distribución Bootstrap del Intercepto", xlab = expression(beta[0]))
hist(beta_boot[, 2], col = "lightcoral", main = "Distribución Bootstrap de Beta_1", xlab = expression(beta[1]))

81.48225/27.26370

## [1] 2.988672

0.56018/0.03493

## [1] 16.03722

# Intervalos de confianza al 25% y 95% para los coeficientes
#IC_25_beta0 <- quantile(beta_boot[, 1], probs = c(0.125, 0.875))  # 25% = 12.5% y 87.5%
IC_95_beta0 <- quantile(beta_boot[, 1], probs = c(0.025, 0.975))  # 95% = 2.5% y 97.5%

#IC_25_beta1 <- quantile(beta_boot[, 2], probs = c(0.125, 0.875))  
IC_95_beta1 <- quantile(beta_boot[, 2], probs = c(0.025, 0.975))  

# Imprimir los resultados
#cat("\nIntervalo de confianza empírico (25%) para Beta_0:", IC_25_beta0, "\n")
cat("Intervalo de confianza empírico (95%) para Beta_0:", IC_95_beta0, "\n")

## Intervalo de confianza empírico (95%) para Beta_0: 41.2005 150.5735

#cat("\nIntervalo de confianza empírico (25%) para Beta_1:", IC_25_beta1, "\n")
cat("Intervalo de confianza empírico (95%) para Beta_1:", IC_95_beta1, "\n")

## Intervalo de confianza empírico (95%) para Beta_1: 0.4702404 0.6137283

summary.rq(modelo_original)

## 
## Call: rq(formula = foodexp ~ income, tau = 0.5, data = engel)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd  upper bd 
## (Intercept)  81.48225     53.25915 114.01156
## income        0.56018      0.48702   0.60199