# --- 0. Cargar datos y librerías ---
load("smoke.RData")   # carga los objetos: data, desc
 
library(stargazer)
## 
## Please cite as:
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(sandwich)
library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
print(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)
head(data)
##   educ cigpric white age income cigs restaurn   lincome agesq lcigpric
## 1 16.0  60.506     1  46  20000    0        0  9.903487  2116 4.102743
## 2 16.0  57.883     1  40  30000    0        0 10.308952  1600 4.058424
## 3 12.0  57.664     1  58  30000    3        0 10.308952  3364 4.054633
## 4 13.5  57.883     1  30  20000    0        0  9.903487   900 4.058424
## 5 10.0  58.320     1  17  20000    0        0  9.903487   289 4.065945
## 6  6.0  59.340     1  86   6500    0        0  8.779557  7396 4.083283

ESTIMAR MODELO

#    cigs = f(cigpric, lcigpric, income, lincome, age, agesq, educ, white, restaurn)
# ============================================================
 
modelo_original <- lm(cigs ~ cigpric + lcigpric + income + lincome +
                         age + agesq + educ + white + restaurn,
                       data = data)

RESUMEN DEL MODELO CON STARGAZER

stargazer(modelo_original,
          type    = "text",
          title   = "Modelo Original: Determinantes del consumo de cigarrillos",
          digits  = 4)
## 
## Modelo Original: Determinantes del consumo de cigarrillos
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                cigs            
## -----------------------------------------------
## cigpric                       2.0023           
##                              (1.4928)          
##                                                
## lcigpric                     -115.2735         
##                              (85.4243)         
##                                                
## income                       -0.00005          
##                              (0.0001)          
##                                                
## lincome                       1.4041           
##                              (1.7082)          
##                                                
## age                          0.7784***         
##                              (0.1606)          
##                                                
## agesq                       -0.0092***         
##                              (0.0017)          
##                                                
## educ                        -0.4948***         
##                              (0.1682)          
##                                                
## white                         -0.5311          
##                              (1.4607)          
##                                                
## restaurn                     -2.6442**         
##                              (1.1300)          
##                                                
## Constant                     340.8044          
##                             (260.0156)         
##                                                
## -----------------------------------------------
## Observations                    807            
## R2                            0.0552           
## Adjusted R2                   0.0445           
## Residual Std. Error     13.4128 (df = 797)     
## F Statistic           5.1693*** (df = 9; 797)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

VERIFICACION MATRIZ-COVARIANZA

vcov_original <- vcov(modelo_original)
vcov_original
##               (Intercept)       cigpric      lcigpric        income
## (Intercept)  6.760811e+04  3.861092e+02 -2.217382e+04  4.361954e-04
## cigpric      3.861092e+02  2.228545e+00 -1.272304e+02 -6.112577e-06
## lcigpric    -2.217382e+04 -1.272304e+02  7.297314e+03  3.973538e-04
## income       4.361954e-04 -6.112577e-06  3.973538e-04  1.781986e-08
## lincome     -8.178466e+00  7.124578e-02 -4.860425e+00 -2.062280e-04
## age          4.955739e-01  3.684842e-03 -2.196985e-01 -2.490054e-07
## agesq       -1.237179e-02 -8.003706e-05  4.619639e-03  4.523652e-09
## educ         9.421945e-02  1.549304e-03 -1.048653e-01 -2.445251e-06
## white       -2.254192e+00  1.177281e-02 -3.119456e-01 -7.733768e-06
## restaurn     4.064836e+01  2.186489e-01 -1.340589e+01 -1.150501e-05
##                   lincome           age         agesq          educ
## (Intercept) -8.1784660668  4.955739e-01 -1.237179e-02  9.421945e-02
## cigpric      0.0712457770  3.684842e-03 -8.003706e-05  1.549304e-03
## lcigpric    -4.8604250078 -2.196985e-01  4.619639e-03 -1.048653e-01
## income      -0.0002062280 -2.490054e-07  4.523652e-09 -2.445251e-06
## lincome      2.9178305410 -2.864418e-02  2.930256e-04 -2.617669e-03
## age         -0.0286441758  2.577810e-02 -2.764371e-04 -2.949036e-03
## agesq        0.0002930256 -2.764371e-04  3.060024e-06  4.144658e-05
## educ        -0.0026176694 -2.949036e-03  4.144658e-05  2.828458e-02
## white        0.1323738372 -1.438531e-02  1.721940e-04  2.035627e-03
## restaurn     0.0828693504  3.811336e-03 -3.412736e-05 -3.821198e-03
##                     white      restaurn
## (Intercept) -2.254192e+00  4.064836e+01
## cigpric      1.177281e-02  2.186489e-01
## lcigpric    -3.119456e-01 -1.340589e+01
## income      -7.733768e-06 -1.150501e-05
## lincome      1.323738e-01  8.286935e-02
## age         -1.438531e-02  3.811336e-03
## agesq        1.721940e-04 -3.412736e-05
## educ         2.035627e-03 -3.821198e-03
## white        2.133708e+00  1.602894e-01
## restaurn     1.602894e-01  1.276897e+00
is.scalar <- is.vector(vcov_original) || length(vcov_original) == 1
is.scalar
## [1] FALSE

