## Carga de librerías
library(stargazer)
library(sandwich)
library(lmtest)
## Cargar entorno de trabajo
load("C:/Users/mendo/Downloads/smoke.RData")
## Verificar objetos disponibles
ls()
## [1] "data" "desc" "self"
## Revisar descripción de variables
desc
## variable label
## 1 educ years of schooling
## 2 cigpric state cig. price, cents/pack
## 3 white =1 if white
## 4 age in years
## 5 income annual income, $
## 6 cigs cigs. smoked per day
## 7 restaurn =1 if rest. smk. restrictions
## 8 lincome log(income)
## 9 agesq age^2
## 10 lcigpric log(cigprice)
## Estimar el modelo lineal
modelo_smoke <- lm(cigs ~ cigpric + lcigpric + income + lincome +
age + agesq + educ + white + restaurn,
data = data)
## Resumen del modelo
summary(modelo_smoke)
##
## Call:
## lm(formula = cigs ~ cigpric + lcigpric + income + lincome + age +
## agesq + educ + white + restaurn, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.169 -9.357 -5.915 7.851 70.744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 340.80437460 260.01558727 1.311 0.19033
## cigpric 2.00226767 1.49283119 1.341 0.18022
## lcigpric -115.27346445 85.42431520 -1.349 0.17758
## income -0.00004619 0.00013349 -0.346 0.72940
## lincome 1.40406118 1.70816584 0.822 0.41134
## age 0.77835901 0.16055561 4.848 0.000001500 ***
## agesq -0.00915035 0.00174929 -5.231 0.000000216 ***
## educ -0.49478062 0.16818020 -2.942 0.00336 **
## white -0.53105164 1.46072181 -0.364 0.71629
## restaurn -2.64424135 1.12999869 -2.340 0.01953 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.41 on 797 degrees of freedom
## Multiple R-squared: 0.05515, Adjusted R-squared: 0.04448
## F-statistic: 5.169 on 9 and 797 DF, p-value: 0.0000007735
## Presentar resumen con stargazer
stargazer(modelo_smoke,
type = "text",
title = "Modelo de Consumo de Cigarrillos")
##
## Modelo de Consumo de Cigarrillos
## ===============================================
## 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.01
########################################################
## Verificación de la matriz Varianza-Covarianza
########################################################
## Matriz Varianza-Covarianza clásica
vcov(modelo_smoke)
## (Intercept) cigpric lcigpric
## (Intercept) 67608.1056227387 386.109182633052 -22173.8151710673
## cigpric 386.1091826331 2.228544957786 -127.2304071208
## lcigpric -22173.8151710673 -127.230407120808 7297.3136265717
## income 0.0004361954 -0.000006112577 0.0003973538
## lincome -8.1784660668 0.071245776994 -4.8604250078
## age 0.4955738865 0.003684842214 -0.2196984501
## agesq -0.0123717935 -0.000080037057 0.0046196387
## educ 0.0942194523 0.001549304349 -0.1048653273
## white -2.2541916142 0.011772809943 -0.3119455870
## restaurn 40.6483588334 0.218648856595 -13.4058899507
## income lincome age
## (Intercept) 0.000436195430145 -8.1784660668 0.4955738865339
## cigpric -0.000006112576584 0.0712457770 0.0036848422138
## lcigpric 0.000397353804854 -4.8604250078 -0.2196984500534
## income 0.000000017819865 -0.0002062280 -0.0000002490054
## lincome -0.000206228025158 2.9178305410 -0.0286441757735
## age -0.000000249005395 -0.0286441758 0.0257781044805
## agesq 0.000000004523652 0.0002930256 -0.0002764370727
## educ -0.000002445251202 -0.0026176694 -0.0029490359616
## white -0.000007733767820 0.1323738372 -0.0143853120956
## restaurn -0.000011505007929 0.0828693504 0.0038113359919
## agesq educ white restaurn
## (Intercept) -0.012371793547420 0.094219452281 -2.254191614183 40.64835883343
## cigpric -0.000080037057102 0.001549304349 0.011772809943 0.21864885659
## lcigpric 0.004619638679171 -0.104865327267 -0.311945587001 -13.40588995075
## income 0.000000004523652 -0.000002445251 -0.000007733768 -0.00001150501
## lincome 0.000293025623241 -0.002617669405 0.132373837244 0.08286935041
## age -0.000276437072676 -0.002949035962 -0.014385312096 0.00381133599
## agesq 0.000003060023723 0.000041446583 0.000172193969 -0.00003412736
## educ 0.000041446582558 0.028284579054 0.002035626571 -0.00382119824
## white 0.000172193969192 0.002035626571 2.133708195634 0.16028942895
## restaurn -0.000034127361210 -0.003821198244 0.160289428946 1.27689703950
## Prueba de heterocedasticidad de Breusch-Pagan
bptest(modelo_smoke)
##
## studentized Breusch-Pagan test
##
## data: modelo_smoke
## BP = 33.525, df = 9, p-value = 0.0001082
## Prueba de autocorrelación de Breusch-Godfrey
bgtest(modelo_smoke)
##
## Breusch-Godfrey test for serial correlation of order up to 1
##
## data: modelo_smoke
## LM test = 0.069737, df = 1, p-value = 0.7917
## Interpretación:
## Si se rechaza H0 en cualquiera de las pruebas,
## la matriz Varianza-Covarianza NO es escalar.
##
## Fuente del problema:
## - Heterocedasticidad -> varianza no constante
## - Autocorrelación -> correlación serial en errores
########################################################
## Estimador HAC
########################################################
## Obtener matriz HAC (Newey-West)
vcov_HAC <- NeweyWest(modelo_smoke)
## Ver matriz HAC
vcov_HAC
## (Intercept) cigpric lcigpric
## (Intercept) 63053.186344524 363.813491868048 -20603.7464784922
## cigpric 363.813491851 2.118848167023 -119.3298061539
## lcigpric -20603.746478883 -119.329806161888 6748.3447957320
## income 0.002515851 0.000008743147 -0.0004662915
## lincome -65.720800311 -0.311098481079 17.2850173544
## age 3.184441336 0.021392544662 -1.1586666560
## agesq -0.045243548 -0.000295620468 0.0161362007
## educ -0.698134869 -0.000977406153 0.1253344765
## white -66.242349492 -0.370628671344 21.4171780907
## restaurn 69.498842477 0.400208669334 -23.2368765344
## income lincome age agesq
## (Intercept) 0.00251585077770 -65.72080033044 3.184441336054 -0.04524354761560
## cigpric 0.00000874314663 -0.31109848118 0.021392544665 -0.00029562046819
## lcigpric -0.00046629152991 17.28501736121 -1.158666656136 0.01613620071441
## income 0.00000001442512 -0.00014019299 -0.000001104724 0.00000001560916
## lincome -0.00014019298530 1.72813522629 -0.007812650741 0.00003301937404
## age -0.00000110472363 -0.00781265074 0.018840993978 -0.00019842600165
## agesq 0.00000001560916 0.00003301937 -0.000198426002 0.00000214548533
## educ -0.00000311098538 -0.00670648461 0.000689809477 0.00000040889072
## white 0.00001226785287 -0.06860415511 -0.000494530317 0.00002014964883
## restaurn -0.00002430139734 0.17085614484 0.007364021897 -0.00008832958382
## educ white restaurn
## (Intercept) -0.6981348700802 -66.24234945776 69.49884247084
## cigpric -0.0009774061614 -0.37062867115 0.40020866929
## lcigpric 0.1253344769420 21.41717807945 -23.23687653272
## income -0.0000031109854 0.00001226785 -0.00002430140
## lincome -0.0067064846149 -0.06860415511 0.17085614485
## age 0.0006898094774 -0.00049453032 0.00736402190
## agesq 0.0000004088907 0.00002014965 -0.00008832958
## educ 0.0285014160639 -0.01281858241 -0.01369174956
## white -0.0128185824073 1.65254142882 0.05461606050
## restaurn -0.0136917495583 0.05461606052 1.08598283251
## Errores estándar HAC
se_HAC <- sqrt(diag(vcov_HAC))
########################################################
## Tabla comparativa con stargazer
########################################################
stargazer(modelo_smoke, modelo_smoke,
se = list(NULL, se_HAC),
column.labels = c("Original", "Corregido HAC"),
type = "text",
title = "Modelo Original vs Modelo Corregido HAC",
align = TRUE,
no.space = TRUE)
##
## Modelo Original vs Modelo Corregido HAC
## ===========================================================
## Dependent variable:
## ----------------------------
## cigs
## Original Corregido HAC
## (1) (2)
## -----------------------------------------------------------
## cigpric 2.002 2.002
## (1.493) (1.456)
## lcigpric -115.273 -115.273
## (85.424) (82.148)
## income -0.00005 -0.00005
## (0.0001) (0.0001)
## lincome 1.404 1.404
## (1.708) (1.315)
## age 0.778*** 0.778***
## (0.161) (0.137)
## 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.286)
## restaurn -2.644** -2.644**
## (1.130) (1.042)
## Constant 340.804 340.804
## (260.016) (251.104)
## -----------------------------------------------------------
## Observations 807 807
## R2 0.055 0.055
## Adjusted R2 0.044 0.044
## Residual Std. Error (df = 797) 13.413 13.413
## F Statistic (df = 9; 797) 5.169*** 5.169***
## ===========================================================
## Note: *p<0.1; **p<0.05; ***p<0.01