log using international03, replace

* Program INTERNATIONAL03.DO
*
* Refines workhorse program INTERNATIONAL02.DO for publicatoin on RIO website.  
*
* Version log: INTERNATIONAL02.DO builds on INTERNATIONAL01.DO. Re-expresses 
* price as fraction of per-capita disease burden. INTERNATIONAL01.DO performs 
* calibrations for intermational HIV pharmaceutical demand paper. Assumes all 
* consumers are homogeneous in a country.

set more 1


*** Read in population data from World Bank "Population, total"
*** Downloaded May 10, 2017 from http://data.worldbank.org/indicator/SP.POP.TOTL

import excel using POP.xlsx, sheet("Data") firstrow
sort ccode
tempfile pop
save "`pop'" 


*** Read in per capita GPD data from World Bank "GDP per capital (current US$)"
*** Downloaded May 10, 2017 from http://data.worldbank.org/indicator/NY.GDP.PCAP.CD

clear
import excel using GDP.xlsx, sheet("Data") firstrow
sort ccode
merge 1:1 ccode using "`pop'"
drop _merge
sort ccode
save "`pop'", replace 


*** Read in HIV prevalence data from World Bank "Prevalence of HIV, total (% of population ages 15-49)"
*** Downloaded May 10, 2017 from http://data.worldbank.org/indicator/SH.DYN.AIDS.ZS

clear
import excel using HIV.xlsx, sheet("Data") firstrow
sort ccode
merge 1:1 ccode using "`pop'"


*** Data cleaning

keep cname ccode pop2003 gdppc2003 hiv2003 hiv2003r
rename pop2003 pop
rename gdppc2003 gdpc
rename hiv2003 hiv
rename hiv2003r hivr
replace hiv = hivr if hiv == .
replace hivr = hiv if hivr == .

drop if pop == .
drop if gdpc == .
drop if hiv == .

*** May not need these if hand-collected HIV database excludes these regional aggregations

drop if ccode == "ARB" | ccode == "CEB" | ccode == "CSS" | ccode == "EAP"
drop if ccode == "EAR" | ccode == "EAS" | ccode == "ECA" | ccode == "ECS"
drop if ccode == "EMU" | ccode == "EUU" | ccode == "FCS" | ccode == "HIC"
drop if ccode == "HPC" | ccode == "IBD" | ccode == "IBT" | ccode == "IDA"
drop if ccode == "IDB" | ccode == "IDX" | ccode == "LAC" | ccode == "LCN"
drop if ccode == "LDC" | ccode == "LIC" | ccode == "LMC" | ccode == "LMY"
drop if ccode == "LTE" | ccode == "MEA" | ccode == "MNA" | ccode == "NAC"
drop if ccode == "OED" | ccode == "PRE" | ccode == "PST" | ccode == "SSA"
drop if ccode == "SSF" | ccode == "SST" | ccode == "TEA" | ccode == "TEC" 
drop if ccode == "TLA" | ccode == "TMN" | ccode == "TSA" | ccode ==  "TSS" 
drop if ccode == "UMC" | ccode == "WLD" | ccode == "MIC"
drop if ccode == "OSS" | ccode == "PSS"

gen x = hiv/pop
gen xr = hivr/pop
gen y = gdpc
gen n = pop


* Table 2 descriptive statistics

egen maxn = max(n)
format n maxn %20.0fc
summ x y
summ n, format
list n if n == maxn

* Table 3 correlations

pwcorr x y n, sig


* Toggle commented out lines and re-run to have just disease risk or income 
* variation reported in Table 6
* Replacing x and y with constant needed to make TS** equal to baseline

egen world_gdp = total(n*y)
egen world_hiv = total(hiv)
egen world_surplus = total(n*x*y)
* replace y = world_surplus / world_hiv
* replace x = world_surplus / world_gdp



*** Parameters (vary to get graphs, table entries, or text statements)

