Indicación: Aplicar la librería stargazer para presentar modelos corregidos con estimadores HAC (Heterocedastisticity And Autocorrelation).
Retomando el dataset visto en clase:
options(scipen = 999999)
library(foreign)
datos_regresion <- read.dta("https://stats.idre.ucla.edu/stat/data/crime.dta")
modelo_estimado<-lm(crime~poverty+single,data=datos_regresion)
print(modelo_estimado)##
## Call:
## lm(formula = crime ~ poverty + single, data = datos_regresion)
##
## Coefficients:
## (Intercept) poverty single
## -1368.189 6.787 166.373Una vez cargado el modelo deben realizarse las pruebas para realizar la estimación robusta adecuada.
Nota: Incluir los términos cruzados. Se trata de un modelo de 2 variables explicativas.
Hipótesis
\(H_o\): La varianza de los residuos es homocedástica. Todos los paramétros del modelo son simultáneamente cero.
\(H_1\): La varianza de los residuos es heterocedástica. No todos los coeficientes del modelo son cero.
Regla de rechazo
\(LM_W ≥ VC\) rechazar \(H_o\)
\(p-value ≤ α\) rechazar \(H_o\)
library(lmtest)
white_test<-bptest(modelo_estimado,~I(poverty^2)+I(single^2)+poverty*single,
data = datos_regresion)
print(white_test)##
## studentized Breusch-Pagan test
##
## data: modelo_estimado
## BP = 10.73, df = 5, p-value = 0.057Se rechaza la hipótesis nula dado que p-value = 0.057 ≤ 0.05. Hay evidencia de heterocedisticidad en el modelo. No todos los coeficientes son simultáneamente cero.
Hipótesis
\(H_o: ρ = 0\) No hay evidencia de autocorrelación de primer orden en los residuos del modelo.
\(H_1: ρ ≠ 0\) Hay evidencia de autocorrelación de primer orden en los residuos del modelo.
Regla de rechazo
\(p-value ≤ α\) rechazar \(H_o\)
library(car)
durbinWatsonTest(model = modelo_estimado)## lag Autocorrelation D-W Statistic p-value
## 1 -0.07014421 2.040007 0.988
## Alternative hypothesis: rho != 0No hay evidencia de autocorrelación de primer orden dado que el p-value es mayor al nivel de significancia.
Hipótesis
\(H_o: ρ_m = 0\) No hay evidencia de autocorrelación de orden “m” en los residuos del modelo.
\(H_1: ρ_m ≠ 0\) Hay evidencia de autocorrelación de orden “m” en los residuos del modelo.
Regla de rechazo
\(LM_B > VC\) rechazar \(H_o\)
\(p-value ≤ α\) rechazar \(H_o\)
Para la detección de autocorrelación de segundo orden se tiene que:
library(lmtest)
prueba_LM <- bgtest(modelo_estimado,order = 2)
print(prueba_LM)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_estimado
## LM test = 0.27165, df = 2, p-value = 0.873El p-value = 0.873 > 0.05, no hay evidencia de autocorrelación de orden 2.
Dado que las pruebas únicamente detectaron heterocedisticad se hará uso de la librería lmtest y sandwich con el argumento HCO.
Modelo sin corregir.
options(scipen = 99999)
library(lmtest)
# Mostrando el modelo original
coeftest(modelo_estimado)##
## 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 ' ' 1Corrección del modelo, HCO
options(scipen = 99999)
library(lmtest)
library(sandwich)
# Matriz VarCov corregida
estimacion_omega <- vcovHC(modelo_estimado, type = "HC0")
coeftest_corregido <- coeftest(modelo_estimado, vcov. = estimacion_omega)
print(coeftest_corregido)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1368.1887 276.4111 -4.9498 0.00000956181 ***
## poverty 6.7874 10.6010 0.6403 0.5251
## single 166.3727 25.4510 6.5370 0.00000003774 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Comparación de modelos con Stargazer
options(scipen = 99999)
library(stargazer)
# Retomando
estimacion_omega <- vcovHC(modelo_estimado, type = "HC0")
se_fix <- sqrt(diag(estimacion_omega)) # Para obtener los errores estándar
# Comparación de modelos
stargazer(modelo_estimado, modelo_estimado,
se=list(NULL, se_fix), # Agregando la correción HAC
type = "html",
title = "Comparación de modelos: corrección HAC.",
column.labels = c("Original", "Corregido"),
model.names = FALSE,
dep.var.labels = "Modelo",
dep.var.caption = "Dependent variable: crime",
model.numbers = FALSE)| Dependent variable: crime | ||
| Modelo | ||
| Original | Corregido | |
| poverty | 6.787 | 6.787 |
| (8.989) | (10.601) | |
| single | 166.373*** | 166.373*** |
| (19.423) | (25.451) | |
| Constant | -1,368.189*** | -1,368.189*** |
| (187.205) | (276.411) | |
| 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 | |
Estimación robusta del modelo.
A través de la estimación robusta de estimadores MM de Yohai se calculan los coeficientes de regresión y la matriz de covarianzas de los errores de forma simultánea.
options(scipen = 999999)
library(robustbase)
library(stargazer)
modelo_robusto <- lmrob(crime~poverty+single, data=datos_regresion)
stargazer(modelo_estimado,modelo_robusto,
type = "html",
title = "Comparación de modelos: Estimación robusta.",
column.labels = c("Original", "Robusto"),
model.names = FALSE,
model.numbers = FALSE,
dep.var.caption = "Dependent variable: crime",
dep.var.labels = "Modelo")| Dependent variable: crime | ||
| Modelo | ||
| Original | Robusto | |
| poverty | 6.787 | 11.466 |
| (8.989) | (9.263) | |
| single | 166.373*** | 176.569*** |
| (19.423) | (23.223) | |
| Constant | -1,368.189*** | -1,539.640*** |
| (187.205) | (235.765) | |
| Observations | 51 | 51 |
| R2 | 0.707 | 0.795 |
| Adjusted R2 | 0.695 | 0.787 |
| Residual Std. Error (df = 48) | 243.610 | 191.864 |
| F Statistic | 57.964*** (df = 2; 48) | |
| Note: | p<0.1; p<0.05; p<0.01 | |
Extra: Uso de la librería NeweyTest.
Para corregir correlación de segundo orden en un modelo el procedimiento a seguir es el siguiente:
library(lmtest)
library(sandwich)
# Correlación de orden 2, lag = 2
estimacion_omega2<-NeweyWest(modelo_estimado,lag = 2)
coeftest(modelo_estimado,vcov. = estimacion_omega2)##
## 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 ' ' 1Comparación de modelos con Stargazer
options(scipen = 99999)
library(stargazer)
# Retomando
estimacion_omega2<-NeweyWest(modelo_estimado,lag = 2)
se_fix2 <- sqrt(diag(estimacion_omega2)) # Para obtener los errores estándar
# Comparación de modelos
stargazer(modelo_estimado, modelo_estimado,
se=list(NULL, se_fix2), # Agregando la correción HAC
type = "html",
title = "Comparación de modelos: corrección HAC de NeweyTest.",
column.labels = c("Original", "Corregido"),
model.names = FALSE,
dep.var.labels = "Modelo",
dep.var.caption = "Dependent variable: crime",
model.numbers = FALSE)| Dependent variable: crime | ||
| Modelo | ||
| Original | Corregido | |
| poverty | 6.787 | 6.787 |
| (8.989) | (10.594) | |
| single | 166.373*** | 166.373*** |
| (19.423) | (25.915) | |
| Constant | -1,368.189*** | -1,368.189*** |
| (187.205) | (303.847) | |
| 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 | |