* do e:\DHS\programs\Bongaarts_proximate_determinants\proximate_determinants_do_27Sep2019.txt

/*

This is a modification of a program sent to Tom Pullum on Oct. 31, 2018; the
 original program was from Ivan Mejia Guevara, written for Jocelyn Finlay, HSPH.

It clearly comes from an earlier version, author unknown, but possibly German Rodriguez.

Modifications by Tom Pullum

This modification is intended for a single survey.

The 3-year age-specific fertility rates are not calculated within the program. 
It is assumed that you have downloaded an Excel file from STATcompiler and 
put it into the same folder as this program. That excel file must include the rates 
needed by this program

Alternatives would be to enter the rates manually, link to a fertility rates program, or link
to an API retrieval program.

To run the program, enter the name of the country, the date of the survey, and the name of the 
IR file (without "FL.dta"), as in this example:

scalar scountry="Malawi"
scalar ssurvey_year="2015-16"
scalar sIRfilename="MWIR7H"

In the Guevara version, i is obtained from the most recent birth and is not a current status measure.
That is a mistake.

*/

cd e:\DHS\programs\Bongaarts_proximate_determinants

log using proximate_determinants_log_27Sep2019, replace


clear all
clear matrix
set more off
set maxvar 32000

********************************************************************************
* TFR = Cm*Ci*Ca*Cp*TF

* 	Cm: Index of proportion marriage
* 	Ci: Index of lactational infecundability
*	Ca: Index of abortion
*	Cp: Index pf pathological sterility
*	Cc: Index of COntraception
*	TF: Total fecundity = 15
 
***************************************************************************

program define get_fertility_rates

import excel STATcompilerExport2018112_74855.xlsx

keep if regexm(A,scountry)==1 & regexm(B,ssurvey_year)==1 
drop A B C
rename D f_1
rename E f_2
rename F f_3
rename G f_4
rename H f_5
rename I f_6
rename J f_7
destring f*, replace

*egen TFR=rowtotal(f*)
*scalar sTFR=5*TFR[1]/1000
*drop TFR

gen dummy=1
reshape long f_, i(dummy) j(age5)
rename *_ *

replace f=f/1000
rename f fa
drop dummy

list, table clean
sort age5
save workf.dta, replace

end   

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

program define make_Cm

  
use IR.dta, clear

* Include women who are married (married=1) or pregnant or using contraception or had sex in past month


collapse (mean) married nonmarital_exposure [iweight=wt], by(age5)
gen Cma=married+nonmarital_exposure
list, table clean

sort age5
save workm.dta, replace

end   

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

program define make_Cc
  
use IR.dta, clear

* reduce to women who are married or living as married
keep if (married==1 | recentsex==1) & insus==0

  
* Stover (1998) & Bongaarts (1982)

/*
(m0) 0=none, 
(m1) 1=pill,
(m2) 2=IUD, 
(m3) 4,5,14=Barrier (Condom[F/M], diaphragm)
(m4) 8,9=Withdrawal or rhythm,
(m5) 11=Implants, 
(m6) 3=Injectable, 
(m7) 6,7=Sterilization,
(m8) 15=Foam or jelly, 
(m9) 10,16,18=Other - traditional,
(m10) 17=Other - modern or specific, 
(m11) 12=Abstinence, 
(m12) 13=Lactational amenorrhea (

m11 and m12 are excluded because included in Ci
*/


*  gen u_m0  = (v312 == 0)
  gen u_m1  = (v312 == 1)
  gen u_m2  = (v312 == 2)
  gen u_m3  = (inlist(v312,4,5,14))
  gen u_m4  = (inlist(v312,8,9))
  gen u_m5  = (v312 == 11)
  gen u_m6  = (v312 == 3)
  gen u_m7  = (inlist(v312,6,7))
  gen u_m8  = (v312 == 15)
  gen u_m9  = (inlist(v312,10,16,18))
  gen u_m10 = (inlist(v312,17))

  * Overlap with postpartum infecundability
