*    -------------------------------------------------- 
*                                                       
*         $$$$$ $$$$$$$$    $    $$$$$$$$    $          
*        $$  $$ $  $$  $   $$    $  $$  $   $$          
*        $$   $    $$      $$$      $$      $$$         
*        $$$       $$     $ $$      $$     $ $$         
*          $$$     $$     $  $$     $$     $  $$        
*           $$$    $$     $$$$$     $$     $$$$$        
*        $   $$    $$    $   $$     $$    $   $$        
*        $$  $$    $$    $    $$    $$    $    $$       
*        $$$$$    $$$$  $$$  $$$$  $$$$  $$$  $$$$      
*                                                       

*---------------------------------------------------------------
*-复现论文:
*	Corno, Lucia, N. Hildebrandt, and A. Voena. 
*   "Age of Marriage, Weather Shocks, 
*	and the Direction of Marriage Payments."
*   Econometrica 88(2020). 
*---------------------------------------------------------------

*----------
*-主要模型和方法
  *- 固定面板效应
  
*-注意：执行后续命令之前，请先执行如下命令
  global path "D:\Stata16\FinalWork\Paper_Corno_Lucia_2020_replication"
  cd $path
  cap mkdir data
  cap mkdir refs
  cap mkdir out
  global D    "$path\Data"    //范例数据
  global R    "$path\refs"    //参考文献
  global Output  "$path\Output"     //结果：图形和表格
  cd "$path"                     //设定工作路径
  set scheme s2color
  
clear
set more off

*目录：
* 一、降雨与农作产出、收入和家庭消费
* 		Figure A1 撒哈拉以南非洲和印度历史干旱和洪水的直方图
*		Table A2 回归分析：【撒哈拉以南非洲】 干旱 与 农作物产量、收入和家庭消费
*		Table A3 回归分析：【印度】 干旱 与 农作物产量、收入和家庭消费
*		Figure A2: 农作物产出与降水分位数分布
* 二、干旱对童婚可能性的影响
*		Table 1: 干旱对结婚时间影响
*		Table 2: 干旱对结婚时间影响
*		Table 3: 干旱对生育率影响
* 三、稳健性检验
*		Figure 2: 干旱的不同定义：百分位截取的不同
*		Table A5: 干旱的不同定义：总体干旱水平
*		Table B9: 洪水（高降水概率对童婚的影响）
*		Table B11: 过去和未来的时间效应
*		Table B12: 以更大的地理单元格单元格聚类
* 四、异质性分析 
*		Table4: 撒哈拉以南地区异质性分析
*		Table5: 印度地区异质性分析


********************************************************************************
*
*
* 一、降雨与农作产出、收入和家庭消费
*
*
********************************************************************************


********************************************************************************
* Figure A1 撒哈拉以南非洲和印度历史干旱和洪水的直方图
********************************************************************************
clear
*Important globals for plotting
global yearStart 	1960
global yearEnd 		2010
global yearStep 	5

* Programs
* Make bar plots based on the type of weather shock input.
cap program drop makeShockBarPlot
program define makeShockBarPlot
	args shockVar shockVarDisp region regionDisp
	
	summ lr_mean_pct_`shockVar'_cy, meanonly
	global lr_mean_`shockVar'_disp = trim("`: di %-10.2f `r(mean)''") + "%"

	#delimit ;
		graph twoway 
			(bar pct_`shockVar'_cy calendar_year, color(gs11))
			(line lr_mean_pct_`shockVar'_cy calendar_year, lcolor(black) lpattern(dash)),
			title("Prevalence of `shockVarDisp': `regionDisp'")
			ytitle("Percent with `shockVarDisp'") xtitle("Year")
			xscale(range(${yearStart} ${yearEnd}))
			xlabel(${yearStart}(${yearStep})${yearEnd})
			scheme(s1mono)
			graphregion(color(white)) bgcolor(white)
			legend(order(2 "Mean ${yearStart}-${yearEnd}: ${lr_mean_`shockVar'_disp}") 
				ring(0) pos(1) col(4) region(lstyle(none) color(none)) symxsize(6))
			name("bar_`region'_`shockVar'", replace)
			note("Figure shows the percent of grid cells in `regionDisp' with `shockVar' each year."
				"The dashed black line indicates the long run mean from ${yearStart} to ${yearEnd}.", span)
		;
		gr export "Output/bar_`region'_`shockVar'.pdf",
			replace name("bar_`region'_`shockVar'");
	#delimit cr
end

* SSA
use "Data/droughts.dta" if data=="SSA"
* Call program and make plots.
makeShockBarPlot drought Drought SSA SSA

* India
use "Data/droughts.dta" if data=="India"
* Call program and make plots.
makeShockBarPlot drought Drought IN India

window manage close graph _all



********************************************************************************
* Table A2 回归分析：【撒哈拉以南非洲】 干旱 与 农作物产量、收入和家庭消费
********************************************************************************
global cropsSSA Maize Sorghum Millet Rice Wheat MainCrops
use "Data/SSA_crops_income", clear
lab var drought_cy "Drought"

**Crop regressions**

foreach crop of global cropsSSA {
	lab var log_yield`crop' "`crop'"
	areg log_yield`crop' drought_cy i.calendar_year ///
		, absorb(country_num) vce(cluster country_num)
	if "`crop'"=="Maize" {
	outreg2 drought_cy using "Output/tableA2.xlsx", /// 
		keep(drought_cy) replace excel label adjr2 nocons
		}
	else {
	outreg2 drought_cy using "Output/tableA2.xlsx", /// 
		keep(drought_cy) append excel label adjr2 nocons
		}
	}

gen gdp_pc = real(gdp_pc_str)
gen loggdp_pc = log(gdp_pc)
lab var loggdp_pc "GDP per capita"
gen hh_consumption_agg = real(hh_consumption_agg_str)
gen loghh_consumption_agg = log(hh_consumption_agg)
lab var loghh_consumption_agg "HH consumption"

