CoCalc -- PREGUNTA_1

GitHub Repository: robertopucp/1eco35_2022_2
Path: blob/main/Trabajo_final/grupo7/PREGUNTA_1_R.R
²⁷¹⁴ views
1
rm(list=ls(all=TRUE))
2

3

4
# cargando librer�as ----
5

6
librarian::shelf(
7
  tidyverse      # dplyr, tidyr, stringr, ggplot2, etc in unique library
8
  , haven        # to read datset: .dta (stata), .spss, .dbf
9
  , fastDummies  # for dummies
10
  , stargazer    # for summary and econometrics tables
11
  , sandwich     # for linear models
12
  , lmtest       # for robust standar error
13
  , estimatr     # for iv, cluster, robust standar error (LM_robust)
14
  , lfe          # for fixed effects, cluster standar error
15
  , caret        # for easy machine learning workflow (mse, rmse)
16
  , texreg       # library for export table
17
  , mfx          # probit , logit model marginal effects
18
  , xtable
19
)
20

21

22
# directorio y base de datos ----
23
user <- Sys.getenv("USERNAME")  # username
24
setwd( paste0("C:/Users/",user,"/Documents/GitHub/Jose_Pastor_r_py_jl/Python & R & Julia/Trabajo Final") ) # set directorio
25
repdata <- read_dta("mss_repdata.dta")
26

27

28
#############################################
29
# exportar en latex una tabla de estad�sticas (media, desviaci�n est�ndar y total de observaciones) ----
30
table1 <- repdata %>% dplyr::select(NDVI_g, tot_100, trade_pGDP, 
31
                                    pop_den_rur, land_crop, va_agr, 
32
                                    va_ind_manf) %>% as.data.frame()
33

34
list_vars <- c("NDVI_g", "tot_100", "trade_pGDP", "pop_den_rur", "land_crop", "va_agr", "va_ind_manf")
35

36
stargazer(table1, title = "Descriptive Statistics", digits = 2, # decimales con 2 digitos
37
          covariate.labels = list_vars,                         # Lista de etiquetas
38
          summary.stat = c("mean", "sd", "n"),                  # se especifica el orden de los estadísticos
39
          min.max = F,                                          # borrar el estad�stico de maximo y minimo
40
          notes = "Source: World Bank's World Development Indicators (WDI)."
41
          ,notes.append = FALSE,                                # TRUE append the significance levels
42
          notes.align = 'l')                                    # notes.align = 'l' : notas a la izquierda
43

44

45
#####################################
46
# Replicar la tabla 3 (pag 737). ---- 
47

48

49
# creando dummy por pa�s para efectos fijos
50
repdata  <- dummy_cols(repdata, select_columns = 'ccode')
51

52

53
# creando variable temporal para que vaya de 1, 2, as�...
54
repdata$time_year <- repdata$year - 1978  
55

56

57
# creando variable trend effect
58
index <- grep("ccode_", colnames(repdata))
59
list_vars <- names(repdata)[index]
60
list_vars[1]
61

62
i = 1
63
while(i < 42){
64
  var <- paste0(list_vars[i],"_","time")
65
  repdata[var]  <- repdata[list_vars[i]]*repdata["time_year"]
66
  i = i + 1
67
}
68

69
repdata$ccode_factor <- as.factor(repdata$ccode)         # variable categ�rica
70

71
index_country_time <- grep("_time$", colnames(repdata))  # grep() devuelve las posiciones
72
country_time_trend <-names(repdata)[index_country_time]
73

74

75

76
#####################
77
####  Modelo 1  #####
78

79
# Usando as.formula pues la f�rmula del modelo se hace extensa, luego de paste lo convertir� en f�rmula
80
model1 <- as.formula(
81
  paste("any_prio", "~", paste("GPCP_g", "GPCP_g_l", "ccode_factor", paste(country_time_trend, collapse = "+"), 
82
                               sep="+")
83
  )
84
)
85

86

87
# Usando LM y luego coestest para error standar robusto
88
m1 <- lm(model1, data = repdata)
89
lm_rmse1 <- round(RMSE(m1$fitted.values, repdata$gdp_g ), 2 )
90
robust_model1 <- coeftest(m1,
91
                          vcov = vcovCL,
92
                          type = "HC1",
93
                          cluster = ~ ccode)
94

95
sd_robust_model1 <- robust_model1[,2]
96

97
coefci(m1, df = Inf, vcov. = vcovCL, cluster = ~ ccode, type = "HC1")
98

99

100

101

102

103
#####################
104
####  Modelo 2  #####
105

106
# Usando as.formula pues la f�rmula del modelo se hace extensa, luego con paste se convertir� en f�rmula
107
model2 <- as.formula(
108
  paste("war_prio", "~", paste("GPCP_g", "GPCP_g_l", "ccode_factor", paste(country_time_trend, collapse = "+"), 
109
                               sep="+")
110
  )
111
)
112

113
# Usando LM y luego coestest para error standar robusto
114
m2 <- lm(model2, data = repdata)
115
lm_rmse2 <- round(RMSE(m2$fitted.values, repdata$gdp_g ), 2 )
116
robust_model2 <- coeftest(m2,
117
                          vcov = vcovCL,
118
                          type = "HC1",
119
                          cluster = ~ ccode)
