Aplicacion de libreria “stargazer”.

library(lmtest)
library(car)
library(sandwich)
library(stargazer)

# Data
load("smoke.RData")

Estimacion de modelo lineal.

modelo <- lm(cigs ~ cigpric + lcigpric + income + lincome +
               age + agesq + educ + white + restaurn,
             data = data)

# Resumen de modelo con stargazer. 

stargazer(
  modelo,
  type = "html",
  title = "Resumen del Modelo Lineal",
  digits = 4
)
Resumen del Modelo Lineal
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.01

Verificar si la matriz Varianza-Covarianza es escalar.

# Matriz de varianza-covarianza usual
vcov_original <- vcov(modelo)
diag_vals <- diag(vcov_original)

if(all(abs(vcov_original - diag(diag_vals)) < 1e-8)){
  print("La matriz es diagonal (aprox. escalar)")
} else {
  print("La matriz NO es escalar")
}

[1] “La matriz NO es escalar”

# Verificación formal de heterocedasticidad
bp_test <- bptest(modelo)
bp_test
studentized Breusch-Pagan test

data: modelo BP = 33.525, df = 9, p-value = 0.0001082

# Verificación de autocorrelación 

bg_test<-bgtest(modelo)
bg_test
Breusch-Godfrey test for serial correlation of order up to 1

data: modelo LM test = 0.069737, df = 1, p-value = 0.7917

# Interpretación:
## - Si Breusch-Pagan rechaza H0 -> heterocedasticidad
## - Si Durbin-Watson rechaza H0 -> autocorrelación
## Estas son las fuentes por las que la matriz no sería escalar

Obtener estimador HAC.

vcov_HAC <- NeweyWest(modelo, prewhite = FALSE)
knitr::kable(
  data.frame(
    Metodo = "Newey-West HAC",
    Variables = ncol(model.matrix(modelo)),
    Observaciones = nobs(modelo)
  ),
  caption = "Resumen del estimador HAC"
)
Resumen del estimador HAC
Metodo Variables Observaciones
Newey-West HAC 10 807

Tabla comparativa modelo original vs corregido

# Errores estándar originales
se_original <- sqrt(diag(vcov(modelo)))

# Errores estándar HAC
se_HAC <- sqrt(diag(vcov_HAC))

# Tabla comparativa.

stargazer(
  modelo,
  modelo,
  type = "html",
  title = "Modelo Original vs Modelo Corregido HAC",

  column.labels = c("Original", "Corregido HAC"),

  se = list(
    sqrt(diag(vcov(modelo))),
    sqrt(diag(vcov_HAC))
  ), digits = 4)
Modelo Original vs Modelo Corregido HAC
Dependent variable:
cigs
Original Corregido HAC
(1) (2)
cigpric 2.0023 2.0023
(1.4928) (1.4770)
lcigpric -115.2735 -115.2735
(85.4243) (83.6700)
income -0.00005 -0.00005
(0.0001) (0.0001)
lincome 1.4041 1.4041
(1.7082) (1.3277)
age 0.7784*** 0.7784***
(0.1606) (0.1348)
agesq -0.0092*** -0.0092***
(0.0017) (0.0014)
educ -0.4948*** -0.4948***
(0.1682) (0.1698)
white -0.5311 -0.5311
(1.4607) (1.2912)
restaurn -2.6442** -2.6442**
(1.1300) (1.0501)
Constant 340.8044 340.8044
(260.0156) (256.4787)
Observations 807 807
R2 0.0552 0.0552
Adjusted R2 0.0445 0.0445
Residual Std. Error (df = 797) 13.4128 13.4128
F Statistic (df = 9; 797) 5.1693*** 5.1693***
Note: p<0.1; p<0.05; p<0.01