################################################################################
####                                                                        ####
####                       Replication script for                           ####
#### 'Heterogeneous Effects of Women’s Schooling on Fertility, Literacy and #### 
####      Work: Evidence from Burundi’s Free Primary Education Policy'      ####
####                                                                        ####
####                             Published in                               ####
####                    Journal of African Economies                        ####
####                                                                        ####
####                              Authors                                   ####
####                   Frederik Wild & David Stadelmann                     ####
####                                                                        ####
####  The script reproduces all tables & figures displayed in the main body ####
####        of the manuscript in the order they appear in the paper.        ####
####                                                                        ####
####                           Script created                               ####
####                            February 2023                               ####
####                 using RStudio 2022.07.1 Build 554                      ####
####                                                                        ####
####                               Datasets                                 ####
####    The enclosed replication datasets are constructed from the female   ####
#### recodes of the Burundi Demographic and Health Surveys (DHS) of 2010/-  #### 
####                    2011 and 2016/2017 available at                     ####
#### https://dhsprogram.com/data/dataset/Burundi_Standard-DHS_2010.cfm and  ####
#### https://dhsprogram.com/data/dataset/Burundi_Standard-DHS_2016.cfm, as  ####
#### well as the 2005 and 2010 Household & Individual Women Modules of the  #### 
####          Multiple Indicator Cluster Surveys (MICS) available at        ####
####                    https://mics.unicef.org/surveys                     ####
####                                                                        ####
################################################################################


#######################################
########## Clean Environment ##########
#######################################
rm(list=ls())

#######################################
########## Install Libraries ##########
#######################################
# install.packages("data.table")
# install.packages("fixest")
# install.packages("ggplot2")
# install.packages("multcomp")
# install.packages("stringr")

#######################################
########### Load Libraries ############
#######################################
library(data.table) #### using data.table_1.14.2
library(fixest)     #### using fixest_0.11.1
library(ggplot2)    #### using ggplot2_3.3.6
library(multcomp)   #### using multcomp_1.4-20
library(stringr)    #### using stringr_1.4.1 

#######################################
###### Load Replication Datasets ######
#######################################
setwd(".../ReplicationFiles/Data")

dhs_fig_sample=fread("Replication_Dataset_DHS_Figs.csv")
dhs_reg_sample=fread("Replication_Dataset_DHS_Regs.csv")
mics_dhs_fig_sample=fread("Replication_Dataset_MICS_DHS_Figs.csv")

#######################################
############# Figure 1 ################
#######################################

#### Creating plotting storage ####
Fig1_left=as.data.frame(matrix(nrow=31,ncol=3))
colnames(Fig1_left)[1]="edu"
colnames(Fig1_left)[2]="yearofbirth_centered"
colnames(Fig1_left)[3]="Z"

#### Fill Storage ####
for (i in 1:31){
  yearofbirth=1960
  subset=subset.matrix(dhs_fig_sample, dhs_fig_sample$yearofbirth==yearofbirth+i)
  summary(subset)
  Fig1_left[i,1]=weighted.mean(subset$edu, subset$pweight, na.rm=T)
  Fig1_left[i,2]=-32+i
  Fig1_left[i,3]="linear regression with slope change"
}

Fig1_right=as.data.frame(matrix(nrow=6,ncol=3))
colnames(Fig1_right)[1]="edu"
colnames(Fig1_right)[2]="yearofbirth_centered"
colnames(Fig1_right)[3]="Z"

for (i in 1:6){
  yearofbirth=1991
  subset=subset.matrix(dhs_fig_sample, dhs_fig_sample$yearofbirth==yearofbirth+i)
  Fig1_right[i,1]=weighted.mean(subset$edu, subset$pweight, na.rm=T)
  Fig1_right[i,2]=-1+i
  Fig1_right[i,3]="linear regression with slope change"
}

#### Fitting data with line ####
Fig1_reg=lm(edu~Z*yearofbirth_centered, data=dhs_fig_sample, weights=pweight)
Fig1_full = cbind(dhs_fig_sample, pred = predict(Fig1_reg))

