Main Equation (the model)
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} w_{t}
\left( \lambda_{1, \gamma} + \frac{p \lambda_{2, \gamma}}{R} \right) -
K \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\label{eq:1}
\tag{1}
\end{equation}\]
# Gamma is used to index gender.
npv <- function(n_male=1/2, n_female=1/2,
interest_r=interst_r_val,
wage=wage_t_val,
lambda1_male=lambda1_vals[1],
lambda1_female=lambda1_vals[2],
tax=tax_val,
saturation=saturation_val,
coverage=coverage_val,
cost_of_schooling=cost_per_student,
delta_ed_male=delta_ed_vals[,1],
delta_ed_female=delta_ed_vals[,1],
lambda2_male=lambda2_vals[1],
lambda2_female=lambda2_vals[2],
s1=0, q1=0, s2=s2_val, q2=q2_val,
periods=periods_val) {
ns <- c(n_male, n_female)
lambda1s <- c(lambda1_male, lambda1_female)
lambda2s <- c(lambda2_male, lambda2_female)
index_t <- 0:periods
delta_ed_s <- cbind(delta_ed_male, delta_ed_female)
delta_ed_s <- rbind(c(0,0), delta_ed_s, matrix(0,41, 2) )
benef <- matrix(NA, 51,2)
for (i in 1:2){
benef[,i] <- ( 1 / (1 + interest_r) )^index_t * wage *
( lambda1s[i] + saturation * lambda2s[i] / coverage )
}
res1 <- sum( ns * ( tax * apply(benef, 2, sum) -
apply( ( 1 / (1 + interest_r) )^index_t *
delta_ed_s * cost_of_schooling, 2, sum) )
) - (s2 * q2 - s1 * q1)
# browser()
return(res1)
}
Sub components:
1 - “\(r\)”
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50} \color{blue}{\left( \frac{1}{1 + r}\right)}^{t} w_{t}
\left( \lambda_{1, \gamma} + \frac{p \lambda_{2, \gamma}}{R} \right) -
K \sum_{t=0}^{50} \color{blue}{\left( \frac{1}{1 + r}\right)}^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\end{equation}\]
The real interest rate \(r\) is obtained from the interest rate on betterment bonds (0.118) minus the inflation rate (0.02).
interst_r_val <- gov_bonds - inflation
2 - “\(w_{t}\)”
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50}\left( \frac{1}{1 + r}\right)^{t} \color{blue}{ w_{t} }
\left( \lambda_{1, \gamma} + \frac{p \lambda_{2, \gamma}}{R} \right) -
K \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\end{equation}\]
\[\begin{equation}
w_t = \text{#weeks} \times w_0 (1 + g)^{Xp}(1 + \hat{\beta_1} Xp + \hat{\beta_2} Xp^2) \quad \text{for } t=10, \dots, 50
\end{equation}\]
individual in the data are assumed to enter the labor force 10 years after the (data) present day (\(w_t = 0\) for \(t<10\)). Wage at time \(t\) is the weekly starting wage in USD (\(w_0\)) that has a base growth rate equal to the per capita GDP growth (\(g\)) applied to however many years of work (\(Xp\)). In addition to this growth, the salaries are adjusted to represent a (concave) wage life cycle profile (\(1 + \hat{\beta_1} Xp + \hat{\beta_2} Xp^2\)).
2.1 - “\(w_0\)”
\[\begin{equation}
w_t = \text{#weeks} \times \color{blue}{w_0} (1 + g)^{Xp}(1 + \hat{\beta_1} Xp + \hat{\beta_2} Xp^2)
\end{equation}\]
\[\begin{equation}
w_0 = \frac{1}{ex} \sum_{l \in \{ag, ww, se\}}w_{l}\alpha_{l} \\ \quad \text{with: } \alpha_{l}= \frac{ h_{l}}{h_{ag} + h_{ww} + h_{se}}
\end{equation}\]
The initial wage in dollars (\(w_{0}\)) is a weighted average of wages for control group in agriculture, working wage, and self-employed sectors (\(ag, ww, se\)). The weights correspond to the average number of hours in each sector (\(h_l\)) relative to the sum of the average number of hours in each sector.
The wage in agriculture comes from research (Suri, 2011), the working wage comes from the data and its defined as hourly wage for the control group for those who reported more than 10 hrs of work per week. The self-employed wage (\(w_{se}\)) was constructed as follows:
\[\begin{equation}
w_{se} = \frac{ \text{Monthly self-employed profits} }{4.5 \times E[h_{se}|h_{se}>0] }
\end{equation}\]
Where both parameters (Monthly self-employed profits and self-employed hours for the control group, conditional on hrs >0 - \(E[h_{se}|h_{se}>0]\) -) come from the data (ww paper). The measure of hours in self employment used to compute wages is (\(E[h_{se}|h_{se}>0]\)) is different from the one is to compute the weights \(\alpha_l\) above. The first one captures hours of work among those actively employed in the self-employed sector, and the second one captures the average hours of work in self-employed among all the population of workin age in the sample (hence capturing the relative inportance of the self employed sector in the economy)
wage_0_f <- function(wage_ag = wage_ag_val,
wage_ww = wage_ww_val,
profits_se = profits_se_val,
hours_se_cond = hours_se_cond_val,
hours_ag = hours_ag_val,
hours_ww = hours_ww_val,
hours_se = hours_se_val,
ex_rate = ex_rate_val) {
wage_se <- profits_se / (4.5 * hours_se_cond)
wage_ls <- c(wage_ag, wage_ww, wage_se)
alpha_ls <- c(hours_ag, hours_ww, hours_se) / sum( c(hours_ag, hours_ww, hours_se) )
res1 <- 1/ex_rate * sum( wage_ls * alpha_ls )
return(res1)
}
#close to value from spreadsheet (Assumps&Panel A Calcs!B137 = 0.1481084),
#but I suspect diff due to computational precision
wage_0_val <- wage_0_f()
experience_val <- 0:periods_val - time_to_jm_val
wage_t <- function(wage_0 = wage_0_val,
growth_rate = growth_rate_val,
experience = experience_val,
coef_exp1 = coef_exp_val[1],
coef_exp2 = coef_exp_val[2]) {
res1 <- 52 * wage_0 *( ( 1 + growth_rate )^experience ) *
( 1 + coef_exp1 * experience + coef_exp2 * experience^2 ) *
ifelse(0:periods_val >= time_to_jm_val, 1, 0)
return(res1)
}
#close to value from spreadsheet (Calcs-Table 5!N21.. = 7.701634678),
#but I suspect diff due to computational precision
wage_t_val <- wage_t()
3 - “\(\lambda_{1,\gamma}\)” and “\(\lambda_{2,\gamma}\)”
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50}\left( \frac{1}{1 + r}\right)^{t} w_{t}
\left(\color{blue}{ \lambda_{1, \gamma} } + \frac{p \color{blue}{\lambda_{2, \gamma}}}{R} \right) -
K \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\end{equation}\]
\(\lambda_{1,\gamma}\) represents the estimated impact of deworming on hours of work for men a women. This two parameter are combined with a underweighted mean:
\[\begin{equation}
\lambda_{1} = \frac{1}{2} \lambda_{1,male} + \frac{1}{2} \lambda_{1,female}
\end{equation}\] Its components come from research (W@W).
\(\lambda_{2,\gamma}\) the estimated externality effect (EXPLAIN) and comes from research (W@W). Note that this parameter in not estimated by gender, so we repeat its value two times.
lambda1_vals <- rep(0.5 * lambda1_vals[1] + 0.5 *lambda1_vals[2], 2)
lambda2_vals <- rep(lambda2_val, 2)
4 - \(R\) and \(p\)
(until clarify will be using \(R\) and \(p\) as the primitive values)
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50}\left( \frac{1}{1 + r}\right)^{t} w_{t}
\left( \lambda_{1, \gamma} + \frac{\color{blue}{p} \lambda_{2, \gamma}}{\color{blue}{R}} \right) -
K \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\end{equation}\]
The coverage, \(R\), is defined as the fraction, among all neighboring schools (within 6 km), that belongs to the treatment group (last paragraph of page 9(1645) of paper). It is computed from the data as:
\[\begin{equation}
R \equiv \frac{\text{Overall Saturation} }{Q(full)}
\end{equation}\]
Overall Saturation is set at 0.511 [ASK TED ABOUT A DEFINITION FOR THIS PARAM]. And \(Q(full)\) is the take-up of the treatment with full subsidy, taken from the literature to be at 0.75 (Miguel and Kremer [add page, table, col, row]).
The saturation of the intervention, \(p\), measures the fraction of the population that is effectively usign the treatment and is defined as:
\[\begin{equation}
p = R \times Q(full) + (1 - R) \times Q(0)
\end{equation}\]
Here R is taken as given [CHECK WITH TED] and \(Q(0)\) is assinged the value of 0 per the same research source.
Note: there is a circularity on how parameters are defined here. Everything depends on Full Saturation (\(Q(full)\)) at the end of the day.
#R
#coverage_val <- saturation_val / q_full_val
#p
#saturation_val <- coverage_val * q_full_val + ( 1 - coverage_val ) * q_zero_val
5 - \(K\) and \(\Delta \overline{E}_{\gamma t}(S1,S2)\)
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50}\left( \frac{1}{1 + r}\right)^{t} w_{t}
\left( \lambda_{1, \gamma} + \frac{p \lambda_{2, \gamma}}{R} \right) -
\color{blue}{K} \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \color{blue}{ \Delta \overline{E}_{\gamma t}(S1,S2) }
\right] - \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)
\end{equation}\]
\(K\) represents the cost per student. This is calculated as the salary of the teacher plus benefits, divided by the average number of students per teacher.
\[\begin{equation}
K = \frac{\text{teacher salary} + \text{teacher benefits}}{\text{# Students}}
\end{equation}\]
cost_per_student <- (teach_sal_val + teach_ben_val) / n_students_val
For \(\Delta \overline{E}_{\gamma t}(S1,S2)\) we use a series of estimated effects the additional direct increase in secondary schooling from 1999 to 2007 obtained from [need to define the source “from Joan” in Assumps&Panel A Calcs!A93].
This series does not take into account the externality effects. To incorporate the we need another series (same source) that estimates the additional secondary schooling increase due to the externality and add it to the original series.
# Nothing here yet with delta_ed_vals, but would like to incorporate model from Joan
# delta_ed_vals <-
delta_ed_ext_total <- delta_ed_ext_vals[,1] + delta_ed_vals[,1]
Note: need to understand better the date of each component (of the model, not only this section).
6 - \(\left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right)\)
\[\begin{equation}
NPV = \sum_{\gamma} N_{\gamma} \left[
\tau \sum_{t=0}^{50}\left( \frac{1}{1 + r}\right)^{t} w_{t}
\left( \lambda_{1, \gamma} + \frac{p \lambda_{2, \gamma}}{R} \right) -
K \sum_{t=0}^{50} \left( \frac{1}{1 + r}\right)^{t} \Delta \overline{E}_{\gamma t}(S1,S2)
\right] - \color{blue}{ \left( S_{2}Q(S_{2}) - S_{1}Q(S_{1}) \right) }
\end{equation}\]
6.1 - \(S_{1}Q(S_{1}) = 0\)
There is no subsidy for deworming under the status quo.
6.2 - \(S_{2}\): complete subsidy to per capita costs of deworming.
With complete subsidy, \(S_2\) represents the total direct costs of deworming in USD. Calculated as follows
\[\begin{equation}
S_{2} = \frac{\text{Cost per person per year (KSH)} }{ex}\times \text{Additional years of treatment} \\
\end{equation}\]
6.3 - \(Q_{2}\)
The take-up with full subsidy (\(Q_2\)) comes from a previous study (Miguel and Kremer 2007) and takes the value of 0.75.
s2_val <- ( unit_cost_local_val / ex_rate_val ) * years_of_treat_val
q2_val <- q_full_val
Main results
#no externality NPV
res_npv_no_ext <- npv(lambda2_male = 0, lambda2_female = 0)
#yes externality NPV
res_npv_yes_ext <- npv(delta_ed_male = delta_ed_ext_total,
delta_ed_female = delta_ed_ext_total )
Montecarlo simulations
Describe approach to MC:
- Everything normal with sd= 0.1 * mean
- Need to work with experts to add more realistic parameters and distributions.
#34.8401
sim.data1 <- function(nsims = 1e4,
gov_bonds_vari, #Data
gov_bonds_sd,
inflation_vari,
inflation_sd,
wage_ag_vari,
wage_ag_sd,
wage_ww_vari,
wage_ww_sd,
profits_se_vari,
profits_se_sd,
hours_se_cond_vari,
hours_se_cond_sd,
hours_ag_vari,
hours_ag_sd,
hours_ww_vari,
hours_ww_sd,
hours_se_vari,
hours_se_sd,
ex_rate_vari,
ex_rate_sd,
growth_rate_vari,
growth_rate_sd,
coverage_vari,
coverage_sd,
full_saturation_vari,
full_saturation_sd,
saturation_vari,
saturation_sd,
tax_vari,
tax_sd,
unit_cost_local_vari,
unit_cost_local_sd,
years_of_treat_vari,
years_of_treat_sd,
lambda1_vari, #Research
lambda1_sd,
lambda2_vari,
lambda2_sd,
q_full_vari,
q_full_sd,
q_zero_vari,
q_zero_sd,
delta_ed_par1,
delta_ed_sd1,
delta_ed_par2,
delta_ed_sd2,
coef_exp_vari, #Guesswork
coef_exp_sd,
teach_sal_vari,
teach_sal_sd,
teach_ben_vari,
teach_ben_sd,
n_students_vari,
n_students_sd,
include_ext_vari=TRUE
) {
set.seed(1234)
#Defaoult dist: normal, default sd: 0.1* mean
## Data
gov_bonds_sim <- rnorm(n = nsims, mean = gov_bonds_vari, sd = gov_bonds_sd)
inflation_sim <- rnorm(nsims, inflation_vari, inflation_sd)
wage_ag_val_sim <- rnorm(nsims, wage_ag_vari, wage_ag_sd)
wage_ww_val_sim <- rnorm(nsims, wage_ww_vari, wage_ww_sd)
profits_se_val_sim <- rnorm(nsims, profits_se_vari, profits_se_sd)
hours_se_cond_val_sim <- rnorm(nsims, hours_se_cond_vari, hours_se_cond_sd)
hours_ag_val_sim <- rnorm(nsims, hours_ag_vari, hours_ag_sd)
hours_ww_val_sim <- rnorm(nsims, hours_ww_vari, hours_ww_sd)
hours_se_val_sim <- rnorm(nsims, hours_se_vari, hours_se_sd)
ex_rate_val_sim <- rnorm(nsims, ex_rate_vari, ex_rate_sd)
growth_rate_val_sim <- rnorm(nsims, growth_rate_vari, growth_rate_sd)
coverage_val_sim <- rnorm(nsims, coverage_vari, coverage_sd)
saturation_val_sim <- rnorm(nsims, saturation_vari, saturation_sd)
full_saturation_val_sim <- rnorm(nsims, full_saturation_vari, full_saturation_sd) ###Check here later
tax_val_sim <- rnorm(nsims, tax_vari, tax_sd)
unit_cost_local_val_sim <- rnorm(nsims, unit_cost_local_vari, unit_cost_local_sd)
years_of_treat_val_sim <- rnorm(nsims, years_of_treat_vari, years_of_treat_sd)
## Research
aux1 <- lapply(1:2,function(x) c(lambda1_vari[x],lambda1_sd[x]) )
lambda1_vals_sim <- sapply(aux1, function(x) rnorm(nsims, mean = x[1], sd = x[2]) )
lambda2_val_sim <- rnorm(nsims, lambda2_vari, lambda2_sd)
q_full_val_sim <- rnorm(nsims, q_full_vari, q_full_sd)
q_zero_val_sim <- rnorm(nsims, q_zero_vari, q_zero_sd)
## Guess work
periods_val <- 50 #Total number of periods to forecast wages
time_to_jm_val <- 10 #Time from intial period until individual join the labor force
aux2 <- lapply(1:2,function(x) c(coef_exp_vari[x],coef_exp_sd[x]) )
coef_exp_val_sim <- sapply(aux2, function(x) rnorm(nsims, mean = x[1], sd = x[2]) )
teach_sal_val_sim <- rnorm(nsims, teach_sal_vari, teach_sal_sd)
teach_ben_val_sim <- rnorm(nsims, teach_ben_vari, teach_ben_sd)
n_students_val_sim <- rnorm(nsims, n_students_vari, n_students_sd)
delta_ed_vals_sim <- sapply(delta_ed_vals[,1], function(x) rnorm(nsims, mean =
x * delta_ed_par1,
sd = delta_ed_sd1 * sd(delta_ed_vals[,1]) ) )
colnames(delta_ed_vals_sim) <- 1999:2007
delta_ed_ext_vals_sim <- sapply(delta_ed_ext_vals[,1], function(x) rnorm(nsims, mean =
x * delta_ed_par2,
sd = delta_ed_sd2 * sd(delta_ed_ext_vals[,1])))
colnames(delta_ed_ext_vals_sim) <- 1999:2007
npv_sim <- rep(NA, nsims)
#yes externality NPV
for (i in 1:nsims) {
interst_r_val <- gov_bonds_sim[i] - inflation_sim[i]
wage_0_val <- wage_0_f(wage_ag = wage_ag_val_sim[i],
wage_ww = wage_ww_val_sim[i],
profits_se = profits_se_val_sim[i],
hours_se_cond = hours_se_cond_val_sim[i],
hours_ag = hours_ag_val_sim[i],
hours_ww = hours_ww_val_sim[i],
hours_se = hours_se_val_sim[i],
ex_rate = ex_rate_val_sim[i])
experience_val <- 0:periods_val - time_to_jm_val
wage_t_val <- wage_t(wage_0 = wage_0_val,
growth_rate = growth_rate_val_sim[i],
experience = experience_val,
coef_exp1 = coef_exp_val_sim[i,1],
coef_exp2 = coef_exp_val_sim[i,2])
lambda1_vals_aux <- rep(0.5 * lambda1_vals_sim[i,1] + 0.5 * lambda1_vals_sim[i,2], 2)
lambda2_vals <- rep(lambda2_val_sim[i], 2)
#coverage_val_aux <- saturation_val_sim[i] / full_saturation_val_sim[i]
#saturation_val_aux <- full_saturation_val_sim[i] * coverage_val_sim[i]
cost_per_student <- (teach_sal_val_sim[i] + teach_ben_val_sim[i]) / n_students_val_sim[i]
q2_val_aux <- q_full_val_sim[i]
s2_val_aux <- ( unit_cost_local_val_sim[i] / ex_rate_val_sim[i] ) * years_of_treat_val_sim[i]
delta_ed_ext_total_sim <- delta_ed_vals_sim[i,] + delta_ed_ext_vals_sim[i,]
if (include_ext_vari==TRUE){
delta_ed_final <- delta_ed_ext_total_sim
}else{
delta_ed_final <- delta_ed_vals_sim[i,]
}
npv_sim[i] <- npv(interest_r = interst_r_val,
wage = wage_t_val,
lambda1_male = lambda1_vals_aux[1],
lambda1_female = lambda1_vals_aux[2],
lambda2_male = lambda2_vals[1],
lambda2_female = lambda2_vals[2],
coverage = coverage_val_sim[i],
saturation = saturation_val_sim[i],
tax = tax_val_sim[i],
cost_of_schooling=cost_per_student,
delta_ed_male = delta_ed_final,
delta_ed_female = delta_ed_final,
q2 = q2_val_aux,
s2 = s2_val_aux)
}
return(npv_sim)
}
npv_sim <- sim.data1(nsims = 1e4,
gov_bonds_vari = gov_bonds, #Data
gov_bonds_sd = 0.1 * gov_bonds,
inflation_vari = inflation,
inflation_sd = 0.1 * inflation,
wage_ag_vari = wage_ag_val,
wage_ag_sd = 0.1 * wage_ag_val,
wage_ww_vari = wage_ww_val,
wage_ww_sd = 0.1 * wage_ww_val,
profits_se_vari = profits_se_val,
profits_se_sd = 0.1 * profits_se_val,
hours_se_cond_vari = hours_se_cond_val,
hours_se_cond_sd = 0.1 * hours_se_cond_val,
hours_ag_vari = hours_ag_val,
hours_ag_sd = 0.1 * hours_ag_val,
hours_ww_vari = hours_ww_val,
hours_ww_sd = 0.1 * hours_ww_val,
hours_se_vari = hours_se_val,
hours_se_sd = 0.1 * hours_se_val,
ex_rate_vari = ex_rate_val,
ex_rate_sd = 0.1 * ex_rate_val,
growth_rate_vari = growth_rate_val,
growth_rate_sd = 0.1 * growth_rate_val,
coverage_vari = coverage_val,
coverage_sd = 0.1 * coverage_val,
saturation_vari = saturation_val,
saturation_sd = 0.1 * saturation_val,
tax_vari = tax_val,
tax_sd = 0.1 * tax_val,
unit_cost_local_vari = unit_cost_local_val,
unit_cost_local_sd = 0.1 * unit_cost_local_val,
years_of_treat_vari = years_of_treat_val,
years_of_treat_sd = 0.1 * years_of_treat_val,
lambda1_vari = c(lambda1_vals[1], lambda1_vals[2]),
lambda1_sd = c(0.17, 0.17),
lambda2_vari = lambda2_val,
lambda2_sd = 0.1 * lambda2_val,
q_full_vari = q_full_val,
q_full_sd = 0.1 * q_full_val,
q_zero_vari = q_zero_val,
q_zero_sd = 0.1 * q_zero_val,
coef_exp_vari = c(coef_exp_val[1], coef_exp_val[2]),
coef_exp_sd = c(0.001 , 0.001),
teach_sal_vari = teach_sal_val,
teach_sal_sd = 0.1 * teach_sal_val,
teach_ben_vari = teach_ben_val,
teach_ben_sd = 0.1 * teach_ben_val,
n_students_vari = n_students_val,
n_students_sd = 0.1 * n_students_val,
include_ext_vari = TRUE,
full_saturation_vari = full_saturation_val,
full_saturation_sd = 0.1 * full_saturation_val,
delta_ed_par1 = 1,
delta_ed_sd1 = 1,
delta_ed_par2 = 1,
delta_ed_sd2 = 1
)
# unit test
if (abs(sd(npv_sim) - 28.38155)>0.0001 ) {
print("Output has change")
}
npv_for_text <- paste("Median NPV:\n ", round(median(npv_sim), 2))
npv_for_text2 <- paste("SD NPV:\n ", round(sd(npv_sim), 2))
ggplot() +
geom_density(aes(x = npv_sim,
alpha = 1/2), kernel = "gau") +
geom_vline(xintercept = c(0, median(npv_sim)), col="blue") +
coord_cartesian(xlim = c(-30,100)) +
guides(alpha = "none", colour="none") +
labs(y = NULL,
x = "NPV" ,
title = "Distribution NPV of Fiscal Impacts of Deworming",
subtitle = "With Externalities")+
annotate("text", x = 2.2 * median(npv_sim), y = 0.012, label = npv_for_text, size = 6)+
annotate("text", x = 80, y = 0.004, label = npv_for_text2, size = 6)+
theme(axis.ticks = element_blank(), axis.text.y = element_blank())