120

121
sd_robust_model2 <- robust_model2[,2]
122

123
coefci(m2, df = Inf, 
124
       vcov. = vcovCL, cluster = ~ ccode, type = "HC1")
125

126

127

128

129
# Exportando tabla a latex
130

131
stargazer( m1, m2,                           # los 5 modelos anterioremente calculados
132
           se=list(sd_robust_model1, sd_robust_model2),
133
           dep.var.labels = c(""),
134
           title = "Rainfall and Civil Conflict (Reduced-Form)",
135
           keep = c("GPCP_g","GPCP_g_l"),   # las variables que mostrar� en la tabla
136
           covariate.labels=c("Growth in rainfall, t","Growth in rainfall, t-1",
137
                              "Growth in rainfall, t+1"),   # label de las variables que mostrar� en la tabla
138
           align = T, no.space = T,
139
           add.lines=list(c("Country fixed effects", "yes","yes"),
140
                          c("Country-specific time trends","yes","yes"),
141
                          c("Root mean square error",lm_rmse1,lm_rmse2)),
142
           keep.stat = c("rsq","n"),
143
           notes.append = FALSE, notes.align = "l",
144
           notes ="Huber robust standard errors are in parentheses", 
145
           style = "qje"         # estilo de tabla como los journals of economics
146
)
147

148

149

150

151
###################################################################################################
152
# Graficar el Coeft plot del coeficiente estimado que corresponde a la variable GPCP g.
153
# El gr�fico debe incluir el coeficiente e intervalo de confianza respectivo de cada modelo.
154

155

156
######### Coefplot Modelo 1 ####################
157
model1 <- lm(model1, data = repdata) 
158

159
model1.tab <- coeftest(model1, vcov=vcovHC(model1, type='HC1')) # vcovHC : heterocedasticidad
160
# type='HC1': matriz var y cov de Huber-White
161
model1.tab
162

163
model1_coef <- model1.tab[2,1]
164
model1_coef_se = model1.tab[2,2]
165

166
model1_lower = coefci(model1, df = Inf, vcov. = vcovHC, type = "HC1")[2,1]   # fila 2, columna 1
167
model1_upper = coefci(model1, df = Inf, vcov. = vcovHC, type = "HC1")[2,2]
168

169

170

171
######## Coefplot Modelo 2 ###################
172
model2 <- lm(model2, data = repdata) 
173

174
# lmtest: coeftest (errores standar robustos para heterocedasticidad)    
175
model2.tab <- coeftest(model2, vcov=vcovHC(model2, type='HC1')) # vcovHC : heterocedasticidad
176
# type='HC1': matriz var y cov de Huber-White
177
model2.tab
178

179
model2_coef <- model2.tab[2,1]
180
model2_coef_se = model2.tab[2,2]
181

182
model2_lower = coefci(model2, df = Inf, vcov. = vcovHC, type = "HC1")[2,1]   # fila 2, columna 1
183
model2_upper = coefci(model2, df = Inf, vcov. = vcovHC, type = "HC1")[2,2]
184

185

186
# Tabla de resultados
187
table<- matrix(0, 2, 4) # matriz llena de ceros con 3 filas y 4 columnas
188

189
table[1,1]<- model1_coef
190
table[1,2]<- model1_coef_se
191

192
table[2,1]<- model2_coef
193
table[2,2]<- model2_coef_se
194

195
table[1,3]<- model1_lower
196
table[1,4]<- model1_upper
197

198
table[2,3]<- model2_lower
199
table[2,4]<- model2_upper
200

201
colnames(table)<- c("Estimate","se","lower_bound","upper_bound")
202
rownames(table)<- c("OLS any_prio", "OLS war_prio")
203

204
# de matriz a dataframe (tab), para poder graficar
205
tab = as.data.frame(table)
206

207
# Exportaci�n a Latex (overleaf)
208
xtable(table)
209

210

211

212
###  Coef-plot  ###
213

214
theme_set(theme_bw(20)) # tema con fondo blanco y cuadro con bordes en negro
215

216
# aes: ejes
217
options(repr.plot.width = 8, repr.plot.height =5)  # tama�o del gr�fico
218

219
tab %>% ggplot(aes(x=rownames(tab), y=Estimate)) +
220
  geom_point(size=2, color = 'red') +
221
  geom_errorbar(aes(ymin=lower_bound, ymax=upper_bound), width = 0.05, color="black", size = 0.8) +
222
  labs(x="", y="") + ggtitle("Global Precipitation Climatology Project coefficient (95% CI)")  +
223
  theme(text=element_text(size =15), plot.title = element_text(hjust = 0.5)) +
224
  geom_hline(yintercept=0, linetype="dashed", color = "black", size=1) +
225
  scale_x_discrete(limits = c("OLS any_prio", "OLS war_prio")) + # order x - axis
226
  theme(panel.grid.major = element_blank(),   # borra las cuadr�culas en el fondo
227
        panel.grid.minor = element_blank())
228

229

230
ggsave("Coef_plot.png"
231
       , height = 8  # alto
232
       , width = 12  # ancho
233
       , dpi = 320   # resoluci�n
234
)
235

236

237

238

239

240

241
Product

Resources

Company
