load("C:/Users/MINEDUCYT/Downloads/smoke.RData")
summary(data)## educ cigpric white age
## Min. : 6.00 Min. :44.00 Min. :0.0000 Min. :17.00
## 1st Qu.:10.00 1st Qu.:58.14 1st Qu.:1.0000 1st Qu.:28.00
## Median :12.00 Median :61.05 Median :1.0000 Median :38.00
## Mean :12.47 Mean :60.30 Mean :0.8786 Mean :41.24
## 3rd Qu.:13.50 3rd Qu.:63.18 3rd Qu.:1.0000 3rd Qu.:54.00
## Max. :18.00 Max. :70.13 Max. :1.0000 Max. :88.00
## income cigs restaurn lincome
## Min. : 500 Min. : 0.000 Min. :0.0000 Min. : 6.215
## 1st Qu.:12500 1st Qu.: 0.000 1st Qu.:0.0000 1st Qu.: 9.433
## Median :20000 Median : 0.000 Median :0.0000 Median : 9.903
## Mean :19305 Mean : 8.686 Mean :0.2466 Mean : 9.687
## 3rd Qu.:30000 3rd Qu.:20.000 3rd Qu.:0.0000 3rd Qu.:10.309
## Max. :30000 Max. :80.000 Max. :1.0000 Max. :10.309
## agesq lcigpric
## Min. : 289 Min. :3.784
## 1st Qu.: 784 1st Qu.:4.063
## Median :1444 Median :4.112
## Mean :1990 Mean :4.096
## 3rd Qu.:2916 3rd Qu.:4.146
## Max. :7744 Max. :4.250##Estimar modelo
library(stargazer)
library(htmltools)
modelo_cigs <- lm(cigs ~ cigpric + lcigpric + income + lincome + age + agesq + educ + white + restaurn, data = data)
stargazer(modelo_cigs,title = "modelo estimado",type = "text")##
## modelo estimado
## ===============================================
## 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.01Verifique si la matriz de Varianza-Covarianza del modelo es escalar, de no ser escalar indique la fuente del problema.
library(skedastic)## Warning: package 'skedastic' was built under R version 4.5.3skedastic::white(modelo_cigs, interactions = FALSE)## # A tibble: 1 × 5
## statistic p.value parameter method alternative
## <dbl> <dbl> <dbl> <chr> <chr>
## 1 36.8 0.00556 18 White's Test greaterconclusion. Como p-value(0.00556) < nivel de signficancia (0.05) se rechaza la hipótesis nula, por lo tanto hay evidencia de que la varianza de los residuos es heterocedasticidad .
library(lmtest)## Cargando paquete requerido: zoo##
## Adjuntando el paquete: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericbgtest(modelo_cigs,order = 2)##
## Breusch-Godfrey test for serial correlation of order up to 2
##
## data: modelo_cigs
## LM test = 0.26889, df = 2, p-value = 0.8742Conclusion. Como p value (0.8742) > nivel de significancia(0.05) No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden 2.
library(lmtest)
bgtest(modelo_cigs,order = 1)##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_cigs
## LM test = 0.069737, df = 1, p-value = 0.7917Conclusion. Como p value (0.7917) > nivel de significancia(0.05) No se rechaza H0, por lo tanto puede concluirse que los residuos del modelo, no siguen autocorrelación de orden 1.
options(scipen = 999) # Evitar notación científica para números grandes
## Obtenga el Estimador HAC apropiado y denominelo vcov_HAC
library(sandwich)## Warning: package 'sandwich' was built under R version 4.5.3library(lmtest)
# 1. Obtención de la matriz de varianza-covarianza HAC
vcov_HAC <- vcovHAC(modelo_cigs)
# 2. Ver la matriz calculada en la consola (opcional)
print(vcov_HAC)## (Intercept) cigpric lcigpric
## (Intercept) 65509.511681212 377.865370642112 -21441.4139125673
## cigpric 377.865370602 2.200149654624 -124.1301211483
## lcigpric -21441.413912023 -124.130121158464 7033.3007081212
## income 0.001988193 0.000005937686 -0.0003050133
## lincome -52.407362232 -0.238225247439 13.1133847030
## age 3.712707114 0.024528749228 -1.3269068765
## agesq -0.050908032 -0.000328365720 0.0179211770
## educ -1.068797713 -0.003399136596 0.2576712546
## white -58.139871067 -0.315878529387 18.5772182742
## restaurn 62.897848883 0.360592090221 -21.0374640659
## income lincome age
## (Intercept) 0.00198819305689 -52.40736223225 3.7127071139277
## cigpric 0.00000593768576 -0.23822524744 0.0245287492279
## lcigpric -0.00030501329732 13.11338470279 -1.3269068765233
## income 0.00000001423620 -0.00013739187 -0.0000008489929
## lincome -0.00013739187421 1.67123889246 -0.0117895204474
## age -0.00000084899293 -0.01178952045 0.0193352098220
## agesq 0.00000001288789 0.00007583843 -0.0002029808634
## educ -0.00000255933602 -0.00970931814 0.0001927667421
## white 0.00001088341401 -0.05795566258 0.0002904116891
## restaurn -0.00002380434876 0.16932400471 0.0060377063173
## agesq educ white restaurn
## (Intercept) -0.05090803162829 -1.068797713156 -58.13987106643 62.89784888385
## cigpric -0.00032836572001 -0.003399136596 -0.31587852938 0.36059209020
## lcigpric 0.01792117702885 0.257671254615 18.57721827400 -21.03746406585
## income 0.00000001288789 -0.000002559336 0.00001088341 -0.00002380435
## lincome 0.00007583842972 -0.009709318138 -0.05795566258 0.16932400471
## age -0.00020298086337 0.000192766742 0.00029041169 0.00603770632
## agesq 0.00000218856115 0.000006056762 0.00001138464 -0.00007430878
## educ 0.00000605676164 0.028439901792 -0.01016281211 -0.01528975880
## white 0.00001138463777 -0.010162812112 1.74983747635 0.09385922100
## restaurn -0.00007430878193 -0.015289758800 0.09385922100 1.08706587907# 3. Aplicación del estimador robusto sobre los coeficientes del modelo
# Esto te mostrará los nuevos errores estándar, estadísticos t y p-valores robustos
coeftest(modelo_cigs, vcov. = vcov_HAC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 340.804374604 255.948259774 1.3315 0.183393
## cigpric 2.002267667 1.483290145 1.3499 0.177437
## lcigpric -115.273464445 83.864776325 -1.3745 0.169668
## income -0.000046194 0.000119316 -0.3872 0.698745
## lincome 1.404061178 1.292764051 1.0861 0.277767
## age 0.778359013 0.139051105 5.5976 0.0000000298971 ***
## agesq -0.009150353 0.001479379 -6.1853 0.0000000009901 ***
## educ -0.494780616 0.168641341 -2.9339 0.003443 **
## white -0.531051635 1.322814226 -0.4015 0.688192
## restaurn -2.644241351 1.042624515 -2.5361 0.011398 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(stargazer)
# 1. Extraemos la raíz cuadrada de la diagonal de vcov_HAC para obtener los errores estándar robustos
errores_se_HAC <- sqrt(diag(vcov_HAC))
# 2. Generamos la tabla conjunta siguiendo la documentación de stargazer
stargazer(modelo_cigs, modelo_cigs,
se = list(NULL, errores_se_HAC),
type = "text",
title = "Resultados del Consumo de Cigarrillos: Inferencia Estándar vs Robustecida (HAC)",
column.labels = c("Original (MCO)", "Corregido (HAC)"),
model.numbers = FALSE, # Opcional: quita los números (1) y (2) sobre las columnas
dep.var.labels = "Cigarrillos consumidos (cigs)",
keep.stat = c("n", "rsq", "adj.rsq", "f")) # Opcional: muestra solo estadísticas clave al final##
## Resultados del Consumo de Cigarrillos: Inferencia Estándar vs Robustecida (HAC)
## ========================================================
## Dependent variable:
## ------------------------------
## Cigarrillos consumidos (cigs)
## Original (MCO) Corregido (HAC)
## --------------------------------------------------------
## 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
## F Statistic (df = 9; 797) 5.169*** 5.169***
## ========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01