options(scipen = 999999)
library(foreign)
datos_regresion <- read.dta("C:/Users/Abigail/Desktop/econometria/ejercicios_de_clase/Estimadores HAC/crime.dta")
modelo_crime<-lm(crime~poverty+single,data=datos_regresion)
print(modelo_crime)
## 
## Call:
## lm(formula = crime ~ poverty + single, data = datos_regresion)
## 
## Coefficients:
## (Intercept)      poverty       single  
##   -1368.189        6.787      166.373

Pruebas de heterocesdasticidad y autocorrelación

Heterocedasticidad

#Prueba de White (prueba de Breusch Pagan)
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

Autocorrelación de 2º Orden:

#Prueba del Multiplicador de Lagrange (Breusch Godfrey)
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 p-value>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.944
##  Alternative hypothesis: rho != 0

No hay evidencia de Autocorrelación de 1º orden ya que p-value>0.05

Estimación Robusta (uso del estimador HAC)

Sin corregir

options(scipen = 99999)
library(lmtest)

coeftest(modelo_crime)   #Sin corregir
## 
## 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

options(scipen = 99999)
library(lmtest)
library(sandwich)
estimacion_omg <- vcovHC(modelo_crime,type="HC1") 
coeftest(modelo_crime, vcov = estimacion_omg)
## 
## 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
options(scipen = 999999)
library(stargazer)
library(lmtest)
library(sandwich)

robust.se <- sqrt(diag(estimacion_omg))

stargazer(modelo_crime, modelo_crime, se=list(NULL, robust.se),
column.labels=c("Sin corregir","Corregido HAC"), align=TRUE,type = "html",title = "comparativa")
comparativa
Dependent variable:
crime
Sin corregir Corregido HAC
(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