Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
robertopucp
GitHub Repository: robertopucp/1eco35_2022_2
 Path: blob/main/Trabajo_final/grupo6/grupo6_pregunta1.R
2714 views
1
################  Trabajo Final ############################

2
## Curso: Laboratorio de R y Python #################################

3
## Grupo 6

4


5
# clear environment

6


7
rm(list=ls(all=TRUE))

8


9
# Load libraries ----

10


11


12
librarian::shelf(

13
  tidyverse  # dplyr, tidyr, stringr, ggplot2, etc in unique library

14
  , haven   # to read datset: .dta (stata), .spss, .dbf

15
  , fastDummies  # for dummies

16
  , stargazer  # for summary and econometrics tables

17
  , sandwich  # for linear models

18
  , lmtest # for robust standar error

19
  , estimatr # for iv, cluster, robust standar error (LM_robust)

20
  , lfe  # for fixed effects, cluster standar error

21
  , caret # for easy machine learning workflow (mse, rmse)

22
  , texreg # library for export table

23
  , mfx # probit , logit model marginal effects

24
)

25


26
user <- Sys.getenv("USERNAME")  # username

27


28
setwd( paste0("C:/Users/",user,"/Documents/GitHub/1ECO35_2022_2/Lab10") ) # set directorio

29


30
repdata <- read_dta("../data/dataverse_files/mss_repdata.dta")

31


32


33
str(repdata)

34


35
lapply(repdata, class)

36


37
names(repdata)

38


39
attributes(repdata)

40


41
summary(repdata)

42


43
################# T A B L A 1: ESTADISTICOS DESCRIPTIVOS  ########################

44


45
#construyendo la tabla de las variables solicitadas

46


47
repdata <- read_dta("../data/dataverse_files/mss_repdata.dta")

48


49
table1 <- repdata %>% dplyr::select(NDVI_g, tot_100, trade_pGDP, pop_den_rur, land_crop, va_agr, va_ind_manf)

50


51


52
stargazer(table1)

53


54
# No obtenemos resultado pues la libreria exige que la base de datos sea DataFrame

55


56
table1 <- repdata %>% dplyr::select(NDVI_g, tot_100, trade_pGDP, pop_den_rur, land_crop, va_agr, va_ind_manf) %>% as.data.frame()

57


58


59
stargazer(table1)

60


61
#etiquetas

62


63
stargazer(table1, title = "Estadísticos Descriptivos", digits = 2,

64
          covariate.labels = c("Tasa de variación del índice de vegetación","Términos de intercambio","Porcentajes de las exportaciones respecto al PBI","Densidad poblacional rural",

65
                               "Porcentaje de tierra cultivable en uso","Valor agregado del sector agrícola respecto al PBI","Valor agregado del secctor manufacturero respecto al PBI"))

66
stargazer(table1, title = "Estadísticos Descriptivos", digits = 2,

67
          covariate.labels = c("Tasa de variación del índice de vegetación","Términos de intercambio","Porcentajes de las exportaciones respecto al PBI","Densidad poblacional rural",

68
                               "Porcentaje de tierra cultivable en uso","Valor agregado del sector agrícola respecto al PBI","Valor agregado del secctor manufacturero respecto al PBI"),

69
                              min.max = F)

70


71
#sum, ordenar las columnas de los estadisticos

72


73
list_vars <- c("Tasa de variación del índice de vegetación (NDVI_g)","Términos de intercambio (tot_100)","Porcentajes de las exportaciones respecto al PBI (trade_pGDP)","Densidad poblacional rural (pop

74
               _den_rur)",

75
               "Porcentaje de tierra cultivable en uso (land_crop)","Valor agregado del sector agrícola respecto al PBI (va_agr)","Valor agregado del secctor manufacturero respecto al PBI (va_ind_manf)")

76


77


78
stargazer(table1, title = "Estadísticos descriptivos", digits = 2, # decimales con 2 digitos

79
          covariate.labels = list_vars,  # Lista de etiquetas

80
          summary.stat = c("mean", "sd", "n"), # se especifica el orden de los estadÃ?sticos

81
          min.max = F, # borrar el estadÃ?stico de maximo y minimo

82
          notes = "Note.— Se tomó las variables indicadas."

83
          , notes.append = FALSE, # TRUE append the significance levels

84
          notes.align = 'l')