foreach y in gdp_pc hh_consumption_agg {
	areg log`y' drought_cy i.calendar_year, absorb(country_num) vce(cluster country_num)
	outreg2 drought_cy using "Output/tableA2.xlsx", ///
		keep(drought_cy) append excel label adjr2 nocons
	}

	
********************************************************************************
* Table A3 回归分析：【印度】 干旱 与 农作物产量、收入和家庭消费
********************************************************************************
global cropsIN Rice Wheat Jowar Maize Bajra MainCrops
use "Data/India_crops", clear
lab var drought_cy "Drought"

foreach crop of global cropsIN {
	lab var log_yield`crop' "`crop'"
	areg log_yield`crop' drought_cy i.calendar_year ///
		, absorb(district_num) vce(cluster district_num)
	if "`crop'"=="Rice" {
		outreg2 drought_cy using "Output/tableA3.xls", ///
			keep(drought_cy) replace excel label adjr2 nocons
		}
	else {
		outreg2 drought_cy using "Output/tableA3.xls", ///
			keep(drought_cy) append excel label adjr2 nocons
		}
	}
	
**NSDP p.c. portion**
use "Data/India_nsdp", clear
gen log_nsdp_pc = log(nsdp_pc)
lab var log_nsdp_pc "NSDP per cap."
replace drought_cy = drought_cy_lag1
lab var drought_cy "Drought"
areg log_nsdp_pc drought_cy i.year, absorb(state) vce(cluster state)
outreg2 drought_cy using "Output/tableA3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**NSS portion**	
use "Data/India_NSS", clear

**Generate real consumption, i.e. adjusted using deflator**
**Note that this generates consumption in 2011 INR**
gen real_cons = (cons / deflator) * 100
gen log_real_cons = log(real_cons)
lab var real_cons "Consumption (2011 INR)"
lab var log_real_cons "HH consumption"

**Generate drought variable**

gen drought_cy = (total_cy_xtile_lag1 < 0.15) if total_cy_xtile_lag1 != . // First lag

lab var drought_cy "Drought"

areg log_real_cons drought_cy i.round if round <= 54, ///
	absorb(district_num) vce(cluster district_num)

outreg2 drought_cy using "Output/tableA3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	

********************************************************************************	
* Figure A2: 农作物产出与降水分位数分布
********************************************************************************
* Globals
global cropsSSA MainCrops 
global cropsIN MainCrops Rice 

* Number of bootstraps.
global numSim 500

* Make rain quantile labels for LaTeX tables.
* Remember to specify the reference distribution for rain quantiles.
* Leave blank for the "raw" quantiles (based on empirical uniform dist)
local quantDist Gam_
global lbldecile
global lblquintile
global lblvingtile
forval iter = 1/20 {
	if `iter' <= 5 {
		global lblquintile $lblquintile `iter'.rain_`quantDist'quintile \hspace{16pt}`iter'
	}
	
	if `iter' <= 10 {
		global lbldecile $lbldecile `iter'.rain_`quantDist'decile \hspace{16pt}`iter'
	}
	
	global lblvingtile $lblvingtile `iter'.rain_`quantDist'vingtile \hspace{16pt}`iter'
}

di "$lblvingtile"
di "$lbldecile"
di "$lblquintile"

global vingtileShort vin
global decileShort dec
global quintileShort qui

global vingtileRef 10 // Actually 10 and 11 but I will redefine 11 to be 10 later.
global decileRef 5
global quintileRef 3

global vingtileRefText The base categories (omitted) are the 10\textsuperscript{th} and 11\textsuperscript{th} rain vingtiles.
global decileRefText The base category (omitted) is the 5\textsuperscript{th} rain decile.
global quintileRefText The base category (omitted) is the 3\textsuperscript{rd} rain quintile.

************* Custom programs *************

