Aplicación de la librería stargazer para presentar modelos corregidos con estimadores HAC

1. Estimación del Modelo

options(scipen = 999999)
library (stargazer)
library(foreign)
datos_regresion <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
modelo_crime<-lm(crime~poverty+single,data=datos_regresion)
stargazer(modelo_crime,type = "html",title = "Modelo Crimen")
Modelo Crimen
Dependent variable:
crime
poverty 6.787
(8.989)
single 166.373***
(19.423)
Constant -1,368.189***
(187.205)
Observations 51
R2 0.707
Adjusted R2 0.695
Residual Std. Error 243.610 (df = 48)
F Statistic 57.964*** (df = 2; 48)
Note: p<0.1; p<0.05; p<0.01

2. Pruebas de Heterocedasticidad y Autocorrelación

Prueba de Heterocedasticidad

Usando la libreria “lmtest”

library(lmtest)
white_test<-bptest(modelo_crime,~I(poverty^2)+I(single^2)+poverty*single,data = datos_regresion)
print(white_test)
studentized Breusch-Pagan test

data: modelo_crime BP = 10.73, df = 5, p-value = 0.057

Hay evidencia de heterocedasticidad ya que pvalue < 0.05

Prueba del Multiplicador de Lagrange (Breusch Godfrey)

Autocorrelación de 2º Orden:

library(lmtest)
prueba_LM<-bgtest(modelo_crime,order = 2)
print(prueba_LM)
Breusch-Godfrey test for serial correlation of order up to 2

data: modelo_crime LM test = 0.27165, df = 2, p-value = 0.873

No hay evidencia de Autocorrelación de 2º orden ya que pvalue > 0.05

Autocorrelación de 1º orden (prueba de Durbin Watson)

library(car)
durbinWatsonTest(model = modelo_crime)

lag Autocorrelation D-W Statistic p-value 1 -0.07014421 2.040007 0.918 Alternative hypothesis: rho != 0

No hay evidencia de Autocorrelación de 1º orden ya que pvalue > 0.05

3. Estimación Robusta (uso del estimador HAC)

Sin Corregir

options(scipen = 99999)
library(lmtest)
#Sin corregir:
coeftest(modelo_crime)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value         Pr(>|t|)    
## (Intercept) -1368.1887   187.2052 -7.3085 0.00000000247861 ***
## poverty         6.7874     8.9885  0.7551           0.4539    
## single        166.3727    19.4229  8.5658 0.00000000003117 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Corregido usando HAC

library(lmtest)
library(sandwich)

#Corregido:
# Para Autocorrelacion de Orden 1
estimacion_omega1<-vcovHC(modelo_crime,type = "HC1")
coeftest(modelo_crime,vcov. = estimacion_omega1) 
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value      Pr(>|t|)    
## (Intercept) -1368.1887   284.9180 -4.8020 0.00001576624 ***
## poverty         6.7874    10.9273  0.6211        0.5374    
## single        166.3727    26.2343  6.3418 0.00000007519 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
robust.se<-sqrt(diag(estimacion_omega1))


#Para Autocorrelacion de Orden 2
estimacion_omega<-NeweyWest(modelo_crime,lag = 2)
coeftest(modelo_crime,vcov. = estimacion_omega)
## 
## t test of coefficients:
## 
##               Estimate Std. Error t value      Pr(>|t|)    
## (Intercept) -1368.1887   303.8466 -4.5029 0.00004279768 ***
## poverty         6.7874    10.5943  0.6407        0.5248    
## single        166.3727    25.9154  6.4198 0.00000005708 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Estimación Robusta usando libreria Stargazer

Nota: Modelo corregido para autocorrelación de primer orden

options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_crime_corregido<-lmrob(crime~poverty+single,data=datos_regresion,)

#Para ver el modelo original y el modelo corregido en dos columnas
stargazer(modelo_crime,modelo_crime,se=list(NULL, robust.se),type = "html",title = "Tabla Comparativa", column.labels = c("Modelo Original", "Modelo Corregido") )
Tabla Comparativa
Dependent variable:
crime
Modelo Original Modelo Corregido
(1) (2)
poverty 6.787 6.787
(8.989) (10.927)
single 166.373*** 166.373***
(19.423) (26.234)
Constant -1,368.189*** -1,368.189***
(187.205) (284.918)
Observations 51 51
R2 0.707 0.707
Adjusted R2 0.695 0.695
Residual Std. Error (df = 48) 243.610 243.610
F Statistic (df = 2; 48) 57.964*** 57.964***
Note: p<0.1; p<0.05; p<0.01