gen epsilon = 1
gen subsidy_v = 0
gen subsidy_d = 0
gen lambda = .7


*** Vaccine 

gen p_v = x * y^epsilon
gsort -p_v
gen q_v = sum(n)
gen s_v = p_v * n
gen ps_v = (p_v + subsidy_v) * q_v
egen PS_v = max(ps_v)
egen P_v = max(p_v*(ps_v == PS_v)) 
egen p_v_max = max(p_v)
gen P_v_pmax = P_v / p_v_max
egen Q_v = max(q_v*(ps_v == PS_v))
egen q_v_max = max(q_v)
gen Q_v_pct = Q_v / q_v_max
egen TS_v = total(s_v)
gen B_v_pcap = TS_v/q_v_max
gen P_v_pmean = P_v / B_v_pcap
gen served_v = p_v >= P_v
gen rho_v = PS_v / TS_v
egen GCS_v = total(s_v*served_v)  
gen CS_v = GCS_v - (P_v * Q_v)
gen GOV_v = subsidy_v * Q_v
gen g_v = GOV_v / TS_v
gen gamma_v = CS_v / TS_v
egen SDWL_v_num = total(s_v*(1-served_v))
gen SDWL_v = SDWL_v_num / TS_v
gen GDWL_v = lambda * g_v
gen DWL_v = 1 - rho_v
gen mu_v = TS_v / (q_v_max * p_v_max)

gen USA = ccode == "USA"
egen s_USA = max(s_v * USA)
gen pivotal = s_v / s_USA
gsort -pivotal
list cname pivotal

summ P_v P_v_pmax P_v_pmean Q_v Q_v_pct gamma_v rho_v g_v SDWL_v GDWL_v DWL_v mu_v
export excel p_v q_v using demands.xlsx, sheet("Vaccines") firstrow(variables) sheetreplace  


*** Drug 

gen p_d = y^epsilon
* Toggle off for calibration with just disease-risk heterogeneity for nice graphs
gsort -p_d
gen q_d = sum(n*x)
gen s_d = p_d * n * x
gen ps_d = (p_d + subsidy_d) * q_d
egen PS_d = max(ps_d)
egen P_d = max(p_d*(ps_d == PS_d)) 
egen p_d_max = max(p_d)
gen P_d_pmax = P_d / p_d_max
egen Q_d = max(q_d*(ps_d == PS_d))
egen q_d_max = max(q_d)
gen Q_d_pct = Q_d / q_d_max
egen TS_d = total(s_d)
gen B_d_pcap = TS_d/q_d_max
gen P_d_pmean = P_d / B_d_pcap
gen served_d = p_d >= P_d
gen rho_d = PS_d / TS_d
egen GCS_d = total(s_d*served_d) 
gen CS_d = GCS_d - (Q_d * P_d)
gen GOV_d = subsidy_d * Q_d
gen g_d = GOV_d / TS_d
gen gamma_d = CS_d / TS_d
egen SDWL_d_num = total(s_d*(1-served_d))
gen SDWL_d = SDWL_d_num / TS_d
gen GDWL_d = lambda * g_d
gen DWL_d = 1 - rho_d
gen mu_d = TS_d / (q_d_max * p_d_max)

summ P_d P_d_pmax P_d_pmean Q_d Q_d_pct q_d_max gamma_d rho_d g_d SDWL_d GDWL_d DWL_d mu_d
export excel p_d q_d using demands.xlsx, sheet("Drugs") firstrow(variables) sheetreplace


*** Both products

gen rho_min = min(rho_v,rho_d)
gen rho_max = max(rho_v,rho_d)
gen DWLC = rho_max - rho_min
gen DWL = 1 - rho_min
summ DWLC DWL


***************************************************************
*** Use revised disease-risk data instead (reported in Table 6)
***************************************************************

*** Vaccine 