#### Plot using ggplot2 ####
Fig1_notes="Notes: Data from the 2010/2011 and the 2016/2017 Burundi DHS Female Surveys. Respondents are restricted to age 20 and older.\nWeighted using DHS sample weights."

Fig1=ggplot(data = Fig1_full, mapping = aes(x = yearofbirth_centered, y = edu), color = "none") +
  geom_line(mapping = aes(y = pred, group = factor(Z), linetype = " Local Linear Fit"), size = .9)+
  geom_vline(aes(xintercept = 0), linetype = "dotted")+
  
  stat_summary(data = Fig1_left, fun = "mean", aes(shape = " Yearly Bins"), size = 4, color = "darkgrey", geom = "point")+
  stat_summary(data = Fig1_right, fun = "mean", aes(shape = " Yearly Bins"), size = 4, color = "black", geom = "point")+
  
  theme(text = element_text(size = 14),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(size = 0.6, color = "black"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(size = 10.4, hjust = 0, vjust = 0),
        legend.position = c(0, 1), 
        legend.justification = c(0, 1),
        legend.box.margin=margin(c(10,0,0,10)),
        legend.direction = "vertical",  
        legend.box = "horizontal",
        legend.title = element_text(size = 0.01, color = "blue", face = "bold"),
        legend.text = element_text( size = 12, color = "black"),
        legend.background = element_rect(fill = FALSE, linetype = "solid", size = 0.3, color = "black"),
        legend.key = element_rect(fill = NA, size = 0.0, color = NA))+
  
  labs(x = "Year of Birth", 
       y = "Years of Schooling",
       caption = Fig1_notes)+
  
  scale_x_continuous(breaks = seq(-32,7, 2), expand = c(0.00 ,1.5),labels = seq("1960","1998", 2))+
  scale_y_continuous(breaks = seq(0, 8, 1), expand = c(0.15,-.02)) +
  scale_shape_manual(values = c(16,1))+
  scale_linetype_manual(values = c(rep("dashed", 1), rep("solid", 2)))+
  
  guides(shape = guide_legend(override.aes = list(shape = 1, size = 3, color = "black", width = 3)))

#### Export plot ####
setwd(".../ReplicationFiles/Output")
ggsave(filename = "Fig1.tiff", 
       plot = Fig1, 
       device = "tiff", 
       dpi = 200, 
       width = 225,
       height = 150, 
       units = "mm")


##########################
####### Figure 2 #########
##########################

#### Using 'poor' subsample first ####
dhs_fig_sample_poor=subset.matrix(dhs_fig_sample, dhs_fig_sample$poor==1)

#### Creating plotting storage ####
Fig2_left_poor=as.data.frame(matrix(nrow=31,ncol=3))
colnames(Fig2_left_poor)[1]="edu"
colnames(Fig2_left_poor)[2]="yearofbirth_centered"
colnames(Fig2_left_poor)[3]="Z"

#### Fill Storage ####
for (i in 1:31){
  yearofbirth=1960
  subset=subset.matrix(dhs_fig_sample_poor, dhs_fig_sample_poor$yearofbirth==yearofbirth+i)
  summary(subset)
  Fig2_left_poor[i,1]=weighted.mean(subset$edu,  subset$pweight, na.rm=T)
  Fig2_left_poor[i,2]=-32+i
  Fig2_left_poor[i,3]="linear regression with slope change"
}

Fig2_right_poor=as.data.frame(matrix(nrow=6,ncol=3))
colnames(Fig2_right_poor)[1]="edu"
colnames(Fig2_right_poor)[2]="yearofbirth_centered"
colnames(Fig2_right_poor)[3]="Z"

for (i in 1:6){
  yearofbirth=1991
  subset=subset.matrix(dhs_fig_sample_poor, dhs_fig_sample_poor$yearofbirth==yearofbirth+i)
  Fig2_right_poor[i,1]=weighted.mean(subset$edu, subset$pweight, na.rm=T)
  Fig2_right_poor[i,2]=-1+i
  Fig2_right_poor[i,3]="linear regression with slope change"
  
}

#### Fitting data with line ####
Fig2_reg_poor=lm(edu~Z*yearofbirth_centered, data=dhs_fig_sample_poor, weights=pweight)
Fig2_full_poor = cbind(dhs_fig_sample_poor, pred = predict(Fig2_reg_poor))


#### Using 'wealthy' subsample ####
dhs_fig_sample_wealthy=subset.matrix(dhs_fig_sample, dhs_fig_sample$poor==0)

#### Creating plotting storage ####
Fig2_left_wealthy=as.data.frame(matrix(nrow=31,ncol=4))
colnames(Fig2_left_wealthy)[1]="edu"
colnames(Fig2_left_wealthy)[2]="yearofbirth_centered"
colnames(Fig2_left_wealthy)[3]="Z"

for (i in 1:31){
  yearofbirth=1960
  subset=subset.matrix(dhs_fig_sample_wealthy, dhs_fig_sample_wealthy$yearofbirth==yearofbirth+i)
  summary(subset)
  Fig2_left_wealthy[i,1]=weighted.mean(subset$edu, subset$pweight, na.rm=T)
  Fig2_left_wealthy[i,2]=-32+i
  Fig2_left_wealthy[i,3]="linear regression with slope change"
}

Fig2_right_wealthy=as.data.frame(matrix(nrow=6,ncol=4))
colnames(Fig2_right_wealthy)[1]="edu"
colnames(Fig2_right_wealthy)[2]="yearofbirth_centered"
colnames(Fig2_right_wealthy)[3]="Z"

for (i in 1:6){
  yearofbirth=1991
  subset=subset.matrix(dhs_fig_sample_wealthy, dhs_fig_sample_wealthy$yearofbirth==yearofbirth+i)
  Fig2_right_wealthy[i,1]=weighted.mean(subset$edu, subset$pweight, na.rm=T)
  Fig2_right_wealthy[i,2]=-1+i
  Fig2_right_wealthy[i,3]="linear regression with slope change"
}

#### Fitting data with line ####
Fig2_reg_wealthy=lm(edu~Z*yearofbirth_centered, data=dhs_fig_sample_wealthy , weights=pweight)
Fig2_full_wealthy = cbind(dhs_fig_sample_wealthy, pred = predict(Fig2_reg_wealthy))


#### Bind datasets ####
Fig2_full_poor=cbind(Fig2_full_poor, ID=seq(1,1,1))
Fig2_full_wealthy=cbind(Fig2_full_wealthy, ID=seq(0,0,0))
Fig2_full=merge(Fig2_full_wealthy, Fig2_full_wealthy, all=TRUE)


#### Plot using ggplot2 ####
Fig2_notes=str_wrap("Notes: Data from the 2010/2011 and the 2016/2017 Burundi DHS Female Surveys. Respondents are restricted to age 20 and older. Weighted using DHS sample weights. Figures for the 'Poor' and 'Wealthy' come from seperate regressions using subsamples split  at the median household wealth level.",
                    width=127)

Fig2=ggplot(data = Fig2_full, mapping = aes(x = yearofbirth_centered, y = edu), color = "grey") +
  geom_line(data = Fig2_full_poor, mapping = aes(y = pred, group = factor(Z), linetype = " Local Linear Fit"), size = .9)+
  geom_line(data = Fig2_full_wealthy, mapping = aes(y = pred, group = factor(Z),linetype = " Local Linear Fit"), size = .9)+
  geom_vline(aes(xintercept = 0), linetype = "dotted")+
  
  stat_summary(data = Fig2_left_poor, fun = "mean", aes(shape = " Poor"), size = 4, color = "darkgrey", geom = "point")+
  stat_summary(data = Fig2_right_poor, fun = "mean", aes(shape = " Poor"), size = 4, color = "black", geom = "point")+
  stat_summary(data = Fig2_left_wealthy, fun = "mean", aes(shape = " Wealthy"), size = 4.5, color = "darkgrey", geom = "point")+
  stat_summary(data = Fig2_right_wealthy, fun = "mean", aes(shape = " Wealthy"), size = 4.5, color = "black", geom = "point")+
  
  theme(text = element_text(size = 14),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(size = 0.6, color = "black"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(size = 10.4, hjust = 0, vjust=0),
        legend.position = c(0, 1), 
        legend.justification = c(0, 1),
        legend.box.margin = margin(c(10,0,0,10)),
        legend.direction = "vertical",  
        legend.box = "horizontal",
        legend.title = element_text(size = 0.01, color = "blue", face = "bold"),
        legend.text = element_text(size = 12, color = "black"),
        legend.background = element_rect(fill = FALSE, linetype = "solid", size = 0.3, color = "black"),
        legend.key = element_rect(fill = NA, size = 0.0, color = NA))+
  
  labs(x = "Year of Birth", 
       y = "Years of Schooling", 
       caption = Fig2_notes)+
  
  scale_x_continuous(breaks = seq(-32,7, 2), expand = c(0.00 ,1.5),labels = seq("1960","1998", 2))+
  scale_y_continuous(breaks=seq(0, 8, 1), expand = c(0.03, 0.05))+
  scale_linetype_manual(values = c(rep("dashed", 1)))+
  scale_shape_manual(values = c(19,18))+
  
  guides(shape = guide_legend(override.aes = list(shape = c(1,5), size = c(3,2.5), fill = NA)))

#### Export plot ####
setwd(".../ReplicationFiles/Output")
ggsave(filename = "Fig2.tiff", 
       plot = Fig2, 
       device = "tiff", 
       dpi = 200,
       width = 225,
       height = 150, 
       units = "mm")


##########################
######## Table 1 #########
##########################

#### Column (1) ####
table1_col1 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z |factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample)
summary(table1_col1)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)

