program define aidsproj
*!   Program -aidsproj.ado- projects AIDS treatment burden
*!   Version 5.31   by Mead Over, Center for Global Development 11Nov2020

*   Purpose of this calling program is to permit the -if- condition and the
*   -using filename- specification to co-exist.
version 12
global apversion 5.31   // Starting in version 5.30, record the aidsproj version

*   Print header
	di _n(2) as txt "AIDSPROJ: Projections of AIDS Treatment Costs"
	di       as txt "              Version ${apversion}"
	di       as txt "{hline 45}"

	local noutstr = 1
	local nowgetfile = 0
	gettoken token 0 : 0 , bind parse(" ,")
	while `"`token'"'~="" {
		if `"`token'"'=="using" | `nowgetfile'==1 {
			if `nowgetfile'==0 {
				local nowgetfile = 1
			}
			else {
				local filename `token'

				di _n as txt "Specifications for input data file: "
				capture qui use `"`filename'"'  //  Force user to clear memory before specifying the using option
				if _rc==4 {
				    di _n as error "There is a problem loading the data file: " as result `"`filename'"'
				    di    as error "There is data in memory which the user must either {stata clear} or {help save}." _n
				    exit 603
				}
				else {
				    qui use `"`filename'"'
				    des , short
				    local nowgetfile = 0
				}
			}
		}
		else {
			local strlen1 : length local outstr`noutstr'
			local strlen2 : length local token
			if `strlen1'+`strlen2' > 240  {
				local noutstr = `noutstr' + 1
			}
			local outstr`noutstr' = `"`outstr`noutstr'' "' + `"`token'"'
		}
		gettoken token 0 : 0 , bind parse(" ,")
	}

	aidsproj_main `outstr1' `outstr2' `outstr3' `outstr4' `outstr5' `outstr6' `outstr7' `outstr8' `outstr9'

end

program define aidsproj_main

syntax [if]  [using/], UPTake(string) [TAkeoff(passthru) HORizon(passthru) adrate1(passthru) adrate2(passthru) BDrate(passthru) /*
			*/ Erate(passthru) NDrate(passthru) INCmult(passthru) COMplete IMPute WEIghts(string) /*
			*/ STRTCov2(passthru) TRGTCov2(passthru) STRTYr(passthru) TRGTYr(passthru) noRESCale /*
			*/ Countries(passthru) Regions(passthru) WOrld GR7 GRAph GRCombine LEGend(passthru) /*
			*/ TITLEN(passthru) TITLEC(passthru) TITLEP(passthru) TITLEE(passthru) SUBTitle(passthru) saving(string) /*
			*/ WEO  LOGOUT Vintage(real 2018) debug CD4(passthru) cost(string) DRate(passthru) REFYear(passthru) /*
			*/ postyears(passthru) posttotal postvars(passthru) postfile(passthru) postrunspecs(passthru) postopen postclose postname(passthru)  *]
*	A syntax command with all the options would contain more than the maximum number of characters
*	Adding new options requires another syntax command.

local 0 , `options'

syntax [, noSMOOTH rachet linesize(integer 100) putexcel fudge(passthru) outpath(string)  /*
            */ rxc1lb(passthru) rxc1ub(passthru) rxc2lb(passthru) rxc2ub(passthru) /*
            */ hsbedn(passthru) hsvstn(passthru) ndaf(passthru) scale(passthru)  /*
Endogenize i*/ iexog ir0(passthru) /*
            */ mincp(passthru) maxcp(passthru) cstrtyr(passthru) ctrgtyr(passthru) csp(passthru)  /*   Circumcision
            */ minep(passthru) maxep(passthru) pstrtyr(passthru) ptrgtyr(passthru)  /*   Prevention "effort"
			*/ maxvp(passthru) vstrtyr(passthru) vtrgtyr(passthru) ve(passthru) /* Vaccination
            */ fp(passthru) gp(passthru) maxdt(passthru) maxdu(passthru) ]

local old_linesize = c(linesize)
set linesize `linesize'

if "`fudge'"~="" & "`world'"=="" {
	di as err "The option -fudge(integer)- only applies to the aggregate projections invoked by the -world- option."
	di as err "Thus the -fudge(integer)- can only be specified with the -world- option."
	exit 198
}

if "`putexcel'" ~= "" {
	capture putexcel
	if _rc==199 {
		di as err "The putexcel option requires Stata version 13.0 or later"
		exit 198
	}
}

if "`debug'"=="" {
	local shhh quietly
}

if `vintage'<2007 | `vintage'>2018 {
    di _n as error "Vintage for the data must be in the range 2007 - 2018"
    exit 198
}
if `vintage'==2009 { /* Projection is the same for vintage 2009 as for vintage 2008 */
    local vintage = 2008
    local vdisplay = 2009
    local a2008 a2008
	local vwdi = `vintage'-2
}
if `vintage'==2007 | `vintage'==2008 {
    local vweo = `vintage'-2
    local vminus1  = `vintage'-1
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
if `vintage' == 2010 {
    local vweo = `vintage'-2
    local vminus1  = `vintage'-1
	local vwdi = `vintage'-3
    local vdisplay = `vintage'
}
if `vintage' == 2011 {				//  This code is a placeholder.  What offsets should we impose here?
    local vweo = `vintage'-2
    local vminus1  = `vintage'-1
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}

if `vintage' == 2012 {				//  This code is a placeholder.  What offsets should we impose here?
    local vweo = `vintage'-1
    local vminus1  = `vintage'-1
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}

if `vintage' == 2013 {				//  This code updates aidsproj to work with aidscgd2014, vintage 2013
    local vweo = `vintage'-2
    local vminus1  = `vintage'-1
    local vminus2  = `vintage'-2
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
if `vintage' == 2014 {				//  This code updates aidsproj to work with aidscgd2014, vintage 2014
    local vweo = `vintage'-2
    local vminus1  = `vintage'-1
    local vminus2  = `vintage'-2
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
if `vintage' == 2015 {				//  This code updates aidsproj to work with aidscgd2015, vintage 2015
    local vweo = `vintage'-1
    local vminus1  = `vintage'-1
    local vminus2  = `vintage'-2
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
if `vintage' == 2016 {				//  This code updates aidsproj to work with aidscgd2016, vintage 2016
    local vweo = `vintage'-1
    local vminus1  = `vintage'-1
    local vminus2  = `vintage'-2
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
if `vintage' == 2018 {				//  This code updates aidsproj to work with aidscgd2018, vintage 2018
    local vweo = `vintage'-1
    local vminus1  = `vintage'-1
    local vminus2  = `vintage'-2
	local vwdi = `vintage'-2
    local vdisplay = `vintage'
}
*	Confirm presence of requisite eight variables with the correct date suffix
local varyes = 0
foreach X in p a n i d rxc1 rxc2 nonrxc {
    capture confirm numeric variable `X'`vminus1'
	local varyes = (!_rc==1) + `varyes'
}
if `varyes' < 8 {
	di _n as error "Vintage of the data is " as result `vintage' as error ".  Therefore, the variables "
	di    as error "p`vminus1', a`vminus1', n`vminus1', i`vminus1', d`vminus1', rxc1`vminus1', rxc2`vminus1' and nonrxc`vminus1' are required"
	exit 198
}
else {
	di _n as txt "Vintage of the data is " as result `vdisplay'
}

*	Verify the data is unchanged from CGD original data
capture datasignature confirm
if _rc==0 {
	di as txt "Data is intact"
}
if _rc==459 {
	di as txt "no data signature was previously set"
}
if _rc~=0 & _rc~=459 {
	datasignature report
	di _n as txt "*** Warning! Data has been changed from CGD's original ***"
	di    as txt "Execution proceeding...."
	if "`debug'" ~="" {
		local oldcodes : char _dta[countrycodes]
		levelsof code
		local newcodes `r(levels)'
		local intersection : list newcodes & oldcodes
		local union : list newcodes | oldcodes
		local gone : list oldcodes - intersection
		local added : list newcodes - intersection

		if `"`gone'"'~="" {
			local first="yes"
			foreach i of numlist `gone' {
				if "`first'" == "yes" {
					di _n as txt "The following countries are missing from the CGD vintage `vintage' data"
					di as txt _n _col(7) "Numeric  Alphanumeric"
					di as txt    _col(7) "  Code      Code"
				}
				local ltrcode : label code `i'
				di as res _col(10) `i' _col(20) "`ltrcode'"
				local first="no"
			}
		}

		if `"`added'"'~="" {
			local first="yes"
			foreach i of numlist `added' {
				if "`first'" == "yes" {
					di _n as txt "The following countries are added to the CGD vintage `vintage' data"
					di as txt _n _col(7) "Numeric  Alphanumeric"
					di as txt    _col(7) "  Code      Code"
				}
				local ltrcode : label code `i'
				di as res _col(10) `i' _col(20) "`ltrcode'"
				local first="no"
			}
		}
	}
}