drop p_v q_v s_v ps_v PS_v TS_v rho_v DWL_v
gen p_v = xr * y^epsilon
gsort -p_v
gen q_v = sum(pop)
gen s_v = p_v * pop
gen ps_v = p_v * q_v
egen PS_v = max(ps_v)
egen TS_v = total(s_v)
gen rho_v = PS_v / TS_v
gen DWL_v = 1 - rho_v
summ rho_v

*** Drug 

drop p_d q_d s_d ps_d PS_d TS_d rho_d DWL_d
gen p_d = y^epsilon
gsort -p_d
gen q_d = sum(pop*xr)
gen s_d = p_d * pop * xr
gen ps_d = p_d * q_d
egen PS_d = max(ps_d)
egen TS_d = total(s_d)
gen rho_d = PS_d / TS_d
gen DWL_d = 1 - rho_d
summ rho_d


*****************
*** US only payer
*****************

gen usa = ccode == "USA"

*** Vaccine 

drop p_v q_v s_v ps_v PS_v TS_v rho_v DWL_v
gen p_v = x * y^epsilon
gsort -p_v
gen q_v = sum(pop)
gen s_v = p_v * pop
gen ps_v = s_v * usa
egen PS_v = max(ps_v)
egen TS_v = total(s_v)
gen rho_v = PS_v / TS_v
gen DWL_v = 1 - rho_v
summ rho_v

*** Drug 

drop p_d q_d s_d ps_d PS_d TS_d rho_d DWL_d
gen p_d = y^epsilon
gsort -p_d
gen q_d = sum(pop*x)
gen s_d = p_d * pop * x
gen ps_d = s_d * usa
egen PS_d = max(ps_d)
egen TS_d = total(s_d)
gen rho_d = PS_d / TS_d
gen DWL_d = 1 - rho_d
summ rho_d


***************************************************************
*** Price discrimination with ceiling at US service requirement
***************************************************************

* Vary this parameter over range to get graph and other findings
gen upsilon = 2

*** Vaccine 

drop p_v q_v s_v ps_v PS_v TS_v rho_v DWL_v
gen p_v = x * y^epsilon
gen p_usa = p_v * usa * upsilon 
egen p_v_usa = max(p_usa)
drop p_usa
gen p_v_ceil = min(p_v,p_v_usa)
replace p_v_ceil = p_v if usa == 1
gen s_v = p_v * pop
gen ps_v = p_v_ceil * pop
egen PS_v = total(ps_v)
egen TS_v = total(s_v)
gen rho_v = PS_v / TS_v
gen DWL_v = 1 - rho_v
summ rho_v

*** Drug 

drop p_d q_d s_d ps_d PS_d TS_d rho_d DWL_d
gen p_d = y^epsilon
gen p_usa = p_d * usa * upsilon
egen p_d_usa = max(p_usa)
gen p_d_ceil = min(p_d,p_d_usa)
replace p_d_ceil = p_d if usa == 1
gen s_d = p_d * pop * x
gen ps_d = p_d_ceil * pop * x
egen PS_d = total(ps_d)
egen TS_d = total(s_d)
gen rho_d = PS_d / TS_d
gen DWL_d = 1 - rho_d
summ rho_d


*******************************************************
*** Price discrimination with ceiling at US dropping US
*******************************************************

*** Vaccine 

drop p_v s_v ps_v PS_v TS_v rho_v DWL_v
gen p_v = x * y^epsilon
gen s_v = p_v * pop
gen ps_v = s_v * (1 - usa)
egen PS_v = total(ps_v)
egen TS_v = total(s_v)
gen rho_v = PS_v / TS_v
gen DWL_v = 1 - rho_v
summ rho_v

*** Drug 

drop p_d s_d ps_d PS_d TS_d rho_d DWL_d
gen p_d = y^epsilon
gen s_d = p_d * pop * x
gen ps_d = s_d * (1 - usa)
egen PS_d = total(ps_d)
egen TS_d = total(s_d)
gen rho_d = PS_d / TS_d
gen DWL_d = 1 - rho_d
summ rho_d