Ejercicio de Pronóstico Simulación
Kevin Alexander Neira García
2026-05-27
# --- 0. Cargar datos y librerías ---
load("smoke.RData") # carga los objetos: data, desc
library(stargazer)##
## Please cite as:## Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.## R package version 5.2.3. https://CRAN.R-project.org/package=stargazerlibrary(sandwich)
library(lmtest)## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericprint(desc)## variable label
## 1 educ years of schooling
## 2 cigpric state cig. price, cents/pack
## 3 white =1 if white
## 4 age in years
## 5 income annual income, $
## 6 cigs cigs. smoked per day
## 7 restaurn =1 if rest. smk. restrictions
## 8 lincome log(income)
## 9 agesq age^2
## 10 lcigpric log(cigprice)head(data)## educ cigpric white age income cigs restaurn lincome agesq lcigpric
## 1 16.0 60.506 1 46 20000 0 0 9.903487 2116 4.102743
## 2 16.0 57.883 1 40 30000 0 0 10.308952 1600 4.058424
## 3 12.0 57.664 1 58 30000 3 0 10.308952 3364 4.054633
## 4 13.5 57.883 1 30 20000 0 0 9.903487 900 4.058424
## 5 10.0 58.320 1 17 20000 0 0 9.903487 289 4.065945
## 6 6.0 59.340 1 86 6500 0 0 8.779557 7396 4.083283ESTIMAR MODELO
# cigs = f(cigpric, lcigpric, income, lincome, age, agesq, educ, white, restaurn)
# ============================================================
modelo_original <- lm(cigs ~ cigpric + lcigpric + income + lincome +
age + agesq + educ + white + restaurn,
data = data)RESUMEN DEL MODELO CON STARGAZER
stargazer(modelo_original,
type = "text",
title = "Modelo Original: Determinantes del consumo de cigarrillos",
digits = 4)##
## Modelo Original: Determinantes del consumo de cigarrillos
## ===============================================
## Dependent variable:
## ---------------------------
## cigs
## -----------------------------------------------
## cigpric 2.0023
## (1.4928)
##
## lcigpric -115.2735
## (85.4243)
##
## income -0.00005
## (0.0001)
##
## lincome 1.4041
## (1.7082)
##
## age 0.7784***
## (0.1606)
##
## agesq -0.0092***
## (0.0017)
##
## educ -0.4948***
## (0.1682)
##
## white -0.5311
## (1.4607)
##
## restaurn -2.6442**
## (1.1300)
##
## Constant 340.8044
## (260.0156)
##
## -----------------------------------------------
## Observations 807
## R2 0.0552
## Adjusted R2 0.0445
## Residual Std. Error 13.4128 (df = 797)
## F Statistic 5.1693*** (df = 9; 797)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01VERIFICACION MATRIZ-COVARIANZA
vcov_original <- vcov(modelo_original)
vcov_original## (Intercept) cigpric lcigpric income
## (Intercept) 6.760811e+04 3.861092e+02 -2.217382e+04 4.361954e-04
## cigpric 3.861092e+02 2.228545e+00 -1.272304e+02 -6.112577e-06
## lcigpric -2.217382e+04 -1.272304e+02 7.297314e+03 3.973538e-04
## income 4.361954e-04 -6.112577e-06 3.973538e-04 1.781986e-08
## lincome -8.178466e+00 7.124578e-02 -4.860425e+00 -2.062280e-04
## age 4.955739e-01 3.684842e-03 -2.196985e-01 -2.490054e-07
## agesq -1.237179e-02 -8.003706e-05 4.619639e-03 4.523652e-09
## educ 9.421945e-02 1.549304e-03 -1.048653e-01 -2.445251e-06
## white -2.254192e+00 1.177281e-02 -3.119456e-01 -7.733768e-06
## restaurn 4.064836e+01 2.186489e-01 -1.340589e+01 -1.150501e-05
## lincome age agesq educ
## (Intercept) -8.1784660668 4.955739e-01 -1.237179e-02 9.421945e-02
## cigpric 0.0712457770 3.684842e-03 -8.003706e-05 1.549304e-03
## lcigpric -4.8604250078 -2.196985e-01 4.619639e-03 -1.048653e-01
## income -0.0002062280 -2.490054e-07 4.523652e-09 -2.445251e-06
## lincome 2.9178305410 -2.864418e-02 2.930256e-04 -2.617669e-03
## age -0.0286441758 2.577810e-02 -2.764371e-04 -2.949036e-03
## agesq 0.0002930256 -2.764371e-04 3.060024e-06 4.144658e-05
## educ -0.0026176694 -2.949036e-03 4.144658e-05 2.828458e-02
## white 0.1323738372 -1.438531e-02 1.721940e-04 2.035627e-03
## restaurn 0.0828693504 3.811336e-03 -3.412736e-05 -3.821198e-03
## white restaurn
## (Intercept) -2.254192e+00 4.064836e+01
## cigpric 1.177281e-02 2.186489e-01
## lcigpric -3.119456e-01 -1.340589e+01
## income -7.733768e-06 -1.150501e-05
## lincome 1.323738e-01 8.286935e-02
## age -1.438531e-02 3.811336e-03
## agesq 1.721940e-04 -3.412736e-05
## educ 2.035627e-03 -3.821198e-03
## white 2.133708e+00 1.602894e-01
## restaurn 1.602894e-01 1.276897e+00is.scalar <- is.vector(vcov_original) || length(vcov_original) == 1
is.scalar## [1] FALSEESTIMADOR HAC
vcov_HAC <- vcovHC(modelo_original, type = "HC1")
# Errores estandar HAC
se_HAC <- sqrt(diag(vcov_HAC))
cat("\n=== Errores Estandar HAC (HC1) ===\n")##
## === Errores Estandar HAC (HC1) ===print(se_HAC)## (Intercept) cigpric lcigpric income lincome age
## 2.803072e+02 1.612751e+00 9.191569e+01 1.163156e-04 1.236657e+00 1.378011e-01
## agesq educ white restaurn
## 1.460626e-03 1.639876e-01 1.370425e+00 1.044748e+00TABLA COMPARATIVA: MODELO ORIGINAL vs MODELO CORREGIDO
# Errores estandar originales (OLS)
se_OLS <- sqrt(diag(vcov(modelo_original)))
stargazer(modelo_original, modelo_original,
type = "text",
title = "Comparacion: Modelo Original vs Corregido con HAC",
column.labels = c("Original (OLS)", "Corregido (HAC)"),
se = list(se_OLS, se_HAC),
dep.var.labels = "Cigarrillos consumidos al dia (cigs)",
covariate.labels = c("Precio cigarrillos (cigpric)",
"Log precio cigarrillos (lcigpric)",
"Ingreso (income)",
"Log ingreso (lincome)",
"Edad (age)",
"Edad al cuadrado (agesq)",
"Educacion (educ)",
"Raza blanca (white)",
"Restaurante libre de humo (restaurn)"),
omit.stat = c("LL", "ser"),
digits = 4,
notes = c("Errores estandar en parentesis.",
"Columna (2): errores robustos HC1 (Huber-White).",
"*p<0.1; **p<0.05; ***p<0.01"))##
## Comparacion: Modelo Original vs Corregido con HAC
## ======================================================================================
## Dependent variable:
## -------------------------------------------------
## Cigarrillos consumidos al dia (cigs)
## Original (OLS) Corregido (HAC)
## (1) (2)
## --------------------------------------------------------------------------------------
## Precio cigarrillos (cigpric) 2.0023 2.0023
## (1.4928) (1.6128)
##
## Log precio cigarrillos (lcigpric) -115.2735 -115.2735
## (85.4243) (91.9157)
##
## Ingreso (income) -0.00005 -0.00005
## (0.0001) (0.0001)
##
## Log ingreso (lincome) 1.4041 1.4041
## (1.7082) (1.2367)
##
## Edad (age) 0.7784*** 0.7784***
## (0.1606) (0.1378)
##
## Edad al cuadrado (agesq) -0.0092*** -0.0092***
## (0.0017) (0.0015)
##
## Educacion (educ) -0.4948*** -0.4948***
## (0.1682) (0.1640)
##
## Raza blanca (white) -0.5311 -0.5311
## (1.4607) (1.3704)
##
## Restaurante libre de humo (restaurn) -2.6442** -2.6442**
## (1.1300) (1.0447)
##
## Constant 340.8044 340.8044
## (260.0156) (280.3072)
##
## --------------------------------------------------------------------------------------
## Observations 807 807
## R2 0.0552 0.0552
## Adjusted R2 0.0445 0.0445
## F Statistic (df = 9; 797) 5.1693*** 5.1693***
## ======================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01
## Errores estandar en parentesis.
## Columna (2): errores robustos HC1 (Huber-White).
## *p<0.1; **p<0.05; ***p<0.01PARTE 2: SIMULACIÓN DE PRONÓSTICO CON DATOS SMOKE
Descripción de la simulación
Se realizan 500 réplicas con semilla
50. En cada réplica se divide aleatoriamente el dataset
en:
- 75% datos de entrenamiento (training): se estima el modelo OLS
- 25% datos de prueba (testing / “desconocidos”): se evalúa el desempeño predictivo
Las métricas de desempeño calculadas en cada réplica son:
- RMSE (Root Mean Squared Error): raíz del error cuadrático medio
- MAE (Mean Absolute Error): error absoluto medio
- R² out-of-sample: coeficiente de determinación en datos de prueba
# ============================================================
# SIMULACIÓN DE PRONÓSTICO - DATOS SMOKE
# 75% entrenamiento | 25% prueba | 500 réplicas | semilla 50
# ============================================================
set.seed(50)
n <- nrow(data) # total de observaciones
n_train <- floor(0.75 * n) # tamaño del conjunto de entrenamiento
replicas <- 500 # número de réplicas
# Vectores para almacenar métricas de cada réplica
rmse_vec <- numeric(replicas)
mae_vec <- numeric(replicas)
r2_vec <- numeric(replicas)
for (i in 1:replicas) {
# --- Dividir datos aleatoriamente ---
idx_train <- sample(1:n, size = n_train, replace = FALSE)
datos_train <- data[idx_train, ] # datos de entrenamiento (conocidos)
datos_test <- data[-idx_train, ] # datos de prueba (desconocidos)
# --- Estimar modelo con datos de entrenamiento ---
modelo_sim <- lm(cigs ~ cigpric + lcigpric + income + lincome +
age + agesq + educ + white + restaurn,
data = datos_train)
# --- Predecir en datos desconocidos ---
y_pred <- predict(modelo_sim, newdata = datos_test)
y_real <- datos_test$cigs
# --- Calcular métricas de desempeño ---
residuos <- y_real - y_pred
rmse_vec[i] <- sqrt(mean(residuos^2))
mae_vec[i] <- mean(abs(residuos))
# R² out-of-sample: 1 - SS_res / SS_tot
ss_res <- sum(residuos^2)
ss_tot <- sum((y_real - mean(y_real))^2)
r2_vec[i] <- 1 - ss_res / ss_tot
}Resultados del desempeño en datos “desconocidos”
# ============================================================
# RESUMEN DE MÉTRICAS (promedio de las 500 réplicas)
# ============================================================
cat("========================================================\n")## ========================================================cat(" DESEMPEÑO DEL MODELO EN DATOS DESCONOCIDOS (TEST)\n")## DESEMPEÑO DEL MODELO EN DATOS DESCONOCIDOS (TEST)cat(" 500 réplicas | 75% entrenamiento | 25% prueba\n")## 500 réplicas | 75% entrenamiento | 25% pruebacat("========================================================\n\n")## ========================================================cat(sprintf(" RMSE promedio : %.4f cigarrillos/día\n", mean(rmse_vec)))## RMSE promedio : 13.4892 cigarrillos/díacat(sprintf(" MAE promedio : %.4f cigarrillos/día\n", mean(mae_vec)))## MAE promedio : 10.7230 cigarrillos/díacat(sprintf(" R² promedio : %.4f\n\n", mean(r2_vec)))## R² promedio : 0.0279cat(sprintf(" RMSE (mediana): %.4f\n", median(rmse_vec)))## RMSE (mediana): 13.4844cat(sprintf(" MAE (mediana): %.4f\n", median(mae_vec)))## MAE (mediana): 10.7136cat(sprintf(" R² (mediana): %.4f\n", median(r2_vec)))## R² (mediana): 0.0308cat("========================================================\n")## ========================================================Distribución de las métricas (histogramas)
# ============================================================
# VISUALIZACIÓN: distribución de RMSE, MAE y R² en 500 réplicas
# ============================================================
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
# RMSE
hist(rmse_vec,
main = "Distribución del RMSE\n(500 réplicas)",
xlab = "RMSE (cigarrillos/día)",
col = "#4472C4", border = "white",
breaks = 30)
abline(v = mean(rmse_vec), col = "red", lwd = 2, lty = 2)
legend("topright", legend = paste("Media =", round(mean(rmse_vec), 3)),
col = "red", lty = 2, lwd = 2, bty = "n", cex = 0.85)
# MAE
hist(mae_vec,
main = "Distribución del MAE\n(500 réplicas)",
xlab = "MAE (cigarrillos/día)",
col = "#ED7D31", border = "white",
breaks = 30)
abline(v = mean(mae_vec), col = "red", lwd = 2, lty = 2)
legend("topright", legend = paste("Media =", round(mean(mae_vec), 3)),
col = "red", lty = 2, lwd = 2, bty = "n", cex = 0.85)
# R²
hist(r2_vec,
main = "Distribución del R² out-of-sample\n(500 réplicas)",
xlab = "R²",
col = "#70AD47", border = "white",
breaks = 30)
abline(v = mean(r2_vec), col = "red", lwd = 2, lty = 2)
legend("topright", legend = paste("Media =", round(mean(r2_vec), 3)),
col = "red", lty = 2, lwd = 2, bty = "n", cex = 0.85)
par(mfrow = c(1, 1))Interpretación de resultados
cat("**Análisis del desempeño predictivo:**\n\n")Análisis del desempeño predictivo:
cat("- El **RMSE promedio** indica el error típico de predicción del modelo en datos no vistos.\n")- El RMSE promedio indica el error típico de predicción del modelo en datos no vistos.
cat(" Un valor bajo sugiere buena capacidad predictiva; un valor alto indica que el modelo\n")Un valor bajo sugiere buena capacidad predictiva; un valor alto indica que el modelo
cat(" tiene dificultades para generalizar.\n\n")tiene dificultades para generalizar.
cat("- El **MAE promedio** es más robusto ante valores atípicos que el RMSE. Si MAE << RMSE,\n")- El MAE promedio es más robusto ante valores atípicos que el RMSE. Si MAE << RMSE,
cat(" hay presencia de observaciones con errores de predicción muy grandes.\n\n")hay presencia de observaciones con errores de predicción muy grandes.
cat("- El **R² out-of-sample** mide qué tan bien explica el modelo la variación en datos\n")- El R² out-of-sample mide qué tan bien explica el modelo la variación en datos
cat(" desconocidos. Valores cercanos a 1 indican buen ajuste predictivo;\n")desconocidos. Valores cercanos a 1 indican buen ajuste predictivo;
cat(" valores negativos indican que el modelo predice peor que la media simple.\n\n")valores negativos indican que el modelo predice peor que la media simple.
cat("- La **variabilidad entre réplicas** (histogramas) refleja la sensibilidad del modelo\n")- La variabilidad entre réplicas (histogramas) refleja la sensibilidad del modelo
cat(" a la partición específica de los datos. Distribuciones más estrechas indican\n")a la partición específica de los datos. Distribuciones más estrechas indican
cat(" mayor estabilidad del modelo.\n")mayor estabilidad del modelo.