ESTIMADOR HAC

vcov_HAC <- vcovHC(modelo_original, type = "HC1")
 
# Errores estandar HAC
se_HAC <- sqrt(diag(vcov_HAC))
 
cat("\n=== Errores Estandar HAC (HC1) ===\n")
## 
## === Errores Estandar HAC (HC1) ===
print(se_HAC)
##  (Intercept)      cigpric     lcigpric       income      lincome          age 
## 2.803072e+02 1.612751e+00 9.191569e+01 1.163156e-04 1.236657e+00 1.378011e-01 
##        agesq         educ        white     restaurn 
## 1.460626e-03 1.639876e-01 1.370425e+00 1.044748e+00

TABLA COMPARATIVA: MODELO ORIGINAL vs MODELO CORREGIDO

# Errores estandar originales (OLS)
se_OLS <- sqrt(diag(vcov(modelo_original)))
 
stargazer(modelo_original, modelo_original,
          type           = "text",
          title          = "Comparacion: Modelo Original vs Corregido con HAC",
          column.labels  = c("Original (OLS)", "Corregido (HAC)"),
          se             = list(se_OLS, se_HAC),
          dep.var.labels = "Cigarrillos consumidos al dia (cigs)",
          covariate.labels = c("Precio cigarrillos (cigpric)",
                               "Log precio cigarrillos (lcigpric)",
                               "Ingreso (income)",
                               "Log ingreso (lincome)",
                               "Edad (age)",
                               "Edad al cuadrado (agesq)",
                               "Educacion (educ)",
                               "Raza blanca (white)",
                               "Restaurante libre de humo (restaurn)"),
          omit.stat      = c("LL", "ser"),
          digits         = 4,
          notes          = c("Errores estandar en parentesis.",
                             "Columna (2): errores robustos HC1 (Huber-White).",
                             "*p<0.1; **p<0.05; ***p<0.01"))
## 
## Comparacion: Modelo Original vs Corregido con HAC
## ======================================================================================
##                                                     Dependent variable:               
##                                      -------------------------------------------------
##                                            Cigarrillos consumidos al dia (cigs)       
##                                           Original (OLS)          Corregido (HAC)     
##                                                (1)                      (2)           
## --------------------------------------------------------------------------------------
## Precio cigarrillos (cigpric)                  2.0023                   2.0023         
##                                              (1.4928)                 (1.6128)        
##                                                                                       
## Log precio cigarrillos (lcigpric)           -115.2735                -115.2735        
##                                             (85.4243)                (91.9157)        
##                                                                                       
## Ingreso (income)                             -0.00005                 -0.00005        
##                                              (0.0001)                 (0.0001)        
##                                                                                       
## Log ingreso (lincome)                         1.4041                   1.4041         
##                                              (1.7082)                 (1.2367)        
##                                                                                       
## Edad (age)                                  0.7784***                0.7784***        
##                                              (0.1606)                 (0.1378)        
##                                                                                       
## Edad al cuadrado (agesq)                    -0.0092***               -0.0092***       
##                                              (0.0017)                 (0.0015)        
##                                                                                       
## Educacion (educ)                            -0.4948***               -0.4948***       
##                                              (0.1682)                 (0.1640)        
##                                                                                       
## Raza blanca (white)                          -0.5311                  -0.5311         
##                                              (1.4607)                 (1.3704)        
##                                                                                       
## Restaurante libre de humo (restaurn)        -2.6442**                -2.6442**        
##                                              (1.1300)                 (1.0447)        
##                                                                                       
## Constant                                     340.8044                 340.8044        
##                                             (260.0156)               (280.3072)       
##                                                                                       
## --------------------------------------------------------------------------------------
## Observations                                   807                      807           
## R2                                            0.0552                   0.0552         
## Adjusted R2                                   0.0445                   0.0445         
## F Statistic (df = 9; 797)                   5.1693***                5.1693***        
## ======================================================================================
## Note:                                                      *p<0.1; **p<0.05; ***p<0.01
##                                                        Errores estandar en parentesis.
##                                       Columna (2): errores robustos HC1 (Huber-White).
##                                                            *p<0.1; **p<0.05; ***p<0.01