85


86
################# T A B L A 2: REGRESIÓN MCO SIMPLE  ########################

87


88
#modelo1

89


90
index_country_time <- grep("_time$", colnames(repdata))

91


92
country_time_trend <-names(repdata)[index_country_time]

93


94
index_country <- grep("^ccode_\\d+$", colnames(repdata))

95


96
country_fe <-names(repdata)[index_country]

97


98
ols_model <- lm( any_prio~ GPCP_g + GPCP_g_l,clusters = ccode, data = repdata, se_type = "stata", fixed_effects = ~ ccode)

99


100
attributes(ols_model)

101


102
ols_model$coefficients  # coeficientes

103


104
ols_model$fitted.values  # Y estimado

105


106
# summary table como en stata

107


108
summary(ols_model)

109


110
summary(ols_model)$call

111


112
summary(ols_model)$coef

113


114
# Test de significancia individual usando Huber robust standard errors

115


116
coeftest(ols_model, vcov = vcovHC(ols_model, "HC")) # Clasical white robust

117
coeftest(ols_model, vcov = vcovHC(ols_model, "HC0"))

118


119
coeftest(ols_model, vcov = vcovHC(ols_model, "HC1")) # Huber-White robust (STATA)

120
coeftest(ols_model, vcov = vcovHC(ols_model, "HC2")) # Eicker-Huber-White robust

121
coeftest(ols_model, vcov = vcovHC(ols_model, "HC3"))

122
coeftest(ols_model, vcov = vcovHC(ols_model, "HC4"))

123


124


125
# vcovHC Matrix varianza-covarianza robusta ante heterocedasticidad

126


127
# robust se and cluster

128


129
robust_model <- coeftest(ols_model,

130
                         vcov = vcovCL,  # Matrix varianza-covarianza cluster (CL)

131
                         type = "HC1",

132
                         cluster = ~ ccode)

133


134
sd_robust <- robust_model[,2]

135


136
#rms

137


138
rmse1 <- RMSE(ols_model$fitted.values, repdata$any_prio ) # root mean squeare error

139


140
RMSE(ols_model$fitted.values, repdata$any_prio )^2 # mean square error

141


142
# using tdyr functions

143


144


145
glance(ols_model)

146


147
tidy(ols_model)

148


149
tidy(ols_model)

150


151
htmlreg(ols_model)

152


153
tidy(ols_model)[2,2] # coeficiente

154
tidy(ols_model)[2,6] # limite inferior

155
tidy(ols_model)[2,7] # limite superior 

156


157


158
"Usando LM y luego coestest para los error standar robusto"

159


160
m1 <- lm(any_prio ~ GPCP_g + GPCP_g_l, data = repdata)

161


162
lm_rmse1 <- round(RMSE(m1$fitted.values, repdata$any_prio ) ,2 )

163


164
robust_model <- coeftest(m1,

165
                          vcov = vcovCL,  # Matrix varianza-covarianza cluster (CL)

166
                          type = "HC1",

167
                          cluster = ~ ccode)

168


169
sd_robust_model <- robust_model[,2]

170


171
# intervalo de confianza 

172


173
model1_lower = coefci(m1, df = Inf, 

174
                      vcov. = vcovCL, cluster = ~ ccode, type = "HC1")[2,1]

175


176
model1_upper = coefci(m1, df = Inf,

177
                      vcov. = vcovCL, cluster = ~ ccode, type = "HC1")[2,2]

178


179
tidy(m1, conf.int = TRUE)

180


181


182
# segundo modelo 

183


184
index_country_time <- grep("_time$", colnames(repdata))

185


186
country_time_trend <-names(repdata)[index_country_time]

187


188
index_country <- grep("^ccode_\\d+$", colnames(repdata))

189


190
country_fe <-names(repdata)[index_country]

191


192
ols_model2 <- lm( war_prio~ GPCP_g + GPCP_g_l,clusters = ccode, data = repdata, se_type = "stata", fixed_effects = ~ ccode)