* Generate quantiles defined based on the gamma rainfall distribution
cap program drop genGammaQuantiles
program define genGammaQuantiles
	gen total_cy_ptile = total_cy_xtile*100

	* Divide percentile by bin size (e.g. vingtile has bin size of 5).
	* Take the floor of the above and add 1 if there is a remainder in the division.
	local binSize 5
	gen rain_Gam_vingtile = floor(total_cy_ptile/`binSize') + 1*((total_cy_ptile - floor(total_cy_ptile/`binSize'))>0)
	tab rain_Gam_vingtile
	scatter rain_Gam_vingtile total_cy_ptile
	replace rain_Gam_vingtile = 10 if rain_Gam_vingtile==11
	tab rain_Gam_vingtile

	local binSize 10
	gen rain_Gam_decile = floor(total_cy_ptile/`binSize') + 1*((total_cy_ptile - floor(total_cy_ptile/`binSize'))>0)
	tab rain_Gam_decile

	local binSize 20
	gen rain_Gam_quintile = floor(total_cy_ptile/`binSize') + 1*((total_cy_ptile - floor(total_cy_ptile/`binSize'))>0)
	tab rain_Gam_quintile
end

* This program plots the regression coefficients of log annual crop yields on rain quantiles
cap program drop plotCropRainCoef
program define plotCropRainCoef
	args modName crop regionShort clustVar
	preserve
		parmest, nores
		gen xvar = real(regexs(1)) if regexm(parm, "^([0-9]*)(b*)\.rain_(.*ile)$")
		keep if xvar != .
		gen temp_base = regexs(2) if regexm(parm, "^([0-9]*)(b*)\.rain_(.*ile)$")
		gen temp_quantile = regexs(3) if regexm(parm, "^([0-9]*)(b*)\.rain_(.*ile)$")
		* drop the "Gam_" modifer (to indicate gamma quantiles) if it exists
		replace temp_quantile = subinstr(temp_quantile, "Gam_", "", .)
		
		summ xvar, meanonly
		local maxCat `r(max)'
		summ xvar if ~missing(temp_base), meanonly
		local baseCat `r(max)'
		local baseCatText `baseCat' `"(Omitted)"'
		if `baseCat'==10 {
			local baseCatText 10 "      (Omitted)" 11 " " 
		}
		
		local quantile = temp_quantile + "s"
		local Quantile = proper("`quantile'")
		
		#delimit ;
		twoway
			(bar estimate xvar, base(0) barwidth(0.7) color(gs10) finten(inten100))
			(rcap max95 min95 xvar, lcolor(blue))
			,
			title("Coefficient for rain `quantile' and 95% CI in `regionShort': `crop'" " ")
			ytitle("Coefficient") xtitle("Rain `Quantile'")
			ylabel(, format(%4.1f))
			xlabel(1(1)`maxCat' `baseCatText')
			legend(off) 
			yline(0, lcolor(black) lpattern(dash))
			note("Regression includes year and `clustVar' FEs. SEs clustered at `clustVar' level.", span)
			name("`modName'", replace)
		;
		#delimit cr
		gr export "Output/`modName'.pdf", replace name("`modName'")
	restore
end

****** SSA rain quantile plots
use "Data/SSA_crops_income", clear
xtset countrycode calendar_year

/* SSA: crop yields vs. rain quantiles (vingtiles, deciles, and quintiles)*/
* Generate quantiles defined based on the gamma rainfall distribution
genGammaQuantiles

* Remember to specify the reference distribution for rain quantiles.
* Leave blank for the "raw" quantiles (based on empirical uniform dist)
local quantDist Gam_
foreach crop of global cropsSSA {
	*** All SSA
	* Vingtiles
	local modName SSAvin1`crop'
	eststo `modName': areg log_yield`crop' ib10.rain_`quantDist'vingtile i.calendar_year ///
		, absorb(country_num) vce(cluster country_num)
	plotCropRainCoef `modName' `crop' SSA country
	
	*** DHS SSA sample
	* Vingtiles
	local modName SSAvin2`crop'
	eststo `modName': areg log_yield`crop' ib10.rain_`quantDist'vingtile i.calendar_year ///
		if in_DHS==1, absorb(country_num) vce(cluster country_num)
	plotCropRainCoef `modName' `crop' SSA country

}

****** India regressions and rain quantile plots
use "Data/India_crops", clear
xtset district_num calendar_year

/*INDIA: crop yields vs. rain quantiles (vingtiles, deciles, and quintiles)*/
* Generate quantiles defined based on the gamma rainfall distribution
genGammaQuantiles

* Remember to specify the reference distribution for rain quantiles.
* Leave blank for the "raw" quantiles (based on empirical uniform dist)
local quantDist Gam_
foreach crop of global cropsIN {
	* Vingtiles
	local modName INvin`crop'
	eststo `modName': areg log_yield`crop' ib10.rain_`quantDist'vingtile i.calendar_year ///
		, absorb(district_num) vce(cluster district_num)
	plotCropRainCoef `modName' `crop' IN district

}

window manage close graph _all
log close



********************************************************************************
*
*
* 二、干旱对童婚可能性的影响
*
*
********************************************************************************	
	
	
	
********************************************************************************	
* Table 1: 干旱对结婚时间影响
********************************************************************************
* 样本选择：<=24岁

**SSA**

use "Data/SSA_duration_final" ///
	if inrange(year,12,24), clear
keep if war_upto24==0
lab var drought_cy "Drought"
areg married_survival drought_cy i.birth_year i.year if year<=24 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table1.xls", ///
	keep(drought_cy) replace excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country if year<=24 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table1.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country i.country#i.cohort_num  if year<=24 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table1.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**India**

use "Data/India_duration_final" if inrange(year,12,24), clear
lab var drought_cy "Drought"
areg married_survival drought_cy i.birth_year i.year ///
	if year<=24, ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy using "Output/table1.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.state_gec_num#i.cohort_num ///
	if year<=24, ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy using "Output/table1.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
	
********************************************************************************	
* Table 2: 干旱对结婚时间影响
********************************************************************************
* 样本选择：<=17岁

use "Data/SSA_duration_final" ///
	if inrange(year,12,17), clear
keep if war_upto17==0
lab var drought_cy "Drought"
areg married_survival drought_cy i.birth_year i.year if year<=17 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table2.xls", ///
	keep(drought_cy) replace excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country if year<=17 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table2.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country i.country#i.cohort_num if year<=17 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table2.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**India**
use "Data/India_duration_final" if inrange(year,12,17), clear
lab var drought_cy "Drought"
areg married_survival drought_cy i.birth_year i.year ///
	if year<=17 & data=="DHS_98", ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy using "Output/table2.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.state_gec_num#i.cohort_num ///
	if year<=17 & data=="DHS_98", ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy using "Output/table2.xls", ///
	keep(drought_cy) append excel label adjr2 nocons	

	
********************************************************************************	
* Table 3: 干旱对生育率影响
********************************************************************************
use "Data/SSA_duration_final" ///
	if year>=12 & year<18, clear
keep if war_upto17==0
lab var drought_cy "Drought"
gen haschild=0 
replace haschild=1 if year==v212
replace haschild=. if year>v212
areg haschild drought_cy i.birth_year i.year ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table3.xls", ///
	keep(drought_cy) replace excel label adjr2 nocons
areg haschild drought_cy i.birth_year i.year i.country ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg haschild drought_cy i.birth_year i.year i.country i.country#i.cohort_num ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**India**
use "Data/India_duration_final" ///
	if data=="DHS_98" & year>=12 & year<18, clear
lab var drought_cy "Drought"
gen haschild=0 
replace haschild=1 if year==v212
replace haschild=. if year>v212
areg haschild drought_cy i.birth_year i.year , absorb(district) vce(cluster district)
outreg2 drought_cy using "Output/table3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg haschild drought_cy i.birth_year i.year i.state_gec_num#i.cohort_num, absorb(district) vce(cluster district)
outreg2 drought_cy using "Output/table3.xls", ///
	keep(drought_cy) append excel label adjr2 nocons

	
********************************************************************************
*
*
* 三、稳健性检验
*
*
********************************************************************************

********************************************************************************
* Figure 2: 干旱的不同定义：百分位截取的不同
********************************************************************************

clear
************ IMPORTANT: Remember to choose age cutoff ************
global ageLowerLim 12
* Whether to limit age up till 17 or 24
global ageUpperLim 17
*** Whether to use weights in India regressions
global INregWeight 
*global INregWeight [pweight = wgt_new]
*** Whether to use SSA high BP
global ssaHighBp 0
*** Whether to include flood (1 means yes, otherwise no)
global includeFlood 0

*-----------------  Robustness to Drought and Flood Cutoffs -----------------*
**********************************  SSA  ***************************************
**** Use surveys conducted >= 1994 
use "Data/SSA_duration_final" ///
	if inrange(year,${ageLowerLim},${ageUpperLim}) & war_upto${ageUpperLim}==0, clear

di " "
di "*****************************************"
di "* Robustness to cutoffs: SSA (High BP)"
di "*****************************************"
di " "

if $ssaHighBp == 1 {
	global highBpControl "& highbp_ge50pct==1"
	global plotTitle "SSA, BP {&ge} 50%"
	global outputStub SSAhighbp
}
else {
	global highBpControl 
	global plotTitle "SSA"
	global outputStub SSAfull
}

* Initialize beta and SE matrices for the new definitions of drought_cy and flood_cy.
matrix matBetas = J(300,2,1)
matrix matSE = J(300,2,1)
matrix mat_drought_cutoff = J(300,1,1)
matrix mat_flood_cutoff = J(300,1,1)
global droughtCutoffMin 1
global droughtCutoffMax 45
global droughtCutoffNum = $droughtCutoffMax - $droughtCutoffMin + 1

**For coefplot**
matrix coef_SSA = J(1,$droughtCutoffNum ,.)
matrix ci_SSA = J(2,$droughtCutoffNum ,.)

* Loop over new cutoffs, run regressions, and save coefs + SEs
local counter = 0
local cutoffRng ${droughtCutoffMin}(1)${droughtCutoffMax}
forvalues cutoff = `cutoffRng' {
	local ++counter
	di "------------------------------------------"
	di "The cutoff is " %-10.2f `cutoff'/100
	di "------------------------------------------"
	di " "
	
	cap drop new_drought
	gen new_drought = (total_cy_xtile <= round(`cutoff'/100, 10^-2))
	
	if $includeFlood == 1 {
		cap drop new_flood
		gen new_flood = (total_cy_xtile > 1 - round(`cutoff'/100, 10^-2))
		
		* Main regression specification: drought + flood
		areg married_survival new_drought new_flood i.birth_year i.year i.country ///
			if year<=$ageUpperLim ${highBpControl} [pweight = wgt_new], ///
			absorb(grid_cell) vce(cluster grid_cell)
			
		matrix matBetas[`counter',2] = _b[new_flood]
		matrix matSE[`counter',2] = _se[new_flood]
		matrix mat_flood_cutoff[`counter',1] = 1 - round(`cutoff'/100, 10^-2)	
	}
	else {
		* Main regression specification: drought
		areg married_survival new_drought i.birth_year i.year i.country ///
			if year<=$ageUpperLim ${highBpControl} [pweight = wgt_new], ///
			absorb(grid_cell) vce(cluster grid_cell)
	}
	
	* Drought is always included.
	matrix matBetas[`counter',1] = _b[new_drought]
	matrix matSE[`counter',1] = _se[new_drought]
	matrix mat_drought_cutoff[`counter',1] = round(`cutoff'/100, 10^-2)	
	matrix coef_SSA[1,`counter'] = _b[new_drought]
	matrix ci_SSA[1,`counter'] = (_se[new_drought] * (-1.96)) + _b[new_drought]
	matrix ci_SSA[2,`counter'] = (_se[new_drought] * 1.96) + _b[new_drought]
}

cap drop new_drought 
cap drop new_flood

preserve
	drop _all
	svmat matBetas
	svmat matSE
	svmat mat_drought_cutoff
	svmat mat_flood_cutoff
	
	drop if matBetas1==1

	***** Drought
	* Rename vars, convert to %, then create 95% upper and lower CI
	rename matBetas1 drought
	rename matSE1 drought_se
	rename mat_drought_cutoff1 drought_cutoff
	replace drought_cutoff = drought_cutoff*100	
	gen drought_CIL = drought - 1.96*drought_se
	gen drought_CIU = drought + 1.96*drought_se
	
	* Make plots of coef (and 95% CI) vs. cutoff.
	#delimit ;
		twoway				
			(rarea drought_CIL drought_CIU drought_cutoff, color(gs13))
			(connected drought drought_cutoff, mcolor(black) lcolor(black) lpattern(solid) msymbol(S) msize(0.8))
			,
			title("Robustness to Rainfall Cutoff: ${plotTitle}")
			ytitle("Effect of Drought") xtitle("Cutoff Percentile for Drought Definition") xscale(range($droughtCutoffMin $droughtCutoffMax))
			xlabel(0(5)$droughtCutoffMax)
			ylabel(, format(%6.3f))
			scheme(s1mono)
			graphregion(color(white)) bgcolor(white)
			legend(off)
			yline(0, lcolor(black) lpattern(dash))
			name("rob_${outputStub}_drought_cutoff", replace)
		;
		gr export "Output/rob_${outputStub}_drought_cutoff.pdf", replace name("rob_${outputStub}_drought_cutoff");
	#delimit cr
	
	if $includeFlood == 1 {
		***** Flood
		* Rename vars, convert to %, then create 95% upper and lower CI
		rename matBetas2 flood
		rename matSE2 flood_se
		rename mat_flood_cutoff1 flood_cutoff
		replace flood_cutoff = flood_cutoff*100
		gen flood_CIL = flood - 1.96*flood_se
		gen flood_CIU = flood + 1.96*flood_se
		
		* Make plots of coef (and 95% CI) vs. cutoff.
		#delimit ;
			twoway				
				(rarea flood_CIL flood_CIU flood_cutoff, color(gs13))
				(connected flood flood_cutoff, mcolor(black) lcolor(black) lpattern(solid) msymbol(S) msize(0.8))
				,
				title("Robustness to Rainfall Shocks Cutoff: ${plotTitle}")
				ytitle("Effect of Flood") xtitle("Cutoff Percentile for Flood Definition") xscale(range(`=100-$droughtCutoffMax' `=100-$droughtCutoffMin'))
				xlabel(`=100-$droughtCutoffMax'(5)`=100-$droughtCutoffMin+1')
				ylabel(, format(%6.3f))
				scheme(s1mono)
				graphregion(color(white)) bgcolor(white)
				legend(off)
				yline(0, lcolor(black) lpattern(dash))
				name("rob_${outputStub}_flood_cutoff", replace)
			;
			gr export "Output/rob_${outputStub}_flood_cutoff.pdf", replace name("rob_${outputStub}_flood_cutoff");
		#delimit cr
	}
	
	mat drop matBetas matSE
	cap mat drop mat_drought_cutoff 
	cap mat drop mat_flood_cutoff
restore

*********************************  India  **************************************

use "Data/India_duration_final" ///
	if inrange(year,${ageLowerLim},${ageUpperLim}), clear

di " "
di "*****************************************"
di "* Robustness to cutoffs: INDIA"
di "*****************************************"
di " "

* Initialize beta and SE matrices for the new definitions of drought_cy and flood_cy.
matrix matBetas = J(300,2,1)
matrix matSE = J(300,2,1)
matrix mat_drought_cutoff = J(300,1,1)
matrix mat_flood_cutoff = J(300,1,1)
global droughtCutoffMin 1
global droughtCutoffMax 45
global droughtCutoffNum = $droughtCutoffMax - $droughtCutoffMin + 1

**For coefplot**
matrix coef_IN = J(1, $droughtCutoffNum , .)
matrix ci_IN = J(2, $droughtCutoffNum , .)

* Loop over new cutoffs, run regressions, and save coefs + SEs
local counter = 0
local cutoffRng ${droughtCutoffMin}(1)${droughtCutoffMax}
forvalues cutoff = `cutoffRng' {
	local ++counter
	di "------------------------------------------"
	di "The cutoff is " %-10.2f `cutoff'/100
	di "------------------------------------------"
	di " "
	
	cap drop new_drought
	gen new_drought = (total_cy_xtile <= round(`cutoff'/100, 10^-2))
	
	if $includeFlood == 1 {
		cap drop new_flood
		gen new_flood = (total_cy_xtile > 1 - round(`cutoff'/100, 10^-2))
	
		* Main regression specification: drought + flood
		areg married_survival new_drought new_flood i.birth_year i.year ///
			if year<=$ageUpperLim & data=="DHS_98" ${INregWeight}, ///
			absorb(district) vce(cluster district)
	
		matrix matBetas[`counter',2] = _b[new_flood]
		matrix matSE[`counter',2] = _se[new_flood]
		matrix mat_flood_cutoff[`counter',1] = 1 - round(`cutoff'/100, 10^-2)
	}
	else {
		* Main regression specification: drought
		areg married_survival new_drought i.birth_year i.year ///
			if year<=$ageUpperLim & data=="DHS_98" ${INregWeight}, ///
			absorb(district) vce(cluster district)
	}
	
	* Drought is always included.
	matrix matBetas[`counter',1] = _b[new_drought]
	matrix matSE[`counter',1] = _se[new_drought]
	matrix mat_drought_cutoff[`counter',1] = round(`cutoff'/100, 10^-2)
	matrix coef_IN[1,`counter'] = _b[new_drought]
	matrix ci_IN[1,`counter'] = (_se[new_drought] * (-1.96)) + _b[new_drought]
	matrix ci_IN[2,`counter'] = (_se[new_drought] * 1.96) + _b[new_drought]
}

cap drop new_drought 
cap drop new_flood

preserve
	drop _all	
	svmat matBetas
	svmat matSE
	svmat mat_drought_cutoff
	svmat mat_flood_cutoff
	
	drop if matBetas1==1

	***** Drought
	* Rename vars, convert to %, then create 95% upper and lower CI
	rename matBetas1 drought
	rename matSE1 drought_se
	rename mat_drought_cutoff1 drought_cutoff
	replace drought_cutoff = drought_cutoff*100	
	gen drought_CIL = drought - 1.96*drought_se
	gen drought_CIU = drought + 1.96*drought_se
	
	* Make plots of coef (and 95% CI) vs. cutoff.
	#delimit ;
		twoway				
			(rarea drought_CIL drought_CIU drought_cutoff, color(gs13))
			(connected drought drought_cutoff, mcolor(black) lcolor(black) lpattern(solid) msymbol(S) msize(0.8))
			,
			title("Robustness to Rainfall Cutoff: India")
			ytitle("Effect of Drought") xtitle("Cutoff Percentile for Drought Definition") xscale(range($droughtCutoffMin $droughtCutoffMax))
			xlabel(0(5)$droughtCutoffMax)
			ylabel(, format(%6.3f))
			scheme(s1mono)
			graphregion(color(white)) bgcolor(white)
			legend(off)
			yline(0, lcolor(black) lpattern(dash))
			name("rob_India_drought_cutoff", replace)
		;
		gr export "Output/rob_India_drought_cutoff.pdf", replace name("rob_India_drought_cutoff");
	#delimit cr
	
	if $includeFlood == 1 {
		***** Flood
		* Rename vars, convert to %, then create 95% upper and lower CI
		rename matBetas2 flood
		rename matSE2 flood_se
		rename mat_flood_cutoff1 flood_cutoff
		replace flood_cutoff = flood_cutoff*100
		gen flood_CIL = flood - 1.96*flood_se
		gen flood_CIU = flood + 1.96*flood_se
		
		* Make plots of coef (and 95% CI) vs. cutoff.
		#delimit ;
			twoway				
				(rarea flood_CIL flood_CIU flood_cutoff, color(gs13))
				(connected flood flood_cutoff, mcolor(black) lcolor(black) lpattern(solid) msymbol(S) msize(0.8))
				,
				title("Robustness to Rainfall Shocks Cutoff: India")
				ytitle("Effect of Flood") xtitle("Cutoff Percentile for Flood Definition") xscale(range(`=100-$droughtCutoffMax' `=100-$droughtCutoffMin'))
				xlabel(`=100-$droughtCutoffMax'(5)`=100-$droughtCutoffMin+1')
				ylabel(, format(%6.3f))
				scheme(s1mono)
				graphregion(color(white)) bgcolor(white)
				legend(off)
				yline(0, lcolor(black) lpattern(dash))
				name("rob_India_flood_cutoff", replace)
			;
			gr export "Output/rob_India_flood_cutoff.pdf", replace name("rob_India_flood_cutoff");
		#delimit cr
		}
		
		mat drop matBetas matSE
		cap mat drop mat_drought_cutoff 
		cap mat drop mat_flood_cutoff
restore

**Combined coefplot**

#delimit ;
	coefplot (matrix(coef_SSA), recast(connected) ci(ci_SSA) label("SSA") color(blue) ciopts(recast(rcap) color(blue)))
		(matrix(coef_IN), recast(connected) lp(solid) m(circle) ci(ci_IN) label("India") color(red) ciopts(recast(rcap) color(red))),
		vertical title()
		ytitle(Effect of Drought) xtitle(Cutoff Percentile for Drought Definition)
		yline(0, lcolor(black)) xlabel(0(5)45) nooffset
	;
#delimit cr

graph export "Output/drought_robustness.pdf", replace

#delimit ;
	coefplot (matrix(coef_SSA), recast(connected) ci(ci_SSA) label("SSA") color(blue) ciopts(recast(rcap) color(blue)))
		(matrix(coef_IN), recast(connected) lp(solid) m(circle) ci(ci_IN) label("India") color(red%0) ciopts(recast(rcap) color(red%0))),
		vertical title()
		ytitle(Effect of Drought) xtitle(Cutoff Percentile for Drought Definition)
		yline(0, lcolor(black)) xlabel(0(5)45) nooffset legend(order(2))
	;
#delimit cr

graph export "Output/drought_robustness_SSA.pdf", replace

#delimit ;
	coefplot (matrix(coef_SSA), recast(connected) ci(ci_SSA) label("SSA") color(blue%0) ciopts(recast(rcap) color(blue%0)))
		(matrix(coef_IN), recast(connected) lp(solid) m(circle) ci(ci_IN) label("India") color(red) ciopts(recast(rcap) color(red))),
		vertical title()
		ytitle(Effect of Drought) xtitle(Cutoff Percentile for Drought Definition)
		yline(0, lcolor(black)) xlabel(0(5)45) nooffset legend(order(4))
	;
#delimit cr

graph export "Output/drought_robustness_IN.pdf", replace

******************************************************************************************
window manage close graph _all
cap log close

******************************************************************************************


********************************************************************************
* Table A5: 干旱的不同定义：总体干旱水平
********************************************************************************
*-----------------------------------------------   SSA    ---------------------------------
clear
************ IMPORTANT: Remember to choose age cutoff ************
global ageLowerLim 12
* Whether to limit age up till 17 or 24
global ageUpperLim 17
******** Globals for the effect of drought by age plots for both SSA and India
global shock_controls_1		drought_cy_${ageLowerLim} - drought_cy_${ageUpperLim} i.birth_year i.year
global step 	1
global xvar 	"Age"
use "Data/SSA_duration_final" ///
	if dhsyear >= 1994 & inrange(year,${ageLowerLim},${ageUpperLim}), clear
	*** Drop war years
keep if war_upto${ageUpperLim}==0

*GENERATE NEW VARIABLES 
*Generate # grids in the country which experienced drought in the same year
**By country**
gen year_cy=calendar_year
bys  year_cy grid_cell: gen countgrid=_n
sort  year_cy  grid

gen aaa=(countgrid==1)
gen bbb=(countgrid==1 & drought_cy==1)

bys country_code year_cy: egen numerator=sum(bbb) if countgrid==1
bys country_code year_cy: egen denominator=sum(aaa) if countgrid==1

gen ndrought_neigh_cty_aux = (numerator - drought_cy) / (denominator - 1) if drought_cy != .
gen ndrought_inclown_cty_aux = numerator/denominator if drought_cy != .

bys grid_cell year_cy: egen ndrought_neigh_cty = max(ndrought_neigh_cty_aux)
bys grid_cell year_cy: egen ndrought_inclown_cty = max(ndrought_inclown_cty_aux)

drop aaa bbb *aux numerator denominator countgrid

gen drought_ndroughtneigh_cty = drought_cy * ndrought_neigh_cty
gen drought_ndroughtinclown_cty = drought_cy * ndrought_inclown_cty

lab var ndrought_neigh_cty "Ratio of in-country droughts (excl. own)"
lab var ndrought_inclown_cty "Ratio of in-division droughts"
lab var drought_ndroughtneigh_cty "Local drought $\times$ ratio of in-country droughts (excl. own)"
lab var drought_ndroughtinclown_cty "Local drought $\times ratio of in-country droughts"
gen ndrought_neigh_clone = ndrought_neigh
gen ndrought_neigh_cty_clone = ndrought_neigh_cty

*REGRESSIONS
// model names: [modelID]_[sampleID]_{Birth Year FE}{Age FE}{Country FE}{Country FE X Cohort FE}, where {} is either y (yes), n (no), or t (not applicable).
	
**Aggregating at the country level**
areg married_survival ndrought_inclown_cty i.birth_year i.year ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster country_code)
areg married_survival ndrought_inclown_cty i.birth_year i.year i.country ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster country_code)
areg married_survival ndrought_inclown_cty i.birth_year i.year  i.country#i.cohort_num ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster country_code)

*------------------------------------    INDIA  ----------------------------------------------------------*
clear
snapshot erase _all
set matsize 5000, perm
************ IMPORTANT: Remember to choose age cutoffs ************
global ageLowerLim 12
* Whether to limit age up till 17 or 24
global ageUpperLim 17

*** Whether to use weights in India regressions
global INregWeight 
*global INregWeight [pweight = wgt_new]

******** Globals for the effect of drought by age plots for both SSA and India
global shock_controls_1		drought_cy_${ageLowerLim} - drought_cy_${ageUpperLim} i.birth_year i.year
global step 	1
global xvar 	"Age"

************** INDIA ************
use "Data/India_duration_final" ///
	if inrange(year,${ageLowerLim},${ageUpperLim}) & data=="DHS_98", clear
tab data, m // DHS_98

**Generating variables**
gen year_cy=calendar_year
bys year_cy district: gen countgrid=_n

gen distarea = total_area if countgrid==1
gen droughtarea = total_area * drought_cy

bys state_gec_num year_cy: egen a_numerator=sum(droughtarea) if countgrid==1
bys state_gec_num year_cy: egen a_denominator=sum(distarea) if countgrid==1

gen ndrought_neigh_aux = (a_numerator - droughtarea) / (a_denominator - distarea) // Leaving out own district
gen ndrought_inclown_aux = a_numerator / a_denominator // Incl. own district

bys district year_cy: egen ndrought_neigh = max(ndrought_neigh_aux)
bys district year_cy: egen ndrought_inclown = max(ndrought_inclown_aux)

drop distarea droughtarea *aux

gen drought_ndrought_neigh = drought_cy * ndrought_neigh
gen drought_ndrought_inclown = drought_cy * ndrought_inclown

lab var ndrought_neigh "Ratio of neighboring droughts"
lab var ndrought_inclown "Ratio of local droughts"
lab var drought_ndrought_neigh "Local drought $\times$ ratio of neighboring droughts"
lab var drought_ndrought_inclown "Local drought $\times$ ratio of local droughts"

gen ndrought_inclown_cty = ndrought_neigh
gen ndrought_neigh_cty = ndrought_neigh

lab var ndrought_inclown_cty "Ratio of local droughts"
lab var ndrought_neigh_cty "Ratio of neighboring droughts"

gen ndrought_neigh_clone = ndrought_neigh

**Regressions**
areg married_survival ndrought_inclown_cty i.birth_year i.year ///
	if year<=$ageUpperLim  ${INregWeight}, ///
	absorb(district) vce(cluster state_gec_num)   
	
areg married_survival ndrought_inclown_cty i.birth_year i.year i.state_gec_num#i.cohort_num ///
	if year<=$ageUpperLim ${INregWeight}, ///
	absorb(district) vce(cluster state_gec_num)  
	
	
********************************************************************************
* Table B9: 洪水（高降水概率对童婚的影响）
********************************************************************************
**SSA**
use "Data/SSA_duration_final" ///
	if inrange(year,12,17), clear
keep if war_upto17==0
areg married_survival drought_cy flood_cy i.birth_year i.year if year<=17 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy flood_cy using "Output/tableB9.xls", ///
	keep(drought_cy flood_cy) replace excel label adjr2 nocons
areg married_survival drought_cy flood_cy i.birth_year i.year i.country if year<=17 ///
	[pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy flood_cy using "Output/tableB9.xls", ///
	keep(drought_cy flood_cy) append excel label adjr2 nocons

	
**India**
use "Data/India_duration_final" if inrange(year,12,17), clear
gen flood_cy_rice=flood_cy*meanshare_RICE
gen drought_cy_rice=drought_cy*meanshare_RICE
areg married_survival drought_cy flood_cy i.birth_year i.year ///
	if year<=17, ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy flood_cy using "Output/tableB9.xls", ///
	keep(drought_cy flood_cy) append excel label adjr2 nocons
areg married_survival drought_cy flood_cy drought_cy_rice flood_cy_rice i.birth_year i.year  ///
	if year<=17, ///
	absorb(district) vce(cluster district)   
outreg2 drought_cy flood_cy using "Output/tableB9.xls", ///
	keep(drought_cy flood_cy drought_cy_rice flood_cy_rice) append excel label adjr2 nocons
	
	

********************************************************************************	
* Table B11: 过去和未来的时间效应
********************************************************************************
************ IMPORTANT: Remember to choose age cutoff ************
global ageLowerLim 12
* Whether to limit age up till 17 or 24
global ageUpperLim 17

**SSA**
use "Data/SSA_duration_final" ///
	if inrange(year,12,17), clear

keep if war_upto17==0

************ Drought Shocks Lags ************

*col 1
areg married_survival drought_cy drought_cy_lag1 i.birth_year i.year i.country ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
*col 2
areg married_survival drought_cy drought_cy_lag1 drought_cy_lead1  ///
	i.birth_year i.year i.country if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
*col 3
	areg married_survival max_drought_cy i.birth_year i.year i.country ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
*col 7
areg married_survival drought_cy_lag1 i.birth_year i.year i.country ///
	if year<=$ageUpperLim [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
	
	
**India**

use "Data/India_duration_final" if inrange(year,12,17), clear

* col 4
areg married_survival drought_cy drought_cy_lag1 ///
	i.birth_year i.year  ///
	if  year<=17, absorb(district) vce(cluster district)
* col 5
areg married_survival drought_cy drought_cy_lag1 drought_cy_lead1 ///
	i.birth_year i.year if year<=17, ///
	absorb(district) vce(cluster district)
*col 6
areg married_survival max_drought_cy i.birth_year i.year ///
	if year<=17, ///
	absorb(district) vce(cluster district)
* col 8
areg married_survival drought_cy_lead1 i.birth_year i.year  ///
	if year<=17, ///
	absorb(district) vce(cluster district)



********************************************************************************
* Table B12: 以更大的地理单元格单元格聚类
********************************************************************************

**SSA**
use "Data/SSA_duration_final" ///
	if inrange(year,12,17), clear
	areg married_survival drought_cy i.birth_year i.year [pweight=wgt_new], absorb(grid_cell) cluster(grid_cell) 
areg married_survival drought_cy i.birth_year i.year i.country [pweight=wgt_new], absorb(grid_cell) cluster(grid_cell) 
areg married_survival drought_cy i.birth_year i.year i.country#i.cohort_num [pweight=wgt_new], absorb(grid_cell) cluster(grid_cell) 
areg married_survival drought_cy i.birth_year i.year  [pweight=wgt_new], absorb(grid_cell) cluster(country) 
areg married_survival drought_cy i.birth_year i.year i.country  [pweight=wgt_new], absorb(grid_cell) cluster(country) 
areg married_survival drought_cy i.birth_year i.year i.country#i.cohort_num  [pweight=wgt_new], absorb(grid_cell) cluster(country) 

**India**

use "Data/India_duration_final" if inrange(year,12,17), clear
tab birth_year,gen(byd)
tab year,gen(yrd)
tab district,gen(distd)
xi: reg         married_survival drought_cy byd* yrd* distd*, cluster(district) 
xi: reg         married_survival drought_cy byd* yrd* distd*, cluster(state_gec_num) 
cap tab birth_year,gen(byd)
cap tab year,gen(yrd)
cap tab district,gen(distd)
egen aaa=group(cohort_num state_gec_num)
clustse regress married_survival drought_cy byd* yrd*, fe(district) cluster(state_gec_num) reps(500) seed(999) method(wild)
tab aaa, gen(aaad)
xi: reg         married_survival drought_cy byd* yrd* distd* aaad*, cluster(district) 
xi: reg         married_survival drought_cy byd* yrd* distd* aaad*, cluster(state_gec_num) 
clustse regress  married_survival drought_cy byd* yrd*  aaad*, fe(district) cluster(state_gec_num) reps(500) seed(999) method(wild)



********************************************************************************
*
*
* 四、异质性分析
*
*
********************************************************************************

********************************************************************************
* Table4: 撒哈拉以南地区异质性分析
********************************************************************************
**SSA**

use "Data/SSA_duration_final" ///
	if inrange(year,12,17), clear
keep if war_upto17==0
lab var drought_cy "Drought"

**High BP**
areg married_survival drought_cy i.birth_year i.year if year<=17 ///
	& highbp_ge50pct==1 [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table4.xls", ///
	keep(drought_cy) replace excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country if year<=17 ///
	& highbp_ge50pct==1 [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table4.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**Low BP**
areg married_survival drought_cy i.birth_year i.year if year<=17 ///
	& lowbp_lt50pct==1 [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table4.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
areg married_survival drought_cy i.birth_year i.year i.country if year<=17 ///
	& lowbp_lt50pct==1 [pweight = wgt_new], absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table4.xls", ///
	keep(drought_cy) append excel label adjr2 nocons
	
**Dowry (Eritrea)**
areg married_survival drought_cy i.birth_year i.year if year<=17 ///
	& country_code=="ERI", absorb(grid_cell) vce(cluster grid_cell)
outreg2 drought_cy using "Output/table4.xls", ///
	keep(drought_cy) append excel label adjr2 nocons




********************************************************************************
* Table5: 印度地区异质性分析
********************************************************************************
use "Data/India_duration_final" if inrange(year,12,17) & data=="DHS_98", clear

gen nohindu=(religion!=1) if (religion!=.)
gen low_dowry=dowry_pop<0.5 if dowry_pop!=. 
gen low_dowry_nohindu=low_dowry*nohindu

gen drought_nohindu=drought_cy*nohindu 
gen drought_low_dowry=drought_cy*low_dowry
gen drought_low_dowry_nohindu=drought_cy*low_dowry*nohindu

lab var drought_cy "Drought"
la var drought_nohindu "Drought $\times$ non-Hindu"
la var drought_low_dowry "Drought $\times$ Low Dowry"
la var drought_low_dowry_nohindu "Drought $\times$ non-Hindu $\times$ Low Dowry"

areg married_survival drought_cy drought_nohindu i.birth_year#i.nohindu i.year#i.nohindu ///
	if year<=17, absorb(district) vce(cluster district)
outreg2 drought_cy drought_nohindu using "Output/table5.xls", ///
	keep(drought_cy drought_nohindu) replace excel label adjr2 nocons
lincom drought_cy + drought_nohindu
areg married_survival drought_cy drought_low_dowry i.birth_year#i.low_dowry i.year#i.low_dowry ///
	if year<=17, absorb(district) vce(cluster district)   
outreg2 drought_cy drought_low_dowry using "Output/table5.xls", ///
	keep(drought_cy drought_low_dowry) append excel label adjr2 nocons
lincom drought_cy + drought_low_dowry
areg married_survival drought_cy drought_low_dowry_nohindu i.birth_year#i.low_dowry_nohindu ///
	i.year#i.low_dowry_nohindu  if year<=17, absorb(district) vce(cluster district)
outreg2 drought_cy drought_low_dowry_nohindu using "Output/table5.xls", ///
	keep(drought_cy drought_low_dowry_nohindu) append excel label adjr2 nocons
lincom drought_cy + drought_low_dowry_nohindu

rename cohort_50t59 cohort1
rename cohort_60t69 cohort2
rename cohort_70t79 cohort3

gen south_cohort1=south_i*cohort1
gen south_cohort2=south_i*cohort2
gen south_cohort3=south_i*cohort3

gen drought_cy1=drought_cy*cohort1
gen drought_cy2=drought_cy*cohort2
gen drought_cy3=drought_cy*cohort3

gen drought_south=drought_cy*south_india
gen drought_south1=drought_cy1*south_india
gen drought_south2=drought_cy2*south_india
gen drought_south3=drought_cy3*south_india

la var drought_south "Drought $\times$ South India"
la var drought_cy1 "Drought $\times$ born 1950-59"
la var drought_cy2 "Drought $\times$ born 1960-69"
la var drought_cy3 "Drought $\times$ born 1970-79"
la var drought_south1 "Drought $\times$ South $\times$ born 1950-59"
la var drought_south2 "Drought $\times$ South $\times$ born 1960-69"
la var drought_south3 "Drought $\times$ South $\times$ born 1970-79"

areg married_survival drought_cy drought_south i.birth_year#i.south_i i.year#i.south_i ///
	if year<=17, absorb(district) vce(cluster district)
outreg2 drought_cy drought_south using "Output/table5.xls", ///
	keep(drought_cy drought_south) append excel label adjr2 nocons
	
lincom drought_cy + drought_south
areg married_survival drought_cy drought_cy2 drought_cy3 drought_south drought_south2 ///
	drought_south3 i.birth_year#i.cohort_num i.year#i.cohort_num i.birth_year#i.south_i ///
	i.year#i.south_i i.birth_year#i.south_cohort2 i.year#i.south_cohort2 i.birth_year#i.south_cohort3 /// 
	i.year#i.south_cohort3 if year<=17, absorb(district) vce(cluster district)
outreg2 drought_cy drought_cy2 drought_cy3 drought_south drought_south2 drought_south3 ///
	using "Output/table5.xls", ///
	keep(drought_cy drought_cy2 drought_cy3 drought_south drought_south2 drought_south3) ///
	append excel label adjr2 nocons
	
lincom drought_cy + drought_cy2
lincom drought_cy + drought_cy3
lincom drought_cy + drought_south
lincom drought_cy + drought_cy2 + drought_south + drought_south2
lincom drought_cy + drought_cy3 + drought_south + drought_south3