*	Confirm presence of expenditure & WEO data and save a file for later merging
*if "`weo'"~="" {     // In version 4.x we must retain these variables whether or not the WEO output is requested
	tempfile tmpfl
	confirm numeric variable hepctgdp`vwdi' pubhepctgdp`vwdi' pubhepcttot`vwdi' oophepctpriv`vwdi' exthepcttot`vwdi' gnipc`vweo'
	unab weovars : cggdp* gdpuskp* pop15_64_*
	confirm numeric variable `weovars'
	local rweovars cggdp gdpuskp pop15_64_
	`shhh'	di _n _n as txt "Reshape the WEO input data to panel data structure:"
	preserve
	qui keep code nation `weovars' rgn
	`shhh' reshape long `rweovars', i(code) j(year)
		lab var  cggdp     "Central government expenditures plus net lending as % GDP"
		lab var  pop15_64_ "Population age 15 - 64"
		lab var  gdpuskp   "GNI in 2009 US $, Billions"
		lab var year "Year"
	qui keep if year >= `vminus2'
	qui save `"`tmpfl'"', replace
	restore
	drop `weovars'
*}

*	SCALING UP OPTIONS: Constant or historical
*
*	First option: Global average uptake of ART as a constant proportion sigma = (1 - lambda) of unmet need
*	To use a specific constant value of sigma, set the macro `uptake' to that constant.
if real("`uptake'")<. {
	if abs(real("`uptake'")) < 1 & abs(real("`uptake'")) >= 0 {
		local sigmaC  = real("`uptake'")
		local lambdaC = 1 - real("`uptake'")
		local uptkmeth CONSTANT
		local rlblarg
	}
	else {
		di _n as error "Uptake option must specify either a constant or an existing numeric variable."
		di _n as error "When specified as a constant, the uptake option must have an absolute value >= 0 & < 1.0 "
		di    as error "Here the option has the value " as result `uptake'
		exit 198
	}
}

*	Second option: Country-specific uptake rates for scaling up ART
*	Use estimated uptake rates based on past performance
if real("`uptake'")>=.  {
	capture confirm numeric variable `uptake'
	if _rc==0 {
		local uptkmeth VARIABLE
*		local rlblarg ", sigma(`uptake') "   //  Version 4_59 moves this to line 299
		qui {
			gen byte uptcat = 1 if `uptake'<=0
			lab var uptcat "Description of uptake variable in input data"
			replace uptcat = 2 if `uptake' >0 & `uptake' <1
			replace uptcat = 3 if `uptake' >=1 & `uptake' <.
			replace uptcat = 4 if `uptake' >=.
		}
		label define uptcat 1 "`uptake'<=0" 2 "0<`uptake'<1" 3 "1<=`uptake'<." 4 "`uptake' missing"
		label value uptcat uptcat
		di _n as txt "Projection uses historical ART uptake data in the input variable: " as res "`uptake'"
		di    as txt "In the input data, the variable " as res "`uptake'" as txt " is distributed as follows."
		tab uptcat rgn
		di _n as txt "For projection purposes, missing values of " as res "`uptake'"
		di    as txt " are set to the regional average and non-missing values of " as res "`uptake'"
		di    as txt " outside the unit interval are set to .01 if <= 0 or to .99 if >= 1."
		tempvar uptmp
		`shhh' gen `uptmp' = `uptake'
		`shhh' replace `uptmp' = .99 if `uptmp' >=1 & `uptmp'<.
		`shhh' replace `uptmp' = .01 if `uptmp' <=0
		`shhh' levelsof rgn
		foreach rgn of numlist `r(levels)' {
			sum `uptmp' if rgn==`rgn', meanonly
			local rgnlbl : label rgn `rgn'
			if "`debug'"=="" {
				`shhh' replace `uptmp' = `r(mean)' if `uptmp'>=. & rgn==`rgn'
			}
			else {
				di _n as txt "For region: " as res "`rgnlbl'"
				repshow `uptmp' = `r(mean)' if `uptmp'>=. & rgn==`rgn', id(nation) l(a20*)
			}
			local sigmaV  "`uptmp'"
			local lambdaV  "(1 - `uptmp' )"
		}
		local rlblarg ", sigma(`sigmaV') "   //  Version 4_59
	}
	else {
		di _n as error "Uptake option must specify either a constant or an existing numeric variable."
		di    as error "Here the option specifes the variable name " as result "`uptake'"
		di    as error "which is not an existing numeric variable."
		exit 198
	}
}

*   Check validity of cost option
if `"`cost'"'~="" & `"`cost'"'~="revert" & `"`cost'"'~="function" {
    di _n as error "The parameter -cost- must be blank or contain one of the two strings: -revert- or -function-"
    exit 198
}

*	Print projection parameters
    di _n(2) as txt _col(15) "Parameters for projections of AIDS Treatment Costs"
    di       as txt _col(15) "=================================================="
    di _n    as txt _col(60)"Parameter"  _col(70) "Default"  _col(80) "Selected"
    di       as txt _col(60)"  Name   "  _col(70) " Value"   _col(80) " Value
    di       as txt _col(60)"---------"  _col(70) "--------" _col(80) "--------"

di as txt "Scale-up of 1st line modeled as constant proportion," _col(70) "   no "
di as txt "  sigma, of unmet need, where sigma is constant"   _col(70) "default"
if "`uptkmeth'"=="CONSTANT" {
	di as txt "  across all countries and equal to:" _col(60) "sigma" _col(70) "  here"  as res _col(80) `sigmaC'
}
else {
	di as txt "    within each country and given by the variable:" _col(70) "  here" as result _col(80)  as res "`uptake'"
}

parmsbycd4 , vintage(`vintage') `adrate1' `adrate2' `bdrate' `erate' `ndrate' `cd4' `takeoff' `horizon' `strtcov2' `trgtcov2' `strtyr' `trgtyr' `incmult' /*
Endogenize i*/ `iexog' `mincp' `maxcp' `csp' `cstrtyr' `ctrgtyr' `minep' `maxep' `pstrtyr' `ptrgtyr' `fp' `gp' `ir0' `maxdt' `maxdu' `vstrtyr' `vtrgtyr' `maxvp' `ve'

tempname adrate1 adrate2 bdrate erate ndrate takeoff horizon incmult strtcov2 trgtcov2 strtyr trgtyr  // Version 3.5 parameters
tempname mincp maxcp csp cstrtyr ctrgtyr minep maxep pstrtyr ptrgtyr fp gp ir0 maxdt maxdu vstrtyr vtrgtyr maxvp ve  // To endogenize incidence Version 4.0

*   Define formatted string macros to appear in graphs and _dta notes.
foreach pname in horizon strtyr trgtyr strtcov2 trgtcov2 cstrtyr ctrgtyr pstrtyr ptrgtyr vstrtyr vtrgtyr {
    scalar ``pname''  = `r(`pname')'
    local `pname'stg = strofreal(scalar(``pname''), "%6.0f")
}

foreach pname in mincp maxcp csp minep maxep fp gp ir0 maxdt maxdu maxvp ve {
    scalar ``pname''  = `r(`pname')'
    if abs(scalar(``pname''))<.0001 {
        local `pname'stg = strofreal(scalar(``pname''), "%3.0f")
    }
    else {
        local `pname'stg = strofreal(scalar(``pname''), "%5.2f")
    }
}
    local cd4strng = strofreal(scalar(`r(cd4)'), "%9.0f")
scalar `adrate1'  = `r(adrate1)'
    local adr1strng = strofreal(scalar(`adrate1'), "%9.3f")
scalar `adrate2'  = `r(adrate2)'
    local adr2strng  = strofreal(`adrate2', "%9.3f")
scalar `bdrate'   = `r(bdrate)'
    local bdrstrng   = strofreal(`bdrate', "%9.3f")
scalar `erate'    = `r(erate)'
    local estrng    = strofreal(`erate', "%9.3f")
scalar `ndrate'   = `r(ndrate)'
    local nstrng   = strofreal(`ndrate', "%9.3f")
local tkoff = `r(takeoff)'
	local takeoff takeoff(`tkoff')
scalar `incmult'  = `r(incmult)'
    local imstrng = strofreal(`incmult', "%9.3f")

if `"`cost'"'=="revert" {
    di as txt _n "Cost computations revert to 2008 assumptions" _n
}
else {
	*  The costmodel program is added in version 5.10, 14Feb2018
	costmodel , `rxc1lb' `rxc1ub' `rxc2lb' `rxc2ub' `hsbedn' `hsvstn' `scale' `ndaf'
	    foreach cname in rxc1lb rxc1ub rxc2lb rxc2ub hsbedn hsvstn ndaf scale {
			local `cname' `r(`cname')'
		}
}

if `"`postyears'`postvars'`postfile'`postrunspecs'`postopen'`postclose'`postname'"'~="" {
    tokenize `postfile'
    di as txt _n "Results of this run are posted to the file " as res "`3'"
    di as txt "  Discount rate for TOTALs is:"                                          as result _col(60) "`drate'"
    di as txt "  Option -postyears- is:"                                                as result _col(60) "`postyears'"
    di as txt "  Option -postvars- is:"                                                 as result _col(60) "`postvars'"
    di as txt "  Option -postrunspecs- is:"                                             as result _col(60) `"`postrunspecs'"'
    di as txt "  Option -posttotal- is:"                                                as result _col(60) "`posttotal'"
    di as txt "  Option -postname- is:"                                                 as result _col(60) "`postname'"
    di as txt "  Option -postopen- is:"                                                 as result _col(60) "`postopen'"
    di as txt "  Option -postclose- is:"                                                as result _col(60) "`postclose'"
}

if "`weights'"=="" {
    local weights "art"
}

*	Check the takeoff year & planning horizon
if `horizon'< `vintage' | `horizon' > 2100 {
	di _n as error "Projection horizon must be between `vintage' and 2100"
	exit 198
}
if `tkoff'< `vminus1' | `tkoff' > 2100 | `tkoff' >= `horizon' {
	di _n as error "Projection -takeoff- year must be between `vminus1' and 2100 and less than the projection horizon"
	exit 198
}
if `horizon'< `trgtyr' & "`debug'" ~= "" {
	di _n as error "Warning! Projection horizon " as result `horizon' as error " is less than -trgtyr- " as result `trgtyr'
	di    as error "Execution continues under this assumption."
}

*	Check the starting year and target year for 2nd line coverage
if `trgtyr'< `vintage' | `trgtyr' < `strtyr' {
	di _n as error "Parameter -strtyr-, the year 2nd line coverage starts to scale up, must precede"
	di    as error "parameter -trgtyr-, the target year for for 2nd line coverage to attain its target coverage."
	di    as error "Parameter -trgtyr- must also be greater than the vintage year of the data."
	di    as error "Here -strtyr- is: " as res `strtyr' as err " and -trgtyr- is: " as res `trgtyr' as err " and -vintage- is: " as res `vintage'
	exit 198
}

*	Check the legitimacy of country and/or region identifiers
if `"`countries'"' ~="" | `"`regions'"' ~="" {
	chkdta , `countries' `regions'
}
local output `r(countries)' `r(regions)' `world'

*	Convert to panel data
`shhh' di _n _n as txt "Reshape the input data to panel data structure:"
`shhh' keep code cntry ccode nation region rgn complete uptake `sigmaV' ///
	p`vminus1' a`vminus1' `a2008' n`vminus1' i`vminus1' d`vminus1' rxc1`vminus1'  ///
	rxc2`vminus1' nonrxc`vminus1' pepfar hepctgdp`vwdi' pubhepctgdp`vwdi'  ///
	pubhepcttot`vwdi' oophepctpriv`vwdi' exthepcttot`vwdi' gnipc`vweo' /* pop15_64_`vweo' */
`shhh' reshape long p a n i d rxc1 rxc2 nonrxc /* pop15_64_ */, i(code) j(year)
	qui keep if year >= `vminus2'
	lab var year "Year"
pause

*	Convert -if- conditions about variables like x`vminus1' into -if- conditions about the variable x, where x =, p, a, n, etc.
	local if = subinstr(`"`if'"',"`vminus1'","",99)

*	Select set of countries for projections
capture confirm numeric variable complete
if !_rc~=1 {
	di _n as error "Existing data must contain a variable called -complete- which designates imputed observations"
	exit 198
}

select a p n i d rxc1 rxc2 nonrxc `if', `shhh' `impute' `complete'

*	Relabel variables and disaggregate ART into stages/lines
lab var rxc1 "Estimated 1st line cost per patient"
lab var rxc2 "Estimated 2nd line cost per patient"
lab var nonrxc "Estimated non-ARV cost per patient"
lab var p "Number living with HIV/AIDS"
lab var n "Number adults needing ART"
lab var d "Number deaths"
lab var i "Number of new HIV infections"
lab var pepfar "PEPFAR Focus Country"
if "`uptkmeth'"=="VARIABLE" {
	lab var `sigmaV' "Country specific adjustment parameter"
}

qui gen  long   a2 = cond((11/18)*a >=1, round((11/18)*a), 0)
	lab var a2   "Number adults receiving 1st line ART, Years 2+"
	assert a2 >= 0
qui gen  long   b  = cond((1/18)*a >=1, round((1/18)*a), 0)
	lab var b   "Number adults receiving 2nd line ART"
	assert b >= 0
qui gen long a1 = cond((a-a2-b) >=1, round(a-a2-b), 1)
	lab var a1   "Number adults receiving 1st line ART, Year 1"
	assert a1 >= 0
	drop a
qui gen  long art =  a1 + a2 + b
	lab var art "Total patients on ART"
	assert art >= 0
/*
*	Version 5.13 code for un
qui gen long un = cond(n-a1-a2-b >=0, n-a1-a2-b, 0 )
*	Version 5.2 code follows
*/
`shhh' repshow n = p  if n>=., id(cntry) listvars(year p n )
qui gen long un = cond(n-a1-a2-b>=0 & n-a1-a2-b<., n-a1-a2-b, . )
	lab var un "Unmet need for ART"
qui gen long es = cond(n-a1-a2-b >=0, 0, -(n-a1-a2-b) )
	lab var es "Early ART starters, included in -art-"
	qui count if (p - n)<0
	if `r(N)'> 0 {
		`shhh' di _n as txt "In the following countries, our estimate of numbers infected (-p-) is less than" ///
				_n "the WHO's estimate of those needing ART (-n-)." ///
				_n "Therefore we adjust -p- upwards to the largest of the variables -n- and -art-."
		`shhh' repshow p = max(n , p, art) if (p - n)<0 , id(nation) l(n art)
	}
qui gen  long h =  p - n  // Early starters are included in -art-
	lab var h    "HIV+, not eligible for ART"  /* i.e. "Healthy" HIV+ */
	`shhh' list nation h i p un art n a1 a2 b es if h < 0
	assert h >= 0
qui recast  long n a1 a2 b es, force

capture confirm numeric imputed
if _rc==0 {
	local impishere imputed
}
order code year cntry ccode nation region rgn complete `impishere' uptake `sigmaV' p i d h n a1 a2 b

format %12.0fc h n h p un art d i a1 a2 b es

*	Expand the sample in order to project to future years
sum year, meanonly
    local maxyr = r(max)
qui expand `horizon' - `maxyr' + 1 if year == `maxyr'
sort code year
qui replace year = year[_n-1]+1 if code==code[_n-1]
qui tsset code year

foreach vname of varlist b p h n un i d a1 a2 art es {
	qui replace `vname'  = . if year>`vminus1'
}

`shhh' {

*	Generate assumed future coverage/scale-up of 2nd line ART
*		If -strtcov2- is zero,
*		set -cvrg2- to be the 2009 region average before and including the -strtyr-
*		with missing values thereafter.
*di as err "Scalar strtyr is: " scalar(`strtyr')
*di as err "Scalar strtcov2 is: " scalar(`strtcov2')

if `strtcov2' < 0 {  // From 2009 UNAIDS/WHO Progress Report on ART
	by code : gen     cvrg2  = 0.020  + ((`trgtcov2' - 0.020 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==1 // EAP
	by code : replace cvrg2  = 0.091  + ((`trgtcov2' - 0.091 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==2 // ECA
	by code : replace cvrg2  = 0.300  + ((`trgtcov2' - 0.300 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==3 // LAC
	by code : replace cvrg2  = 0.027  + ((`trgtcov2' - 0.027 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==4 // MNA  assumed same as SSA
	by code : replace cvrg2  = 0.002  + ((`trgtcov2' - 0.002 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==5 // SAS
	by code : replace cvrg2  = 0.027  + ((`trgtcov2' - 0.027 ) / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr' & rgn==6 // SSA
/*
    qui gen      cvrg2 = 0.020 if rgn==1 & year <= max(`strtyr',`vminus1')  // EAP
    qui replace  cvrg2 = 0.091 if rgn==2 & year <= max(`strtyr',`vminus1')  // ECA
    qui replace  cvrg2 = 0.300 if rgn==3 & year <= max(`strtyr',`vminus1')  // LAC
    qui replace  cvrg2 = 0.027 if rgn==4 & year <= max(`strtyr',`vminus1')  // MNA assumed same as SSA
    qui replace  cvrg2 = 0.002 if rgn==5 & year <= max(`strtyr',`vminus1')  // SAS
    qui replace  cvrg2 = 0.027 if rgn==6 & year <= max(`strtyr',`vminus1')  // SSA
*/

}

*	But if -strtcov2- is greater than zero,
*		set -cvrg2- to zero before and including the -strtyr-
*		and to -strtcov2- thereafter ... to potentially be modified below.
else {

	di as txt _n  ///
		"by code : replace cvrg2  = `strtcov2' +  ((`trgtcov2' - `strtcov2') / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr'"
	by code : gen cvrg2  = `strtcov2' + ((`trgtcov2' - `strtcov2') / ( `trgtyr' - `strtyr')) * ( year - `strtyr' )  if year >= `strtyr'

*	qui gen cvrg2 = cond(year < `strtyr', 0, `strtcov2' )  // Modified in version 5.25
}
/*
noisily{
	di as err "Tab of cvrg2 by region and year"
	tab year region if rgn>2, sum(cvrg2)
}
*/

	replace cvrg2 = `trgtcov2' if year >= `trgtyr'
		lab var cvrg2 "Assumed coverage of 2nd line ART"
		format cvrg2 %9.3f

// For Version 4.00 this code is moved here from down below and is no longer conditional on -weo-
	tempvar mrg3
	merge code year using `"`tmpfl'"', unique
	sort code year
	by code : egen `mrg3' = max(_merge==3)  /* 1 if country merges in any year */
	qui keep if `mrg3'==1                                   /* Keep all annual observations > `vminus1' for any country that merged */
	extweo gdpuskp, rate cnvrgyr(2040) cnsryr(2015)		/* Extrapolate data to end of horizon */
	extweo pop15_64_, rate cnvrgyr(2040)				/* Extrapolate data to end of horizon */
	drop _merge
	drop if year > `horizon'
	tempvar gdp2pop
	gen `gdp2pop' = gdpuskp/pop15_64_

	levelsof year if `gdp2pop' <. & year < `vintage' + 5, local(gdp2popyrs)
	levelsof year if gnipc`vweo' <. & year < `vintage' + 5, local(gnipcyrs)  //  Bug fix in 4.52: Completed the name of the gnipc variable
		local overlap : list gdp2popyrs & gnipcyrs
		local overlap = subinstr("`overlap'"," ",", ",.)
	reg gnipc`vweo' `gdp2pop' if year==min(`overlap')
	predict gnipc_hat, xb
		lab var gnipc_hat "Projected GNI per capita, 2009 US Dollars"
		char define gnipc_hat[scale] 1
	drop if year<`vminus1'  //  Version 4.51
}


`shhh' {
    if "`iexog'"~="" {                                  // Generate incidence independent of AIDS treatment
        replace i = round( `incmult' * L.i )                      if year>`vminus1'
    	    assert i >= 0
    }
    else {                                              // Model incidence endogenously (Version >= 4.0)
*	Generate assumed coverage/scale-up of male circumcision and effort
		local mcstrtyr = cond("${apversion}"~="", `cstrtyr',`pstrtyr')
		local mctrgtyr = cond("${apversion}"~="", `ctrgtyr',`ptrgtyr')
		local mcmin = cond("${apversion}"~="", `mincp',0)
        gen ccvrg =  `mcmin'  if year <= `mcstrtyr'
        by code : replace ccvrg  = (`maxcp' - `mcmin') * ( year - `mcstrtyr')  / ( `mctrgtyr' - `mcstrtyr') + `mcmin'   if year >= `mcstrtyr'
    	    replace ccvrg = `maxcp' if year >= `mctrgtyr'
    	    lab var ccvrg "Proportion of males circumcised"
    	    format ccvrg %9.3f

		local emin = cond("${apversion}"~="", `minep',0)
        gen effort =  `emin' if year <= `pstrtyr'
            by code : replace effort  = (`maxep' - `emin') * ( year - `pstrtyr')  / ( `ptrgtyr' - `pstrtyr') + `emin'   if year >= `pstrtyr'
    	    replace effort = `maxep' if year >= `ptrgtyr'
    	    lab var effort "HIV prevention -effort-"
    	    format effort %9.3f

        gen vr =  0 if year <= `vstrtyr'
            by code : replace vr  = `maxvp' * ( year - `vstrtyr')  / ( `vtrgtyr' - `vstrtyr')   if year >= `vstrtyr'
    	    replace vr = `maxvp' if year >= `vtrgtyr'
    	    lab var vr "Fraction of population vaccinated"
    	    format vr %9.3f

        gen  prevprop = .
            lab var prevprop "Prevalence as proporortion of pop 15 to 64"

		// Add to subsequent output the generated variables ccvrg, vr
		local ccvrg_vr ccvrg vr
    }
}

*   Generate uptake variable -lambda-
if "`uptkmeth'"=="VARIABLE" {
 	local lambda "`lambdaV'"
}
// When uptake is a constant parameter and -takeoff- is later than `vintage',
// adjust each country to the target uptake rate from its historical rate
// over three years after -takeoff-.  The -rachet- option forces instant adjustment
// of the historical uptake to the new one in the year -takeoff-.
// This calculation requires that the historical uptake be in the data
// as the variable named "uptake".
else {
    if `tkoff'< `vintage' {
    	local lambda = `lambdaC'
    }
    else {
        tempvar lmbtmp
        confirm numeric variable uptake
        qui gen `lmbtmp' = 1-uptake                     if year < `tkoff'
			qui replace `lmbtmp' = .01 if `lmbtmp' <=0 & year < `tkoff'
			qui replace `lmbtmp' = .99 if `lmbtmp' >=1 & `lmbtmp'<. & year < `tkoff'
			qui replace `lmbtmp' = .25 if `lmbtmp' >=.& year < `tkoff'
		if "`rachet'"=="" & "`smooth'"=="" {
			qui replace `lmbtmp' =  ///
				(1-uptake) + (( year-`tkoff') / ( 3 )) * (`lambdaC'-(1-uptake))   ///
														if `lmbtmp'>=.&year-`tkoff'<=3
		}
		qui replace `lmbtmp' = `lambdaC'                if `lmbtmp'>=.
    	local lambda `lmbtmp'
    }
}

`shhh' {
    gen cvrgpct = .
	    lab var cvrgpct "ART patients as percent of need"
	gen ir = `ir0'                                              if year==`vminus1'
        lab var ir "Incidence rate in the core group"          // Applies only to core group.  Noncore is immune.
	// Assume i = ir0·(k·pop15_64_ - p).  Then solve for k yielding:
    gen k = (i + `ir0'*p)/(`ir0'*pop15_64_)                    if year==`vminus1'
        by code: replace k = k[_n-1] if k >=.               // Incidence rate is on this fraction of pop only
        lab var k "Fraction of population in the core group"  // Constant within a country over time
    gen kmod = .
        lab var kmod "k modified by disinhibition of the untreated"
    gen prevprop0 = p / pop15_64_                               if year==`vminus1'
        by code: replace prevprop0 = prevprop0[_n-1] if prevprop0 >=.
        lab var prevprop0 "Adult prevalence rate as prop. of 15-64"
    gen artf0 = min(1,max(0,(a1 + a2 + b)/p ))    if year==`vminus1'  // Version 4.06
        lab var artf0 "Fraction of HIV+ on ART, Year 0"
    gen artf = .
        lab var artf "Fraction of HIV+ on ART"
    gen tau0 = 1  -  `fp' * `gp' * artf0                        if year==`vminus1'
        by code: replace tau0 = tau0[_n-1] if tau0 >=.
        lab var tau0 "Treatment attenuation of incidence rate, Year 0"
    gen tau = .
        lab var tau "Treatment attenuation of incidence rate"
	tempvar ndratev ndrate_dfltv eratev
	gen `ndrate_dfltv' =           ///   Set the default rate to be consistent with deaths in vminus1
		(d - (a1 * scalar(`adrate1') + a2 * scalar(`adrate2') * (1 - cvrg2) + b * scalar(`bdrate'))) ///
		 /     ((`lambda') * un)											if year==`vminus1'
		replace `ndrate_dfltv'  =   min(1, max(0, `ndrate_dfltv' ) ) 		if year==`vminus1'
		replace `ndrate_dfltv'  =   L.`ndrate_dfltv'						if year>`vminus1'

	gen `ndratev' = max(0,min(1,(year -`vminus1'-1)*`erate'))  * `ndrate'      ///   Converge to the new ndrate in 1/erate years
		+ (1 - max(0,min(1,(year -`vminus1' -1)*`erate')))     * `ndrate_dfltv'   if year > `vminus1'

	local erate_dflt   = 0.11         //  Value from parmsbycd4
	gen `eratev'  = max(0,min(1,(year -`vminus1'-1)*`erate'))   * `erate'      ///   Converge to the new erate in 1/erate years
		+ (1 - max(0,min(1,(year -`vminus1' -1)*`erate')))      * `erate_dflt'    if year > `vminus1'

}


*	Loop to simulate AIDS epidemics
********************************************
`shhh' {
di _n as txt "Loop to simulate AIDS epidemics"

local hrzn=scalar(`horizon')
forvalues t = `vintage'/`hrzn' {  /* If system is recursive, the loop is unnecessary */

	di _n as txt "Year = " as res `t'
	replace cvrgpct = 100*(L.a1 + L.a2 + L.b)/(L.un + L.a1 + L.a2 + L.b) 	if year == `t'

*	Endogenize HIV incidence as function of previous year's ART coverage                       Version 4.00...
if "`iexog'"== "" {
    replace p   = L.h + L.un + L.a1 + L.a2 + L.b                 		if year > `vminus1' & year== `t'
    replace prevprop   =  p/pop15_64_ 		                     		if year==`t' // Current adult prevalence proportion
    replace artf   =  min(1,max(0,(L.a1 + L.a2 + L.b)/p  ))      		if year==`t' // Art as fraction of all HIV+
    replace tau = 1  -  `fp' * `gp' * artf * (1 - `maxdt' * artf)		if year==`t' // Treatment reduction of incidence rate  Version 4.03

/*	Versions before 4.2 used this expression for the incidence rate in the core:
    replace ir = `ir0' *                       /// Starting incidence rate modified by...
        (1 -  ccvrg * 0.5 * `csp' ) *          /// Circumcision benefit (assumes males are .5 of pop)
        (1 -  effort ) *                       /// Prevention effort (non-circumcision)
        (1 +  `maxdt' * artf ) *               /// Disinhibition of treated
        (prevprop*tau)/(prevprop0*tau0)        /// Treatment effect
                                                                          if year == `t'
*/
    replace ir = `ir0' *                       /// Starting incidence rate modified by...
        (1 -  ccvrg * 0.5 * `csp' ) *          /// Circumcision benefit (assumes males are .5 of pop)
        (1 -  effort ) *                       /// Prevention effort (non-circumcision)
        (1 +  `maxdt' * artf ) *               /// Disinhibition of treated
        (k-prevprop)/(k-prevprop0) *           /// Saturation of the core  (New on 9/20/11)
        (prevprop*tau)/(prevprop0*tau0)        /// Treatment effect
                                                                          if year == `t'

    replace kmod = (1-artf)*k + artf * (min(0.5*(1-k), `maxdu'*k) + k ) if year == `t'

    replace i = round(                         /// Absolute incidence in Version 4.00
        ir * pop15_64_ *                       /// Incidence rate times population size
        (kmod - prevprop0) *                   /// Susceptible fraction = Core fraction - Initial prevalence
        (1 - `ve'*vr)   )                      /// Vaccination effect (assuming one were to exist)
                                                                            if year==`t'
}   //  End of code to endogenize incidence

*	"Healthy" HIV+
	replace h  = L.h  + L.i    - (`eratev' * L.h )		if year == `t'

*	Unmet need for treatment
	replace un  = L.un                          /*   Those present last year
	*/  + (`eratev' * L.h )                     /*   Plus those newly in need (v4.13 uses eratev)
	*/	- (1 - `lambda') * L.un                 /*   Less those recruited to ART
	*/	- `ndratev' * `lambda' * L.un           /*   Less deaths  (v4.13 uses ndratev)
	*/																if year == `t'

*	First year of first line ART
	replace a1  = L.a1                          /*   Those present last year
	*/  + (1 - `lambda') * L.un                 /*   Plus those recruited to ART
	*/	- (1-`adrate1')  * L.a1                 /*   Survivors move to second & subsequent years 1st line
	*/	- `adrate1'      * L.a1                 /*   Less deaths
	*/																if year == `t'

*	Assure that a1 is never negative
	replace a1  = 0 												if a1 < 0 & year==`t'

*	Second and subsequent years of first line ART
	replace a2  = L.a2                          /*   Those present last year
	*/  + (1-`adrate1') * L.a1                  /*   Plus survivors of 1st year
	*/	- `adrate2' *    cvrg2    * L.a2        /*   Less those who move to 2nd line
	*/	- `adrate2' * (1 - cvrg2) * L.a2        /*   Less deaths
	*/																if year == `t'

*	Second line ART
	replace b  = L.b                            /*   Those present last year
	*/  + `adrate2' *    cvrg2    * L.a2        /*   Plus those coming from 1st line ART
	*/	- `bdrate'  * L.b                       /*   Less deaths & loss to follow-up
	*/																if year == `t'

	di as txt _n  ///   Added in version 5.24
	"sum p h un a1 a2 b cvrg2 if year == `t'"
	sum p h un a1 a2 b cvrg2 if year == `t'		// Added in version 5.24

}    /* End of projection loop */
}    //  End of quietly braces for the loop

*	Annual deaths
/*	Version 4.12 and prior code
qui replace d =   `ndrate' * `lambda' * L.un     /* Deaths of those waiting for treatment
	*/  + `adrate1' * L.a1                       /* Deaths during 1st year of first line ART
	*/  + `adrate2' * (1 - cvrg2) * L.a2         /* Deaths during 2 + years of first line ART
	*/  + `bdrate'  * L.b                        /* Deaths during 2nd line ART
	*/																if year > `vminus1'
*/
*	Version 4.13 replaces ndrate with ndratev
qui replace d =   `ndratev' * `lambda' * L.un    /* Deaths of those waiting for treatment
	*/  + `adrate1' * L.a1                       /* Deaths during 1st year of first line ART
	*/  + `adrate2' * (1 - cvrg2) * L.a2         /* Deaths during 2 + years of first line ART
	*/  + `bdrate'  * L.b                        /* Deaths during 2nd line ART
	*/																if year > `vminus1'

*	Version 4.14 smooths the death variable
if "`smooth'" ~="nosmooth" & `horizon' >= `tkoff' +6 {
	tempvar dsmooth dunsmooth
	qui gen `dunsmooth' = d
	`shhh' di _n as txt "Smooth the death projections:"
	`shhh' tssmooth ma `dsmooth' =  d, window(4,1,2)
	qui replace d = `dsmooth'											if year > `vminus1'
}

if "`debug'"~="" {
	if "`smooth'" ~="nosmooth" {
		gen dbg_dunsmooth = `dunsmooth'             if year > `vminus1'
		label var dbg_dunsmooth "Annual deaths prior to smoothing"
	}
	gen dbg_L_d = L.d                               if year > `vminus1'
	gen dbg_ndrate_dfltv = `ndrate_dfltv'
	gen dbg_ndrate = `ndratev'                       if year > `vminus1'
	gen dbg_lambda = `lambda'                       if year > `vminus1'
	gen dbg_sigma = 1 - `lambda'                    if year > `vminus1'
	gen dbg_L_un = L.un                             if year > `vminus1'
	gen dbg_dA=dbg_ndrate * dbg_lambda * dbg_L_un	if year > `vminus1'  //  Part A of deaths

	gen dbg_adrate1 = `adrate1'                     if year > `vminus1'
	gen dbg_L_a1 = L.a1                             if year > `vminus1'
	gen dbg_dB=dbg_adrate1 * dbg_L_a1				if year > `vminus1'  //  Part B of deaths

	gen dbg_adrate2 = `adrate2'                     if year > `vminus1'
	gen dbg_cvrg2 = cvrg2                           if year > `vminus1'
	gen dbg_L_a2 = L.a2                             if year > `vminus1'
	gen dbg_dC=dbg_adrate2*(1-dbg_cvrg2)*dbg_L_a2	if year > `vminus1'  //  Part C of deaths

	gen dbg_bdrate = `bdrate'                       if year > `vminus1'
	gen dbg_L_b = L.b                               if year > `vminus1'
	gen dbg_dD=dbg_bdrate * dbg_L_b					if year > `vminus1'  //  Part D of deaths

	gen dbg_erate = `eratev'                        if year > `vminus1'
	gen dbg_erate_dflt = `erate_dflt'               if year > `vminus1'
}

*	Recast all values as longs
foreach X in h un a1 a2 b d i {
	qui replace `X' = round(`X')
	qui recast long `X'
}

*	Total need for ART
qui replace n   = un + a1 + a2 + b  									if year > `vminus1'
qui replace art =      a1 + a2 + b  									if year > `vminus1'

*	Prevalence of HIV+							// Replace lagged value of p with current value
qui replace p   = h + un + a1 + a2 + b		    						if year > `vminus1'
												// And assure that p is legitimate
capture assert p < pop15_64_  | pop15_64_>=. 							if year > `vminus1'
if _rc~=0 {
	noi di as err "For the following countries and for years between " as res "`vminus1'" as err " and " as res "`hrzn'"  ///
		_n as err "the variable p, the number of HIV infected, is missing or exceeds the adult population." ///
		_n as err "Execution terminating. Specify one of the two options: "  ///
		as res "-complete-" as err " or " as res "-impute-"
	tab nation if ~(p < pop15_64_  | pop15_64_>=.) & year >= `vminus1'
	error 9
}

*	Compute cost of ART
if `"`cost'"'~="revert" {
    qui gen  lgtrxc1 = -6.9 + 1.02 * ln(gnipc_hat)                     // logit bounded between rxc1lb & rxc1ub
        lab var lgtrxc1 "Predicted logit of 1st line drug cost"
        qui replace rxc1 = (`rxc1ub' * exp(lgtrxc1) + `rxc1lb')/(exp(lgtrxc1) + 1)
        lab var rxc1 "1st line drug cost, annual"
    qui gen  lgtrxc2 = -10.4 + 1.34 * ln(gnipc_hat)                    // logit bounded between rxc2lb & rxc2ub
        lab var lgtrxc2 "Predicted logit of 2nd line drug cost"
        qui replace rxc2 = (`rxc2ub' * exp(lgtrxc2) + `rxc2lb')/(exp(lgtrxc2) + 1)
        lab var rxc2 "2nd line drug cost, annual"
        tempvar ln_gnipcppp
    qui gen  `ln_gnipcppp' = 2.405 + 0.781 * ln(gnipc_hat)
        lab var `ln_gnipcppp' "Fitted values from a regression of gnipcppp on gnipc_hat in 2006"
    qui gen  hsbedc = exp(-2.75 + 0.77 * `ln_gnipcppp')
        lab var hsbedc "Cost of hospital bed day, Secondary"
        notes hsbedc : "Relationship "
        notes hsbedc : by Taghreed Adam and David B Evans, forthcoming
    qui gen  hsvstc = exp(-4.8 + 0.89 * `ln_gnipcppp')
        lab var hsvstc "Cost of hospital visit, Secondary"
		notes hsvstc : "Determinants of variation in the cost of inpatient stays versus outpatient visits in hospitals: a multi-country analysis."
        notes hsvstc : "Soc Sci Med. 2006 Oct;63(7):1700-10."
        notes hsvstc : by Taghreed Adam and David B Evans
    qui gen  hcvstc = exp(-1.7 + 0.49 * `ln_gnipcppp')             // Not currently used
        lab var hcvstc "Predicted unit cost of health center visit"
        notes hcvstc : "Capacity utilization and the cost of primary care visits: Implications for the costs of scaling up health interventions"
        notes hcvstc : by Taghreed Adam, Steve Ebener, Benjamin Johns and David B. Evans
        notes hcvstc : Cost Effectiveness and Resource Allocation 2008, 6:22doi:10.1186/1478-7547-6-22

*	DEBUG start
if "`debug'" ~="" {
	di _n(5) in re "DEBUG Cost computation"
	des gnipc_hat lgtrxc1 rxc1 lgtrxc2 rxc2 `ln_gnipcppp' hsbedc hsvstc
	di _n(2) in re "Cost elements for all countries and years"
	sum gnipc_hat lgtrxc1 rxc1 lgtrxc2 rxc2 `ln_gnipcppp' hsbedc hsvstc
	di _n(2) in re "Cost elements for Zambia in all years"
	sum gnipc_hat lgtrxc1 rxc1 lgtrxc2 rxc2 `ln_gnipcppp' hsbedc hsvstc if ccode=="ZMB"
	foreach yr of numlist 2015(15)2050 {
		di _n as txt "Cost elements for Zambia for the year: " as res "`yr'"
		sum gnipc_hat lgtrxc1 rxc1 lgtrxc2 rxc2 `ln_gnipcppp' hsbedc hsvstc if ccode=="ZMB" & year==`yr'
	}
	di in re "End of DEBUG Cost computation"
}
*	DEBUG end

    if `"`cost'"'=="function" {
*   Followiing coefficients are from -avg_trtmnt_cost.do- 10/27/2009
*                          -smear-                                                  -yr23-   -yr3-   (Subtract off estimated first-line drug cost)
        qui gen  nonrxcaf = 1.757 * exp( 3.224 + 0.706 * ln(gnipc_hat) - 0.142 * ln(art) -.801 -.186 ) - rxc1
            qui replace nonrxcaf = 0 if nonrxcaf<0
            qui replace nonrxcaf = exp((ln(`ndaf') - `scale'*6.91) + `scale'*ln(art))  if nonrxcaf>=.
         lab var nonrxcaf "Average fixed non-drug cost"
         notes nonrxcaf : Calculated from regression of ln(total ART financing/art) on ln(art) ln(gnipc_hat) and year dummies
    }
    else {
        qui gen  nonrxcaf = exp((ln(`ndaf') - `scale'*6.91) + `scale'*ln(art))  //  `scale' :=  elasticity <=0
         lab var nonrxcaf "Average fixed non-drug cost"
         notes nonrxcaf : Calculated from regression of ln(total ART financing/art) on ln(art) less $250
         notes nonrxcaf : Choice of a scale elasticity other than -0.146 pivots the line around the point [ln(art),ln(nonrxcaf)]=[ln(1000 patients),ln($750)]
    }

    qui gen  nonrxcav = `hsbedn'*hsbedc + `hsvstn'*hsvstc
        lab var nonrxcav "Average variable non-drug cost"
        format %9.0fc nonrxcav nonrxcaf
    drop lgtrxc1 lgtrxc2
    qui replace nonrxc = nonrxcav + nonrxcaf

    local costbrkdn costbrkdn(nonrxcav nonrxcaf)
}

qui gen long cost1st     = ((rxc1 + nonrxc) * (a1 + a2) ) / 1000
	lab var cost1st "Cost of 1st Line ART ($1,000's)"
	char define cost1st[scale] 1000
qui gen long cost2nd     = ((rxc2 + nonrxc) * b) / 1000
	lab var cost2nd "Cost of 2nd Line ART ($1,000's)"
	char define cost2nd[scale] 1000
qui gen long costtot     = cost1st + cost2nd
	lab var costtot "Current total fiscal burden of ART ($1,000's)"
	char define costtot[scale] 1000
qui gen fundgap = ((rxc1 + nonrxc) * ( un ) + (rxc2 + nonrxc) * `adrate2' * (1 - cvrg2) * L.a2) / 1000
	qui replace fundgap = ((rxc1 + nonrxc) * ( un ) + (rxc2 + nonrxc) * `adrate2' * (1 - cvrg2) * a2) / 1000 if fundgap>=.
	lab var fundgap "Funding gap ($1,000's)"
	char define fundgap[scale] 1000
	notes fundgap : "Funding gap is additional amount to serve those who need care in current year"
qui gen long title1st =  int( ((rxc1 + nonrxc) * (a1 + a2 - (1 - `lambda') * L.un ) ) / 1000 )   if year > `vminus1'
	qui replace         title1st = cost1st if year==`vminus1'
	lab var title1st "Entitlement, 1st line ($1,000's)"
	char define title1st[scale] 1000
qui gen long titletot =  title1st + cost2nd
	lab var titletot "Entitlement, Total ($1,000's)"
	char define titletot[scale] 1000
qui gen titlepct = 100*titletot/costtot
	lab var titlepct "Entitlement as percent of costs"
*	assert titlepct<=100 if year > `vminus1'
*	assert titlepct>=0

qui gen long disctot =  costtot - titletot  // Tautologically = 0 the first year
	lab var disctot "Discretionary ($1,000's)"
	char define disctot[scale] 1000
qui gen discpct = 100*disctot/costtot
	lab var discpct "Discretionary as percent of costs"
*	assert discpct<=100 if year > `vminus1'
*	assert discpct>=0

format  %15.0fc cost1st cost2nd costtot fundgap title1st titletot disctot
format  %9.0fc rxc1 rxc2 nonrxc
format  %9.1f  titlepct discpct cvrgpct

qui gen _uptake = 1 - (`lambda')
	lab var _uptake "Version of uptake variable used in projections"

if "`weo'"~="" {    /* N.B. check whether these projections should cross one another  */
	qui {
	tempvar mrg4
	merge code year using `"`tmpfl'"', unique
	sort code year
	by code : egen `mrg4' = max(_merge==3)    /* 1 if country merges in any year */
	qui keep if `mrg4'==1                         /* Keep all annual observations > `vminus1' for any country that merged */

	qui extweo cggdp, level cnvrgyr(2025)     /* Extrapolate data to end of horizon */
*	qui extweo gdpuskp, rate cnvrgyr(2040)    /* Extrapolate data to end of horizon */ // For Version 4.10 this is moved up top
*	qui extweo pop15_64_, rate cnvrgyr(2040)  /* Extrapolate data to end of horizon */ // For Version 4.x this is moved up top

*	Question arises whether to model HE and Public HE as a constant
*	or a rising share of GDP.  Cross country regression shows
*	little evidence of these shares rising with GDP inside Africa.
*		ln(pubhepctgdp[it]) = C[i] +  Elast * ln(gnipc_hat[it])
*	expr: reg ln(hepctgdp2009)  ln(gnipc_hat) if year==2010 & rgn==6
*	yields an Elast estimate of zero in Africa.
*	Should we decide to model this, we can use the following parameters:
	local hepctelast = 0	//	Elast of public share of total HE
	local pubpctelast = 0	//	Elast of public share of total HE
	local ooppctelast = 0	//	Elast of uninsured share of private HE
	local extpctelast = 0	//	Elast of external share of total HE
	gen double he =  ///
		cond(hepctgdp`vwdi'/100>0 & hepctgdp`vwdi'/100<1, hepctgdp`vwdi'/100, .) * gdpuskp*1000
		lab var he "Total health expenditure ($1,000,000's)"
	gen double pubhe =  ///
		cond(pubhepctgdp`vwdi'/100>0 & pubhepctgdp`vwdi'/100<1, pubhepctgdp`vwdi'/100, .) * gdpuskp*1000
		lab var pubhe "Public health expenditure ($1,000,000's)"
	gen double exthe =  ///
		cond(exthepcttot`vwdi'/100>0 & exthepcttot`vwdi'/100<1, exthepcttot`vwdi'/100, .) * he
		lab var exthe "External health expenditure ($1,000,000's)"
	gen double privhe = he - pubhe
		lab var privhe "Private health expenditure ($1,000,000's)"
	gen double oophe = ///
		cond(oophepctpriv`vwdi'/100>0 & oophepctpriv`vwdi'/100<1, oophepctpriv`vwdi'/100, .) * privhe
		lab var oophe "Out-of-pocket health expenditure ($1,000,000's)"

	format  %9.1f  he pubhe privhe oophe *hepct*

*	gen artpcthe = costtot/(hepctgdp`vwdi'*gdpuskp*100)
	gen artpcthe = 100*costtot/(he*1000)
		lab var artpcthe "ART as % total health spending"

*	gen artpctpub = costtot/(pubhepctgdp`vwdi'*gdpuskp*100)
	gen artpctpub = 100*costtot/(pubhe*1000)
		lab var artpctpub "ART as % public health spending"

*	gen artpctprv = costtot/((hepctgdp`vwdi' - pubhepctgdp`vwdi') * gdpuskp*100)
	gen artpctprv = 100*costtot/(privhe*1000)
		lab var artpctprv "ART as % private health spending"

*	gen artpctoop = costtot/((oophepctpriv`vwdi'*(hepctgdp`vwdi' - pubhepctgdp`vwdi')/100) * gdpuskp*100)
	gen artpctoop = 100*costtot/(oophe*1000)
		lab var artpctoop "ART as % OOP health spending"

	gen artpctcg = 100*costtot/(max(cggdp/100,0)* gdpuskp*1000)
		lab var artpctcg "ART as % central govt spending"

	gen artpctpop = 100*art/pop15_64_
		lab var artpctpop "Adult ART patients as % population 15-64"

	format  %9.1f  artpcthe artpctpub artpctprv artpctoop
	format  %9.1f  artpctcg artpctpop

	drop _merge
	drop if year > `horizon'
	}
}

label data "AIDS Projections with uptake parameter = `uptake'"
notes : Data created and saved by -aidsproj-, version ${apversion} for projections on the following country(ies): `countries' and the following regions: `regions'
notes : ART parameters are: uptake = `uptake'; cd4 = `cd4strng'; adrate1 = `adr1strng'; adrate2 = `adr2strng'; bdrate = `bdrstrng'; erate = `estrng'; ndrate = `nstrng'
if "`iexog'"~="" {
    notes : Incidence is assumed exogenous to ART, prevention effort or vaccination with the annual rate of incidence reduction being:  incmult = `imstrng'
}
else {
    notes : Incidence is assumed endogenous to ART: Fp = `fpstg'; Gp = `gpstg'; maxDt = `maxdtstg'; maxDu = `maxdustg'.
    notes : Incidence is potentially affected by vaccination or prevention effort: ir0 = `ir0stg'; minCp = `mincpstg'; maxCp = `maxcpstg'; CSp = `cspstg'; minEp = `minepstg'; maxEp = `maxepstg'; pstrtyr = `pstrtyrstg'; ptrgtyr = `ptrgtyrstg'; vstrtyr = `vstrtyrstg'; vtrgtyr = `vtrgtyrstg'; maxVp = `maxvpstg'; ve = `vestg'
}
char _dta[vintage] `vintage'
char _dta[morebadgood] present

relabel `rlblarg'  /* Add the stockflow and morebadgood characteristics */
if `"`saving'"' ~= "" {  /* Save the projections */
    local starpos = strpos(`"`saving'"',"*")
    if `starpos' {
        di as error _n "Option -saving- may not contain the punctuation mark * "
        di as error    "Use the -saving- option to specify only the stub of a filename."
        error 198
    }
	local rplpos = strpos(`"`saving'"',"replace")
	if `rplpos' > 0 {
		local saving = subinstr(`"`saving'"',"replace","",1)
		local replace "replace"
	}
	local saving = trim(subinstr(`"`saving'"',",","",1))
    if strpos(`"`outpath'"',"\") {
		local outpath = subinstr(`"`outpath'"',"\","/", .)
		local outpath = subinstr(`"`outpath'"',"//","/", .)  // In case there are multiple slashes,
		local outpath = subinstr(`"`outpath'"',"//","/", .)  // repeat this command.
		local outpath = subinstr(`"`outpath'"',"//","/", .)
		if substr("`outpath'",strlen("`outpath'"),1)~="/"  {
			local outpath = "`outpath'" + "/"
		}
	}

	preserve
	qui compress
	qui drop __*
	di _n
	if  c(stata_version) > 12 {
		saveold `"`outpath'`saving'"',`replace'  version(12)
	}
	else {
		save `"`outpath'`saving'"',`replace'
	}
	restore
}

if "`shhh'" =="" {
	di _n as txt "xtdes"
	xtdes
}

if "`output'" ~= "" {
    	local grnote1 "ART parameters: cd4 = `cd4strng'; adrate1 = `adr1strng'; adrate2 = `adr2strng'; bdrate = `bdrstrng'; erate = `estrng'; ndrate = `nstrng'"
	if "`iexog'"~="" {  // If incidence is exogenous; grnote is only one line:
		local grnote2 "Incidence is determined without reference to ART with parameters: `iexog' and 'incmult = `imstrng'"

    }
    else {  // N.B.: Semi-colons instead of commas - or it breaks
		local grnote2 "Incidence is a function of ART coverage and quality, with parameters: ir0 = `ir0stg' and ... "
		local grnote3 grnote3("    ART: gp = `gpstg'; fp = `fpstg'; maxdt = `maxdtstg'; maxdu = `maxdustg'; trgtcov2 = `trgtcov2stg'; trgtyr = `trgtyrstg'")
		local grnote4 grnote4("    Male circumcision: mincp = `mincpstg'; maxcp = `maxcpstg'; csp = `cspstg'; cstrtyr = `cstrtyrstg'; ctrgtyr = `ctrgtyrstg'")
		//  Making following definition of grnote5 conditional on -maxep->0 weirdly also knocks out grnote6. (????)
		local grnote5 grnote5("    Prevention effort: minep = `minepstg'; maxep = `maxepstg'; pstrtyr = `pstrtyrstg'; ptrgtyr = `ptrgtyrstg'")
		if `ve' * `maxvp' > 0 {  //  Both -ve- and -maxvp- must be g> 0 to display this note
			local grnote6 grnote6("    Vaccination: maxvp = `maxvpstg'; ve = `vestg'; vstrtyr = `vstrtyrstg'; vtrgtyr = `vtrgtyrstg'")
		}
    }

	`shhh' di _n as txt "Output requested for: " as result "`output'"
	`shhh' di _n
	`shhh' notes _dta
	if "`uptkmeth'"=="VARIABLE" {
		local uparg "sigmav(`sigmaV')"
	}
	else {
		local uparg "sigmac(`sigmaC')"
	}
	if "`outpath'"~="" {
		local outpathopt outpath(`outpath')
	}
	genout , `uparg'  uptkmeth(`uptkmeth')  `countries' `regions' `world' /*
		*/ `gr7' `graph' `grcombine' `legend' /*
		*/ `titlen' `titlec' `titlep' `titlee' `subtitle' /*
		*/ saving(`"`outpath'`saving'"') `replace' `rescale' `weo' `logout' /*
		*/ weights(`weights') vintage(`vintage') `costbrkdn' `drate' `refyear' /*
		*/ `postyears' `posttotal' `postvars'  `postfile' `postrunspecs'  /*
		*/ `postopen' `postclose' `postname' `takeoff' `putexcel' `fudge'  /*
		*/ `outpathopt' /*
		*/ grnote1(`grnote1') grnote2(`grnote2') `grnote3' `grnote4' `grnote5' `grnote6'
}
else {
	di _n as txt "No output requested."
	if `"`postyears'`posttotal'`postvars'`postfile'`postrunspecs'`postopen'`postclose'`postname'"' ~= "" {
    	di as err "Since no output has been requested, nothing can be harvested for export to Excel either."
		error 198
	}
}

*	Print footer
	di _n _n as txt "AIDSPROJ: Completion of -aidsproj- run"
		local len1 = length("AIDSPROJ: Completion of -aidsproj- run")

	if "`saving'"~="" {
		di as txt "Projections saved in file: " as res `"`outpath'`saving'"'
		local len2 = length("Projections saved in file: ") + length(`"`outpath'`saving'"')
		local len = min(max(`len1',`len2'),c(linesize))
		di  as txt "{hline `len'}" _n
	}
	else {
		di as txt "Projections not saved."
		local len = `len1'
		di  as txt "{hline `len'}" _n
	}

	set linesize `old_linesize'

end /* of main program -aidsproj- */

*************************SUBROUTINES************************************
program define chkdta, rclass
*	Check the country & region lists to make sure they are in the data
	syntax [, countries(namelist) regions(namelist) ]

	if "`countries'"=="" & "`regions'"=="" {
		di as error "This error should never occur"
		exit 198
	}
	if "`countries'"'~="" {
		confirm numeric variable code
		tempvar cntrystr
		qui decode code, gen(`cntrystr')
		foreach X in `countries' {
			qui count if `cntrystr'==upper("`X'")
			if r(N) == 0 {
				di _n as error "Country code `X' is not in the data."
				exit 198
			}
		}
		return local countries = upper("`countries'")
	}
	if "`regions'"~="" {
		confirm numeric variable rgn
		tempvar rgnstr
		qui decode rgn,gen(`rgnstr')
		foreach X in `regions' {
			qui count if `rgnstr'==upper("`X'")
			if r(N) == 0 {
				di _n as error "Region code `X' is not in the data."
				exit 198
			}
		}
		return local regions = upper("`regions'")
	}
end    /* End of program chkdta */

program define select
*	Impute additional information for some countries.
	syntax varlist [if] [, quietly complete impute]

	if "`complete'"~="" & "`impute'"~="" {
		di _n as error "Only one of the options -complete- & -impute- may be selected."
		exit 198
	}

	`quietly' {

	*	Filter the countries on -if-
		tempvar touse
		mark `touse' `if'

	*	Filter the countries on whether data is complete
		if "`complete'"~="" {
			noi di as txt "Option " as result "complete" as txt " selected. Only include countries with complete epidemiological data."
			replace `touse' = 0 if complete==0
		}
		if "`impute'"~="" {
			noi di as txt "Option " as result "impute" as txt " selected. Include additional countries"  ///
			_n "for which the varialbles " as res "`varlist'" as txt " can be imputed."

			//  Alternative imputation code new in version 4.59
			//  First impute the epidemiological data, but only if only one of p, a, i or d is missing
			local paid_list p a i d n
			foreach var of varlist `paid_list' {
				tempvar `var'tmp ln`var'
				local tmp_list `tmp_list' ``var'tmp'
				local ln_list `ln_list' `ln`var''
				qui gen double `ln`var'' = ln(`var')
			}

			foreach var of varlist `paid_list'  {
				local depvar `ln`var''
				local rhsvars : list ln_list - depvar
				local rhsnames : list paid_list - var
				di _n as txt "Interpolate missing values of " as res "`var'"  ///
					as txt " from observed values of " as res "`rhsnames'"
				reg `depvar' `rhsvars'
					tempvar `depvar'hat `var'hat
					qui predict ``depvar'hat' if `depvar'>=.
					// Assure that minimum imputed value is at leat 1.0
					qui gen long ``var'hat' = max(exp(``depvar'hat'),1) if ``depvar'hat'<.
				noi di _n as txt "values of " as res "`var'" as txt " are interpolated for the following country(ies):"
				noi di    as txt "list cntry `var' `var'hat exp(`var'hat) if `var'>=. & `var'hat<."
				noi list cntry `depvar' ``depvar'hat' ``var'hat' if `depvar'>=. & ``depvar'hat'<.
			}

			gen byte imputed = 0
				label var imputed "Country-specific imputation indicator"
				label define imputed 0 "0 vrbls imputed" 1 ">=1 vrbls imputed"
				label values imputed imputed

			di _n as txt "Replace the missing values of the epi variables with the imputed values"
			foreach var of varlist `paid_list'  {
				local rhsnames : list paid_list - var
				qui replace imputed = 1 if `var' >= . & ``var'hat'<.
				noi repshow `var' = ``var'hat' if `var' >= . & ``var'hat'<., id(cntry) listvars(ccode `rhsnames' complete)
				drop ```ln`var'''hat' ``var'hat'
			}

			di _n as txt "Impute the missing values of the other variables with region averages"
			local othvars: list varlist - paid_list
			qui levelsof rgn, local(rgns)
			foreach var in `othvars' {
				tempvar `var'hat
				qui gen double ``var'hat' = .
				foreach r of numlist `rgns' {
					sum `var' if rgn==`r', meanonly
					qui replace ``var'hat' = r(mean) if rgn==`r'
				}
				qui replace imputed = 1 if `var' >= . & ``var'hat'<.
				noi repshow `var' = round(``var'hat') if `var' >=. & ``var'hat'<., id(cntry) listvars(rgn complete)
			}
			*	Define the subsample of data to be kept and used for projections
			markout `touse' `varlist'

			di _n as txt "Classify countries as -complete- or -imputed- before dropping observations"
			tab complete imputed, mi

			di _n as txt "List of countries to be dropped for lack of data on one or more variables"
			list ccode `varlist' if `touse' == 0

			di _n as txt "Classify countries as -complete- or -imputed- after dropping observations"
			tab complete imputed if `touse', mi

		}  //  End of the impute code
		qui keep if `touse'
	}  //  End of `quietly' loop

end  /* End of program select */

program define genout
*	Generate output from aidsproj
	syntax  [using] ,  uptkmeth(string) [ sigmac(real 0.5) sigmav(varname) countries(namelist) /*
		*/ regions(namelist) world gr7 graph grcombine legend(passthru) /*
		*/ grnote1(passthru) grnote2(passthru) grnote3(passthru) grnote4(passthru) grnote5(passthru) grnote6(passthru)  /*
		*/ titlen(passthru) titlec(passthru) titlep(passthru) titlee(passthru) subtitle(passthru)  /*
		*/ saving(string) replace norescale weo logout weights(varname) vintage(real 2008) /*
		*/ costbrkdn(string) drate(string) refyear(real 2008)  takeoff(string) /*
		*/ postyears(passthru) posttotal postvars(varlist) postfile(passthru) /*
		*/ postrunspecs(passthru) postopen postnow postclose postname(passthru) putexcel fudge(int 100) /*
		*/ outpath(passthru)]

	if "`countries'"=="" & "`regions'"=="" & "`world'"=="" {
		di as error "This error should never occur"
		exit 198
	}

	if "`drate'"~="" {
		tempname dfactor
		local rate = real("`drate'")
		if `rate' <0 | `rate'>1 {
			di _n as err "Discount rate must be between zero and one."
			exit 198
		}
		if `refyear'< `vintage' | `refyear'>2011 {
			di _n as err "Reference year (-refyear-) must be between `vintage' and 2011."
			exit 198
		}
		qui gen `dfactor' = 1/(1 + real("`rate'"))^(year - `refyear') if year >= `refyear'
			qui replace `dfactor' = 0 if year < `refyear'
		local discwt [aw=`dfactor']
		local discalert `" "Values in this table have been discounted to the year " as res "`refyear'" as txt " at the rate " as res "`rate'" as txt"." "'
	}


	if "`postopen'"~="" {
		aidspost `postvars', `postyears' `posttotal' `postfile' countries(`countries') regions(`regions') `world' `postrunspecs' postopen `postname'
	}
	if "`postyears'`postvars'`postfile'" ~="" {
		local tabsave save
		local harvest	"aidspost `postvars', `postyears' `posttotal'  \`group'  harvest `postname'"  /* countries(`countries') regions(`regions') `world' */
		local postnow postnow
	}	/* For option -postnow-, see the end of program -genout-  */

	if "`uptkmeth'"=="VARIABLE" {
		local sigma "`sigmav'"
	}
	else {
		tempvar sigma
		gen `sigma' = `sigmac'
	}

	if "`putexcel'"~="" & "`postopen'`postyears'`postvars'`postfile'" ~="" {
		di as err "The option, -putexcel-, is not compatible with the options -postopen-, -postfile-, etc."
		exit 198
	}
	if "`putexcel'"~="" {
		local tabsave save
		local putexcel putexcel A1 = matrix(r(StatAll),names) using `outpath'`saving'_\`Y' , sheet(\`sheetname')
		if "`replace'"~="" {
			local excelreplace replace
			local excelmodify modify
		}
	}

	local vminus1  = `vintage'-1

	if "`weo'"~="" {
		local weovars artpcthe artpctpub artpctprv artpctoop artpctcg artpctpop
		capture confirm numeric variable `weovars'
		if _rc  {
			di _n as error "WEO option requires the variables `weovars'"
			exit 198
		}
	}

	if "`takeoff'"~="" {
		local iftkoff  if year >= real("`takeoff'")
		local andtkoff & year >= real("`takeoff'")
	}

	*	Projections for countries
	if "`countries'"~="" {
		foreach X in `countries' {
			local group countries(`X')
			if "`logout'"~="" {
				log using `"`outpath'`saving'_cntry_`X'"',`replace' name(_aidsproj_logout_cntry_`X')
*				log using _logout_cntry_`X', replace name(aidsproj_cntry_`X')
			}
			local Y = upper("`X'")
			di _n as txt "Projections for Country " as result "`Y'"
			di    as txt "============================"
			des p i h n un d h un a1 a2 b art cvrgpct cvrg2 `ccvrg_vr' rxc1 rxc2  ///
				`costbrkdn' nonrxc cvrg2 `lambda' cost1st cost2nd costtot   ///
				fundgap disctot discpct titletot title1st titlepct `weovars'
			di _n as txt "----Epidemiologic projections for country " as result "`Y'" as txt"----"
				local sheetname Epidemiology
				di as txt  `discalert'
				tabstat2  p i h n un d    if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe  /* Epidemic */
					`harvest'
					`putexcel' `excelreplace'
			di _n as txt "----Treatment projections for country " as result "`Y'" as txt"----"
				local sheetname Treatment
				di as txt  `discalert'
				tabstat2  h un a1 a2 b art  if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe /* ART */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Prevalence and coverage projections for country " as result "`Y'" as txt"----"
				local sheetname Coverage
				tabstat2  prevprop k cvrgpct cvrg2 `ccvrg_vr'  if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe  /* Coverage */
					`harvest'
					`putexcel' `excelreplace'
			di _n as txt "----Unit costs and rates for country " as result "`Y'" as txt"----"
				local sheetname UnitCosts
				di as txt  `discalert'
				tabstat2  rxc1 rxc2 `costbrkdn' nonrxc cvrg2 `lambda' if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe /* Cost, Unit */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Total cost projections (US$1,000's) for country " as result "`Y'" as txt"----"
				local sheetname TotalCosts
				di as txt  `discalert'
				tabstat2  cost1st cost2nd costtot fundgap if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe/* Cost, Amount */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Spending for country (US$1,000's)----"
				local sheetname Disc&NondiscAmts
				di as txt  `discalert'
				tabstat2 disctot discpct titletot title1st titlepct if code=="`Y'":code `andtkoff' `discwt' , by(year) format mi noheader `tabsave' describe /* Cost, Discretionary Amounts */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Percentages for country " as result "`Y'" as txt"----"
				local sheetname Disc&NondiscPct
				tabstat2 cvrgpct discpct titlepct if code=="`Y'":code `andtkoff', by(year) format mi noheader `tabsave' describe /* Cost, Discretionary percentages */
					`harvest'
					`putexcel' `excelmodify'
			if "`weo'"~="" {
				di _n as txt "----AIDS treatment costs as % of health and of central government expenditures for country " as result "`Y'" as txt"----"
					local sheetname RelativeCostBurden
					tabstat2 artpcthe artpctpub artpctprv artpctoop artpctpop if code=="`Y'":code  `andtkoff', by(year) format mi noheader `tabsave' describe /* Cost, Affordability */
						`harvest'
						`putexcel' `excelmodify'
			}
			di _n

			if "`gr7'"~="" | "`graph'"~="" {
				makegph [fw=`weights'],  sigma(`sigma') uptkmeth(`uptkmeth') country(`X') `gr7' `graph' `grcombine' `legend' `titlen' `titlec' `titlep' `titlee' `subtitle' saving(`"`outpath'`saving'"') `replace' `rescale' `weo' discalert(`discalert') discwt(`discwt') takeoff(`takeoff') `outpath' `grnote1'  `grnote2' `grnote3' `grnote4' `grnote5' `grnote6'
			}

			if "`logout'"~="" {
				log close _aidsproj_logout_cntry_`X'
			}

		}
	}

	*	Projections for regions
	if "`regions'"~="" {
		foreach X in `regions' {
			local group regions(`X')
			if "`logout'"~="" {
				log using `"`outpath'`saving'_rgn_`X'"',`replace' name(_aidsproj_logout_rgn_`X')
			}
			local Y = upper("`X'")
			di _n as txt "Projections for Region " as result "`Y'"
			di    as txt "==========================="
			sum code if year==`vminus1' & rgn=="`Y'":rgn, meanonly
			di _n as txt "There are " as result "`r(N)'" as txt " countries included in the projection for this region"
			des p i h n un d h un a1 a2 b art cvrgpct cvrg2 `ccvrg_vr' rxc1 rxc2  ///
				`costbrkdn' nonrxc cvrg2 `lambda' cost1st cost2nd costtot   ///
				fundgap disctot discpct titletot title1st titlepct `weovars'
			di _n as txt "----Epidemiologic projections for region " as result "`Y'" as txt"----"
				local sheetname Epidemiology
				di as txt  `discalert'
				tabstat2 p i h n un d                  if rgn=="`Y'":rgn `andtkoff' `discwt' , by(year) stat(sum) format mi noheader `tabsave' describe /* Epidemic */
					`harvest'
					`putexcel' `excelreplace'
			di _n as txt "----Treatment projections for region " as result "`Y'" as txt"----"
				local sheetname Treatment
				di as txt  `discalert'
				tabstat2 h un a1 a2 b art               if rgn=="`Y'":rgn `andtkoff' `discwt' , by(year) stat(sum) format mi noheader `tabsave' describe /* ART */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Population prevalence and coverage projections for region " as result "`Y'" as txt"----"
			di    as txt "      Average weighted by " as result "population 15-64"
				local sheetname Coverage
				tabstat2  prevprop k cvrgpct `ccvrg_vr'    if rgn=="`Y'":rgn `andtkoff' [aw=pop15_64_], by(year) stat(mean) format mi noheader `tabsave' describe  /* Coverage */
					`harvest'
					`putexcel' `excelreplace'
			di _n as txt "----Unit costs and rates for region " as result "`Y'" as txt"----"
			di    as txt "      Weighted by " as result "`weights'"
				local sheetname UnitCosts
				tabstat2 rxc1 rxc2 `costbrkdn' nonrxc cvrg2 `lambda' if rgn=="`Y'":rgn `andtkoff' [freq=`weights'], by(year) stat(mean) format mi noheader `tabsave' describe /* Cost, Unit  */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Total cost projections (US$1,000's) for region " as result "`Y'" as txt"----"
				local sheetname TotalCosts
				di as txt  `discalert'
				tabstat2 cost1st cost2nd costtot fundgap if rgn=="`Y'":rgn `andtkoff' `discwt' , by(year) stat(sum) format mi noheader `tabsave' describe /* Cost, Amount */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Spending (US$1,000's) for region " as result "`Y'" as txt"----"
				local sheetname Disc&NondiscAmts
				di as txt  `discalert'
				tabstat2 disctot titletot title1st  if rgn=="`Y'":rgn `andtkoff' `discwt' , by(year) stat(sum) format mi noheader `tabsave' describe /* Cost, Discretionary Amounts*/
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Percentages for region " as result "`Y'" as txt"----"
			di    as txt "      Average weighted by " as result "`weights'"
				local sheetname Disc&NondiscPct
				tabstat2 cvrgpct discpct titlepct if rgn=="`Y'":rgn `andtkoff' [freq=`weights'], by(year) stat(mean) format mi noheader `tabsave' describe /* Cost, Discretionary percentages */
					`harvest'
					`putexcel' `excelmodify'
			if "`weo'"~="" {
				di _n as txt "----AIDS treatment costs as % of health and of central government expenditures for region " as result "`Y'" as txt"----"
				di    as txt "      Weighted by " as result "`weights'"
					local sheetname RelativeCostBurden
					tabstat2 artpcthe artpctpub artpctprv artpctoop artpctpop if rgn=="`Y'":rgn `andtkoff' [freq=`weights'], by(year) stat(mean) format mi noheader `tabsave' describe /* Cost, Affordability */
						`harvest'
						`putexcel' `excelmodify'
			}
			di _n

			if "`gr7'"~="" | "`graph'"~="" {
				makegph [fw=`weights'],  sigma(`sigma') uptkmeth(`uptkmeth') region(`X') `gr7' `graph' `grcombine' `legend' `titlen' `titlec' `titlep' `titlee' `subtitle' saving(`"`outpath'`saving'"') `replace' `rescale' `weo' discalert(`discalert') discwt(`discwt') takeoff(`takeoff') `outpath' `grnote1'  `grnote2' `grnote3' `grnote4' `grnote5' `grnote6'
			}

			if "`logout'"~="" {
				log close _aidsproj_logout_rgn_`X'
			}
		}
	}

	*	Projections:  Global totals
	if "`world'" ~= "" {
		local group `world'
		local Y world
		if "`logout'"~="" {
			log using `"`outpath'`saving'_world"',`replace' name(_aidsproj_logout_world)
		}

		if `fudge'~=100  {  // Added version 4.56
			local fudgeprop = `fudge'/100
			foreach var of varlist  art n  ///
				p i h a1 a2 b un d a1 a2 b ///
				cost1st cost2nd costtot fundgap disctot titletot title1st {
					qui replace `var' = `fudgeprop'*`var'
					local fudgetext _n `"(All results have been multiplied by `fudgeprop')"'
					local padleft  "           "
					local padright "====================="
			}
			foreach var of varlist art n  {
				qui recast long `var' , force
			}
		}

		di _n as txt "`padleft'" as result "Global " as txt "projections" `fudgetext'
		di    as txt "=====================" "`padright'"
		sum code if year==`vminus1', meanonly
		di _n as txt "The " as result "`r(N)'" as txt " countries included in the projection"
		di    as txt "are distributed across the regions as follows:"
		tab rgn if year==`vminus1', mi
		des p i h n a1 a2 b art un d cvrgpct disctot discpct titletot title1st titlepct `weovars'
		di _n as txt "----Epidemiologic projections----"
				local sheetname Epidemiology
				di as txt  `discalert'
			tabstat2 p i h n un d     `discwt' `iftkoff' , by(year) stat(sum) format mi noheader casewise `tabsave' describe /* Epidemic */
				`harvest'
				`putexcel' `excelreplace'
		di _n as txt "----Treatment projections----"
				local sheetname Treatment
				di as txt  `discalert'
			tabstat2 h un a1 a2 b art   `discwt' `iftkoff' , by(year) stat(sum) format mi noheader casewise `tabsave' describe /* ART */
				`harvest'
				`putexcel' `excelmodify'
		di _n as txt "----Unit costs and rates (Weighted averages)----"
		di    as txt "      Weighted by " as result "`weights'"
				local sheetname UnitCosts
			tabstat2 rxc1 rxc2 `costbrkdn' nonrxc cvrg2 `lambda' [freq=`weights']`iftkoff' , by(year) stat(mean) format mi noheader `tabsave' describe /* Cost */
				`harvest'
				`putexcel' `excelmodify'
		di _n as txt "----Total cost projections (US$1,000's)----"
				local sheetname TotalCosts
				di as txt  `discalert'
			tabstat2 cost1st cost2nd costtot fundgap `discwt' `iftkoff' , by(year) stat(sum) format mi noheader `tabsave' describe casewise /* Cost */
				`harvest'
				`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Spending (US$1,000's) for all countries----"
				local sheetname Disc&NondiscAmts
				di as txt  `discalert'
				tabstat2 disctot titletot title1st `discwt' `iftkoff' , by(year) stat(sum) format mi noheader `tabsave' describe casewise /* Cost */
					`harvest'
					`putexcel' `excelmodify'
			di _n as txt "----Discretionary & Non-discretionary Percentages for for all countries----"
			di    as txt "      Weighted by " as result "`weights'"
				local sheetname Disc&NondiscPct
				tabstat2 cvrgpct discpct titlepct [freq=`weights']`iftkoff' , by(year) stat(mean) format mi noheader `tabsave' describe casewise /* Cost */
					`harvest'
					`putexcel' `excelmodify'
		if "`weo'"~="" {
			di _n as txt "----AIDS treatment costs as percentages of health and of central government expenditures----"
			di    as txt "      Weighted by " as result "`weights'"
				local sheetname RelativeCostBurden
				tabstat2 artpcthe artpctpub artpctprv artpctoop artpctpop [freq=`weights']`iftkoff' , by(year) stat(mean) format mi noheader `tabsave' describe casewise /* Cost */
					`harvest'
					`putexcel' `excelmodify'
		}
		if "`gr7'"~="" | "`graph'"~="" {
			makegph [fw=`weights'],  sigma(`sigma') uptkmeth(`uptkmeth') world `gr7' `graph' `grcombine' `legend' `titlen' `titlec' `titlep' `titlee' `subtitle' saving(`"`outpath'`saving'"') `replace' `rescale' `weo' discalert(`discalert') discwt(`discwt') takeoff(`takeoff') `outpath' `grnote1'  `grnote2' `grnote3' `grnote4' `grnote5' `grnote6'
		}

		if "`logout'"~="" {
			log close _aidsproj_logout_world
		}

	}

	if "`postnow'" ~= "" {
		aidspost `postvars', `postyears' `postfile' countries(`countries') regions(`regions') `world' `postrunspecs' postnow `postname'
	}
	if "`postclose'" ~= "" {
		aidspost `postvars', `postyears' `postfile' postclose countries(`countries') regions(`regions') `world' `postname'
	}
end  /* End of program genout */

program define makegph
*	Construct, display & save one or more graphs for a specific country, region or the world.
	syntax  [fw/], sigma(varname) uptkmeth(string) [country(string) region(string) world /*
		*/ gr7 graph grcombine legend(string) /*
		*/ grnote1(string asis) grnote2(string asis) grnote3(string /* asis */) grnote4(string /* asis */) grnote5(string /* asis */) grnote6(string /* asis */) /*
		*/ titlen(string) titlec(string) titlep(string) titlee(string) subtitle(string) /*
		*/ saving(string) replace norescale weo discalert(string) discwt(string) takeoff(string) /*
		*/ outpath(string)]

	tempvar tmpordr
		gen long `tmpordr' = _n

	if ("`country'"=="") + ("`region'"=="") + ("`world'"=="") ~= 2  {
		di _n as error "This subroutine is only called once for each unit of analysis, be it country, region or the world."
		di    as error "Thus, this error in aidsproj:makegph should never occur."
		exit 198
	}

	local suffix "`country'`region'"
	if "`suffix'"=="" {
		local suffix "World"
	}

	if `"`saving'"'~="" & "`replace'"~="" {
		local rpl ", replace"
	}

	if "`legend'"=="" {
		local legend  "position(10) col(1) size(small) ring(0)"
		local legendp "position(12) size(small) rowgap(*.1) span"
*		local legendt "position(6) size(small) ring(0)"
		local legendt "position(12) size(small) span"
	}
	else {
		local legendp "`legend'"
		local legendt "`legend'"
	}

	if "`titlen'"=="" {
		local titlen "AIDS Treatment: Projected patients and unmet need"
	}

	if "`titlet'"=="" {
		local titlet "AIDS Transition: Projected prevalence, incidence and deaths"
	}

	if "`titlec'"=="" {
		local titlec "AIDS Treatment: Projected cost & funding gap"
	}

	if "`titlep'"=="" {
		local titlep "AIDS Treatment: Entitlements as Percent Total"
	}

	if "`titlee'"=="" {
		local titlee "ART Expenditure as Percent Other Expenditures"
	}

	local gr7subtitle "   "
	if "`subtitle'"=="" {
		if "`country'" ~= "" {
			local subtitle  "Total for country \`suffix'"
		}
		if "`region'"  ~= "" {
			local subtitle  "Total for region \`suffix'"
		}
		if "`world'"   ~= "" {
			local subtitle  "Total for all countries"
		}
	}
	else {
		local gr7subtitle "`subtitle'"
	}

	if "`weo'"~="" {
		local weovars artpcthe artpctpub artpctprv artpctoop artpctpop
	}

	if "`country'"~="" {
		local country = upper("`country'")
		capture confirm numeric variable code
		if !_rc ~= 1 {
			di _n as error "The labeled numeric variable -code- to designate countries must be included in the input data"
			exit 198
		}
		local clbl : value label code
		if "`clbl'"~="code"  {
			di _n as error "The variable -code- must have an associated value label also named -code-"
			exit 198
		}
		preserve
		qui keep if code=="`country'":code
		local dlbl : data label
	*	Graphing the results for a single country also requires a collapse command
		collapse (mean) code `tmpordr' `sigma' discpct titlepct cvrgpct `weovars'  ///
			(rawsum) p i h n un art b d a2 a1 cost1st cost2nd costtot fundgap disctot titletot title1st   ///
				[fweight=`exp']  if year>=real("`takeoff'"), by(year)
			label data "`dlbl', total for country `country'"
			lab var code   "World Bank country code"
				lab values code code
			relabel `rlblarg' /* Relabel the variables after the collapse */

		sum `sigma' [fweight=n], meanonly
			global glblsigma = string(100*r(mean), "%9.1f")

		local gr7ttl = "Country `country', enrolling $glblsigma % per yr"
	}

	if "`region'"~="" {
		local region = upper("`region'")
		capture confirm numeric variable rgn
		if !_rc ~= 1 {
			di _n as error "The numeric variable -rgn- to designate rgns must be included in the input data"
			exit 198
		}
		local rgnlbl : value label rgn
		if "`rgnlbl'"~="rgn"  {
			di _n as error "The variable -rgn- must have an associated value label also named -rgn-"
			exit 198
		}
		preserve
		qui keep if rgn=="`region'":rgn
		local dlbl : data label
		collapse (mean) rgn `tmpordr' `sigma' discpct titlepct cvrgpct `weovars'   ///
			(rawsum) p i h n un art b d a2 a1 cost1st cost2nd costtot fundgap disctot titletot title1st  ///
				[fweight=`exp']  if year>=real("`takeoff'"), by(year)
			label data "`dlbl', total for region `region'"
			lab var rgn    "Region"
				lab values rgn rgn
			relabel `rlblarg' /* Relabel the variables after the collapse */

		sum `sigma' [fweight=n], meanonly
			global glblsigma = string(100*r(mean), "%9.1f")

		local gr7ttl = "Region `region', enrolling $glblsigma % per yr"

	}

	if "`world'"~="" {
		preserve
		local dlbl : data label
		collapse (mean) `tmpordr' `sigma' discpct titlepct cvrgpct `weovars'   ///
			(rawsum) p i h n un art b d a2 a1 cost1st cost2nd costtot fundgap disctot titletot title1st  ///
				[fweight=`exp']  if year>=real("`takeoff'"), by(year)
			label data "`dlbl', total for the world"
			relabel `rlblarg' /* Relabel the variables after the collapse */

		sum `sigma' [fweight=n], meanonly
			global glblsigma = string(100*r(mean), "%9.1f")

		local gr7ttl = "World, enrolling $glblsigma % per yr"
	}

	if `"`saving'"'~= "" {
		local svngn `"saving(`"`outpath'`saving'_n_`suffix'.gph"' `rpl')"' /* for ART numbers graphs */
		local svngt `"saving(`"`outpath'`saving'_t_`suffix'.gph"' `rpl')"' /* for AIDS transition graphs */
		local svngc `"saving(`"`outpath'`saving'_c_`suffix'.gph"' `rpl')"' /* for cost graphs */
		local svngp `"saving(`"`outpath'`saving'_p_`suffix'.gph"' `rpl')"' /* for percent discretionary graphs */
		local svnge `"saving(`"`outpath'`saving'_e_`suffix'.gph"' `rpl')"' /* for percent expenditure graphs */
		local svng  `"saving(`"`outpath'`saving'_`suffix'.gph"'   `rpl')"' /* for combined graphs */
	}
	else {
		local namen `"name(`"n_`suffix'"' , replace)"' /* for ART numbers graphs */
		local namet `"name(`"t_`suffix'"' , replace)"' /* for AIDS transition graphs */
		local namec `"name(`"c_`suffix'"' , replace)"' /* for cost graphs */
		local namep `"name(`"p_`suffix'"' , replace)"' /* for percent discretionary graphs */
		local namee `"name(`"e_`suffix'"' , replace)"' /* for percent expenditure graphs */
		local name  `"name(`"`suffix'"'   , replace)"' /* for combined graphs */
	}
	if "`rescale'"=="" {
		rescale art b un p, auto
			local unitlblstk  "`r(unitlbl)'"
		rescale d i, auto
			local unitlblflw  "`r(unitlbl)'"
		rescale fundgap costtot cost2nd, auto already(1000)
		local unitlblc  "`r(unitlbl)'"
	}
	sum d, meanonly
	local dmax = r(max)
	local dmin = r(min)
	sum i, meanonly
	local imax = r(max)
	local imin = r(min)
	local rsmax = max(`dmax', `imax')
	local rsmin = min(`dmin', `imin')
	sum year if i < d & d<. & i[_n-1]>=d[_n-1] & i[_n-1]<., meanonly
	local tyear = r(min)
	if `tyear' <. {
		local transition xline(`tyear') ttick(`tyear') tlabel(`tyear')
		local transttl "Transition in the year `tyear'"
	}
	local incidlbl : variable label i

	if "`gr7'"~="" {  // Use Version 7 graphs, fast but ugly
		if "`grcombine'"=="" {
			gr7 art b un year   , xlabel ylabel c(lll) title(`gr7ttl')  b2("`gr7subtitle'") `svngn' l1title("`unitlblstk' of people") gap(3)
			if `"`saving'"'~="" {
				di as txt "(Version 7 graph file saved with option: " as result `"`svngn'"'
			}
			pause
			gr7 p d i year   , xlabel rscale(`rsmin' `rsmax') c(lll) `transition' gap(3) title(`gr7ttl')  b2("`transttl'") `svngt' r2("`incidlbl'") ylabel rlabel
			if `"`saving'"'~="" {
				di as txt "(Version 7 graph file saved with option: " as result `"`svngt'"'
			}
			pause
			gr7 costtot cost2nd fundgap year   , xlabel ylabel c(lll) title(`gr7ttl')  b2("`gr7subtitle'") `svngc' l1title("`unitlblc' of `vminus1' US dollars") gap(3)
			if `"`saving'"'~="" {
				di as txt "(Version 7 graph file saved with option: " as result `"`svngc'"'
			}
			pause
			format  %9.0f  titlepct discpct cvrgpct
			gr7 cvrgpct titlepct discpct year if year > 2007 , xlabel ylabel c(lll) title(`gr7ttl')  b2("`gr7subtitle'") `svngp' l1title("Percentage") gap(3)
			if `"`saving'"'~="" {
				di as txt "(Version 7 graph file saved with option: " as result `"`svngp'"'
			}
			pause
			if "`weo'" ~= "" {
				format  %9.0f  artpctpub artpcthe artpctpop
				gr7 artpctpub artpcthe artpctpop year if year > 2007 , xlabel ylabel c(lll) title(`gr7ttl')  b2("`gr7subtitle'") `svnge' l1title("Percentage") gap(3)
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as result `"`svnge'"'
				}
				pause
			}
		}
		else {
			tempfile nmbrs trnstn costs disc
			local gphset=c(graphics)
			if "`gphset'" == "on" {
				set graph off
			}
			set textsize 150
			gr7 art b un year , xlabel ylabel c(lll) gap(3)  saving(`"`nmbrs'"')  l1title("`unitlblstk' of people")
			gr7 p d i year    , xlabel rscale(`rsmin' `rsmax') c(lll) `transition' gap(3)  saving(`"`trnstn'"') r2("`incidlbl'") ylabel rlabel
			gr7 costtot cost2nd fundgap year   , xlabel ylabel c(lll) gap(3)  saving(`"`costs'"') l1title("`unitlblc' of `vminus1' US dollars")
			format  %9.0f  titlepct discpct cvrgpct
			gr7 cvrgpct titlepct discpct year       if year > 2007  , xlabel ylabel c(lll) gap(3) saving(`"`disc'"') l1title("Percentage")
			if "`weo'" ~= "" {
				tempfile weogr7
				format  %9.0f  artpctpub artpcthe artpctpop
				gr7 artpctpub artpcthe artpctpop year if year > 2007  , xlabel ylabel c(lll) gap(3) saving(`"`weogr7'"') l1title("Percentage")
			}
			if "`gphset'" == "on" {
				set graph on
			}
			set textsize 100
			if "`weo'" == "" {
				gr7 using `nmbrs' `trnstn' `costs' `disc'         , margin(5) `svng' title(`gr7ttl') b2("`gr7subtitle'")
			}
			else {  // Replace graph on discretionary with weogr7 so as to display only four graphs
				gr7 using `nmbrs' `trnstn' `costs'        `weogr7', margin(5) `svng' title(`gr7ttl') b2("`gr7subtitle'")
			}
			if `"`saving'"'~="" {
				di as txt "(Version 7 graph file saved with option: " as result `"`svng'"'
			}
			pause
		}
	}
	if "`graph'"~="" {  // Use modern graphs, somewhat slower, but elegant
		if "`transition'"~="" {
			local transition ": AIDS transition in `tyear'"
			local trmtxt  `"xline(`tyear',  lwidth(medthick) lpattern(dot)) xmtick(`tyear') xmlabel(`tyear')  xlabel(#5) text(`rsmax'  `tyear' "AIDS" "Transition")"'
		}
		if "`grcombine'"=="" {  // Create, display and save one graph at a time
			twoway (line art year, sort lwidth(thick) lpattern(solid) ) /*
				*/ (line b year  ,      lwidth(thick) lpattern(dash) ) /*
				*/ (line un year ,      lwidth(thick) lpattern(dot) ) , /*
			  */ ytitle("`unitlblstk' of people") ylabel(, angle(horizontal)) /*
			  */ title("`titlen'") /*
			  */ subtitle("`subtitle'") /*
			  */ xtitle("Enrolling $glblsigma % of unmet need each year") /*
			  */ legend(order(1 "Total number on ART" 2 "Number on 2nd Line" 3 "Unmet need for ART") `legend') /*
			  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
					/* span */ size(vsmall)) `svngn' `namen'
			pause

			twoway  (line i  year ,      lwidth(thick) lpattern(dash) yaxis(1) )/*
				*/  (line d  year , sort  lwidth(thick) lpattern(solid) yaxis(1) ) /*
				*/  (line p  year ,      lwidth(thick) lpattern(dot) yaxis(2) ) , /*
			  */ `trmtxt'   /*
			  */ ytitle("New infections or deaths" "(`unitlblflw' )" "  " , axis(1)) /*
			  */ ytitle("Number living with HIV/AIDS" "(`unitlblstk')" "  " , axis(2)) /*
			  */ ylabel(, angle(horizontal) format(%9.1f) axis(1))  ylabel(, angle(horizontal) axis(2)) /*
			  */ title("`titlet'") /*
			  */ subtitle("`subtitle'`transition'") /*
			  */ xtitle("Enrolling $glblsigma % of unmet need each year") /*
			  */ legend(order(1 "New HIV infections" 2 "Number AIDS deaths" 3 "Number living with HIV/AIDS" ) `legendt') /*
			  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
					/* span */ size(vsmall)) `svngt' `namet'
			pause

			twoway (line costtot  year , sort lwidth(thick) lpattern(solid) ) /*
				*/ (line cost2nd  year ,      lwidth(thick) lpattern(dash) ) /*
				*/ (line fundgap  year ,      lwidth(thick) lpattern(dot) ) , /*
			  */ ytitle("`unitlblc' of dollars") ylabel(, angle(horizontal)) /*
			  */ title("`titlec'") /*
			  */ subtitle("`subtitle'") /*
			  */ xtitle("Enrolling $glblsigma % of unmet need each year") /*
			  */ legend(order(1 "Total cost of ART" 2 "Cost of 2nd Line" 3 "ART funding gap") `legend') /*
			  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
					/* span */ size(vsmall)) `svngc' `namec'
			pause

			twoway (line cvrgpct  year             ,  sort lwidth(thick) lpattern(dot)   ) /*
				*/ (line titlepct  year if year>2007, sort lwidth(thick) lpattern(solid) ) /*
				*/ (line discpct   year if year>2007,      lwidth(thick) lpattern(dash)  ),/*
			  */ ytitle("Percentages") ylabel(, angle(horizontal)) /*
			  */ title("`titlep'") /*
			  */ subtitle("`subtitle'") /*
			  */ xtitle("Enrolling $glblsigma % of unmet need each year") /*
			  */ legend(order(1 "ART coverage" 2 "Entitlement percent" 3 "Discretionary percent" ) `legendp') /*
			  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
					/* span */ size(vsmall)) `svngp' `namep'
			pause

			if "`weo'" ~= "" { /* artpcthe artpctpub artpctprv artpctoop artpctpop */
				twoway (line artpctpub  year  , sort lwidth(thick) lpattern(dot)   ) /*
				  */ (line artpcthe   year  , sort lwidth(thick) lpattern(solid) ) /*
				  */ (line artpctpop   year ,      lwidth(thick) lpattern(dash)  ),/*
				  */ ytitle("Percentages") ylabel(, angle(horizontal) format(%9.0f)) /*
				  */ title("`titlee'") /*
				  */ subtitle("`subtitle'") /*
				  */ xtitle("Enrolling $glblsigma % of unmet need each year") /*
				  */ legend(order(1 "ART % Pub HE" 2 "ART % Total HE" 3 "Reference: ART patients as % population" ) `legendp') /*
				  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
					  /* span */ size(vsmall)) `svnge' `namee'
				pause
			}
		}
		else {    // Create, display and save four graphs as panels in a single chart
			local gphset=c(graphics)
			if "`gphset'" == "on" {
				set graph off
			}

			local panel_lbla subtitle("Panel a. Numbers of people")
			local panel_lblb subtitle("Panel b. Infections and deaths")
			local panel_lblc subtitle("Panel c. Expenditures")
			local panel_lbld subtitle("Panel d. Affordability")
			local panel_lble subtitle("Panel d. Coverage")

			twoway (line art year, sort lwidth(thick) lpattern(solid) ) /*
			  */ (line b year  ,      lwidth(thick) lpattern(dash) ) /*
			  */ (line un year ,      lwidth(thick) lpattern(dot) ) , /*
			  */ ytitle("`unitlblstk' of people") ylabel(, angle(horizontal)) /*
			  */ xtitle("") xlabel(,grid) scale(0.80) legend(span) /*
			  */ legend(order(1 "Total number on ART" 2 "Number on 2nd Line" 3 "Unmet need for ART") `legendp') `panel_lbla'

			tempname nmbrs
			qui graph rename `nmbrs', replace

			twoway  (line i  year ,      lwidth(thick) lpattern(dash) yaxis(1) )/*
			  */  (line d  year , sort lwidth(thick) lpattern(solid) yaxis(1) ) /*
			  */  (line p  year ,      lwidth(thick) lpattern(dot) yaxis(2) ) , /*
			  */ `trmtxt'   /*
			  */ ytitle("New infections" "or deaths (`unitlblflw' )" "  " , axis(1)) /*
			  */ ytitle("Number living with" "HIV/AIDS (`unitlblstk')" "  " , axis(2)) /*
			  */ ylabel(, angle(horizontal)  format(%9.1f) axis(1))  ylabel(, angle(horizontal) axis(2)) /*
			  */ xtitle("") xlabel(,grid) scale(0.80) legend(span) /*
			  */ legend(order(1 "New HIV infections" 2 "Number AIDS deaths" 3 "Number living with HIV/AIDS" ) `legendp') `panel_lblb'

			tempname trnstn
			qui graph rename `trnstn', replace

			twoway (line costtot  year , sort lwidth(thick) lpattern(solid) ) /*
			  */ (line cost2nd  year ,      lwidth(thick) lpattern(dash) ) /*
			  */ (line fundgap  year ,      lwidth(thick) lpattern(dot) ) , /*
			  */ ytitle("`unitlblc' of `vminus1' US dollars") ylabel(, angle(horizontal)) /*
			  */ xtitle("") xlabel(,grid) scale(0.80) legend(span) /*
			  */ legend(order(1 "Total cost of ART" 2 "Cost of 2nd Line" 3 "ART funding gap") `legendp')  `panel_lblc'

			tempname costs
			qui graph rename `costs', replace

			twoway (line cvrgpct  year              ,  sort lwidth(thick) lpattern(dot)   ) /*
			  */ (line titlepct  year if year>2007, sort lwidth(thick) lpattern(solid) ) /*
			  */ (line discpct   year if year>2007,      lwidth(thick) lpattern(dash)  ),/*
			  */ ytitle("Percentages") ylabel(, angle(horizontal)) /*
			  */ xtitle("") xlabel(,grid) scale(0.80) legend(span) /*
			  */ legend(order(1 "ART coverage" 2 "Entitlement percent" 3 "Discretionary percent" ) `legendp') `panel_lble'

			tempname pcnt
			qui graph rename `pcnt', replace

			if "`weo'" ~= "" { /* artpcthe artpctpub artpctprv artpctoop artpctpop */
				twoway (line artpctpub  year             ,  sort lwidth(thick) lpattern(dot)   ) /*
				  */ (line artpcthe   year , sort lwidth(thick) lpattern(solid) ) /*
				  */ (line artpctpop   year ,      lwidth(thick) lpattern(dash)  ),/*
				  */ ytitle("Percentages") ylabel(, angle(horizontal) format(%9.0f)) /*
				  */ xtitle("") xlabel(,grid) scale(0.80) legend(span) /*
				  */ legend(order(1 "ART % Pub HE" 2 "ART % Total HE" 3 "Reference: ART patients as % population" ) `legendp')  `panel_lbld'

				tempname afford
				qui graph rename `afford', replace
			}
			if "`gphset'" == "on" {
				set graph on
			}
			if "`weo'" == "" {
				gr combine `nmbrs' `trnstn' `costs' `pcnt'         , rows(2) xcommon /*
				  */ title("AIDS treatment: Patients, mortality and costs") /*
				  */ subtitle("`subtitle'`transition'") /*
				  */ b1title("Enrolling $glblsigma % of unmet need each year") /*
				  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
						/* span */ size(vsmall)) `svng' `name' imargin(small)
			}
			else { // Replace graph on discretionary with -afford- so as to display only four graphs
				gr combine `nmbrs' `trnstn' `costs'        `afford', rows(2) xcommon /*
				  */ title("AIDS treatment: Patients, mortality and costs") /*
				  */ subtitle("`subtitle'`transition'") /*
				  */ b1title("Enrolling $glblsigma % of unmet need each year") /*
				  */ note(`"`grnote1'"' `"`grnote2'"' `"`grnote3'"' `"`grnote4'"' `"`grnote5'"' `"`grnote6'"',  ///
						/* span */ size(vsmall)) `svng' `name' imargin(small)
				 }
			pause
		}
	}

	restore
end		/* End of program makegph */

program define relabel
*	Variables that are either a stock or a flow receive a "characteristic"
*   Variables are assigned a value of the characteristic [morebadgood]
*       When [morebadgood] is set to "moreisbad" ("moreisgood"), an increase in the variable
*       harms (improves) human welfare in the analyzed country.
*       For some variables like HIV prevalence (p), an increase is
*       both "bad" (because the country has more sick people to eventually need treatment)
*       and "good" (because people who might otherwise have died are alive).
*       The [morebadgood] characteristic is assigned the value "both" for these variables.
	version 9.2
	syntax [, sigma(varname)]
	lab var p      "Number living with HIV/AIDS"
		char def p[stockflow] stock
		char def p[morebadgood] both
	lab var i      "Number of new HIV infections"
		char def i[stockflow] flow
		char def i[morebadgood] moreisbad
	lab var h      "HIV+, not eligible for ART"  /* i.e. "Healthy" HIV+ */
		char def h[stockflow] stock
		char def h[morebadgood] both
	lab var n      "Number adults needing ART"
		char def n[stockflow] stock
		char def n[morebadgood] both
	lab var un     "Unmet need for ART"
		char def un[stockflow] stock
		char def un[morebadgood] moreisbad
	lab var art    "Total patients on ART"
		char def art[stockflow] stock
		char def art[morebadgood] moreisgood
	lab var b      "Number adults receiving 2nd line ART"
		char def b[stockflow] stock
		char def b[morebadgood] moreisgood
	lab var d      "Number deaths"
		char def d[stockflow] flow
		char def d[morebadgood] moreisbad
	lab var a2     "Number adults receiving 1st line ART, Years 2+"
		char def a2[stockflow] stock
		char def a2[morebadgood] moreisgood
	lab var a1     "Number adults receiving 1st line ART, Year 1"
		char def a1[stockflow] stock
		char def a1[morebadgood] moreisgood
	lab var fundgap "Funding gap for ART"
		char define fundgap[scale] 1000
		char def fundgap[stockflow] flow
		char def fundgap[morebadgood] moreisbad
	lab var costtot "Total cost of ART ($1,000's)"
		char define costtot[scale] 1000
		char def costtot[stockflow] flow
		char def costtot[morebadgood] moreisbad
	lab var cost2nd "Cost of 2nd line ART ($1,000's)"
		char define cost2nd[scale] 1000
		char def cost2nd[stockflow] flow
		char def cost2nd[morebadgood] moreisbad
	lab var disctot "Discretionary ($1,000's)"
		char define disctot[scale] 1000
		char def disctot[stockflow] flow
		char def disctot[morebadgood] moreisgood
	lab var discpct "Discretionary as % of costs"
	lab var titletot "Entitlement, Total ($1,000's)"
		char define titletot[scale] 1000
		char def titletot[stockflow] flow
		char def titletot[morebadgood] both
	lab var title1st "Entitlement, 1st line ($1,000's)"
		char define title1st[scale] 1000
		char def title1st[stockflow] flow
		char def title1st[morebadgood] both
	lab var titlepct "Entitlement as % of costs"
	lab var cvrgpct "ART patients as % of need"

	capture confirm numeric variable artpcthe artpctpub artpctprv artpctoop artpctpop
	if _rc==0  {
		lab var artpcthe "ART as % total health spending"
		lab var artpctpub "ART as % public health spending"
		lab var artpctprv "ART as % private health spending"
		lab var artpctoop "ART as % OOP health spending"
		lab var artpctpop "ART patients as % population"
	}

	if "`sigma'"~="" {
		lab var `sigma' "Country specific adjustment parameter"
	}
end    /* End of program relabel */

program define parmsbycd4, rclass
	syntax , vintage(real) [adrate1(string) adrate2(string) bdrate(string) erate(string) ndrate(string) cd4(string) ///
		takeoff(string) horizon(string) strtcov2(string) trgtcov2(string) trgtyr(string) strtyr(string) incmult(string) ///
		iexog ir0(string)  ///
		mincp(string) maxcp(string) cstrtyr(string) ctrgtyr(string) csp(string)   ///  Circumcision
		minep(string) maxep(string) pstrtyr(string) ptrgtyr(string)   ///  Prevention "effort"
		maxvp(string) vstrtyr(string) vtrgtyr(string) ve(string)  ///Vaccination
		fp(string) gp(string) maxdt(string) maxdu(string) ]

    foreach pname in adrate1  adrate2  bdrate  erate  ndrate strtcov2 trgtcov2 {
        // Check to make sure user has not entered a string or an out of bounds value into a parameter
        if "``pname''"~="" & (real("``pname''")>=. | real("``pname''")>1 | real("``pname''")<0) {
            noi di _n as error "Parameter `pname' must be numeric and between zero and one.
            noi di    as error "Instead it is: ``pname''"
            error 198
        }
    }
    foreach pname in strtyr takeoff horizon trgtyr ctrgtyr ptrgtyr {
        // Check to make sure user has not entered a year that is out of bounds
		if "`pname'"=="strtyr" {
*			local lowbnd = `vintage' - 1  // Lower bound
			local lowbnd = 2007			// Lower bound
			local uprbnd = `vintage' + 5  // Upper bound
		}
		if "`pname'"=="takeoff" {
			local lowbnd = `vintage' - 1  // Lower bound
			local uprbnd = `vintage' + 5  // Upper bound
		}
		if "`pname'"=="horizon" | "`pname'"=="trgtyr" | "`pname'"=="ctrgtyr" | "`pname'"=="prgtyr" {
			local lowbnd = 2010  // Lower bound
			local uprbnd = 2100  // Upper bound
		}
        if "``pname''"~="" & (real("``pname''")>=. | real("``pname''")>`uprbnd' | real("``pname''")<`lowbnd') {
            noi di _n as error "Parameter " as res "`pname'" as err " must be a valid year between "  ///
					as res "`lowbnd'" as err " and " as res "`uprbnd'."
            di    as error "Instead it is: " as res "``pname''"
            error 198
        }
    }

    foreach pname in incmult ir0 minep maxep mincp maxcp maxvp ve fp gp {
		if `"``pname''"'~="" {
			if ``pname'' < . & (``pname'' > 1.0 | ``pname'' < 0 ) {
				noi di _n as err "The parameter " as res "`pname'" as err " must be between zero and 1.0, inclusive."  ///
					_n "Here it is set to : " as res `"``pname''"'
				exit 198
			}
		}
	}
    foreach pname in maxdt maxdu {
		if `"``pname''"'~="" {
			if ``pname'' < . & (``pname'' > 0.5 | ``pname'' < -0.5 ) {
				noi di _n as err "The parameter " as res "`pname'" as err " must be between -0.5 and 0.5, inclusive."  ///
					_n "Here it is set to : " as res `"``pname''"'
				exit 198
			}
		}
	}

	    // Set default values of incidence parameters
        local incmult_dflt = 0.97  //  Year t incidence as fraction of year t-1 (For exogenous incidence model only)
        local ir0_dflt     =  0.2  //  Incidence rate in the population at most risk for sexual tranmission
        local minEp_dflt   = 0     //  Minimum HIV prevention policy Effort year<=pstrtyr: 0 == No effort, 1.0 = Effort sufficient to eliminate incidence
        local maxEp_dflt   = 0     //  Maximum HIV prevention policy Effort year<=pstrtyr: 0 == No effort, 1.0 = Effort sufficient to eliminate incidence
        local pstrtyr_dflt = 2020  //  Start year of prevention efforts (In version 5.30, separate pstrtyr from cstrtyr)
        local ptrgtyr_dflt = 2030  //  Target year of prevention efforts (In version 5.30, separate ptrgtyr from ctrgtyr)
        local minCp_dflt   = .1    //  Minimum proportion circumcised year<=cstrtyr
        local maxCp_dflt   = .1    //  Maximum proportion circumcised year>=ctrgtyr (In version 5.30, this was .8)
        local cstrtyr_dflt = 2020  //  Start year of circumcision program  (In version 5.30, separate pstrtyr from cstrtyr)
        local ctrgtyr_dflt = 2030  //  Target year of circumcision program  (In version 5.30, separate ptrgtyr from ctrgtyr)
        local CSp_dflt     = .6    //  Proportion susceptible of circumcised to uncircumcised male
        local maxVp_dflt   =  0    //  Maximum proportion of population to be vaccinated year>=vtrgtyr: 0 < maxVp < 1
        local ve_dflt      =  .6   //  Degree of protection provided by vaccine: ve = 1 is a perfect vaccine
        local vstrtyr_dflt = 2025  //  Start year of vaccination scale up
        local vtrgtyr_dflt = 2030  //  Target year of vaccination scale-up
        local Fp_dflt      =  0.7  //  Proportion of transmission after primary infection: (i.e. after the first 2.2 months)
        local Gp_dflt      =  0.8  //  Infectiousness of person on treatment as a proportion of infectiousness of untreated person
        local maxDt_dflt   =  0    //  Maximum disinhibition effect of ART coverage on treated: maxDt < 1
        local maxDu_dflt   =  0    //  Maximum disinhibition effect of ART coverage on UNtreated:
                                   //           maxDu < infinity but effect is bounded to half of (1-k)


    if "`iexog'"~="" {  // If incidence is to be modelled as exogenous
        // Check to make sure user has not specified an implausible precentage change for incidence
        if "`incmult'"~="" & (real("`incmult'")>=. | real("`incmult'")>1.5 | real("`incmult'")<=0) {
            di _n as error "Parameter incmult must be numeric and between zero and 1.5.
            di    as error "Instead it is: `incmult'"
            error 198
        }
        di _n
        di as txt "Incidence is modeled as exogenous to ART"
        di as txt "  Ratio that HIV incidence is of last year's incidence" _col(60) "incmult" as result _col(72) `incmult_dflt'  as res _col(82) `incmult'
    }
    else {              // If incidence is to be modelled as endogenous
        if "`incmult'"~="" {
            di _n as error "Parameter -incmult- cannot be specified unless option -iexog- has also been specified."
            error 198
        }
        di _n
        di as txt "Incidence is modeled as endogenous to ART and affected by prevention scale-up:"
        di as txt "	The following list gives in each row the name of an"
        di as txt "	-aidsproj- incidence parameter, its default value and"
        di as txt "	the user-chosen value, if there is one, which overrides"
        di as txt "	that default.  In each row, only the right-most value"
        di as txt "	is used in the subsequent projections." _n
        di as txt "    Incidence rate in core group, Year 0"           _col(60) "ir0"      as result _col(72) `ir0_dflt'        as res _col(82) `ir0'
        di as txt `"    Start year of prevention "effort""'            _col(60) "pstrtyr"  as result _col(72) `pstrtyr_dflt'    as res _col(82) `pstrtyr'
        di as txt `"    Target year of maximum prevention "effort""'   _col(60) "ptrgtyr"  as result _col(72) `ptrgtyr_dflt'    as res _col(82) `ptrgtyr'
        di as txt "    Minimum HIV prevention effort before pstrtyr"   _col(60) "minep"    as result _col(72) `minEp_dflt'      as res _col(82) `minep'
        di as txt "    Maximum HIV prevention effort after ptrgtyr"    _col(60) "maxep"    as result _col(72) `maxEp_dflt'      as res _col(82) `maxep'
        di as txt "  Circumcision"
        di as txt "    Start year of circumcision program"             _col(60) "cstrtyr"  as result _col(72) `cstrtyr_dflt'    as res _col(82) `cstrtyr'
        di as txt "    Target year of maximum circumcision coverage"   _col(60) "ctrgtyr"  as result _col(72) `ctrgtyr_dflt'    as res _col(82) `ctrgtyr'
        di as txt "    Minimum proportion circumcised before cstrtyr"  _col(60) "mincp"    as result _col(72) `minCp_dflt'      as res _col(82) `mincp'
        di as txt "    Maximum proportion circumcised after cstrtyr"   _col(60) "maxcp"    as result _col(72) `maxCp_dflt'      as res _col(82) `maxcp'
        di as txt "  Vaccination"
        di as txt "    Start year of vaccination scale-up"             _col(60) "vstrtyr"  as result _col(72) `vstrtyr_dflt'    as res _col(82) `vstrtyr'
        di as txt "    Target year of maximum vaccination scale-up"    _col(60) "vtrgtyr"  as result _col(72) `vtrgtyr_dflt'    as res _col(82) `vtrgtyr'
        di as txt "    Maximum proportion vaccinated after vtrgtyr"    _col(60) "maxvp"    as result _col(72) `maxVp_dflt'      as res _col(82) `maxvp'
        di as txt "    Protection of vaccine: ve=1 perfect"            _col(60) "ve"       as result _col(72) `ve_dflt'         as res _col(82) `ve'
        di as txt "  Treatment effects"
        di as txt "    Fraction transmission after primary infection"  _col(60) "fp"       as result _col(72) `Fp_dflt'         as res _col(82) `fp'
        di as txt "    Protection of treatment: gp=1 perfect"          _col(60) "gp"       as result _col(72) `Gp_dflt'         as res _col(82) `gp'
        di as txt "    Max disinhibition of ART on treated < 1"        _col(60) "maxdt"    as result _col(72) `maxDt_dflt'      as res _col(82) `maxdt'
        di as txt "    Max disinhibition of ART on untreated < inf"    _col(60) "maxdu"    as result _col(72) `maxDu_dflt'      as res _col(82) `maxdu'
    }

	    // Set default values of projection parameters
        local takeoff_dflt  = real("`vintage'")
        local horizon_dflt  = 2030
        local strtcov2_dflt = -9999  // Flags that user left strtcov2 blank
        local trgtcov2_dflt = .50
        local strtyr_dflt   = real("`vintage'") - 1
        local trgtyr_dflt   = 2050

    di _n
    di as txt "Projection parameters:"
    di as txt "  First year of projections (projection takeoff)" _col(60) "takeoff"  as result _col(72) `takeoff_dflt'    as res _col(82) `takeoff'
    di as txt "  Last year of projections (projection horizon)"  _col(60) "horizon"  as result _col(72) `horizon_dflt'    as res _col(82) `horizon'
    di as txt "  2nd line ART to start in year"                  _col(60) "strtyr"   as result _col(72) `strtyr_dflt'     as res _col(82) `strtyr'
    di as txt "  2nd line ART to reach target in year"           _col(60) "trgtyr"   as result _col(72) `trgtyr_dflt'     as res _col(82) `trgtyr'
    di as txt "  Starting coverage rate for 2nd line ART"        _col(60) "strtcov2" as result _col(72) "regional"        as res _col(82) `strtcov2'
    di as txt "  Target coverage rate for 2nd line ART"          _col(60) "trgtcov2" as result _col(72) `trgtcov2_dflt'   as res _col(82) `trgtcov2'

    if "`cd4'"=="revert" {
        di _n as txt "The CD4 parameter has been set to: "  as res _col(80) "`cd4'"
        di    as txt "	With this selection, the user requests that -aidsproj-"
        di    as txt "	revert to its version 2 behavior with the following"
        di    as txt "	default values of the parameters.  Any user-chosen"
        di    as txt "	values override those defaults.  In each row, only "
        di    as txt "	the right-most value is used in the subsequent"
        di    as txt "	projections." _n

	    // Revert default values to version 2 levels
        local adrate1_dflt = 0.13
        local adrate2_dflt = 0.04   // Equivalent to a life-expectancy of 25 years from initiation
        local bdrate_dflt  = 0.01
        local erate_dflt   = 0.11
        local ndrate_dflt  = 0.325  // The value of ndrate at CD4=130
        local cd4 = .
    }
    else {
       local dflt_cd4 = 130
       if "`cd4'"=="" {
			local usr_cd4
			local cd4 = `dflt_cd4'
        }
		else {
			local usr_cd4 = `cd4'
			local cd4 = `usr_cd4'
		}
        if real("`cd4'")>=. | real("`cd4'")<50 | real("`cd4'")>800 {
            di _n as error "Parameter cd4 must be numeric and between 50 and 800.
            di    as error "Instead it is: " as res "`cd4'"
            error 198
        }
        di _n as txt "The median CD4 count at ART initiation is: " _col(60) "cd4" as result _col(72) "`dflt_cd4'"  as res _col(82) `usr_cd4'
        di    as txt "	The following list gives in each row the name of an"
        di    as txt "	-aidsproj- treatment parameter, its default value for "
        di    as txt "	this choice of the CD4 initiation parameter and "
        di    as txt "	the user-chosen value, if there is one, which overrides"
        di    as txt "	that default.  In each row, only the right-most value"
        di    as txt "	is used in the subsequent projections." _n

    	// Set default values as functions of user-selected CD4 at initiation

        // Estimated hazard ratio from the WTSC April 2009 paper is:
        local lnhrhat = cond(`cd4' < 450, 11.21569 + 0.0208304*(1000 - `cd4') - 0.9640528*(1000 - `cd4')^0.5, 0.0633927 - 0.0001409*(`cd4' - 450))
            // The line that best fits the data for cd4<450 reaches a minimum about there and then implausibly rises for higher cd4.
            // Therefore we splice that line with a linear trend downward which reaches zero at cd4=900.

        // New version of ndrate is estimated from the ART-LINC data, excluding the Thai cohorts
        // see the file N:\Health\AIDS Research\AIDS_Proj_Model\Garnett\Copenhagen_Consensus\Smpl_CE\HIV_Progression_Data_and_Models_Revised_for_CC_CE_paper.xlsx
        // see the sheet "ndrate estimation"
        local ndrate_dflt  =  round(exp(-0.2241 - 0.00644*`cd4') , .0001)

        // Version 4.18 used the following fitted equation based on ART-LINC including the Thai cohorts:
        // see the files N:\Health\AIDS Research\Book\Analysis\hazrat.do & .log
        // local ndrate_dflt  =  round(exp(-0.427349 - 0.0053641*`cd4') , .0001)
        // Previous version of ndrate_dflt was:
        // local ndrate_dflt  =  round(1.41*(0.993^(`cd4')), .001)

        // Estimate of adrate1 is a splice at CD4=250 of the meta-analysis-derived result
        // and the ART-LINC-derived ndrate deflated by the WTSC hazard ratio
        // The result of the meta-analyisis is given in the file:
        //      "N:\Health\AIDS Research\Book\Analysis\Mistata\mistata5_smr_60collapse.gph"
        // and is graphed in "N:\Health\AIDS Research\Book\Analysis\pepfar_blog.dfw"
        // Smear coefficient for this estimate is 1.31.  Assumes 60% one-year mortality among LTFU
        local adrate1_dflt = round( ///
            cond(`cd4' > 250, `ndrate_dflt'/exp(`lnhrhat'), 1.31*exp(-0.85 - 0.296*ln(`cd4')) ///
            ) , .0001)
        // Previous version of adrate1_dflt is here:
    	// local adrate1_dflt =  round(1.31*exp(-0.85-0.296*(ln(`cd4'))), .001)

		// For version 4.19, allow adrate2 to vary as a function of initating cd4 according to:
        //		N:\Health\AIDS Research\AIDS_Proj_Model\Garnett\Copenhagen_Consensus\life_expect.do
        // Since -aidsproj- does not distinguish ages, set the age range at 25-29 for the typical patient.
        local adrate2_dflt = round( ///
            cond(`cd4' > 317.5, 1/36.8, 1/(-10.97143 - 0.256862 * 25 - 0.2529094 * `cd4' + 7.551944 * sqrt(`cd4')) ///
            ) , .0001)

*        local adrate2_dflt = .04  // Value of adrate2 was a fixed constant in version 4.18

        local bdrate_dflt  = .08
*        local bdrate_dflt  = .04   // Previous value of bdrate was equal to the constant value of adrate2
*        local erate_dflt   =  round(exp(-7.042134334 + 0.960788107*ln(`cd4')) , .0001)  // From Collins and eART-LINC (Version 3.5)
*        local erate_dflt   =  round(exp(0.67*(`cd4')/1000) , .0001)  // Approximation of above relationships, from Collins and eART-LINC (Version 4)
        // Note implication that with a CD4 recruitment threshold of 900, 60% of HIV+ would become eligible each year.
		// Previous version:      local erate_dflt   =  round(1/(11-0.01*(`cd4')),.001)
		// Version 4.13 semilog from N:\Health\AIDS Research\Book\Analysis\HIV_Progression_Data_and_Models_Revised_erate.xlsx   2/17/11
        local erate_dflt   =  round(exp(-2.62192 + 0.003269515*`cd4') , .0001)  // From Collins and eART-LINC (Version 4.13, 2/17/11)
    }

    di as txt "  Proportion of HIV+ newly eligible for ART"           _col(60) "erate"   as result _col(72) `erate_dflt'    as res _col(82) `erate'
    di as txt "  ART Death Rate during first year on 1st line"        _col(60) "adrate1" as result _col(72) `adrate1_dflt'  as res _col(82) `adrate1'
    di as txt "  ART Death Rate during subsequent years  on 1st line" _col(60) "adrate2" as result _col(72) `adrate2_dflt'  as res _col(82) `adrate2'
    di as txt "  ART Death Rate of AIDS patients on 2nd line"         _col(60) "bdrate"  as result _col(72) `bdrate_dflt'   as res _col(82) `bdrate'
    di as txt "  Non-ART Death Rate of AIDS patients"                 _col(60) "ndrate"  as result _col(72) `ndrate_dflt'   as res _col(82) `ndrate'

    // set the parameters to default values or override them with user specified values
    return scalar cd4       = round(`cd4',0.1)
    foreach pname in takeoff horizon strtyr trgtyr strtcov2 trgtcov2 cstrtyr ctrgtyr pstrtyr ptrgtyr vstrtyr vtrgtyr {
        local `pname' = cond("``pname''"=="",``pname'_dflt',real("``pname''"))
        return scalar `pname' = round(``pname'')
    }
    foreach pname in minCp maxCp CSp minEp maxEp maxVp ve Fp Gp maxDt maxDu ir0 {
        local lpname = lower("`pname'")
        local `lpname' = cond("``lpname''"=="",``pname'_dflt',real("``lpname''"))
        return scalar `lpname' = round(``lpname'', 0.01)
    }
    foreach pname in strtcov2 trgtcov2 incmult {
        local `pname' = cond("``pname''"=="",``pname'_dflt',real("``pname''"))
        return scalar `pname' = round(``pname'', 0.001 )
    }
    foreach pname in adrate1  adrate2  bdrate  erate  ndrate  {
        local `pname' = cond("``pname''"=="",``pname'_dflt',real("``pname''"))
        return scalar `pname' = round(``pname'', 0.0001 )
    }
end    /* End of program parmsbycd4 */

prog define aidspost
*	Called four ways.
*	1. With option -postopen-, constructs lists of wide variable names for the postfile & opens file
*	2. With option -harvest-, it collects desired summary statistics from the return macros of a tabstat command
*		and assigns each one to a scalar with the appropriate name
*	3. With the -postnow- option, the harvested values are posted in the order of the names
*	4. With the -postclose- option, the program drops the scalars and closes the postfile
	version 9.2
	syntax varlist , [ postyears(numlist) posttotal postrunspecs(string) postfile(string) ///
		countries(string) regions(string) world postopen harvest postclose postnow postname(string)]
	if `"`postname'"'=="" {
		local postname __postname
	}

	if "`postyears'"=="" | "`posttotal'"~="" {
		local posttotal Total
	}

	if "`countries'`regions'`world'" =="" {
		di _n as err "Postfile options require that output be requested either by specifying "
		di    as err "the -world- option or by choosing at least one region or country."
		exit 198
	}

	if "`world'"~="" {
		local world ALL
	}

	if "`postopen'"~="" {
		global svarlist
		global wvarlist
		foreach var of local varlist {
			foreach nm in `countries' `regions' `world' {
				foreach year in `postyears' `posttotal' {
					local strlen1 : length global svarlist
					local strlen2 = length(`"(s`var'`nm'`year')"')
					if `strlen1'+`strlen2' > 240  {
						di _n as err "Error in postfile options"
						di    as err "User is attempting to harvest too many values from each run of -aidsproj-."
						di    as err "The product of the number of requested years in the -postyears- option"
						di    as err "times the number of requested aggregates with the -countries-, -regions- and -world- options"
						di    as err "times one plus the number of requested variables in the -postvars- option"
						di    as err "should not exceed 12."
						exit 198
					}
				global wvarlist = "${wvarlist} " + "`var'`nm'`year'"
				global svarlist = "${svarlist} " + "(s`var'`nm'`year')"
				}
			}
		}
		di _n as txt "postfile " as res "`postname' runspecs $wvarlist " as txt " using " as res "`postfile' " as txt ", replace"
		capture postfile `postname' str200 runspecs double(${wvarlist}) using `postfile' , replace
		if _rc~=0 {
		    postclose `postname'
		    postfile `postname' str200 runspecs double(${wvarlist}) using `postfile' , replace
		}
	}

	else if "`harvest'"~="" {
*	Harvest the desired tabstat result(s) & place them in approriately named scalars
		foreach var of local varlist {
			global `var'unit : char `var'[scale]
			local thiscoln
			local colnms : colname(r(Stat1))
			local thiscoln : list posof "`var'" in colnms
			if `thiscoln'~=0 {
				foreach nm in `countries' `regions' `world' {
					foreach year in `postyears' `posttotal' {
						scalar s`var'`nm'`year' =  .
						if "`year'"~="Total" {
							local  i = 1
							local yr = 0
							while "`r(name`i')'" ~= "" {
								if r(name`i')=="`year'" {
									local yr = `i'
								}
								local i = `i' + 1
							}
							if `yr'==0 {
								di _n as err "Option -postyears- specifies a year outside the projection window, "
								di _n as err "which is between the -takeoff- and the -horizon- years."
								error 198
							}
						}
						else {
							local yr Total
						}
						scalar s`var'`nm'`year' =  el(r(Stat`yr'),1,`thiscoln')
*						di _n as res "The scalar "  as res  "s`var'`nm'`year'"  as res " is set to "  as res s`var'`nm'`year'
					}
				}
			}
		}
	}

	else if "`postnow'"~="" {
		local postrunspecs = trim(itrim(`"`postrunspecs'"'))
		post `postname' (`"`postrunspecs'"') ${svarlist}
	}
	else if "`postclose'"~="" {
		local svarlist = subinstr(subinstr("${svarlist}", "(", " ",.),")"," ",.)
		di as txt "Drop the following scalars here:"
		di as res "    `svarlist' "
		scalar drop `svarlist'
		global svarlist
		postclose `postname'
		preserve
		use `postfile', clear
		foreach var in `varlist' {
			foreach wvar of varlist ${wvarlist} {
				if strpos("`wvar'","`var'") >0 {
					char def `wvar'[scale] $`var'unit
				}
			}
		}
		gen run_nmbr = _n
		    label var run_nmbr "Sequence number of run"
		local nwords = wordcount(runspecs)
        foreach i of numlist 1/`nwords' {
    		qui gen runspec`i' = real(word(runspecs,`i'))
        }
        order run_nmbr runspec*
        sort run_nmbr
        di _n as txt "Saving the file of harvested results." _n
		if  c(stata_version) > 12 {
			saveold `postfile', replace version(12)
		}
		else {
			save `postfile', replace
		}
		global wvarlist
        restore
	}

end   /* End of program aidspost */


prog define costmodel, rclass
*	Added in version 5.10 and greater, 14Feb2018
	syntax, [rxc1lb(string) rxc1ub(string) rxc2lb(string) rxc2ub(string) hsbedn(string) hsvstn(string) scale(string) ndaf(string) ]

    foreach cname in rxc1lb rxc1ub rxc2lb rxc2ub hsbedn hsvstn  {
        // Check to make sure user has not entered a string or an out of bounds value into a parameter
        if "``cname''"~="" {
			if (real("``cname''")>=. | real("``cname''")<0) {
				di _n as error "Parameter `cname' must be numeric and greater than or equal to zero."
				di    as error "Instead it is: ``cname''"
				error 198
			}
			local usr_`cname' = ``cname''
        }
    }
	if "`scale'"~="" {
		if (real("`scale'")>=. | real("`scale'")<-1 | real("`scale'")>0) {
			di _n as error "Parameter -scale- must be numeric and between -1.0 and zero."
			di    as error "Instead it is: `scale'"
			error 198
		}
		local usr_scale = `scale'
	}
	if "`ndaf'"~="" {
		if real("`ndaf'")<1.0 {
			di _n as error "Average non-drug fixed cost must be at least $1. "
			error 198
		}
		local usr_ndaf = `ndaf'
	}

	local dflt_rxc1lb    88
	local dflt_rxc1ub    261
	local dflt_rxc2lb    819
	local dflt_rxc2ub    2634
	local dflt_hsbedn    1.56
	local dflt_hsvstn    9.50
	local dflt_ndaf  750
	local dflt_scale     -.146

    di as txt _n "Cost computations based on following parameters"
    di as txt "  Lower bound for 1st-line drug costs"           _col(60)  "rxc1lb" as res _col(72) "   $`dflt_rxc1lb'" as res _col(82) `usr_rxc1lb'
    di as txt "  Upper bound for 1st-line drug costs"           _col(60)  "rxc1ub" as res _col(72) "  $`dflt_rxc1ub'"  as res _col(82) `usr_rxc1ub'
    di as txt "  Lower bound for 2nd-line drug costs"           _col(60)  "rxc2lb" as res _col(72) "  $`dflt_rxc2lb'"  as res _col(82) `usr_rxc2lb'
    di as txt "  Upper bound for 2nd-line drug costs"           _col(60)  "rxc2ub" as res _col(72) " $`dflt_rxc2ub'"   as res _col(82) `usr_rxc2ub'
    di as txt "  Number of bed-days per year per patient"       _col(60)  "hsbedn" as res _col(72) "  `dflt_hsbedn'"   as res _col(82) `usr_hsbedn'
    di as txt "  Number of out-patient visits per patient"      _col(60)  "hsvstn" as res _col(72) "  `dflt_hsvstn'"   as res _col(82) `usr_hsvstn'
    if `"`cost'"'=="function" {
        di as txt "  Average fixed non-drug cost is function of scale and GNI per capita"
        di as txt "  Elasticity of average fixed cost w.r.t. ART"   _col(60)  "scale"       as txt _col(72) " -.142" _n
    }
    else {
        di as txt "  Average fixed non-drug cost at ART=1000"       _col(60)  "ndaf"    as result _col(72) "  $`dflt_ndaf'"  as res _col(82) `usr_ndaf'
        di as txt "  Elasticity of average fixed cost w.r.t. ART"   _col(60)  "scale"       as result _col(72) " `dflt_scale'"  as res _col(82) `usr_scale' _n
    }

    foreach cname in rxc1lb rxc1ub rxc2lb rxc2ub hsbedn hsvstn ndaf scale {
		if "`usr_`cname''" ~= "" {
			return local `cname' `usr_`cname''
		}
		else {
			return local `cname' `dflt_`cname''
		}
    }
end  /* End of program costmodel  */

	*	Version 1.6 removes the non-standard programs -expr- and -mmerge-
*	 and adds weo graph to gr7 output
*	Version 1.7 adds a log option, to generate a separate log file for each country, region or for the global totals.
*	Version 1.8 fixes problem that titlepct can exceed 100 % in a few cases.
*	Version 1.9 graph labels when graphs not saved, textsize on combined gr7 graphs
*	Version 2.0 adds the weights option, with a weights characteristic added to the data file
*	Version 2.1 adds the vintage option, designed to operate with vintage 2007 or vintage 2008 data
*       Assumption is that the weo data will be vintage less 2 years
*	Version 2.11 fixes bug in the expand command at line 260
*	Version 2.12/3 fix bug in the save command at line 116 and add completion footer
*	Version 2.16 adds check of datasignature
*	Version 2.17 replaces command -tabstat- with -tabstat2-
*	Version 2.17 defines vweo as vintage-2 instead of vintage-3
*	Version 2.19 fixes bug in the datasignature check
*	Version 2.20 adds variable descriptions to table headers in projections
*	Version 2.21 changes -uptake()- option to expect (1-lambda) instead of -lambda-
*	Version 2.22 fixes another bug in the datasignature
* 	Version 2.23 updates datasignature to reflect dataset changes
*   Version 3.1 adds specification of CD4 at initiation   /* Ver. 3.0 Flag */
*   Version 3.11 adds the "| x < 0" condition for replaicing missing values of d, n, i, p
*   Version 3.12 adds the [morebadgood] characteristic to most of the variables in the saved dta file
*       and checks that the saving option contains only the stub of a filename. Also it assures
*       that cvrg2 hits its target in the target year.
*   Version 3.13 updates cost specification, fixes missing 2007 values for cvrg2, fundgap
*	Version 3.2 includes options to post results.
*	Version 3.21 removes 240 character limit on options allowed by aidsproj
*	Version 3.22 changes all "as input" displays to "as res" displays
*	Version 3.23 fixes parsing of -runspecs- by using -gettoken- instead of -foreach token-
*	Version 3.24 add -drate- option to call for discounted sums of flow variables:
*       	-dfactor- added as an analytic weight to tabstat stat(sum) commands
*	Version 3.25 postclose option now separates out the runspec variables, but does not name them
*	Version 3.26 posttotal option fixed: It's chosen if -postyears- are omitted or if explicitly selected
*	Version 3.27 adds characteristic `var'[unit] to cost variables
*	Version 3.27 fixes bug in -aidspost, harvest- that assigned world total to all sub-aggregates
*	Version 3.28 replaces variable -art- with variable -un- in first tabstat output
*	Version 3.29 adds the -TAkeoff(passthru)- option to allow program to "taxi" for a few years
*			before "taking off" and dropping output in -genout- using the -tabstat2- command.
*	Version 3.30 changes the formulae for computing -ndrate_dflt- and -adrate1_dflt- to conform to
*           the incremental hazard rate estimated in the April, 2009 WTSC paper.
*	Version 3.31 changes the formula for computing -edrate_dflt- to conform to eART-LINC (2008) with Thai cohorts
*	Version 3.32 forces uptake to remain historical through the "takeoff" if the takeoff parameter is specified
*           Requires that a variable named -uptake- be in the input data to capture the historical uptake of each country.
*           Also fix the grcombine option.
*	Version 3.33 changes the length of the -runspecs- variable from str80 to str200 in line 1588
*   Version 3.42 requires the uptake option to be not only < 1, but also >= 0.
*           Also the -incmult- option now reduces incidence only after the year -takeoff-
*           Change the upper bound of the simulation loop to -horizon-.
*           When takepff option is specified, each country goes from historical to selected constant uptake rate
*           over a three year period.  Lines 381-383.
*   Version 3.43 Capture error and postclose -postname- if -postname- exists from a previous run.
*           Keep posted scalars after -postclose-.
*   Version 3.44 Add -rachet- option which forces adjustment to new uptake target immediately after -takeoff-
*           instead of with a three-year smoothed adjustment.  change <= to < in line 380.
*   Version 3.47 Add a graph of new infections and deaths in order to track the AIDS treansition.
*   Version 3.48 Change "aidspct*" vrbls to "artpct*" vrbls
*   Version 3.49 Fix -extweo.ado- so that projections of the artpct variables are plausible.  Make -extweo- external to -aidsproj.ado-
*   Version 3.50 Fix reversed variable labels on AIDS transition graph.  Adjust defualt value of incmult to 0.97.  Tweak graph titles, legends.
*   Version 4.01 Add error trap in -parmsbycd4- for user specifying -incmult- without -iexog-.  Require Stata version 10.1
*   Version 4.02 Assure that the -postnow- and -postclose- commands at the end of -genout- need to function for any call to -genout-, not just -world- calls.
*   Version 4.04 Fix typo in definition of -tau- in line 482.
*   Version 4.05 Add panel labels and narrow margins for four-panel graphs selected with the -grcombine- option.
*   Version 4.08 Add the -vintage- designation to the arguments passed to -parmsbycd4-
*   Version 4.09 Fix bug in scale up of 2nd-line using -cvrg2-
*   Version 4.10 Predict GNI per capita (-gnipc_hat-) in order to use -function- option for cost function
*   Version 4.12 Get logout function to work for -world- option
*   Version 4.13 Smooth transition of ndrate/erate to new values given by CD4 to account for preexistinng stocks of un/h
*   Version 4.14 Smooth the death variable using tssmooth with window (4,1,2).  Subsume the rachet option into new noSMOOTH option
*   Version 4.15 Fix bug in projection of -artpctcg-
*   Version 4.16 Modified so that vintage 2010 expects WDI afforditure data to be from 2007 and fix bug when -strtyr- < 2009.
*   Version 4.17 Fix the vintage display and drop years after the projection horizon.
*   Version 4.19 Update -ndrate- to use the ART-Linc results which exclude Thialand, and thus are more optimistic.
*   Version 4.20 Extend horizon to 2100 & add the multiplicative term (k-prevprop)/(k-prevprop0) to the computation of -ir-
*   Version 4.21 Truncate historical uptake variab, le from input data to lie within the range .01 and .99.  (Lines 210-50)
*      Adjust vintage screen to allow vintage 2011.  Choose NOT to add calculations of -dy- & -ly- to -aidsproj-,
*      since they are cumulated variables and thus raise the issue of discounting.  See -aidsdif- version 2.13.
*      Make the first year of 2nd line coverage, parameter -strtyr-, default to 2007
*   Version 4.22 (20Mar2013) To run with vintage 2011, change the date 2009 to `vminus1' in lines 461-468
*      Assure that the truncation of `lmbtmp' in lines 556-558 is only for the observations for which uptake is defined
*      Assure that a1 is never negative in line 669
*   Version 4.23 (8/5/2013) Changes the third plotted curve in the affordability graph
*      to the percentage that ART patients are of the total population.  Adds memory names to
*      graphs when the saving option is not used.  Introduces intermediate calculations of HE and Public HE
*      Change legendt position from 6 to 0, so that the legend will not obscure the bottom axis.
*	Version 4.50 (1/29/2014) Enable vintage 2013
*   Version 4.51 Drop years < vminus1 from aidsproj output file.
*   Version 4.52 Add the options: rachet, linesize, putexcel, change the gnipc variables to gnipc`vweo' & gnipc_hat
*   Version 4.53 Distribution version of 4.52
*   Version 4.54 Add the "& i[_n-1]>=d[_n-1]" condition to declare an "AIDS Transition"
*   Version 4.55 (2 June 2014) Add the line 705 containing "if year == `t' " within the projection loop.
*   Version 4.56 (16 Dec 2014) Add a multiplicative fudge factor option to calibrate
*      a regional or world total number of HIV infected to include countries
*      missing from UNAIDS published data.  Default is 100%, or no adustment.
*   Version 4.57 (25 Jan 2015) Add the ability to check the aidscgd2014.dta data signature
*   Version 4.58 (27 Jan 2015) Allow vintage 2014
*   Version 4.59 (17 Feb 2015) Change all ==. to >=. and ~=. to <. to accommodate code for censored data.
*       Edit the local rlblarg ", sigma(`uptake') "  command
*       Simplify the subroutine select's imputation code
*   Version 4.60 (20 Feb 2015) Investigate determinants of cost per patient-year of treatment
*		Insert error trap for the paramters incmult ir0 maxep maxcp maxvp ve fp gp dt du
*   Version 4.61 (23 Feb 2016) Update to accept new data aidscgd2015 as "vintage 2015"
*   Version 4.62 (23 Mar 2016) Protect against the tssmooth command failing when horizon < takeoff + 6 (line 763)
*		Add noisily qualifier in the -select- subroutine
*   Version 5.00 (26 Jan 2017) Allow vintage 2016. Change default value of -maxep- from .6 to 0.
*   Version 5.01 (7 Oct 2017) Always use 2 lines to specify incidence assumptions.
*	Version 5.02 (13Feb2018) Convert to UTF 8 encoding.  Add option -outpath(string)-
*	Version 5.03 (14Feb2018) Add default text for incidence as well as treatment paramters.
*		Also change vintage range to 2007 <-> 2017 and default value to 2014
*		Change maxvp default (maxVp_dflt) from 0.7 to zero.
*		Replace note(`grnote') with note(`grnote',span) beginning in line 1811
*	Version 5.10 (14Feb2018) Add the program -costmodel- to prepare for better cost modeling
*		Add the option -outpath- so that user may specify a path for output files.
*	Version 5.11 (15Feb2018) Fix bug in saving graphs by adding "gph" extension in lines 1681-1687
*	Version 5.12 (16Feb2018) Once the gph extension is added explicitly to the `saving' options,
*		the new option outpath is no longer necessary.  Retain the option for backward compatibility with
*		version 5.11.  But change error trap to only catch the asterisk character in the -saving(stub)- option.
*	Version 5.13 (19Feb2018) Change the option nonrxcaf to ndaf so that the name nonrxcaf is reserved for
*		the name of the generated variable
*	Version 5.20 (4Dec2018) Set total needingg treatmment to p if n is missing. Line 481
*	Version 5.21 (5Dec2018) Switch from -save- to -saveold ..., version(11)-
*		Only keep years after vminus2.  Silence most -keep- commands.
*	Version 5.22 (8Feb2019) Display saved data files at run completion
*	Version 5.23 (12Feb2019) Use -saveold, version(12)- only if c(stata_version) > 12
*	Version 5.24 (19Nov2019) Define default -cvrg2- in line 546 as "gen cvrg2 = cond(year <= `strtyr', 0, `strtcov2' ) "
*	Version 5.25 (25Nov2019) Smooth -cvrg2- definition.  Fix -logout-, -excel- options.
*	Version 5.26 (26Nov2019) Precede all instances of `saving' with `outpath'
*	Version 5.27 (31Jan2020) Enable different values of strtcov2 by region.  Lines 532-573
*	Version 5.30 (13Feb2020) Add options for minimum values of ep and cp (mincp and minep)
*		Add cstrtyr and ctrgtyr for the circumcision intervention.
*		Add an output table to log and optionally as an Excel sheet for "Coverage" variables
*		Replace -note(`grnote',span)- with -note(`"`grnote'"', /* span */ size(vsmall))- beginning in line 1875
*		In the syntax command, change the default vintage to 2018
*		When -iexog- selected, remove ccvrg, vr from output
*		Change the cost variable characteristic `var'[unit] to `var'[scale]
*		Set the value of ccvrg, effort and vr to their minimum prior to start year.
*	Version 5.31 (11Nov2020) Under construction
*		option -impute- produces an ugly table. Commenting out the variable -complete- in lines 1263 & 1278 doesn't help
*		Possibly could incorporate in the -repshow- command Stata's option for -list-: -display- to prevent the table format
*		Enhance appearance of the footer and message when projections are not saved. Lines 1163-72