options(scipen = 9999)
library(stargazer)
library(lmtest)
library(sandwich)
library(car)
load("smoke.RData")
### 1. Estimación del Modelo Lineal
modelo_original <- lm(cigs ~ cigpric + lcigpric + income + lincome + age +
agesq + educ + white + restaurn, data = data)
print("--- RESUMEN DEL MODELO ORIGINAL ---")
## [1] "--- RESUMEN DEL MODELO ORIGINAL ---"
stargazer(modelo_original, type = "text", title = "Resumen del Modelo Lineal Original")
##
## Resumen del Modelo Lineal Original
## ===============================================
## Dependent variable:
## ---------------------------
## cigs
## -----------------------------------------------
## cigpric 2.002
## (1.493)
##
## lcigpric -115.273
## (85.424)
##
## income -0.00005
## (0.0001)
##
## lincome 1.404
## (1.708)
##
## age 0.778***
## (0.161)
##
## agesq -0.009***
## (0.002)
##
## educ -0.495***
## (0.168)
##
## white -0.531
## (1.461)
##
## restaurn -2.644**
## (1.130)
##
## Constant 340.804
## (260.016)
##
## -----------------------------------------------
## Observations 807
## R2 0.055
## Adjusted R2 0.044
## Residual Std. Error 13.413 (df = 797)
## F Statistic 5.169*** (df = 9; 797)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01
### 2. Verificación de la Matriz de Varianza-Covarianza (¿Es Escalar?)
print("--- PRUEBAS PARA VERIFICAR SI LA MATRIZ ES ESCALAR ---")
## [1] "--- PRUEBAS PARA VERIFICAR SI LA MATRIZ ES ESCALAR ---"
# Para que la matriz sea escalar no debe haber Heterocedasticidad ni Autocorrelación.
# A) Prueba de White / Breusch-Pagan para Heterocedasticidad
white_test <- bptest(modelo_original, data = data)
print(white_test)
##
## studentized Breusch-Pagan test
##
## data: modelo_original
## BP = 33.525, df = 9, p-value = 0.0001082
# Si p-value < 0.05, hay evidencia de Heterocedasticidad.
# B) Prueba de Autocorrelación de 1er Orden (Durbin-Watson)
dw_test <- durbinWatsonTest(modelo_original)
print(dw_test)
## lag Autocorrelation D-W Statistic p-value
## 1 -0.009243664 2.017442 0.92
## Alternative hypothesis: rho != 0
# C) Prueba del Multiplicador de Lagrange (Breusch-Godfrey) - Autocorrelación de 2° Orden
bg_test <- bgtest(modelo_original, order = 2)
print(bg_test)
##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_original
## LM test = 0.26889, df = 2, p-value = 0.8742
# Si p-value < 0.05 en estas pruebas, existe autocorrelación
# CONCLUSIÓN SOBRE LA MATRIZ ESCALAR:
# "Si cualquiera de las pruebas anteriores rechaza la hipótesis nula (p-value < 0.05),
# la matriz de Varianza-Covarianza NO es escalar.
# La fuente del problema será la Heterocedasticidad (si falla White) y/o la Autocorrelación (si falla BG/DW)."
### 3. Obtención del Estimador Robustecido (HAC)
print("--- CÁLCULO DEL ESTIMADOR ROBUSTO HAC ---")
## [1] "--- CÁLCULO DEL ESTIMADOR ROBUSTO HAC ---"
vcov_HAC <- vcovHAC(modelo_original)
### 4. Presentación de Modelos (Original vs Corregido)
print("--- TABLA COMPARATIVA FINAL (STARGAZER) ---")
## [1] "--- TABLA COMPARATIVA FINAL (STARGAZER) ---"
stargazer(modelo_original, modelo_original,
se = list(NULL, sqrt(diag(vcov_HAC))),
type = "text",
column.labels = c("Modelo Original", "Modelo Corregido (HAC)"),
title = "Comparación de Modelos: Errores Estándar Originales vs Robustos HAC")
##
## Comparación de Modelos: Errores Estándar Originales vs Robustos HAC
## =====================================================================
## Dependent variable:
## --------------------------------------
## cigs
## Modelo Original Modelo Corregido (HAC)
## (1) (2)
## ---------------------------------------------------------------------
## cigpric 2.002 2.002
## (1.493) (1.483)
##
## lcigpric -115.273 -115.273
## (85.424) (83.865)
##
## income -0.00005 -0.00005
## (0.0001) (0.0001)
##
## lincome 1.404 1.404
## (1.708) (1.293)
##
## age 0.778*** 0.778***
## (0.161) (0.139)
##
## agesq -0.009*** -0.009***
## (0.002) (0.001)
##
## educ -0.495*** -0.495***
## (0.168) (0.169)
##
## white -0.531 -0.531
## (1.461) (1.323)
##
## restaurn -2.644** -2.644**
## (1.130) (1.043)
##
## Constant 340.804 340.804
## (260.016) (255.948)
##
## ---------------------------------------------------------------------
## Observations 807 807
## R2 0.055 0.055
## Adjusted R2 0.044 0.044
## Residual Std. Error (df = 797) 13.413 13.413
## F Statistic (df = 9; 797) 5.169*** 5.169***
## =====================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01