193


194
attributes(ols_model2)

195


196
ols_model2$coefficients  # coeficientes

197


198
ols_model2$fitted.values  # Y estimado

199


200
# summary table como en stata

201


202
summary(ols_model2)

203


204
summary(ols_model2)$call

205


206
summary(ols_model2)$coef

207


208
# Test de significancia individual usando Huber robust standard errors

209


210
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC")) # Clasical white robust

211
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC0"))

212


213
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC1")) # Huber-White robust (STATA)

214
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC2")) # Eicker-Huber-White robust

215
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC3"))

216
coeftest(ols_model, vcov = vcovHC(ols_model2, "HC4"))

217


218


219
# vcovHC Matrix varianza-covarianza robusta ante heterocedasticidad

220


221
# robust se and cluster

222


223
robust_model2 <- coeftest(ols_model2,

224
                         vcov = vcovCL,  # Matrix varianza-covarianza cluster (CL)

225
                         type = "HC1",

226
                         cluster = ~ ccode)

227


228
sd_robust2 <- robust_model2[,2]

229


230
#rms

231


232
rmse2 <- RMSE(ols_model2$fitted.values, repdata$war_prio ) # root mean squeare error

233


234
RMSE(ols_model2$fitted.values, repdata$war_prio )^2 # mean square error

235


236
# using tdyr functions

237


238


239
glance(ols_model2)

240


241
tidy(ols_model2)

242


243
tidy(ols_model2)

244


245
htmlreg(ols_model2)

246


247
tidy(ols_model2)[2,2] # coeficiente

248
tidy(ols_model2)[2,6] # limite inferior

249
tidy(ols_model2)[2,7] # limite superior 

250


251


252
"Usando LM y luego coestest para los error standar robusto"

253


254
m2 <- lm(war_prio ~ GPCP_g + GPCP_g_l, data = repdata)

255


256
lm_rmse2 <- round(RMSE(m2$fitted.values, repdata$war_prio ) ,2 )

257


258
robust_model2 <- coeftest(m2,

259
                         vcov = vcovCL,  # Matrix varianza-covarianza cluster (CL)

260
                         type = "HC1",

261
                         cluster = ~ ccode)

262


263
sd_robust_model2 <- robust_model2[,2]

264


265
# intervalo de confianza 

266


267
model2_lower = coefci(m2, df = Inf, 

268
                      vcov. = vcovCL, cluster = ~ ccode, type = "HC1")[2,1]

269


270
model2_upper = coefci(m2, df = Inf,

271
                      vcov. = vcovCL, cluster = ~ ccode, type = "HC1")[2,2]

272


273
tidy(m2, conf.int = TRUE)

274


275


276
# stargazer::stargazer(ols_model5, se = starprep(ols_model5))

277


278
texreg(list(ols_model, ols_model2),

279
       custom.coef.map = list("GPCP_g"="Growth in rainfall, t",

280
                              "GPCP_g_l"="Growth in rainfall, t-1"), digits = 3,

281
       stars = c(0.01, 0.05, 0.1),

282
       custom.gof.rows = list("Country fixed effects" = c("yes", "yes"),

283
                              "Country-specific time trends" = c("yes", "yes"),

284
                              "RMSE" = c(rmse1,rmse2)),

285
       caption = "Dependent Variable: Economic Growth Rate, t")

286


287


288
##################################################

289
#Gráfico de Coefplot del coeficiente para GPCP_g #                               

290
##################################################

291


292
#tabla de resultados

293


294
table<- matrix(0, 3, 4)

295


296
#

297


298


299
library(coefplot)

300
coefplot(ols_model, ols_model2)  + 

301
  theme_minimal() + 

302
  labs(title="Estimación de coeficientes con error estandar", 

303
       x="Estimación", 

304
       y="Variable", 

305
       caption="Elaboración propia con datos de mss_repdata")

306


307


308
ggsave("../plots/Coef_plot.png"

309
       , height = 8  # alto

310
       , width = 12  # ancho

311
       , dpi = 320   # resolución (calidad de la imagen)

312
)

313


314