*  gen u_m11 = (v312 == 12)
*  gen u_m12 = (v312 == 13)
 
recode u_m* (0 1 = .) if missing(v312) | inlist(v312,16,98)

collapse (mean) u_m* [iweight=wt], by(age5)
egen u=rowtotal(u_m*)

list, table clean

* u*e is equivalent to the inner product of the u's and e's 

* Enter use-effectiveness of the methods as scalars 
 
  * Stover(1998)
scalar e_m0 = 0
scalar e_m1 = 0.92
scalar e_m2 = 0.96
scalar e_m3 = 0.81
scalar e_m4 = 0.50
scalar e_m5 = 1
scalar e_m6 = 1
scalar e_m7 = 1
scalar e_m8 = 0.91
scalar e_m9 = 0.10
scalar e_m10 = 0.91
*scalar e_m11 = 1
*scalar e_m12 = 1

* the mean use-effectiveness as the mean of e weighted by u 
* when this is multiplied by u, we just have the inner product of the u's and e's

* The following numbers are age-specific estimates of the fecundity adjustment from Bongaarts (2015) p. 559
scalar sr1=0.62
scalar sr2=0.81
scalar sr3=0.99
scalar sr4=1.08
scalar sr5=1.14
scalar sr6=1.26
scalar sr7=1.62

gen r=1
local lages 1 2 3 4 5 6 7
foreach la of local lages {
replace r = sr`la' if age5 == `la'
}

* alternative to take account of overlap with infecundability is just to sum from 1 to 10
gen ue=0
local li=1
quietly while `li'<=10 {
replace ue=ue+u_m`li'*e_m`li'
local li=`li'+1
}

gen rue=r*ue
gen Cca = 1 - rue 
replace Cca = 0.1 if Cca <= 0.1

list age5 ue r rue Cca, table clean
drop u* r*
  
list, table clean 
sort age5

save workc.dta, replace

end

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

program define make_Ci

use IR.dta, clear

* calculate i, the mean duration of insusceptibility, in months

gen n=1

keep if interval<36
gen newint=int(interval/2)

collapse (sum) insus n [iweight=wt], by(newint age5)
gen width=2
replace width=1 if newint==0

gen p_insus=insus/n

* pr for p revised
gen pr_insus=width*p_insus

rename pr_insus i

collapse (sum) i, by(age5)

*mean i_2

gen Cia = 20 / (18.5 + i)  

* b is a potential adjustment for the abortion index
*gen b = 14 / (18.5 + i_2)

list, table clean

sort age5
save worki.dta, replace    

end   

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

program define make_Ca

clear
*** Ca ***
/*
import excel "C:\Users\ivm761\Dropbox (Personal)\research\harvard\projects\jocelyn\africa\Fertility_Determ2015\data\Sedgh_etal_Lancet2012.xlsx", sheet("table2") firstrow clear
merge 1:m region using "${master_dhs}data\countries_un_dhs.DTA"
  keep if _merge == 3
  drop _merge
  keep total* region
merge 1:m using `lpath1'\tmp2_dhs.DTA
  keep if _merge == 3
  drop _merge
  gen myear1995 = abs(year_end - 1995)
  gen myear2003 = abs(year_end - 2003)
  gen myear2008 = abs(year_end - 2008)
  gen min_difyear = min(myear1995, myear2003, myear2008)
  
  gen TAR = total_1995 if abs(year_end - 1995) == min_difyear
  replace TAR = total_2003 if abs(year_end - 2003) == min_difyear
  replace TAR = total_2008 if abs(year_end - 2008) == min_difyear
  
  replace TAR = (TAR*30) / 1000  				
  gen TAR_a = TAR / 30 *5
  order total* min_difyear year_end TAR
  keep age5 TAR TAR_a
