Ejercicio de Pronóstico Simulación

options(scipen = 9999)
set.seed(50)

library(foreign)
library(lmtest)
library(sandwich)
library(robustbase)
library(stargazer)

library(haven)
crime <- read_dta("C:/Users/Rebeca/Downloads/crime.dta")
View(crime)

cat("--- PARTE 1: SIMULACIÓN DE PRONÓSTICO (MODELO BASE) ---\n")

## --- PARTE 1: SIMULACIÓN DE PRONÓSTICO (MODELO BASE) ---

modelo_crime <- lm(crime ~ poverty + single, data = crime)
print(summary(modelo_crime))

## 
## Call:
## lm(formula = crime ~ poverty + single, data = crime)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -811.14 -114.27  -22.44  121.86  689.82 
## 
## Coefficients:
##              Estimate Std. Error t value        Pr(>|t|)    
## (Intercept) -1368.189    187.205  -7.308 0.0000000024786 ***
## poverty         6.787      8.989   0.755           0.454    
## single        166.373     19.423   8.566 0.0000000000312 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 243.6 on 48 degrees of freedom
## Multiple R-squared:  0.7072, Adjusted R-squared:  0.695 
## F-statistic: 57.96 on 2 and 48 DF,  p-value: 0.0000000000001578

nuevo_escenario <- data.frame(
  poverty = mean(crime$poverty, na.rm = TRUE),
  single  = mean(crime$single, na.rm = TRUE)
)

pronostico_tabla <- predict(modelo_crime, newdata = nuevo_escenario, interval = "prediction")
print(pronostico_tabla)

##        fit      lwr      upr
## 1 612.8431 118.2534 1107.433

PARTE 2: Simulación de Pronóstico con Enfoque de Validación Cruzada (HAC)

cat("\n--- PARTE 2: SIMULACIÓN ENFOQUE TRAIN/TEST (500 RÉPLICAS) ---\n")

## 
## --- PARTE 2: SIMULACIÓN ENFOQUE TRAIN/TEST (500 RÉPLICAS) ---

n_replicas <- 500
proporcion_train <- 0.75
n_total <- nrow(crime)
n_train <- round(n_total * proporcion_train)
# Vectores para almacenar las métricas de desempeño
rmse_mecanico <- numeric(n_replicas)
mae_mecanico  <- numeric(n_replicas)
# Bucle de simulación (Monte Carlo / Validación Cruzada Aleatoria)
for(i in 1:n_replicas) {# 1. División aleatoria de datos (Train 75% / Test 25%)
  indices_train <- sample(seq_len(n_total), size = n_train)
  datos_train   <- crime[indices_train, ]
  datos_test    <- crime[-indices_train, ]
  
  # 2. Ajuste del modelo con datos de entrenamiento
  modelo_train <- lm(crime ~ poverty + single, data = datos_train)
  
  # 3. Predicción sobre los datos "desconocidos" (Test)
  predicciones <- predict(modelo_train, newdata = datos_test)
  
  # 4. Cálculo de errores
  errores <- datos_test$crime - predicciones
  
  # 5. Guardar métricas de desempeño para la réplica i
  rmse_mecanico[i] <- sqrt(mean(errores^2, na.rm = TRUE))
  mae_mecanico[i]  <- mean(abs(errores), na.rm = TRUE)
}

ANÁLISIS DE RESULTADOS DESEMPEÑO DEL MODELO

cat("ANÁLISIS DE DESEMPEÑO EN DATOS DESCONOCIDOS (TEST)\n")

## ANÁLISIS DE DESEMPEÑO EN DATOS DESCONOCIDOS (TEST)

resultados_desempeno <- data.frame(
  Metrica = c("RMSE (Error Cuadrático Medio Raíz)", "MAE (Error Absoluto Medio)"),
  Media   = c(mean(rmse_mecanico), mean(mae_mecanico)),
  Mediana = c(median(rmse_mecanico), median(mae_mecanico)),
  Minimo  = c(min(rmse_mecanico), min(mae_mecanico)),
  Maximo  = c(max(rmse_mecanico), max(mae_mecanico)),
  Desv_Est = c(sd(rmse_mecanico), sd(mae_mecanico))
)

print(resultados_desempeno, row.names = FALSE)

##                             Metrica    Media  Mediana    Minimo   Maximo
##  RMSE (Error Cuadrático Medio Raíz) 275.0489 284.3751 130.25144 429.0787
##          MAE (Error Absoluto Medio) 200.3083 198.3948  97.97605 325.7022
##  Desv_Est
##  76.57143
##  47.24460

# Histograma de los errores de predicción (RMSE) distribuidos en las 500 réplicas
hist(rmse_mecanico, 
     main = "Distribución del RMSE en 500 Réplicas de Simulación",
     xlab = "RMSE (Datos de Validación)", 
     col  = "skyblue", 
     border = "white")
abline(v = mean(rmse_mecanico), col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("RMSE Medio"), col = c("red"), lty = 2, lwd = 2)