#### Column (2) ####
table1_col2 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample)
summary(table1_col2)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
summary(glht(table1_col2, linfct = c("Z + Z:yearofbirth_centered + yearofbirth_centered = 0")))

#### Column (3) ####
table1_col3 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z |factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample,
                    subset = dhs_reg_sample$poor==1)
summary(table1_col3)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)

#### Column (4) ####
table1_col4 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample,
                    subset = dhs_reg_sample$poor==1)
summary(table1_col4)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
summary(glht(table1_col4, linfct = c("Z + Z:yearofbirth_centered + yearofbirth_centered = 0")))

#### Column (5) ####
table1_col5 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z |factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample,
                    subset = dhs_reg_sample$poor==0)
summary(table1_col5)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)

#### Column (6) ####
table1_col6 = feols(edu~Z*yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey),
                    weights = dhs_reg_sample$pweight,
                    cluster = "sample_cluster",
                    data = dhs_reg_sample,
                    subset = dhs_reg_sample$poor==0)
summary(table1_col6)
weighted.mean(dhs_reg_sample$edu [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
summary(glht(table1_col6, linfct = c("Z + Z:yearofbirth_centered + yearofbirth_centered = 0")))


##########################
######## Table 2 #########
######## Panel A #########
##########################

#### Column (1) ####
table2_panelA_col1 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table2_panelA_col1, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table2_panelA_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table2_panelA_col2 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table2_panelA_col2, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table2_panelA_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table2_panelA_col3 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table2_panelA_col3, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table2_panelA_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table2_panelA_col4 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table2_panelA_col4, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table2_panelA_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table2_panelA_col5 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table2_panelA_col5, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table2_panelA_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table2_panelA_col6 = feols(literacy~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table2_panelA_col6, stage = 2)
weighted.mean(dhs_reg_sample$literacy [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table2_panelA_col6, ~ ivwald1, cluster = "sample_cluster")



##########################
######## Table 2 #########
######## Panel B #########
##########################

#### Column (1) ####
table2_panelB_col1 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table2_panelB_col1, stage = 2)
weighted.mean(dhs_reg_sample$des_children [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table2_panelB_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table2_panelB_col2 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table2_panelB_col2, stage = 2)
weighted.mean(dhs_reg_sample$des_children [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table2_panelB_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table2_panelB_col3 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table2_panelB_col3, stage = 2)
weighted.mean(dhs_reg_sample$des_children [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table2_panelB_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table2_panelB_col4 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table2_panelB_col4, stage = 2)
weighted.mean(dhs_reg_sample$des_children[dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table2_panelB_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table2_panelB_col5 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table2_panelB_col5, stage = 2)
weighted.mean(dhs_reg_sample$des_children [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table2_panelB_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table2_panelB_col6 = feols(des_children~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table2_panelB_col6, stage = 2)
weighted.mean(dhs_reg_sample$des_children [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table2_panelB_col6, ~ ivwald1, cluster = "sample_cluster")


##########################
######## Table 3 #########
######## Panel A #########
##########################

#### Column (1) ####
table3_panelA_col1 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelA_col1, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelA_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table3_panelA_col2 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelA_col2, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelA_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table3_panelA_col3 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelA_col3, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelA_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table3_panelA_col4 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelA_col4, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelA_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table3_panelA_col5 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelA_col5, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelA_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table3_panelA_col6 = feols(age1stbirth~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelA_col6, stage = 2)
weighted.mean(dhs_reg_sample$age1stbirth [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelA_col6, ~ ivwald1, cluster = "sample_cluster")


##########################
######## Table 3 #########
######## Panel B #########
##########################

#### Column (1) ####
table3_panelB_col1 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelB_col1, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelB_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table3_panelB_col2 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelB_col2, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelB_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table3_panelB_col3 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelB_col3, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelB_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table3_panelB_col4 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelB_col4, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelB_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table3_panelB_col5 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelB_col5, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelB_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table3_panelB_col6 = feols(birthbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelB_col6, stage = 2)
weighted.mean(dhs_reg_sample$birthbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelB_col6, ~ ivwald1, cluster = "sample_cluster")


##########################
######## Table 3 #########
######## Panel C #########
##########################

#### Column (1) ####
table3_panelC_col1 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelC_col1, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelC_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table3_panelC_col2 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelC_col2, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelC_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table3_panelC_col3 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelC_col3, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelC_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table3_panelC_col4 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelC_col4, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelC_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table3_panelC_col5 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelC_col5, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelC_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table3_panelC_col6 = feols(marbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelC_col6, stage = 2)
weighted.mean(dhs_reg_sample$marbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelC_col6, ~ ivwald1, cluster = "sample_cluster")


##########################
######## Table 3 #########
######## Panel D #########
##########################

#### Column (1) ####
table3_panelD_col1 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelD_col1, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelD_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table3_panelD_col2 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table3_panelD_col2, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table3_panelD_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table3_panelD_col3 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelD_col3, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelD_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table3_panelD_col4 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table3_panelD_col4, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table3_panelD_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table3_panelD_col5 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelD_col5, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelD_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table3_panelD_col6 = feols(sexbefore20~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table3_panelD_col6, stage = 2)
weighted.mean(dhs_reg_sample$sexbefore20 [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table3_panelD_col6, ~ ivwald1, cluster = "sample_cluster")



##########################
######## Table 4 #########
######## Panel A #########
##########################

#### Column (1) ####
table4_panelA_col1 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table4_panelA_col1, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table4_panelA_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table4_panelA_col2 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table4_panelA_col2, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table4_panelA_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table4_panelA_col3 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table4_panelA_col3, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table4_panelA_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table4_panelA_col4 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table4_panelA_col4, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table4_panelA_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table4_panelA_col5 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table4_panelA_col5, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table4_panelA_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table4_panelA_col6 = feols(workedlastyear~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table4_panelA_col6, stage = 2)
weighted.mean(dhs_reg_sample$workedlastyear [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table4_panelA_col6, ~ ivwald1, cluster = "sample_cluster")

mean(dhs_reg_sample$workedlastyear)


##########################
######## Table 4 #########
######## Panel B #########
##########################

#### Column (1) ####
table4_panelB_col1 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table4_panelB_col1, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table4_panelB_col1, ~ ivwald1, cluster = "sample_cluster")

#### Column (2) ####
table4_panelB_col2 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample)
summary(table4_panelB_col2, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0], na.rm=T)
fitstat(table4_panelB_col2, ~ ivwald1, cluster = "sample_cluster")

#### Column (3) ####
table4_panelB_col3 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table4_panelB_col3, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table4_panelB_col3, ~ ivwald1, cluster = "sample_cluster")

#### Column (4) ####
table4_panelB_col4 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==1)
summary(table4_panelB_col4, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==1], na.rm=T)
fitstat(table4_panelB_col4, ~ ivwald1, cluster = "sample_cluster")

#### Column (5) ####
table4_panelB_col5 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z |factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table4_panelB_col5, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table4_panelB_col5, ~ ivwald1, cluster = "sample_cluster")

#### Column (6) ####
table4_panelB_col6 = feols(work_outsideHH_cash~yearofbirth_centered+yearofbirth_centered:Z+rural+protestant+muslim+rel_other+n_siblings |factor(province)+factor(survey)|edu~Z,
                           weights = dhs_reg_sample$pweight,
                           cluster = "sample_cluster",
                           data = dhs_reg_sample,
                           subset = dhs_reg_sample$poor==0)
summary(table4_panelB_col6, stage = 2)
weighted.mean(dhs_reg_sample$work_outsideHH_cash [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0],
              dhs_reg_sample$pweight [dhs_reg_sample$Z==0 & dhs_reg_sample$poor==0], na.rm=T)
fitstat(table4_panelB_col6, ~ ivwald1, cluster = "sample_cluster")



##########################
####### Figure 3 #########
##########################

#### Row (1) ####
Fig3_reg_edulevel = feols(edulevel~treated*post*poor+age+rural |factor(region)+factor(survey),
                          weights = mics_dhs_fig_sample$pweight,
                          cluster = c("sample_cluster", "yearofbirth"),
                          data = mics_dhs_fig_sample)
summary(Fig3_reg_edulevel)

#### Creating plotting storage ####
Fig3_plot_storage=as.data.frame(matrix(nrow=2, ncol=5))
colnames(Fig3_plot_storage)=c("Y", "lwr", "est", "upr", "poor")

#### Fill Storage ####
conf_int_row1_wealthy=confint(glht(Fig3_reg_edulevel, linfct = c("treated:post = 0")),level = 0.95)$confint
Fig3_plot_storage[1,1]=2.0
Fig3_plot_storage[1,2]=conf_int_row1_wealthy[1,2]
Fig3_plot_storage[1,3]=conf_int_row1_wealthy[1,1]
Fig3_plot_storage[1,4]=conf_int_row1_wealthy[1,3]
Fig3_plot_storage[1,5]=0

conf_int_row1_poor=confint(glht(Fig3_reg_edulevel, linfct = c("treated:post + treated:post:poor = 0")),level = 0.95)$confint
Fig3_plot_storage[2,1]=1.8
Fig3_plot_storage[2,2]=conf_int_row1_poor[1,2]
Fig3_plot_storage[2,3]=conf_int_row1_poor[1,1]
Fig3_plot_storage[2,4]=conf_int_row1_poor[1,3]
Fig3_plot_storage[2,5]=1


#### Row (2) ####
Fig3_reg_literacy = feols(literacy~treated*post*poor+age+rural |factor(region)+factor(survey),
                          weights = mics_dhs_fig_sample$pweight,
                          cluster = c("sample_cluster", "yearofbirth"),
                          data = mics_dhs_fig_sample)
summary(Fig3_reg_literacy)

#### Fill Storage ####
conf_int_row1_wealthy=confint(glht(Fig3_reg_literacy, linfct = c("treated:post = 0")),level = 0.95)$confint
Fig3_plot_storage[3,1]=1.4
Fig3_plot_storage[3,2]=conf_int_row1_wealthy[1,2]
Fig3_plot_storage[3,3]=conf_int_row1_wealthy[1,1]
Fig3_plot_storage[3,4]=conf_int_row1_wealthy[1,3]
Fig3_plot_storage[3,5]=0

conf_int_row1_poor=confint(glht(Fig3_reg_literacy, linfct = c("treated:post + treated:post:poor = 0")),level = 0.95)$confint
Fig3_plot_storage[4,1]=1.2
Fig3_plot_storage[4,2]=conf_int_row1_poor[1,2]
Fig3_plot_storage[4,3]=conf_int_row1_poor[1,1]
Fig3_plot_storage[4,4]=conf_int_row1_poor[1,3]
Fig3_plot_storage[4,5]=1


#### Row (3) ####
Fig3_reg_wanted_later = feols(wanted_later~treated*post*poor+age+rural |factor(region)+factor(survey),
                              weights = mics_dhs_fig_sample$pweight,
                              cluster = c("sample_cluster", "yearofbirth"),
                              data = mics_dhs_fig_sample)
summary(Fig3_reg_wanted_later)

#### Fill Storage ####
conf_int_row1_wealthy=confint(glht(Fig3_reg_wanted_later, linfct = c("treated:post = 0")),level = 0.95)$confint
Fig3_plot_storage[5,1]=0.8
Fig3_plot_storage[5,2]=conf_int_row1_wealthy[1,2]
Fig3_plot_storage[5,3]=conf_int_row1_wealthy[1,1]
Fig3_plot_storage[5,4]=conf_int_row1_wealthy[1,3]
Fig3_plot_storage[5,5]=0

conf_int_row1_poor=confint(glht(Fig3_reg_wanted_later, linfct = c("treated:post + treated:post:poor = 0")), level = 0.95)$confint
Fig3_plot_storage[6,1]=0.6
Fig3_plot_storage[6,2]=conf_int_row1_poor[1,2]
Fig3_plot_storage[6,3]=conf_int_row1_poor[1,1]
Fig3_plot_storage[6,4]=conf_int_row1_poor[1,3]
Fig3_plot_storage[6,5]=1


#### Row (4) ####
Fig3_reg_wealthindex = feols(wealthindex~treated*post*poor+age+rural |factor(region)+factor(survey),
                             weights = mics_dhs_fig_sample$pweight,
                             cluster = c("sample_cluster", "yearofbirth"),
                             data = mics_dhs_fig_sample)
summary(Fig3_reg_wealthindex)

#### Fill Storage ####
conf_int_row1_wealthy=confint(glht(Fig3_reg_wealthindex, linfct = c("treated:post = 0")),level = 0.95)$confint
Fig3_plot_storage[7,1]=0.2
Fig3_plot_storage[7,2]=conf_int_row1_wealthy[1,2]
Fig3_plot_storage[7,3]=conf_int_row1_wealthy[1,1]
Fig3_plot_storage[7,4]=conf_int_row1_wealthy[1,3]
Fig3_plot_storage[7,5]=0

conf_int_row1_poor=confint(glht(Fig3_reg_wealthindex, linfct = c("treated:post + treated:post:poor = 0")), level = 0.95)$confint
Fig3_plot_storage[8,1]=0.0
Fig3_plot_storage[8,2]=conf_int_row1_poor[1,2]
Fig3_plot_storage[8,3]=conf_int_row1_poor[1,1]
Fig3_plot_storage[8,4]=conf_int_row1_poor[1,3]
Fig3_plot_storage[8,5]=1


#### Add identifiers #####
Fig3_plot_storage=cbind(Fig3_plot_storage, legend_est=seq(1,8,1))
Fig3_plot_storage$legend_est=c("Point Estimate (Poor)", "Point Estimate (Wealthy)", 
                               "Point Estimate (Poor)", "Point Estimate (Wealthy)", 
                               "Point Estimate (Poor)", "Point Estimate (Wealthy)",
                               "Point Estimate (Poor)", "Point Estimate (Wealthy)")

Fig3_plot_storage=cbind(Fig3_plot_storage, ci=seq(1,8,1))
Fig3_plot_storage$ci=c(" 95% CI (Wealthy)", " 95% CI (Poor)", 
                       " 95% CI (Wealthy)", " 95% CI (Poor)",  
                       " 95% CI (Wealthy)", " 95% CI (Poor)",
                       " 95% CI (Wealthy)", " 95% CI (Poor)")


#### Plot using ggplot2 ####
Fig3_notes="Notes: The 'Point Estimates' represent difference-in-differences estimates, including a 95% confidence\ninterval, for the respective dependent variable indicated on the Y-axis. Data come from the 2000 and\n2005 Multiple Cluster Indicator Surveys (MICS) as well as the 2010/2011 Burundi DHS survey, using \nthe women's recodes. The coefficients are estimated through a 'triple-difference' estimation, comparing\nindividuals aged 15-19 to individuals aged 20-24 in 2000, 2005 and 2010 for both 'wealthy' and for \n'poor' households (splitting at the median household wealth level). The estimates are weighted using\nMICS/DHS sample weights. Standard errors are clustered at the birth-year and enumeration-area level."

Fig3=ggplot(data = Fig3_plot_storage, aes(x = Y, y = est, ymin = lwr, ymax = upr, fill = legend_est))+
  geom_errorbar(data = Fig3_plot_storage, aes(ymin = lwr,ymax = upr, linetype = ci), size = 0.6, color = "black", width = .075, alpha = 1)+
  geom_hline(aes(yintercept = 0), linetype = "solid")+
  geom_point(data = Fig3_plot_storage, aes(x = Y, y = est, size = 0.5, shape = c("23", "16",
                                                                                 "23", "16",
                                                                                 "23", "16",
                                                                                 "23", "16")),
             color = c("black", "black",
                       "black", "black",
                       "black", "black",
                       "black", "black"), size = 4)+
  
  theme(text = element_text(size = 15),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        panel.border = element_rect(size = 0.6, color = "black", fill=NA),
        axis.line = element_line(size = 0, color = "black"),
        axis.title = element_text(size = 13.5, color = "black"),
        axis.text = element_text(size = 13.5, color = "black"),
        plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(size = 10.4, hjust = 0, vjust = 1),
        legend.position = c(1.01, 0.275), 
        legend.direction = "vertical",  legend.box = "vertical",
        legend.key = element_rect(fill = NA,  size = 0.0, color = NA),
        legend.justification = c(1, 1),
        legend.box.margin=margin(c(10,10,0,0)),
        legend.title = element_text(size = 0.01, color = "blue", face = "bold"),
        legend.text = element_text(size = 14, color = "black"),
        legend.background = element_rect(fill = NA, size = 0.3, linetype = "solid", color = "black"))+
  
  labs(x = "",
       y = "Difference-in-Differences Estimate",
       caption = Fig3_notes)+
  
  scale_x_continuous(breaks = seq(0.1, 2.5, .615), expand = c(0.1, -0.05),
                     labels = c("Wealth Index (1-5)", "Wanted Child Later (0/1)", 
                                "Literacy (0/1)","Education Level (0-3)"))+
  scale_y_continuous(breaks = seq(-.25, 0.5, .125), expand = expansion(add = c(.025, .35)),
                     labels = seq("-.25","0.5", .125))+
  
  guides(fill = guide_legend(override.aes = list(size = 3.5, shape = c(16, 17), order = 1, fill = c("Black", "grey"), color = c("Black", "Black"))),
         linetype = guide_legend(override.aes = list(
           linetype = c("solid", "dashed"), size = .65, alpha = 1, color= c("black", "black")), order = 2),
         shape = "none", color = "none", size = "none")+
  
  coord_flip()

#### Export plot ####
setwd(".../ReplicationFiles/Output")
ggsave(filename = "Fig3.tiff", 
       plot = Fig3, 
       device = "tiff", 
       dpi = 200, 
       width = 235,
       height = 200,
       units = "mm")