merge 1:1 age5 using `lpath1'\Cm_Q1.DTA
  keep if _merge == 3
  drop _merge
merge 1:1 age5 using `lpath1'\Ci_Q1.DTA
  keep if _merge == 3
  drop _merge
  
  *gen Ca_a = fa_a / (fa_a + b_a * TAR_a)
  gen Ca_a = TFR / (TFR + b_a * TAR)
*/

set obs 7
gen age5=_n

* Assumed value for Malawi 2010 survey
gen Caa=.92


sort age5
save worka.dta, replace

end

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

program define combine_and_save

*** Agregate Determinants ***

use workf.dta, clear
merge age5 using workm.dta
drop _merge
sort age5

merge age5 using workc.dta
drop _merge
sort age5

merge age5 using worki.dta
drop _merge
sort age5

merge age5 using worka.dta
drop _merge
sort age5

/*

* Calculate r if not using the adjusted fecundity numbers given by Bongaarts
gen fnca = f*1000 / (Cma * Cia * Caa)
gen uoe = u * e
gen r = .
gen mff = .

* within age groups, regress fnc on uoe; estimate what fnc would be if uoe=0

local lages 1 2 3 4 5 6 7
foreach la of local lages {
  reg fnc uoe if age5 == `la'
  replace r = abs(_b[uoe] / _b[_cons]) if age5 == `la'
}
*/

* calculate the implied natural fertility rate within each age group
* suffix a denotes age group (1 through 7)

* Insert the factor of 5 at this point
gen TFRa=5*fa
gen TFa = TFRa / (Cma*Cca*Cia*Caa)

* Calculate the age-specific versions of all the intermediate fertility rates 
* Note that T for total is not really appropriate at this stage, before adding

* The observed rate is fa and does not have to be re-calculated

*gen TFRa = Cma * Cca * Cia * Caa * TFa
gen  TMa = Cca * Cia * Caa * TFa
gen  TNa = Cia * Caa * TFa
gen  TAa = Ca * TFa

list age5 C* T*, table clean

* Now add across age groups
collapse (sum) T*

* Drop the a suffix because we now have totals
rename T*a T*

order TFR TM TN TA TF

* Construct the overal C's on the basis of the T's
gen   Cm = TFR / TM
gen   Cc = TM / TN
gen   Ci = TN / TA
gen   Ca = TA / TF


* Calculate the expected TFR if the TF were 15.3
gen TFRe = Cm * Cc * Ci * Ca * 15.3 

    

gen delta=TFR-TFRe  

list, table clean

save results_27Sep2019.dta, replace 

end

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

program define main

local lIRfilename=sIRfilename
use "C:\Users\26216\ICF\Analysis - Shared Resources\Data\DHSdata\\`lIRfilename'FL.DTA" , clear
rename v013 age5
gen wt=v005/1000000
gen married=0
replace married=1 if v501==1 | v501==2

gen pregnant=0
replace pregnant=1 if v213==1

gen recentsex=0
replace recentsex=1 if v536==1 | v536==2

gen user=0
replace user=1 if v313>0

gen nonmarital_exposure=0
replace nonmarital_exposure=1 if married==0 & (pregnant==1 | recentsex==1 | user==1)

rename v405 amen
rename v406 abst
gen insus=0
replace insus=1 if amen==1 | abst==1
gen interval=v008-b3_01

*keep wt age5 v312 married nonmarital_exposure m4_* m5_* m7_* m9_* insus newint
keep wt age5 v312 married pregnant recentsex user nonmarital_exposure insus interval 

* save with a generic name
save IR.dta, replace
clear

get_fertility_rates
make_Cm
make_Cc
make_Ci
make_Ca

combine_and_save

end

***************************************************************************
***************************************************************************
***************************************************************************
***************************************************************************
* EXECUTION BEGINS HERE

cd e:\DHS\programs\Bongaarts_proximate_determinants

scalar scountry="Malawi"
scalar ssurvey_year="2010"
scalar sIRfilename="MWIR61"

main