program define aidsdif, rclass
*!	Program -aidsdif- by Mead Over  Version 2.51     8Dec2020
*	Compute the benefits and costs of scenario 1 (& 2, 3, 4, 5, 6) compared to the baseline scenario.
*	Calls programs -rescale.ado- and -grc1leg_aidscost- which must be on the -adopath-

	version 12

	syntax namelist(name=avnames id="Impact analysis variables"), bl(string) s1(string) [DRate(real .03) OTHVar(namelist) Weights(string) /*
		*/ bl_name(string) s1_name(string) s2(string) s2_name(string) /*
		*/ s3(string) s3_name(string) s4(string) s4_name(string) /*
		*/ s5(string) s5_name(string) s6(string) s6_name(string) /*
		*/ REFYear(integer 2008) COSTVar(namelist) COSTUnit(string) BENVar(namelist) BENUnit(string) /*
		*/ Countries(passthru) Regions(passthru) World GR7 GRAph GRCOmbine noRESCale LEGend(passthru) /*
		*/ GRScenarios GRDifferences grre grcebyh noGRCE GRNote(string) /*
		*/ TITLES(passthru) TITLED(passthru) TITLER(passthru) TITLEB(passthru) /*
		*/ TITLEV(passthru) TITLEC(passthru) SUBTitle(passthru) saving(string) altshrink noTRUNCATE putexcel /*
		*/ grc1leg_options(passthru) grc1leg outpath(string) inpath(string) gricer DEBUG]

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

*	Print header
	di _n _n as txt "AIDSDIF: Comparison & C/E of AIDS Projection Scenarios"
	di       as txt "               Version 2.51"
	di       as txt "{hline 54}"

	if `"`saving'"'~="" {
		if strpos(`"`saving'"',",")>0 {
			local cposm1 = strpos(`"`saving'"',",")-1
			if strpos(`"`saving'"',"replace")>0 {
				local replace replace
				local rpl ", replace"
			}
			local saving = trim(substr(`"`saving'"',1,`cposm1'))
		}
		else {
			local saving = trim(`"`saving'"')
		}
	}
	if (`"`outpath'"'~="" ) & (strpos(`"`saving'"',"/") | strpos(`"`saving'"',"\")) {
		di as err "Include the output path either as part of the filename in the -saving(filename)- option"  ///
			_n "or in the -outpath(filepath)- option, but not in both.  Here:"  ///
			_n "the -saving(filename)- option is: " `"`saving'"'  ///
			_n "and the -outpath(filepath)- option is: " `"`outpath'"'  ////
			_n "-aidsdif- is trying to save to the file: " `"${outpath}`saving'"'
		exit 198
	}
	if `"`outpath'"'~="" {
		local saving `"`outpath'`saving'"'
	}
	foreach scen in bl s1 s2 s3 s4 s5 s6 {
		if (strpos(`"``scen''"',"/") | strpos(`"``scen''"',"\")) {
			local scenispath scenispath
		}
	}
	if (`"`inpath'"'~="") & "`scenispath'"~="" {
		di as err "Include the input path either as part of the filename"  ///
			_n "in all of the selected scenario options, bl(), s1(), ..., s6(),"  ///
			_n "or in the -inpath(filepath)- option, but not in both.  Here:"  ///
			_n "the -inpath(filepath)- option is: " `"`inpath'"' ///
			_n "and the scenario options are:"
			foreach scen in bl s1 s2 s3 s4 s5 s6 {
				if `"``scen''"'~="" {
					di as err "    The expected file for scenario `scen' is: " `"``scen''"'
				}
			}
		exit 198
	}
	if `"`inpath'"'~="" {
		foreach scen in bl s1 s2 s3 s4 s5 s6 {
			if `"``scen''"'~="" {
				local `scen' = `"`inpath'``scen''"'
			}
		}
	}

	di _n as txt "Impact analysis requested for the following variables: " as res "`avnames'"

	local vnum  : word count "`avnames'"
	if `vnum' > 5 {  /* Set maximum number of variables to be analyzed */
		di _n as error "Error: Number of variables greater than 5, the maximum set by program -aidsdif-"
		exit 198
	}

	di    as txt "Using discount rate " as res `drate'
*		Since the data files saved by -aidsproj.ado- are only annual numbers, and never include results
*		summed across years, there is no danger that discounted sums here in -aidsdif.ado- are of results
*		that have already been discounted

	di _n as txt "Costs and effects discounted to the reference year of " as res "`refyear'"

	di _n as txt "Baseline scenario: " as res "`bl_name'"
	di _col(4) as txt " defined by data file: " as res "`bl'"
	if "`bl_name'"~="" {
		local scenes = `"bl_name(`"`bl_name'"')"'
	}

	if "`costvar'"=="" {
		local costvar "costtot"
	}
	else {
	    if wordcount(`"`costvar'"')>1 | strpos(`"`avnames'"',`"`costvar'"')==0 {
			di _n as error "-costvar()- option specified as: " as res "`costvar'" as err " ,"  ///
				"should specify only one variable from list of variables to be analyzed: " as res "`avnames'"
			error 198
		}
		else {
			di _n as txt "User has defined the cost variable as " as res "`costvar'."  ///
				_n as txt "-aidsdif- assumes this variable is a -flow- variable measured in thousands of current US dollars."
		}
	}
	if "`costunit'"=="" {
		local costunit "Dollars" /* "\$1,000" */
	}

	//	Accommodate using life-years saved or death-years averted as measures of benefits
	//	In version 2_34, switch back to making life-years saved the default benefit variable

	if "`benvar'"=="" {
		local benvar "ly"	// Default is life-years
	}
	else {
	    if wordcount(`"`benvar'"')>1  {
			di _n as error "The -benvar()- option is specified as: " as res "`benvar'" as error " ."
			di    as error "The -benvar()- option should specify only one variable."
			error 198
		}
		if strpos(`"`avnames'"',`"`benvar'"')==0 & `"`benvar'"'~="ly" & `"`benvar'"'~="dy" & `"`benvar'"'~="vsly" {
			di _n as error "The -benvar()- option is specified as: " as res "`benvar'" as error " , but the list of variables"
			di    as error "selected for -aidsdif- to analyze is " as res `"`avnames'"' as err " and does not include this variable"
			error 198
		}
	}
	if "`benunit'"=="" {
		if "`benvar'"=="d" {
			local benunit "death averted"
		}
		if "`benvar'"=="dy" {
			local benunit "death-year averted"
		}
		if "`benvar'"=="ly" {
			local benunit "life-year saved"
		}
		if "`benvar'"=="vsly" {
			local benunit "dollars of life-years saved"
		}
		if "`benvar'"=="p" {
			local benunit "prevalent case averted"
		}
		if "`benvar'"=="i" {
			local benunit "new infection averted"
		}
	}
	if "`grcombine'"~="" & "`gr7'`graph'"=="" {
		di _n as error "Option -grcombine- specified without specifying either option -gr7- or option -graph-"
		exit 198
	}
	if ("`grc1leg'"~="" | "`grc1leg_options'"~="") & "`graph'"=="" {
		di _n as error "Option -grc1leg- specified without specifying option -graph-"
		exit 198
	}

	di _n as txt "Scenario 1:  " as res "`s1_name'"
	di _col(4) as txt " defined by data file: " as res "`s1'"
	local s = 2
	local last_s = 1
	while `s'<=6 {
		if "`s`s''"~= "" {
			if `s'>`last_s'+1 {
				local sm1 = `s'-1
				di _n as error "Define Scenario " as res "`sm1'" as error " before defining Scenario " as res "`s'" as error "."
				error 198
			}
			di _n as txt "Scenario `s':  " as res "`s`s'_name'"
			di _col(4) as txt " defined by data file: " as res "`s`s''"
			local last_s = `last_s' + 1
		}
		if "`s`s'_name'"~="" {
			local scenes = `"`scenes'"' + " " + `"s`s'_name(`s`s'_name')"'
			note _dta : Scenario S`s' is `s`s'_name'
		}
		local s = `s' + 1
	}

*	Merge the simulation data to be used in the analysis of impact
    tempfile bl_vars aidsdif

*   Before proceeding, set aside weight variables & other variables from the baseline not to be renamed
	qui use `"`bl'"'
	if "`weights'"=="" {
	    local weights art
	}
    keepwts , weights("`weights'") confirm
    local wtvlst `r(wtvlst)'

	keep  code year `wtvlst' `othvar'
    qui compress
    sort code year
    qui save `"`bl_vars'"'

*	Now reload the baseline scenario & get started constructing aidsdif.dta
	qui use `"`bl'"', clear

	if "`costvar'"~="" {
			char define `costvar'[stockflow] flow
			char define `costvar'[scale]  1000
			char define `costvar'[morebadgood] moreisbad
			char define `costvar'[unit] Dollar
	}
	if "`benvar'"=="dy" | "`benvar'"=="ly" {
		qui makedyly , refyear(`refyear')
		local dyly dy ly
		local avnames : list avnames | dyly
		local avnames : list uniq avnames
	}
	if "`benvar'"=="vsly" {
		qui makevsly , refyear(`refyear')
		local dyly dy ly vsly
		local avnames : list avnames | dyly
		local avnames : list uniq avnames
	}

	local bl_vntg : char _dta[vintage]

	di _n as txt "Keep selected variables from file: " as res "`bl'"
	keep code year cntry ccode nation region rgn `avnames'

	isflow, costvar(`costvar') benvar(`benvar')

	sum year, meanonly
		local maxyr = r(max)
		local vntg : char _dta[vintage]
		if "`vntg'"=="" {
			local vntg 2007
		}
		local minyr = max(r(min), `vntg')
	if `refyear' < `minyr' |  `refyear' > `maxyr' {
		di as error "Reference year " as res `refyear' as error " out of bounds.  Must be >= " ///
			as res "`minyr'" as err " and < " as res "`maxyr'"
		exit 198
	}

*	Compute discount factor
	gen vrefyr = 1/(1+`drate')^(year-`refyear')
		lab var vrefyr "`refyear' value of a future unit at r=`drate'"

	tokenize `avnames'

	local v = 1
	while "``v''" ~= "" { /* Rename all the variables under analysis */
		local varlbl : variable label ``v''
		rename ``v'' bl_``v''
			lab var bl_``v'' "`varlbl', BL"
			char define bl_``v''[origvname] "``v''"
			notes bl_``v'': "Baseline values of ``v''"
		local v = `v' + 1
	}

	sort code year
	label data "Baseline data for impact analysis"
		note _dta : Baseline scenario is `bl_name'
	qui compress
	qui save `"`aidsdif'"', replace

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

*	Load analysis scenarios & compute impacts and CE relative to the baseline
*	one scenario at a time
	local s = 1  // Index of the scenarios to be compared to the baseline
	local slbl = "S1"
	while `"`s`s''"' ~= "" {
		qui use `"`s`s''"', clear

		if "`costvar'"~="" {
				char define `costvar'[stockflow] flow
				char define `costvar'[scale]  1000
				char define `costvar'[morebadgood] moreisbad
				char define `costvar'[unit] Dollar
		}
		if "`benvar'"=="dy" | "`benvar'"=="ly" {
			qui makedyly , refyear(`refyear')
		}
		if "`benvar'"=="vsly"  {
			qui makevsly , refyear(`refyear')
		}

		keep  code year cntry ccode nation region rgn `avnames'

		local v = 1  /* Rename all the variables under analysis */
		while "``v''" ~= "" {
			local varlbl : variable label ``v''
			rename ``v'' s`s'_``v''
				lab var s`s'_``v'' "`varlbl', S`s'"
				char define s`s'_``v''[origvname] "``v''"
				notes s`s'_``v'': "Scenario `s' values of ``v''"
			local v = `v' + 1
		}

		sort code year

		di _n as txt "Merge the selected variables from file: " as res "`s`s''"
		qui merge  code year using `"`aidsdif'"', unique _merge(_merge`s')
			lab var _merge`s' "Results of merge  of scenario `s' with aidsdif dta file"

		sort code year

		local v = 1  /* Compute impact on all the variables under analysis */
		while "``v''" ~= "" { /* For each variable in -avnames- */
			qui gen dif`s'_``v'' = s`s'_``v'' - bl_``v''
				lab var dif`s'_``v'' "Impact of S`s' on ``v''"

*				format dif`s'_``v'' %15.0fc	       //  On 15Feb2020, this line replaced by:

			sum dif`s'_``v'' , meanonly            //  this line to line 329 in order to
				local range = `r(max)' - `r(min)'  //  allow digits after the decimal
				if `range' > 2 {                   //  when rescaled range is small
					format dif`s'_``v'' %15.0fc
				}
				else {
					if `range' > 1 {
						format dif`s'_``v'' %15.1f
					}
					if `range' > .5 {
						format dif`s'_``v'' %15.2f
					}
					else {
						format dif`s'_``v'' %15.3f
					}
				}                                 //  Lines new on 15Feb2020 end here

			local thisscale : char bl_``v''[scale]
				char define dif`s'_``v''[scale] `thisscale'

			chkmoresign bl_``v''
				char dif`s'_``v''[morebadgood] `r(morebadgood)'
			local v = `v' + 1
		}

*	Use the discount factor -vrefyr- to compute present values in the reference year
*	of impact variables (with -dif- prefix) for all the variables under analysis
		local v = 1
		while "``v''" ~= "" { /* For each country and each variable in -avnames- */
			qui by code: egen pv`s'_``v'' = total(vrefyr * dif`s'_``v'') if year>`refyear'-1
				lab var pv`s'_``v'' "PV of impact of S`s' on ``v'' at r=`drate'"
			qui replace pv`s'_``v'' =. if year > `refyear'
			local thisscale : char dif`s'_``v''[scale]
				char define pv`s'_``v''[scale] `thisscale'
			local thisbadgood : char dif`s'_``v''[morebadgood]
				char define pv`s'_``v''[morebadgood] `thisbadgood'

			sum pv`s'_``v'', meanonly /* Check for zero denominator of CE and RE ratios */
			if abs(r(max) - r(min)) < .001 & abs(r(mean)) < .001 & r(mean)~=. {
				di _n as error "Warning.  Present value of the impact of Scenario " as res "`s'"
				di    as error "on variable " as res "``v''" as error " is smaller than 0.001 for all countries."
				if "``v''"=="`benvar'" {
					di _n as error "Since variable " as res "``v''" as error " is the benefit variable,"
					di    as error "its present value appears in the denominator of the CE calculation."
					di    as error "Either choose another benefit variable or assure that  " as res "``v''"
					di    as error "differs between the scenarios and the baseline."
					exit 198
				}
			}

			sum year, meanonly  // For pvbyhorizon, use the max of the horizons in the compared projections
				local horizon `r(max)'
			pvbyhorizon , var(``v'') scenario(`s') refyear(`refyear') horizon(`horizon') drate(`drate')

			local v = `v' + 1
		}
		capture confirm variable bl_`costvar'
		if _rc~=0 {
			di _n as txt "Cost variable " as res "`costvar'" as txt " not included in -varlist-."  ///
					as txt "No cost-effectiveness calculations can be performed."
			error 198
		}
		else {       /* Compute cost-effectiveness of scenarios vis-a-vis baseline */
			local v = 1
			while "``v''" ~= "" { /* For each variable in -avnames- */
				if "``v''" ~= "`costvar'" {
					//  Numerator is "bad" or "good"
				    local numbg   : char dif`s'_`costvar'[morebadgood]
					//  Denominator is "bad" or "good"
				    local denombg : char dif`s'_``v''[morebadgood]
                    qui gen     re`s'_``v'' = .
                    if "`numbg'"=="`denombg'" | ("`numbg'"=="both" | "`denombg'"=="both") {
                        qui replace re`s'_``v'' = -(dif`s'_`costvar'/dif`s'_``v'')
                    }
                    else {
                        qui replace re`s'_``v'' = (dif`s'_`costvar'/dif`s'_``v'')
					}

					if "``v''"=="`benvar'" {
						if "`costunit'"=="Dollars" {  //  Assume cost units are in $1,000's
							local thsnds `" ($1,000's)"'
						}
						lab var re`s'_``v'' "S`s': Annual ratio of `costunit' to `benunit'`thsnds'"
						qui gen ce`s'_``v'' = .     //  CE over entire horizon
						qui gen cebyh`s'_``v'' = .  //  CE to any given horizon
                        if "`numbg'"=="`denombg'" | ("`numbg'"=="both" | "`denombg'"=="both") {
                            qui replace    ce`s'_``v'' = -(   pv`s'_`costvar'/   pv`s'_``v'') if year == `refyear'
                            qui replace cebyh`s'_``v'' = -(pvbyh`s'_`costvar'/pvbyh`s'_``v'')
                        }
                        else {
                            qui replace    ce`s'_``v'' = (   pv`s'_`costvar'/   pv`s'_``v'') if year == `refyear'
                            qui replace cebyh`s'_``v'' = (pvbyh`s'_`costvar'/pvbyh`s'_``v'')
                        }
						lab var    ce`s'_``v'' "CE of S`s': `costunit' per `benunit'`thsnds'"
						lab var cebyh`s'_``v'' "CE of S`s': `costunit' per `benunit'`thsnds' by horizon"
					}
					else {
						lab var re`s'_``v'' "S`s': Annual ratio of `costunit' to unit of ``v''`thsnds'"
					}
				}
				local v = `v' + 1
			}
		}
		label data "Baseline data plus scenarios `slbl' for impact analysis"
		qui compress
			note _dta : Scenario `s' is `s`s'_name'
		qui save `"`aidsdif'"', replace

		local s = `s' + 1
		local slbl = "`slbl'" + ", S`s'"
		local thisvntg : char _dta[vintage]
		local minvntg = min(`bl_vntg', `thisvntg')
	}
	local ns = `s' - 1  /* Number of scenarios */

*   Merge back in the variables from the baseline that have been set aside and not renamed
	qui merge  code year using `"`bl_vars'"', unique _merge(_mrgoth)
		lab var _mrgoth "Results of merge  of scenarios with other variables"
	order code year cntry ccode nation region rgn vrefyr `othvar' `wtvlst'
	qui keep if year>=`minvntg'
	sort code year
	qui save `"`aidsdif'"', replace

	di _n(3)  as txt "Variables containing the annual ratios of costs to effects: by country, year"
	des re?_`benvar'
	qui xtsum year

	di _n  as txt "Averages of annual ratio of cost to effect over the " as res r(n) as txt " countries by year:"
    di     as txt "  Weighted by " as res "`weights'"
	tabstat2 re?_`benvar' [aweight=`weights'], by(year) format(%9.3f) noheader describe

	di _n(2) as txt "Variables containing the cost-effectiveness over the projection horizon: by country"
	des cebyh?_`benvar' ce?_`benvar'
	sum cebyh?_`benvar' ce?_`benvar'

	di _n  as txt "Averages of cost-effectiveness over the " as res `r(N)' as txt " countries by region:"
    di     as txt "  Weighted by " as res "`weights'"
	tabstat2 ce?_* [aweight=`weights'], by(rgn) format(%9.3f) stats(mean count) noheader describe

	if "`output'" ~= "" {
		di _n as txt "Output requested for: " as res "`output'"
		notes _dta
		genout2 , refyear(`refyear') ns(`ns') costvar(`costvar') costunit(`costunit') vnames(`avnames') /*
			*/ benvar(`benvar') benunit(`benunit') `countries' `regions' `world' weights(`weights') /*
			*/ `gr7' `graph' `grcombine' grnote(`grnote') `legend' /*
			*/ `grscenarios' `grdifferences' `grre' `grcebyh' `grce' drate(`drate') /*
			*/ `titlen' `titlec' `titlep' `titlee' `subtitle' /*
			*/ saving(`"`saving'"') `replace' `rescale' `weo' bl_name(`bl_name')  /*
			*/ s1_name(`s1_name') s2_name(`s2_name') s3_name(`s3_name')  /*
			*/ s4_name(`s4_name') s5_name(`s5_name') s6_name(`s6_name')  /*
			*/ `altshrink' minvntg(`minvntg') horizon(`horizon') `truncate' /*
			*/ `putexcel' `grc1leg' `grc1leg_options' `gricer' `debug'

	}
	else {
		di _n as txt "No output requested."
	}

    if `"`saving'"'~="" {
        qui compress
        save `"`saving'"' `rpl'
    }

	if "`putexcel'"~="" {
		local putexcel putexcel A1 = matrix(\`rmatrix',names) using `saving', sheet(\`matrix')
		if "`replace'"~="" {
			local excelmodify modify
		}
	}

	local matrixlist `r(ce_matrixlist)'  /* `r(pv_matrixlist)'  not used*/
    foreach matrix of local matrixlist  {
		tempname rmatrix
		matrix rename r(`matrix') `rmatrix'
			`putexcel' `excelmodify'
		return matrix `matrix' `rmatrix'
/*
		if "`putexcel'"~="" {
			local excelreplace modify
		}
*/
    }

*	Print footer
	di _n _n as txt "AIDSDIF: Completion of run"
	di       as txt "=========================="

end     /* End of -aidsdif- main program */

*************************SUBROUTINES************************************

program define genout2, rclass
*	Generate output from aidsdif
syntax , refyear(integer) ns(integer) costvar(namelist) costunit(string) vnames(namelist) /*
	*/ benvar(namelist) benunit(string) drate(real)  horizon(integer) [ countries(namelist) regions(namelist) world  weights(string) /*
	*/ gr7 graph grcombine grnote(passthru) legend(passthru) /*
	*/ grscenarios grdifferences grre grcebyh nogrce  /*
	*/ titlen(passthru) titlec(passthru) titlep(passthru) titlee(passthru) subtitle(passthru)  /*
	*/ saving(string) replace norescale weo bl_name(string)  /*
	*/ s1_name(string) s2_name(string) s3_name(string) s4_name(string) /*
	*/ s5_name(string) s6_name(string) altshrink minvntg(integer 2006) noTRUNCATE  /*
	*/ putexcel grc1leg_options(passthru) grc1leg gricer debug]

if "`countries'"=="" & "`regions'"=="" & "`world'"=="" {
	di as error "This error should never occur, since any call to genout2 is specific to one of these three options"
	exit 198
}

if "`weo'"~="" {
	local weovars aidspcthe aidspctpub aidspctprv aidspctoop aidspctcg
	capture confirm numeric variable `weovars'
	if _rc  {
		di _n as error "WEO option NOT YET IMPLEMENTED"
		exit 198
	}
}

local lols
foreach V in `vnames' {
	local vlist`V' bl_`V'
	local s = 1
	while `s' <= `ns' { /* List of baseline, scenario and difference variables for variable V */
		local vlist`V' `vlist`V'' s`s'_`V' dif`s'_`V'
		local s = `s' + 1
	}
	local lols `lols' `vlist`V''   /* Assemble a list of the lists */
}

if "`putexcel'"~="" {
	local tabsave save
	local putexcel putexcel A1 = matrix(r(StatAll),names) using `saving' , sheet(\`sheetname')
	if "`replace'"~="" {
		local excelreplace replace
		local excelmodify modify
	}
}

*	Projections for countries
if "`countries'"~="" {
	foreach X in `countries' {
		local pv_costlist
		local pv_effectlist
		local Y = upper("`X'")
		di _n as txt "Projections for Country " as res "`Y'"
		di    as txt "============================"
		des `lols'

		if "`putexcel'"~="" & "`xlr'"=="" {
			local xlr `excelreplace'
		}
		foreach V in `vnames' {
			local sheetname Impact_on_`V'_for_`Y'
			di _n as txt "----Projected impact on variable " as res "`V'" in gr" for country " as res "`Y'" in gr"----"
			tabstat2 `vlist`V''    if code=="`Y'":code , by(year) stat(sum) format mi noheader describe `tabsave'
				`putexcel' `xlr'
			if "`putexcel'"~="" {
				local xlr `excelmodify'
			}
		}

		di _n as txt "----Present values & cost-effectiveness in `refyear' by the horizon year for country " as res "`Y'" in gr"----"
		di    as txt "        Cost variable: " as res "`costvar'" as txt "        Effect variable: " as res "`benvar'"
		foreach s of numlist 1/`ns' {
*			local sheetname PV_CE_Scen_`s'
			di  _n as txt "        --For scenario: " as res "`s': `s`s'_name'" as txt " --"
			tabstat2 pvbyh`s'_* cebyh`s'_* if code=="`Y'":code , by(year) stat(mean) format(%13.1fc) mi casewise nototal noheader describe save
*				`putexcel' `xlr'
			tempname cebyh`s'_cntry`Y'
			local nrowsofpv = rowsof(r(StatAll))
			local rownmsofpv : rownames r(StatAll)
			matrix rename r(StatAll) `cebyh`s'_cntry`Y''

			local pvcost   = `cebyh`s'_cntry`Y''[`nrowsofpv',colnumb(`cebyh`s'_cntry`Y'',"pvbyh`s'_`costvar'")]
			local pveffect = `cebyh`s'_cntry`Y''[`nrowsofpv',colnumb(`cebyh`s'_cntry`Y'',"pvbyh`s'_`benvar'")]

if "`debug'"~=""  {
	di _n(2) as err "DEBUG: Line 599:  Extract PV of costs and effects for ICER graph "
	di _n as inp "matrix list `cebyh`s'_cntry`Y''"
	matrix list `cebyh`s'_cntry`Y''
	di _n as inp `"local pvcost   = `cebyh`s'_cntry`Y''[`nrowsofpv',colnumb(`cebyh`s'_cntry`Y'',"pvbyh`s'_`costvar'")]"'
	di as txt "The value of local -pvcost- is: " as res "`pvcost'"
	di _n as inp `"local pveffect = `cebyh`s'_cntry`Y''[`nrowsofpv',colnumb(`cebyh`s'_cntry`Y'',"pvbyh`s'_`benvar'")]"'
	di as txt "The value of local -pveffect- is: " as res "`pveffect'"
}

			local pv_costlist  =  "`pv_costlist' "    + strofreal(scalar(`pvcost'),"%20.4g")
			local pv_effectlist  = "`pv_effectlist' " + strofreal(scalar(`pveffect'),"%20.4g")

			return matrix cebyh`s'_cntry`Y' `cebyh`s'_cntry`Y''
			local ce_matrixlist `ce_matrixlist' cebyh`s'_cntry`Y'
		}

		di _n as txt "----Cost-effectiveness in `refyear' over the entire horizon for country " as res "`Y'" as txt "----"
		di    as txt "        For all scenarios"
		qui tabstat2 ce?_* if year==`refyear' & code=="`Y'":code , by(year) stat(mean) format mi casewise nototal noheader describe save

		tempname ce_cntry`Y'
		matrix rename r(Stat1) `ce_cntry`Y''
			local nrows : word count `rownmsofpv'
			local lastrwnm = word("`rownmsofpv'", `nrows')
		matrix rownames `ce_cntry`Y'' = `lastrwnm'
		matrix list `ce_cntry`Y'', noheader
		return matrix ce_cntry`Y' `ce_cntry`Y''
		local ce_matrixlist `ce_matrixlist' ce_cntry`Y'

		di _n

		if "`gr7'"~="" | "`graph'"~="" {
			makegph `vnames'  ,  refyear(`refyear') ns(`ns') country(`Y')   weights(`weights')  ///
				costvar(`costvar') costunit(`costunit') benvar(`benvar') benunit(`benunit')  ///
				`gr7' `graph' `grnote' `grcombine' `grscenarios' `grdifferences' `grre' `grcebyh' `grce' ///
				`legend' `titlen' `titlec' `titlep' `titlee' `subtitle'  drate(`drate') ///
                bl_name(`bl_name') s1_name(`s1_name') s2_name(`s2_name') s3_name(`s3_name') s4_name(`s4_name') s5_name(`s5_name') s6_name(`s6_name') ///
				saving(`"`saving'"') `replace' `rescale' `weo' `scenes' `altshrink'  horizon(`horizon') `truncate' ///
				`grc1leg' `grc1leg_options' `gricer' pv_costlist(`pv_costlist') pv_effectlist(`pv_effectlist') `debug'

		}
	}
}

*	Projections for regions
if "`regions'"~="" {
	foreach X in `regions' {
		local pv_costlist
		local pv_effectlist
		local Y = upper("`X'")
		di _n as txt "Projections for Region " as res "`Y'"
		di    as txt "==========================="
		des `lols'

		if "`putexcel'"~="" & "`xlr'"=="" {
			local xlr `excelreplace'
		}
		foreach V in `vnames' {
			local sheetname Impact_on_`V'_for_`Y'
			di _n as txt "----Projected impact on variable " as res "`V'" in gr" for region " as res "`Y'" in gr"----"
			tabstat2 `vlist`V''    if rgn=="`Y'":rgn , by(year) stat(sum) format mi noheader describe `tabsave'
				`putexcel' `xlr'
			if "`putexcel'"~="" {
				local xlr `excelmodify'
			}
		}

		//	Use tabstat to make the matrices cebyh`s'_rgn`Y' and ce_rgn`Y'
		tempname ce_rgn`Y'
		local strtngce
		foreach s of numlist 1/`ns' {
			qui tabstat2 pvbyh`s'_* if rgn=="`Y'":rgn , by(year) stat(sum) format(%13.1fc) mi casewise nototal noheader describe save
			tempname pvbyh`s'_rgn`Y'
			matrix rename r(StatAll) `pvbyh`s'_rgn`Y''

			qui svmat `pvbyh`s'_rgn`Y'', names(matcol)
				local nrowsofpv = rowsof(`pvbyh`s'_rgn`Y'')
				local rownmsofpv : rownames `pvbyh`s'_rgn`Y''

			sum `pvbyh`s'_rgn`Y''pvbyh`s'_`costvar' if year==`horizon', meanonly
				tempname pvcost`s'_rgn`Y'
			scalar `pvcost`s'_rgn`Y'' = `r(mean)'

			sum `pvbyh`s'_rgn`Y''pvbyh`s'_`benvar' if year==`horizon', meanonly
				tempname pveffect`s'_rgn`Y'
			scalar `pveffect`s'_rgn`Y'' = `r(mean)'

			local pv_costlist    =  "`pv_costlist' "  + strofreal(scalar(`pvcost`s'_rgn`Y''),"%20.4g")
			local pv_effectlist  = "`pv_effectlist' " + strofreal(scalar(`pveffect`s'_rgn`Y''),"%20.4g")

				tempvar vcebyh`s'_rgn`Y'
			qui gen `vcebyh`s'_rgn`Y'' = `pvbyh`s'_rgn`Y''pvbyh`s'_`costvar' / `pvbyh`s'_rgn`Y''pvbyh`s'_`benvar'
				tempname cebyh`s'_rgn`Y'
			mkmat `pvbyh`s'_rgn`Y''pvbyh`s'_`costvar' `pvbyh`s'_rgn`Y''pvbyh`s'_`benvar' `vcebyh`s'_rgn`Y'' if _n <= `nrowsofpv', matrix(`cebyh`s'_rgn`Y'')
				drop `pvbyh`s'_rgn`Y''pvbyh`s'_`costvar' `pvbyh`s'_rgn`Y''pvbyh`s'_`benvar'
				matrix rownames `cebyh`s'_rgn`Y'' = `rownmsofpv'
					local nrows : word count `rownmsofpv'
				matrix colnames `cebyh`s'_rgn`Y'' = pvbyh`s'_`costvar' pvbyh`s'_`benvar' cebyh`s'
					local pv_col_nms : colnames `cebyh`s'_rgn`Y''
				local ce_col_nm = word("`pv_col_nms'",3)
				local ce_rgn_colnms `ce_rgn_colnms' `ce_col_nm'
			matrix define `ce_rgn`Y'' = `strtngce' `cebyh`s'_rgn`Y''[`nrows',3]
			local strtngce \`ce_rgn\`Y'' ,
		}

		di _n as txt "----Present values & cost-effectiveness in `refyear' by the horizon year for region " as res "`Y'" in gr"----"
		di    as txt "        Cost variable: " as res "`costvar'" as txt "        Effect variable: " as res "`benvar'"
*			local sheetname PV_CE_Scen_`s'
		foreach s of numlist 1/`ns' {
			di  _n as txt "        --For scenario: " as res "`s': `s`s'_name'" as txt " --"
			matrix list `cebyh`s'_rgn`Y'' , noheader
			return matrix cebyh`s'_rgn`Y' `cebyh`s'_rgn`Y''
			local ce_matrixlist `ce_matrixlist' cebyh`s'_rgn`Y'
		}

		di _n as txt "----Cost-effectiveness in `refyear' over the entire horizon for region " as res "`Y'" as txt "----"
		di    as txt "        For all scenarios"
		local lastrwnm = word("`rownmsofpv'", `nrows')
		matrix rownames `ce_rgn`Y'' = `lastrwnm'
		matrix colnames `ce_rgn`Y'' = `ce_rgn_colnms'
		matrix list `ce_rgn`Y'' , noheader
		return matrix ce_rgn`Y' `ce_rgn`Y''
		local ce_matrixlist `ce_matrixlist' ce_rgn`Y'

		di _n

		if "`gr7'"~="" | "`graph'"~="" {
			makegph `vnames'  ,  refyear(`refyear') ns(`ns') region(`Y')   weights(`weights') ///
				costvar(`costvar') costunit(`costunit') benvar(`benvar') benunit(`benunit') ///
				`gr7' `graph' `grnote' `grcombine' `grscenarios' `grdifferences' `grre' `grcebyh' `grce' ///
				`legend' `titlen' `titlec' `titlep' `titlee' `subtitle'  drate(`drate') ///
                bl_name(`bl_name') s1_name(`s1_name') s2_name(`s2_name') s3_name(`s3_name') s4_name(`s4_name') s5_name(`s5_name') s6_name(`s6_name') ///
				saving(`"`saving'"') `replace' `rescale' `weo' `scenes' `altshrink'  horizon(`horizon') `truncate' ///
				`grc1leg' `grc1leg_options' `gricer' pv_costlist("`pv_costlist'") pv_effectlist("`pv_effectlist'") `debug'
		}
	}
}

*	Projections:  Global totals
if "`world'" ~= "" {
	di _n as res "Global " as txt "projections"
	di    as txt "====================="
	sum code if year==`minvntg', meanonly
	di _n as txt "The " as res "`r(N)'" as txt " countries included in the projection"
	di    as txt "are distributed across the regions as follows:"
	tab rgn if year==`minvntg', mi

	if "`putexcel'"~="" & "`xlr'"=="" {
		local xlr `excelreplace'
	}
	foreach V in `vnames' {
			local sheetname Impact_on_`V'_for_wrld
		di _n as txt "----Projected impact on variable " as res "`V'" in gr" for all countries----"
		tabstat2 `vlist`V'' , by(year) stat(sum) format mi casewise noheader describe `tabsave'
			`putexcel' `xlr'
		if "`putexcel'"~="" {
			local xlr `excelmodify'
		}
	}

		//	Use tabstat to make the matrices cebyh`s'_wrld
		tempname ce_wrld
		local strtngce
		foreach s of numlist 1/`ns' {
			qui tabstat2 pvbyh`s'_* , by(year) stat(sum) format(%13.1fc) mi casewise nototal noheader describe save
			tempname pvbyh`s'_wrld
			matrix rename r(StatAll) `pvbyh`s'_wrld'

			qui svmat `pvbyh`s'_wrld', names(matcol)
				local nrowsofpv = rowsof(`pvbyh`s'_wrld')
				local rownmsofpv : rownames `pvbyh`s'_wrld'
				tempvar vcebyh`s'_wrld
			qui gen `vcebyh`s'_wrld' = `pvbyh`s'_wrld'pvbyh`s'_`costvar' / `pvbyh`s'_wrld'pvbyh`s'_`benvar'
				tempname cebyh`s'_wrld
			mkmat `pvbyh`s'_wrld'pvbyh`s'_`costvar' `pvbyh`s'_wrld'pvbyh`s'_`benvar' `vcebyh`s'_wrld' if _n <= `nrowsofpv', matrix(`cebyh`s'_wrld')
				drop `pvbyh`s'_wrld'pvbyh`s'_`costvar' `pvbyh`s'_wrld'pvbyh`s'_`benvar'
				matrix rownames `cebyh`s'_wrld' = `rownmsofpv'
					local nrows : word count `rownmsofpv'
				matrix colnames `cebyh`s'_wrld' = pvbyh`s'_`costvar' pvbyh`s'_`benvar' cebyh`s'
					local pv_col_nms : colnames `cebyh`s'_wrld'
				local ce_col_nm = word("`pv_col_nms'",3)
				local ce_wrld_colnms `ce_wrld_colnms' `ce_col_nm'
			matrix define `ce_wrld' = `strtngce' `cebyh`s'_wrld'[`nrows',3]
			local strtngce \`ce_wrld' ,
		}


		di _n as txt "----Present values & cost-effectiveness in `refyear' by the horizon year for all countries" as txt "----"
		di    as txt "        Cost variable: " as res "`costvar'" as txt "        Effect variable: " as res "`benvar'"
		foreach s of numlist 1/`ns' {
			di  _n as txt "        --For scenario: " as res "`s': `s`s'_name'" as txt " --"
			matrix list `cebyh`s'_wrld' , noheader
			return matrix cebyh`s'_wrld `cebyh`s'_wrld'
			local ce_matrixlist `ce_matrixlist' cebyh`s'_wrld
		}

		di _n as txt "----Cost-effectiveness in `refyear' over the entire horizon for all countries" as txt "----"
		di    as txt "        For all scenarios"
		local lastrwnm = word("`rownmsofpv'", `nrows')
		matrix rownames `ce_wrld' = `lastrwnm'
		matrix colnames `ce_wrld' = `ce_wrld_colnms'
		matrix list `ce_wrld' , noheader
		return matrix ce_wrld `ce_wrld'
		local ce_matrixlist `ce_matrixlist' ce_wrld

	di _n

	if "`gr7'"~="" | "`graph'"~="" {
		makegph `vnames'  ,  refyear(`refyear') ns(`ns') world  weights(`weights') ///
			costvar(`costvar') costunit(`costunit') benvar(`benvar') benunit(`benunit') ///
			`gr7' `graph' `grnote' `grcombine' `grscenarios' `grdifferences' `grre' `grcebyh' `grce' ///
			`legend' `titlen' `titlec' `titlep' `titlee' `subtitle' drate(`drate') ///
            bl_name(`bl_name') s1_name(`s1_name') s2_name(`s2_name') s3_name(`s3_name') s4_name(`s4_name') s5_name(`s5_name') s6_name(`s6_name') ///
			saving(`"`saving'"') `replace' `rescale' `weo' `scenes' `altshrink'  horizon(`horizon') `truncate' ///
			`grc1leg'  `grc1leg_options' `gricer' pv_costlist(`pv_costlist') pv_effectlist(`pv_effectlist') `debug'
	}

}

*return local pv_matrixlist `pv_world' `pv_matrixlist'  // These matrix names are not used
return local ce_matrixlist `ce_world' `ce_matrixlist'

end  /* End of program genout2 */


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 " as res "`X'" as error " 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 " as res "`X'" as error " is not in the data."
			exit 198
		}
	}
	return local regions = upper("`regions'")
}
end    /* End of program chkdta */

program define makegph
*	Construct, display & save one or more graphs for a specific country, region or the world.
	syntax  namelist (name=gphvars) , refyear(integer) ns(integer) drate(real) horizon(integer)  /*
		*/ [country(string) region(string) world  weights(string) /*
		*/ COSTVar(namelist) COSTUnit(string) BENVar(namelist) BENUnit(string) /*
		*/ gr7 graph grcombine grnote(string) legend(string) /*
		*/ GRScenarios GRDifferences grre grcebyh nogrce  /*
		*/ titles(string) titled(string) titler(string) titleb(string) titlec(string) titlev(string) subtitle(string) /*
		*/ bl_name(passthru) s1_name(passthru) s2_name(passthru) s3_name(passthru) s4_name(passthru) s5_name(passthru) s6_name(passthru) /*
		*/ saving(string) replace norescale weo altshrink noTRUNCATE /*
		*/ grc1leg_options(string) grc1leg gricer pv_costlist(string) pv_effectlist(string) debug]

	makelists `gphvars', ns(`ns') refyear(`refyear') benvar(`benvar') `bl_name' `s1_name' `s2_name' `s3_name' `s4_name' `s5_name' `s6_name' `debug'

if "`debug'"~=""  {
	di _n(2) as err "DEBUG: Line 874:  What lists are returned by -makelists-?"
	di as inp "return list"
	return list
	di _n(2) "DEBUG: Line 877: What variables are available to -makelists-?"
	di as inp "`gphvars'"
}

	foreach V in `gphvars' {  /* Transfer returned macros to local macros */
		foreach X in blist slist slbls dlist tslist tdlist {
			local `X'`V'  `r(`X'`V')'
		}
	}

	foreach X in blist slist slbls bkey skey skeyce dkey  ///
		border sorder sorderce dorder dlist rlist tslist tdlist trlist tcebyhlist  ///
		cebyhlist celist pvbyhlist cebarlist vtypelists `r(vtypelists)' {
		local `X' `r(`X')'
	}

	tempvar tmpordr
		gen long `tmpordr' = _n

	di _n as txt "Graph output follows:" _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 aidsdif:makegph should never occur."
		exit 198
	}

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

	if ~(lower("`legend'")=="embed")  ///
		& ~("`legend'"=="" | lower("`legend'")=="below" | "`legend'"=="position(6)" | "`legend'"=="pos(6)") {
		local legends  `legend'   /* For scenario graphs */
		local legendd  `legend'   /* For differences graphs */
		local legendr  `legend'   /* For ratio graphs */
		local legendb  `legend'   /* For CE bar graphs */
		local legendc  `legend'   /* For graphs to be combined */
	}
	if lower("`legend'")=="embed" {
		local legends  "position(4)  ring(0) stack textfirst style(column) size(small)"   /* For scenario graphs */
		local legendd  "position(4)  ring(0) stack textfirst style(column) size(small)"   /* For differences graphs */
		local legendr  "position(4)  ring(0) stack textfirst style(column) size(small)"   /* For ratio graphs */
		local legendb  "position(4)  ring(0) stack textfirst style(column) size(vsmall)"   /* For CE bar graphs */
		local legendc  "position(4)  ring(0) stack textfirst style(column) size(small)"   /* For graphs to be combined */
	}
	if "`legend'"=="" | lower("`legend'")=="below" | "`legend'"=="position(6)" | "`legend'"=="pos(6)" {
		local legends  "position(3) style(column) size(small) stack"   /* For scenario graphs */
		local legendd  "position(3) style(column) size(small) stack"   /* For differences graphs */
		local legendr  "position(6) style(column) size(small)"   /* For ratio graphs */
		local legendb  "position(6) style(column) size(vsmall)"   /* For CE bar graphs */
		local legendc  "position(6) style(column) size(small)"   /* For graphs to be combined */
	}

	if `ns'> 1 {
		local plrl "s"
	}

	if "`titles'"=="" {     /* For scenario graphs */
		local titles "Baseline & scenario`plrl'"
	}
	if "`titled'"=="" {     /* For differences graphs */
		local titled "Impact"
	}
	if "`titler'"=="" {     /* For ratio graphs */
		local titler "Annual ratio of costs to effects"
	}
	if "`titlecebyh'"=="" {     /* For CE by horizon graphs */
		local titlecebyh "CE in `refyear' by horizon at r=`drate'"
	}
	if "`titleb'"=="" {     /* For CE bar graphs */
		local titleb "Cost-effectiveness in `refyear' at r=`drate'"
	}
	if "`titlecba'"=="" {     /* For CBA bar graphs */
		local titlecba "Cost-benefit ratio in `refyear' at r=`drate'"
	}
	if "`titlev'"=="" {     /* Overall title for combined variable graphs */
		local titlev "Cost-effectiveness in `refyear' relative to baseline"
	}
	if "`titlei'"=="" {     /* Overall title for impact-only graphs */
		local titlei "Impact relative to baseline"
	}
	if "`titlec'"=="" {     /* Overall title for 4-way combined CE graphs */
		local titlec "Impact & cost-effectiveness relative to baseline"
	}

	local suffix "`country'`region'"
	if "`suffix'"~="" {
		local suffix = upper("`suffix'")
	}
	if "`suffix'"=="" {
		local suffix "all"
	}
	if "`subtitle'"~="" {
		local subtitle `"`subtitle': "'
	}
	if "`country'" ~= "" {
		local subtitle  "`subtitle'For country `suffix'"
		local subsfx
	}
	if "`region'"  ~= "" {
		local subtitle  "`subtitle'For region `suffix'"
		local subsfx
	}
	if "`world'"   ~= "" {
		local subtitle  "`subtitle'For all countries"
		local subsfx
	}

	if "`grnote'"=="" {
		local grnote "Compare scenario variables \`slist\`V'' to baseline"
	}

	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
	}

	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
		}

		//	Beginning with version 2.20, the regional CE ratio is no longer computed as a weighted average
		//	of the country-specific CE ratios.  While this version still uses the collapse command and
		//	computes weighted average CE ratios, it overwrites these weighted averages with the ratio of the
		//	regional "rawsum" of the present value of costs to the regional "rawsum" of the present value
		//	of benefits.  The collapse command here is simply used to reduce the dimensionality of the data.
		preserve

		//	Remove duplicate variables from variables to be collapsed ("vrbls2bc") into -rawsums-
		local vrbls2bc `blist' `slist' `dlist' `pvbyhlist'
		local vrbls2bc  : list uniq vrbls2bc

		qui keep if rgn=="`region'":rgn
		local dlbl : data label
		retainlbl rgn `rlist' `celist' `cebyhlist' `vrbls2bc', store
			gen byte nobs = 1
		collapse (mean) rgn `tmpordr' `rlist' `celist' `cebyhlist' (rawsum) `vrbls2bc' nobs [fweight=`weights'] , by(year)
			label data "`dlbl', total for region `region'"
			lab var nobs "Number of countries"
			retainlbl rgn `vrbls2bc' `rlist' `celist' `cebyhlist', restore  /* Restore the variable labels */
				lab values rgn rgn

			qui xtset rgn year
			local lead = `horizon' - `refyear'
			foreach s of numlist 1/`ns' {
				qui replace re`s'_`benvar' = dif`s'_`costvar'/dif`s'_`benvar'
				qui replace cebyh`s'_`benvar' = pvbyh`s'_`costvar'/pvbyh`s'_`benvar'
				qui replace ce`s'_`benvar' = F`lead'.cebyh`s'_`benvar' if year==`refyear'
			}
	}

	if "`world'"~="" {
		//	Beginning with version 2.20, the global CE ratio is no longer computed as a weighted average
		//	of the country-specific CE ratios.  While this version still uses the collapse command and
		//	computes weighted average CE ratios, it overwrites these computations with the ratio of the
		//	global "rawsum" of the present value of costs to the global "rawsum" of the present value
		//	of benefits.  The collapse command here is simply used to reduce the dimensionality of the data.
		preserve
		local dlbl : data label
		retainlbl rgn `blist' `slist' `dlist'  `rlist' `celist' `cebyhlist' `pvbyhlist', store    /* Added cebyhlist & pvbyhlist */
			gen byte nobs = 1
		collapse (mean) rgn `tmpordr' `rlist' `celist' `cebyhlist' (rawsum) `blist' `slist' `dlist' `pvbyhlist' nobs [fweight=`weights'] , by(year)   /* Added cebyhlist & pvbyhlist */
			label data "`dlbl', total for the world"
			lab var nobs "Number of countries"
			retainlbl rgn `rlist' `celist' `cebyhlist' `blist' `slist' `dlist' `pvbyhlist', restore  /* Restore the variable labels */

			qui tsset year
			local lead = `horizon' - `refyear'
			foreach s of numlist 1/`ns' {
				qui replace re`s'_`benvar' = dif`s'_`costvar'/dif`s'_`benvar'
				qui replace cebyh`s'_`benvar' = pvbyh`s'_`costvar'/pvbyh`s'_`benvar'
				qui replace ce`s'_`benvar' = F`lead'.cebyh`s'_`benvar' if year==`refyear'
			}
	}

	if `"`saving'"'~= "" {
		local svngs  `"saving(`"`saving'_s_`suffix'_\`V'.gph"' `rpl')"' // for scenario graphs
		local svngd  `"saving(`"`saving'_d_`suffix'_\`V'.gph"' `rpl')"' // for difference graphs
		local svngr  `"saving(`"`saving'_r_`suffix'.gph"' `rpl')"'      // for ratio graphs
		local svngbh `"saving(`"`saving'_cebyh_`suffix'.gph"' `rpl')"'  // for CE by horizon graphs
		local svngb  `"saving(`"`saving'_ce_`suffix'.gph"' `rpl')"'     // for CE bar graphs
		local svngv  `"saving(`"`saving'_`suffix'_\`V'.gph"' `rpl')"'   // for combined variable graphs
		local svng   `"saving(`"`saving'_`suffix'.gph"' `rpl')"'        // for combined CE graphs
		local svngi  `"saving(`"`saving'_icer_`suffix'.gph"' `rpl')"'   // for ICER quadrant line graphs
	}

	if "`rescale'"~="norescale" {  //  Start of rescale loop

		foreach list in `vtypelists'  {

			if strpos("`list'","cost") >0 {
				rescale ``list'' , auto already(1000)
				local unitlblc  "`r(unitlbl)'"
*				local already "Thousands"
				local costscale `r(unit)'
			}
			else if substr("`list'",1,2)=="d_"  {
				rescale ``list'' , auto
				local unitlbldth  unitlbldth("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else if substr("`list'",1,3)=="dy_" {
				rescale ``list'' , auto
				local unitlbldy  unitlbldy("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else if substr("`list'",1,3)=="ly_" {
				rescale ``list'' , auto
				local unitlblly  unitlblly("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else if substr("`list'",1,5)=="vsly_" {
				rescale ``list'' , auto
				local unitlblvsly  unitlblvsly("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else if substr("`list'",1,2)=="i_" {
				rescale ``list'' , auto
				local unitlbli  unitlbli("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else if substr("`list'",1,2)=="p_" {
				rescale ``list'' , auto
				local unitlblvsp  unitlblvsp("`r(unitlbl)'")
				local benscale `r(unit)'
			}
			else {
				rescale ``list'' , auto
				local unitlbln  unitlbln("`r(unitlbl)'")
				local benscale `r(unit)'
			}
		}

		local unitce = strofreal(round(real("`costscale'")/real("`benscale'"),0),"%15.0f")
		local unitlblce =  ///
			cond("`unitce'"=="1", " ",  ///
			cond("`unitce'"=="100", "Hundreds of ",  ///
			cond("`unitce'"=="1000", "Thousands of ",  ///
			cond("`unitce'"=="1000000", "Millions of ",  ///
			cond("`unitce'"=="1000000000", "Billions of ",  ///
			"Unknown")))))

	}  //  End of rescale loop
	else {
		local unitlblce "1000"
	}

	if "`gr7'"~="" {  // Faster, but lower quality, legacy graph engine
		if "`grcombine'"=="" { /* Produce individual gr7 graphs */
			foreach V in `gphvars' {
				titlev `V',  unitlblc(`unitlblc') `unitlbldth' `unitlbldy' `unitlblly' `unitlblvsly' `unitlbli' `unitlblvsp'  `unitlbln'
					local lbl_V = r(lbl_V)
					local lbl_Y = r(lbl_Y)
				if "`grscenarios'"~="" {
					gr7 bl_`V' `slist`V'' year if year>`refyear'-1, /*
						*/ xlabel ylabel c(llllll) /*
						*/ b1("`titles' for `lbl_V'") b2("`subtitle'") /*
						*/ `svngs' `bkey' `skey' gap(3) l1(`"`lbl_Y'"')
					if `"`saving'"'~="" {
						di as txt "(Version 7 graph file saved with option: " as res `"`svngs'"' as txt ")"
					}
					pause
				}
				if "`grdifferences'"~="" {
					local dmax = 0
					local dmin = 1
					foreach d of varlist `dlist`V'' {
						sum `d' if year>`horizon'-6 | year<`refyear'+6, meanonly
						local dmax = max(2*r(max), `dmax', 0.5 )
						local dmin = min(2*r(min), `dmin' )
					}
					if `dmax' > 0 & `dmin' <= 0  local dyline yline(0)
					gr7 `dlist`V'' year if year>`refyear'-1,  /*
						*/ xlabel ylabel c(llllll)  /*
						*/ b1("`titled' on `lbl_V'") b2("`subtitle'") /*
						*/ `svngd' `dkey' gap(3) l1(`"`lbl_Y'"') `dyline'
					if `"`saving'"'~="" {
						di as txt "(Version 7 graph file saved with option: " as res `"`svngd'"'  as txt ")"
					}
					pause
				}
			}
			/*  The following three graphs use the variables costvar and benvar */
			if "`grre'"~="" {
				local revmax = 0
				local revmin = 1
				foreach rev of varlist `rlist' {
					sum `rev' if year>`horizon'-6 | year<`refyear'+6, meanonly
					local revmax = max(2*r(max), `revmax', 0.5 )
					local revmin = min(0.5*r(min), `revmin' )
				}
				if `revmax' > 0 & `revmin' <= 0 {
					local ryline yline(0)
				}
				if "`truncate'"~="notruncate" {
					rescale `rlist', maxvalue(`revmax') minvalue(`revmin') none
					local truncnote "Ratios outside the range of the y-axis are hidden"
				}
				gr7 `rlist' year if year>`refyear'-1,  /*
					*/ xlabel ylabel c(llllll)  /*
					*/ b1("`titler'") b2("`subtitle'`subsfx'") gap(3) /*
					*/ l1(`"`unitlblc' of `costunit'"') l2(`"per `benunit'"')  /*
					*/ `svngr' `skeyce' `ryline'
				if "`truncate'"~="notruncate" {
					rescale `rlist', undo
				}
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as res `"`svngr'"'  as txt ")"
				}
				pause
			}
			if "`grcebyh'"~="" {
				local cevmax = 0
				local cevmin = 1
				foreach cev of varlist `cebyhlist' {
					sum `cev' if year>`horizon'-6 | year<`refyear'+6, meanonly
					local cevmax = max(2*r(max), `cevmax', 0.5 )
					local cevmin = min(2*r(min), `cevmin' )
				}
				if `cevmax' > 0 & `cevmin' <= 0  local ceyline yline(0)
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', maxvalue(`cevmax') minvalue(`cevmin') none
					local truncnote "Cost-effectiveness values outside the range of the y-axis are hidden"
				}
				gr7 `cebyhlist' year if year>`refyear'-1,  /*
					*/ xlabel ylabel c(llllll)  /*
					*/ b1("`titlecebyh'") b2("`subtitle'`subsfx'") gap(3)  /*
					*/ l1(`"`unitlblce'`costunit'"') l2(`"per `benunit'"')  /*
					*/ `svngbh' `skeyce' `ceyline'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', undo
				}
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as res `"`svngr'"'  as txt ")"
				}
				pause
			}
			if "`grce'"~="nogrce" {
				gr7 `celist' if year==`refyear', bar mean ylabel  /*
					*/ b1("`titleb'") b2("`subtitle'`subsfx'")  gap(3)  /*
					*/ l1(`"`unitlblce'`costunit'"') l2(`"per `benunit'"') `svngb'
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as res `"`svngb'"' as txt ")"
				}
				pause
			}
		}
		else {  /* Combine gr7 graphs */
			set textsize 150
			if "`grre'"=="" & "`grce'"=="nogrce" & "`grscenarios'"=="" & "`grdifferences'"=="" {
				di _n as error "Option -grcombine- selected without specifying any specific graphs to combine."
				exit 198
			}
			set graph off
			/*  The following three graphs use the variables costvar and benvar */
			if "`grre'"~="" {
				local revmax = 0
				local revmin = 1
				foreach rev of varlist `rlist' {
					sum `rev' if year>`horizon'-6 | year<`refyear'+6, meanonly
					local revmax = max(2*r(max), `revmax', 0.5 )
					local revmin = min(0.5*r(min), `revmin' )
				}
				if `revmax' > 0 & `revmin' <= 0  {
					local ryline yline(0)
				}
				if "`truncate'"~="notruncate" {
					rescale `rlist', maxvalue(`revmax') minvalue(`revmin') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Ratios outside the range of the y-axis are hidden"
					}
				}
				tempfile regph
				gr7 `rlist' year if year>`refyear'-1, xlabel ylabel  /*
					*/ c(llllll) saving(`"`regph'"') /*
					*/ b1("`titler'") gap(3) l1(`"`unitlblc' of `costunit'"') /*
					*/ l2(`"per `benunit'"') `skeyce' `ryline'
				if "`truncate'"~="notruncate" {
					rescale `rlist', undo
				}
				local gphlist = "`gphlist' " + "`regph'"
			}
			if "`grcebyh'"~="" {
				local cevmax = 0
				local cevmin = 1
				foreach cev of varlist `cebyhlist' {
					sum `cev' if year>`horizon'-6 | year<`refyear'+6, meanonly
					local cevmax = max(2*r(max), `cevmax', 0.5 )
					local cevmin = min(2*r(min), `cevmin' )
				}
				if `cevmax' > 0 & `cevmin' <= 0  local ceyline yline(0)
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', maxvalue(`cevmax') minvalue(`cevmin') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Cost-effectiveness values outside the range of the y-axis are hidden"
					}
				}
				tempfile cebyhgph
				gr7 `cebyhlist' year if year>`refyear'-1, xlabel ylabel  /*
					*/ c(llllll) saving(`"`cebyhgph'"')  /*
					*/ b1("`titlecebyh'") gap(3)  /*
					*/ l1(`"`unitlblce'`costunit'"') /*
					*/ l2(`"per `benunit'"') `skeyce'  `ceyline'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', undo
				}
				local gphlist = "`gphlist' " + "`cebyhgph'"
			}
			if "`grce'"~="nogrce" {
				tempfile cegph
				gr7 `celist' if year==`refyear', bar mean ylabel saving(`"`cegph'"') b1(`"`titleb'"') gap(3) l1(`"`unitlblce'`costunit'"') l2(`"per `benunit'"')
				local gphlist = "`gphlist' " + "`cegph'"
			}
			if "`grscenarios'"~="" | "`grdifferences'"~="" {
				foreach V in `gphvars' {
					titlev `V',  unitlblc(`unitlblc')  `unitlbldth' `unitlbldy' `unitlblly' `unitlblvsly' `unitlbli' `unitlblvsp'  `unitlbln'
						local lbl_V = r(lbl_V)
						local lbl_Y = r(lbl_Y)
					set graph off
					set textsize 150
					if "`grscenarios'"~="" {
						tempfile sns`V'
						gr7 bl_`V' `slist`V'' year if year>`refyear'-1, xlabel ylabel c(llllll) saving(`"`sns`V''"') b1(`"`titles' for `lbl_V'"') `bkey' `skey' gap(3) l1(`"`lbl_Y'"')
						local gphlist`V' = "`gphlist`V'' " + "`sns`V''"
					}
					if "`grdifferences'"~="" {
						tempfile dif`V'
						local dmax = 0
						local dmin = 1
						foreach d of varlist `dlist`V'' {
							sum `d' if year>`horizon'-6 | year<`refyear'+6, meanonly
							local dmax = max(2*r(max), `dmax', 0.5 )
							local dmin = min(2*r(min), `dmin' )
						}
						if `dmax' > 0 & `dmin' <= 0  local dyline yline(0)
						gr7 `dlist`V'' year        if year>`refyear'-1,  /*
							*/ xlabel ylabel c(llllll) saving(`"`dif`V''"')  /*
							*/ b1(`"`titled' on `lbl_V'"') `dkey' gap(3) l1(`"`lbl_Y'"') `dyline'
						local gphlist`V' = "`gphlist`V'' " + "`dif`V''"
					}
					set graph on
					set textsize 100
					local grnote2 = substr(`"`grnote'"',1,38)
					gr7 using `gphlist`V'', margin(5) `svngv' b1(`"`grnote2'"')  b2(`"`subtitle'"')
					if `"`saving'"'~="" {
						di as txt "(Version 7 graph file saved with option: " as res `"`svngv'"' as txt ")"
					}
					pause
				}
				if "`gphlist'"=="" {
					local ttl "`titlei'"
				}
				else {
					local ttl "`titlec'"
				}
				set graph on
				gr7 using `dif`costvar'' `dif`benvar'' `gphlist', margin(5) `svng' title(`"`ttl'"') b2("`subtitle'`subsfx'")
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as res `"`svng'"' as txt ")"
				}
				pause
			}
			else {
				set graph on
				set textsize 100
				gr7 using `gphlist', margin(5) `svng' title("`titlev'")  b2("`subtitle'`subsfx'")
				if `"`saving'"'~="" {
					di as txt "(Version 7 graph file saved with option: " as res `"`svng'"' as txt ")"
				}
				pause
			}
		}
	}
if "`debug'"~=""  {
	di _n(2) as err "DEBUG: Line 1340: The list of graph variables for the -graph- option is: " as txt "`gphvars'"
}

	if "`graph'"~="" {  // Slower, but higher quality, modern graph engine
		if "`grcombine'`grc1leg'`grc1leg_options'"=="" {  // Graphs are combined if either option is specified
			foreach V in `gphvars' {
				titlev `V',  unitlblc(`unitlblc') `unitlbldth' `unitlbldy' `unitlblly' `unitlblvsly' `unitlbli' `unitlblvsp'  `unitlbln'
					local lbl_V = r(lbl_V)
					local lbl_Y = r(lbl_Y)
				if "`grscenarios'"~="" {
					ylbl_line bl_`V' `slist`V''
					twoway (connected bl_`V' year if year>(2008-1), sort  lwidth(thick) lpattern(solid) ) `tslist`V'', /*
					  */ ytitle(`"`lbl_Y'"', margin(right)) ylabel(, angle(horizontal) `r(ylblfmt)' )  /*
					  */ title(`"`titles' for `lbl_V'"' ) /*
					  */ subtitle("`subtitle'") /*
					  */ xtitle("Year") /*
					  */ legend(order(`border' `sorder') `legends') /*
					  */ note(`"`grnote'"', span size(vsmall)) `svngs'
					pause
				}
				if "`grdifferences'"~="" {
					ylbl_line `dlist`V''
					twoway `tdlist`V'' , /*
					  */ ytitle(`"`lbl_Y'"', margin(right)) ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
					  */ title("`titled' on `lbl_V'")  /*
					  */ subtitle("`subtitle'") /*
					  */ xtitle("Year") /*
					  */ legend(order(`dorder') `legendd') /*
					  */ note(`"`grnote'"', span size(vsmall)) `svngd' `r(yline)'
					pause
				}
			}

			/*  The following three graphs use the variables costvar and benvar */
			if "`grre'"~="" {
				ylbl_line `rlist'
				if "`truncate'"~="notruncate" {
					rescale `rlist', maxvalue(`r(ymax)') minvalue(`r(ymin)') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Ratios outside the range of the y-axis are hidden"
					}
					ylbl_line `rlist'
				}
				twoway `trlist' , /*
				  */ ytitle(`"`unitlblc' of `costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
				  */ title("`titler'")  /*
				  */ subtitle("`subtitle'`subsfx'") /*
				  */ legend(order(`sorderce') `legendr') /*
				  */ note(`"`grnote'"' "`truncnote'", span size(vsmall)) `svngr' `r(yline)'
				if "`truncate'"~="notruncate" {
					rescale `rlist', undo
				}
				pause
			}
			if "`grcebyh'"~="" {
				ylbl_line `cebyhlist'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', maxvalue(`r(ymax)') minvalue(`r(ymin)') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Cost-effectiveness values outside the range of the y-axis are hidden"
					}
					ylbl_line `cebyhlist'
				}
				twoway `tcebyhlist' , /*
				  */ ytitle(`"`unitlblce'`costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
				  */ xtitle("Planning horizon")  /*
				  */ title("`titlecebyh'")  /*
				  */ subtitle("`subtitle'`subsfx'") /*
				  */ legend(order(`sorderce') `legendr') /*
				  */ note(`"`grnote'"' "`truncnote'", span size(vsmall)) `svngbh' `r(yline)'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', undo
				}
 				pause
			}
			if "`grce'"~="nogrce" {
				if "`benvar'"~="vsly"  {
					graph bar ce?_* ,  bargap(50) blabel(bar, format(%9.2f) size(medlarge)) /*
						*/ `cebarlist' /*
						*/ ytitle(`"`unitlblce'`costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal)) /*
						*/ title("`titleb'") name(ce_`suffix', replace) /*
						*/ subtitle("`subtitle'`subsfx'", margin(bottom))  /*
						*/ legend(order(`sorderce') `legendb') /*
						*/ note(`"`grnote'"', span size(vsmall)) `svngb'
				}
				else {
					graph bar ce?_* ,  bargap(50) blabel(bar, format(%9.2f) size(medlarge)) /*
						*/ `cebarlist' /*
						*/ ytitle("Cost-benefit ratio", margin(right))  ylabel(, angle(horizontal)) /*
						*/ title("`titlecba'") name(vsly_`suffix', replace)   /*
						*/ subtitle("`subtitle'`subsfx'", margin(bottom))  /*
						*/ legend(order(`sorderce') `legendb') /*
						*/ note(`"`grnote'"', span size(vsmall)) `svngb'
				}
				pause
			}
		}
		*****************************************************
		else {  /* Combine Stata's modern -graph- graphs */
			if "`grre'"=="" & "`grce'"=="nogrce" & "`grcebyh'" =="" & "`grscenarios'"=="" & "`grdifferences'"=="" {
				di _n as error "Option -grcombine- or -grc1leg- selected without specifying any specific graphs to combine."
				exit 198
			}


			if "`grc1leg'"~= "" | "`grc1leg_options'"~= "" {
				if "`grc1leg_options'"== "" {
					// Options to be used when grc1leg_options are not specified
					local grc_command grc1leg_aidscost
					local grc1args_for2 ring(1) pos(3) //  For 2-graph combinations
					local grc1args_for3 ring(0) pos(5)  //  For 3-graph combinations
					local grc1args_for4 ring(0) pos(5) holes(5 6) //  For 4-graph combinations
					local grc1args_for5 ring(0) pos(5)  //  For 5-graph combinations
					local grc1args_for6 ring(1) pos(6)  //  For 6-graph combinations
					local grc1args_for7 ring(0) pos(5) rows(3) holes(8 9) //  For 7-graph combinations
					local grc1args_for8 ring(0) pos(5) rows(3) //  For 8-graph combinations
				}
				grc1leg_opts , `grc1leg_options'
					local grc1_lf `r(legendfrom)'
					local grc1_pos `r(position)'
					local grc1_ring `r(ring)'
					local grc1_span `r(span)'
					local grcl_opts `r(options)'  //  The -gr combine- options
*					local grc1args `grc1_pos' `grc1_ring' `grc1_span' `grcl_opts'
			}
			else {
				local grc_command gr combine
				local rows2 rows(2)
			}

			set graph off
			//  The following three graphs use the variables costvar and benvar
			if "`grre'"~="" {
				tempname regph
				ylbl_line `rlist'
				if "`truncate'"~="notruncate" {
					rescale `rlist', maxvalue(`r(ymax)') minvalue(`r(ymin)') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Values outside the range of the y-axis are hidden"
					}
					ylbl_line `rlist'
				}
				twoway `trlist' , /*
				  */ ytitle(`"`unitlblc' of `costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
				  */ title("`titler'")  name(`regph', replace)  /*
				  */ xtitle(`"Year"') /*
				  */ legend(order(`sorderce') `legendc') `r(yline)'
				if "`truncate'"~="notruncate" {
					rescale `rlist', undo
				}
				local gphlist = "`gphlist' " + "`regph'"
			}
			if "`grcebyh'"~="" {
				tempname cebyhgph
				ylbl_line `cebyhlist'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', maxvalue(`r(ymax)') minvalue(`r(ymin)') none
					if `r(nsetmiss)' > 0 {
						local truncnote "Values outside the range of the y-axis are hidden"
					}
					ylbl_line `cebyhlist'
				}
				twoway `tcebyhlist' , /*
				  */ ytitle(`"`unitlblce'`costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
				  */ title("`titlecebyh'") name(`cebyhgph', replace)  /*
				  */ xtitle(`"Year of Horizon"') /*
				  */ legend(order(`sorderce') `legendc') `r(yline)'
				if "`truncate'"~="notruncate" {
					rescale `cebyhlist', undo
				}
				local gphlist = "`gphlist' " + "`cebyhgph'"
			}

			if "`grce'"~="nogrce" {
				tempname cegph
				if "`benvar'"~="vsly"  {
					graph bar ce?_* ,  bargap(50) blabel(bar, format(%9.2f) size(medlarge)) /*
					  */ `cebarlist' /*
					  */ ytitle(`"`unitlblce'`costunit'"' `"per `benunit'"'', margin(right))  ylabel(, angle(horizontal)) /*
					  */ title("`titleb'" ) name(`cegph', replace) /*
					  */ legend(order(`sorderce') `legendc')
				}
				else {
					graph bar ce?_* ,  bargap(50) blabel(bar, format(%9.2f) size(medlarge)) /*
						*/ `cebarlist' /*
						*/ ytitle("Cost-benefit ratio", margin(right))  ylabel(, angle(horizontal)) /*
						*/ title("`titlecba'") name(`cegph', replace)   /*
						*/ legend(order(`sorderce') `legendc')
				}
				local gphlist = "`gphlist' " + "`cegph'"
			}
			if "`grscenarios'"~="" | "`grdifferences'"~="" {
				foreach V in `gphvars' {
					titlev `V',  unitlblc(`unitlblc')  `unitlbldth' `unitlbldy' `unitlblly' `unitlblvsly' `unitlbli' `unitlblvsp'  `unitlbln'
						local lbl_V = r(lbl_V)
						local lbl_Y = r(lbl_Y)
					set graph off
					if "`grscenarios'"~="" {
						tempname sns`V'
						ylbl_line bl_`V' `slist`V''
						twoway (connected bl_`V' year if year>(2008-1), sort lwidth(thick) lpattern(solid) ) `tslist`V'', /*
						  */ ytitle(`"`lbl_Y'"', margin(right)) ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
						  */ title(`"`titles' for `lbl_V'"') /*
						  */ xtitle("Year") name(`sns`V'') /*
						  */ legend(order(`border' `sorder') `legends')  `r(yline)'
						local gphlist`V' = "`gphlist`V'' " + "`sns`V''"
					}

					if "`grdifferences'"~="" {
						tempname dif`V'
						ylbl_line `dlist`V''
						twoway `tdlist`V'' , /*
						  */ ytitle(`"`lbl_Y'"', margin(right)) ylabel(, angle(horizontal) `r(ylblfmt)' ) /*
						  */ title(`"`titled' on `lbl_V'"') /*
						  */ xtitle("Year") name(`dif`V'') /*
						  */ legend(order(`dorder') `legendd') `r(yline)'
						local gphlist`V' = "`gphlist`V'' " + "`dif`V''"
					}
					set graph on
					`grc_command'  ///
						`gphlist`V'', col(1)  ///
						note(`"`grnote'"',span) `svngv' `altshrink' subtitle(`"`subtitle'"', span size(vsmall)) ///
						title("Scenario projections and impact relative to baseline", span) ///
						`grc1args_for2'
					pause
				}
				if "`gphlist'"=="" {
					local ttl "`titlei'"
				}
				else {
					local ttl "`titlec'"
				}
				set graph on

				if "`grc_command'"=="grc1leg_aidscost" {
					local ngphs : word count `dif`costvar'' `dif`benvar'' `gphlist'
					local grc1args `grc1args_for`ngphs''
					//  Legend should always be from the last graph.
						local lastgph : word `ngphs' of `dif`costvar'' `dif`benvar'' `gphlist'
						local grcl_lf legendfrom(`lastgph')
				}

				`grc_command'  /*
				  */ `dif`costvar'' `dif`benvar'' `gphlist', `rows2' /*
				  */ title(`"`ttl'"') subtitle(`"`subtitle'"') /*
				  */ note(`"`grnote'"' "`truncnote'", span size(vsmall)) `svng' `altshrink' /*
				  */ `grcl_lf' `grc1args'
				pause

			}
			else {  // When neither -grscenarios- nor -grdifferences- is requested
				set graph on

				if "`grc_command'"=="grc1leg_aidscost" & "`grcl_lf'"=="" {
					local ngphs : word count `gphlist'
					local lastgph : word `ngphs' of `gphlist'
					local grcl_lf legendfrom(`lastgph')
				}

			`grc_command'  /*
				  */ `gphlist' ,  /*
				  */ title("`titlev'") subtitle(`"`subtitle'"') /*
				  */ note(`"`grnote'"' "`truncnote'", span size(vsmall)) `svng' `altshrink'  /*
				  */ `grcl_lf' `grc1args'
				pause
			}
		}

		//  Added in version 2.50  19Mar2020
		//  -aidsdif- does not combine this -icer- graph with other graphs

		if "`gricer'"~="" {

			local location =   ///
				cond("`world'"~="","world",   ///
				cond("`region'"~="",`"region("`region'")"',  ///
				cond("`country'"~="",`"country("`country'")"',"")))

			graph_icers, pv_costlist(`pv_costlist') pv_effectlist(`pv_effectlist')   ///
				`location' ns(`ns') benunit(`benunit') costunit(`costunit') costscale(`costscale') benscale(`benscale')  ///
				`svngi'
			pause

		}

	}

restore

end		/* End of program makegph */

program define graph_icers
*	Added in version 2.50  19Mar2020
*	-aidsdif- does not combine this -icer- graph with other graphs
	syntax , pv_costlist(string) pv_effectlist(string) ns(integer) [costscale(real 1000) benscale(real 1) costunit(string) benunit(string) world region(string) country(string) saving(passthru)]

	if "`world'`region'`country'"=="" {
		di as error "This error should never occur"
		exit 198
	}
	else {
		local location "`world'`region'`country'"
	}

	di _n(2) as txt "Constructing the ICER figure"  ///
		_n as txt "  The costs are " as res "`pv_costlist'"  ///
		_n as txt "  The effects are: " as res "`pv_effectlist'"  ///
		_n as txt "  For " as res "`ns'" as txt " scenarios for the location: " as res "`location'"

	local ncosts   : word count `pv_costlist'
	local neffects : word count `pv_effectlist'
	if `ncosts'~=`neffects' {
		di as err "ncosts~=neffects"
		exit 198
	}
	if `ncosts'~=`ns' {
		di as err "ncosts==neffects~=ns"
		exit 198
	}

	tempvar mlbl
	qui gen     `mlbl' = "BL"                                        if _n==1
	foreach s of numlist 1/`ns' {
		local sp1 = `s' + 1
		qui replace `mlbl' = "S`s'"                                  if _n==`sp1'
	    tempvar cst_icer_`s' eff_icer_`s'
		qui gen     `cst_icer_`s'' = 0                               if _n==1
		qui replace `cst_icer_`s'' = real(word("`pv_costlist'",`s')) if _n==`sp1'
			local cvarlist `cvarlist' `cst_icer_`s''
			lab var `cst_icer_`s'' "PV of future costs of scenario S`s'"
			char `cst_icer_`s''[varname] c_icer_`s'

		qui gen     `eff_icer_`s'' = 0                                 if _n==1
		qui replace `eff_icer_`s'' = real(word("`pv_effectlist'",`s')) if _n==`sp1'
			local evarlist `evarlist' `eff_icer_`s''
			local svarlist `svarlist' `cst_icer_`s'' `eff_icer_`s''
			lab var `eff_icer_`s'' "PV of future effects of scenario S`s'"
			char `eff_icer_`s''[varname] e_icer_`s'

		local ticerlist  `ticerlist' ///
			(line      `cst_icer_`s'' `eff_icer_`s''              , pstyle(p`sp1') lpattern(dash) lwidth(thick))   ///
			(scatter   `cst_icer_`s'' `eff_icer_`s'' if _n==`sp1' , pstyle(p`sp1') mlabel(`mlbl'))   ///
			(scatter   `cst_icer_`s'' `eff_icer_`s'' if _n==1     , mstyle(p1) mlabel(`mlbl') mlabpos(12) mlabstyle(p1))
		local lo = 2 + 3*(`s' - 1)
		local lgnd_order = `"`lgnd_order'"'  +  " `lo' " + `""S`s' vs BL""'
	}

	capture drop ___*

	//  Rescaling new variables, so previous rescale results don't apply
	rescale `cvarlist', auto already(1000)
		local costunitlbl  =lower("`r(unitlbl)'")

	rescale `evarlist', auto
		local benunitlbl = lower("`r(unitlbl)'")

	graph twoway `ticerlist' ,  ///
		yline(0, lcolor(black)) ylabel(,angle(hor))  ///
		title("Costs, effects and ICERs in `location'" "of `ns' scenarios relative to a baseline")  ///
		ytitle("Present value of future costs" "in `costunitlbl' of dollars")  ///
		xtitle("Present value of future effects" "in `benunitlbl' of `benunit's")  ///
		legend(order(`lgnd_order') size(small))  ///
		name(icer_for_`location', replace) `saving'

	rescale `cvarlist', undo
	rescale `evarlist', undo

end  /* End of program -graph_icers-  */

program define makelists, rclass
*	Conventions used for naming lists:
*		Names of lists for the following types of graphs:
*			Scenario graphs: slist
*			Difference graphs: dlist (NOT diflist)
*			Cost-effectiveness ratio graphs: celist
*			Other ratio graphs: rlist
*			Cost-effectiveness by horizon graphs: cebyhlist
*			Twoway graph definitions for a type of graph: Prefix -t- before list name
*		Lists of various versions of the same variable.
*			Benefit variable lists: dthlist, dylist, lylist, vslylist, ilist plist
*		Lists of death variables (i.e. for the variable -d-), use dthlist, not dlist
	version 9.2
	syntax namelist(name=avnames id="Impact analysis variables"), ns(integer) refyear(integer) benvar(namelist) 	///
	[bl_name(string) s1_name(string) s2_name(string) s3_name(string) s4_name(string) s5_name(string) s6_name(string) debug]

	tokenize `avnames'

    local ptrn1 solid  /* Line patterns */
    local ptrn2 dash
    local ptrn3 dash_dot
    local ptrn4 shortdash
    local ptrn5 shortdash_dot
    local ptrn6 longdash
    local ptrn7 longdash_dot

	local v_is_listed
	local v = 1
	while "``v''" ~= "" { // For each of the impact analysis variables in -avnames-

		local blist`v' = "bl_``v''"

		local slist`v'
		local slbls`v'
		local dlist`v'
		local tslist`v'
		local tdlist`v'

		local s = 1
		while `s'<= `ns' {  /* Loop over each scenario other than the baseline (i.e. S1, S2, ... ) */
			local sp1 = `s' + 1
			local sp2 = `s' + 2

			tempvar  S`s'var`v'     /* Generate symbols for labeling points */
			tempname S`s'lab`v'
			gen byte `S`s'var`v'' = 1
				lab def `S`s'lab`v'' 1 "S`s'"
				lab val `S`s'var`v'' `S`s'lab`v''

			/* List of scenario variables for variable -v- */
			local slist`v' = "`slist`v''" + " " + "s`s'_``v''"
			/* List of temporary variables for scenario labels for variable -v- */
			local slbls`v' = "`slbls`v''" + " " + "[`S`s'var`v'']"
			/* List of legend keys for -gr7- graphs for all variables */
			if `v'==1 & `s' < 4 {
				if `"`s`s'_name'"' == "" {
    				local skey    = `"`skey' "' + `"key`sp1'(s(T) c(l) p(`sp2') "Scenario `s'")"'
    				local skeyce  = `"`skeyce' "' + `"key`s'(s(T) c(l) p(`sp1') "Scenario `s'")"'
    				local dkey    = `"`dkey' "' + `"key`s'(s(T) c(l) p(`sp1') "S`s' - Baseline")"'
				}
				else {
    				local skey    = `"`skey' "' + `"key`sp1'(s(T) c(l) p(`sp2') `"`s`s'_name'"')"'
    				local skeyce  = `"`skeyce' "' + `"key`s'(s(T) c(l) p(`sp1') `"`s`s'_name'"')"'
    				local dkey    = `"`dkey' "' + `"key`s'(s(T) c(l) p(`sp1') `"`s`s'_name' vs. `bl_name'"')"'
				}
			}
			/* List of legend keys for -graph- graphs for all variables */
			if `v'==1 & `s' <= `ns' {
				if `"`s`s'_name'"' == "" {
				    local sorder    = `"`sorder' "'   + `"`sp1' "Scenario `s'""'
				    local sorderce  = `"`sorderce' "' + `"`s' "Scenario `s'""'
				    local dorder    = `"`dorder' "'   + `"`s' "S`s' - Baseline""'
				}
				else {
					chklength `s' `"`sorder' "' `"`sp1' `"`s`s'_name'"'"'
						local sorder    = `"`sorder' "'   + `"`sp1' `"`s`s'_name'"'"'
					chklength `s' `"`sorderce' "'  `"`s' `"`s`s'_name'"' "Relative to" `"`bl_name'"'"'
						local sorderce  = `"`sorderce' "' + `"`s' `"`s`s'_name'"' "Relative to" `"`bl_name'"'"'
					chklength `s' `"`dorder' "' `"`s' `"`s`s'_name'"' "Relative to" `"`bl_name'"'"'
						local dorder    = `"`dorder' "'   + `"`s' `"`s`s'_name'"' "Relative to" `"`bl_name'"'"'
				}
			}
			/* List of present value variables for variable -v- */
			chklength `s' "`pvbyhlist' "  "pvbyh`s'_``v''"
				local pvbyhlist ="`pvbyhlist' " + " " + "pvbyh`s'_``v''"
			/* List of difference variables for variable -v- */
			chklength `s' "`dlist`v''"  "dif`s'_``v''"
				local dlist`v' = "`dlist`v''" + " " + "dif`s'_``v''"
			/* List of twoway graph descriptors for scenario graphs */
			chklength `s' "`tslist`v''"  "(connected s`s'_``v'' year if year>(`refyear'-1), lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
				local tslist`v' = "`tslist`v''" + " " + "(connected s`s'_``v'' year if year>(`refyear'-1), lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
			/* List of twoway graph descriptors for difference graphs */
			chklength `s' "`tdlist`v''"  "(connected dif`s'_``v'' year if year>(`refyear'-1), lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
				local tdlist`v' = "`tdlist`v''" + " " + "(connected dif`s'_``v'' year if year>(`refyear'-1), lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"

			if "``v''"=="`benvar'" {
				chklength `s' "`rlist' "  "re`s'_``v''"
					local rlist = "`rlist' " + " " + "re`s'_``v''"
				chklength `s' "`trlist' "  "(connected re`s'_``v'' year if year>(`refyear'-1), cmissing(n) lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
					local trlist = "`trlist' " + " " + "(connected re`s'_``v'' year if year>(`refyear'-1), cmissing(n) lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
				chklength `s' "`cebyhlist' "  "cebyh`s'_``v''"
					local cebyhlist ="`cebyhlist' " + " " + "cebyh`s'_``v''"
				chklength `s' "`tcebyhlist' "  "(connected cebyh`s'_``v'' year if year>(`refyear'-1), cmissing(n) lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
					local tcebyhlist = "`tcebyhlist' " + " " + "(connected cebyh`s'_``v'' year if year>(`refyear'-1), cmissing(n) lpattern(`ptrn`sp1'') lstyle(p`sp1') mstyle(p`sp1'))"
				chklength `s' "`celist' "  "ce`s'_``v''"
					local celist ="`celist' " + " " + "ce`s'_``v''"
				chklength `s' "`cebarlist' "  "bar(`s', bstyle(p`sp1'bar))"
					local cebarlist ="`cebarlist' " + " " + "bar(`s', bstyle(p`sp1'bar))"
			}
			local s = `s' + 1

		}  //  End of the scenario loop

		return local blist``v''  `blist`v''
		return local slist``v''  `slist`v''
		return local slbls``v''  `slbls`v''
		return local dlist``v''  `dlist`v''
		return local tslist``v''  `tslist`v''
		return local tdlist``v''  `tdlist`v''
		local blist = "`blist' " + " " + "`blist`v''"
		local slist = "`slist' " + " " + "`slist`v''"
		local slbls = "`slbls' " + " " + "`slbls`v''"
		local dlist = "`dlist' " + " " + "`dlist`v''"
		local tslist = "`tslist' " + " " + "`tslist`v''"
		local tdlist = "`tdlist' " + " " + "`tdlist`v''"

		qui des *_``v'', varlist
		if "`r(varlist)'"~="" {
			local ``v''_list = "```v''_list' " + " " + "`blist`v'' " + " " + "`slist`v'' " + " " + "`dlist`v''"
			return local ``v''_list  ```v''_list'

			local vtypelists `vtypelists' ``v''_list
		}

		local v = `v' + 1

	}  // End of loop over the impact analysis variables in -avnames-

	capture des bl*d s*d dif*d, varlist
	local dthlist `r(varlist)'
	local ncostlist : list ncostlist - dthlist

if "`debug'"~=""  {
	di _n(2) as err "DEBUG: Line 1874: -ncostlist- is: " as res `"`ncostlist'"'  ///
		_n as err "and -lt1000list- is: " as res `"`lt1000list'"'
}

	return local blist `blist'
	if `"`bl_name'"'=="" {
	    return local bkey = `"key1(s(o) c(l) p(2) "Baseline")"'
	    return local border = `"1 "Baseline""'
	}
	else {
	    return local bkey = `"key1(s(o) c(l) p(2) `"`bl_name'"')"'
	    return local border = `"1 `"`bl_name'"'"'
	}

	return local slist `slist'
	return local slbls `slbls'
	return local skey  `skey'
	return local skeyce  `skeyce'
	return local dkey  `dkey'
	return local sorder  `sorder'
	return local sorderce  `sorderce'
	return local dorder  `dorder'
	return local dlist `dlist'
	return local rlist `rlist'
	return local tslist `tslist'
	return local tdlist `tdlist'
	return local trlist `trlist'
	return local vtypelists  `vtypelists'
	return local tcebyhlist `tcebyhlist'
	return local pvbyhlist `pvbyhlist'
	return local cebyhlist `cebyhlist'
	return local celist `celist'
	return local cebarlist `cebarlist'
	return local benvarlist `benvarlist'

end  /* End of program makelists */

program define titlev, rclass
	version 9.2
	syntax name (name=V) , unitlblc(string) [unitlbldth(string) unitlbldy(string) unitlblly(string) unitlblvsly(string) unitlbli(string) unitlblp(string) unitlbln(string) ]

	if strpos(lower("`V'"),"cost")  {
	    if "`V'"=="costtot" {
			return local lbl_V = "total ART costs"
		}
		else {
			return local lbl_V = "`V'"
		}
		return local lbl_Y = `"`unitlblc' of dollars"'
		exit
	}
	if "`V'"=="d" {
		return local lbl_V = "deaths"
		return local lbl_Y = `"`unitlbldth' of deaths"'
		exit
	}
	if "`V'"=="dy" {
		return local lbl_V = "death-years"
		return local lbl_Y = `"`unitlbldy' of death-years"'
		exit
	}
	if "`V'"=="ly" {
		return local lbl_V = "life-years"
		return local lbl_Y = `"`unitlblly' of life-years"'
		exit
	}
	if "`V'"=="vsly" {
		return local lbl_V = "dollar value of life-years"
		return local lbl_Y = `"`unitlblvsly' of dollars"'
		exit
	}
	if "`V'"=="i" {
		return local lbl_V = "new infections"
		return local lbl_Y = `"`unitlbli' of new infections"'
		exit
	}
	if "`V'"=="p" {
		return local lbl_V = "HIV+ population 15+"
		return local lbl_Y = `"`unitlblp' of PLWHA"'
		exit
	}
	if "`V'"=="cvrg2" {
		return local lbl_V = "proportion to 2nd line ART"
		return local lbl_Y = `"Proportion to 2nd line"'
		exit
	}
	if "`V'"=="ccvrg" {
		return local lbl_V = "proportion circumcised"
		return local lbl_Y = `"Proportion circumcised"'
		exit
	}
	else {  // Default axis, title labels
		return local lbl_V = "`V'"
		return local lbl_Y = `"`unitlbln' of `V'"'
	}
end  /* End of program titlev */

program define retainlbl
*   version 1.0                         by Mead Over              3/10/08
*	Preserves and restores variable labels, for example before and after execution of -collapse-
*	Is not relevant for value labels.  Requires that the variable name is unchanged.
	version 9.2

	syntax varlist , [store restore]

	if ("`store'"=="store")+("`restore'"=="restore") ~= 1{
		di as error "Program retainlbl requires one and only one of the options -store- and -restore-"
		error 198
	}
	if "`store'"=="store" {
		foreach V of varlist `varlist' {
			global vlbl_`V' : variable lab `V'
		}
	}
	if "`restore'"=="restore" {
		foreach V of varlist `varlist' {
			if `"${vlbl_`V'}"' == "" {
				di _n as error "Label has not been stored for variable: " as res "`V'"
				di    as error "and therefore cannot be -re-stored."
				error 198
			}
			local lbl : variable lab `V'
			local parenplace = strpos("`lbl'",")")
			if `parenplace'~= 0 {
				local prfx = substr("`lbl'",1,`parenplace')
				char define `V'[collapse_type] = "`prfx'"
			}
			label var `V' `"${vlbl_`V'}"'
			global vlbl_`V' /* Clean out the globals when they are no longer needed. */
		}
	}
end /* End of program retainlbl */

prog def isflow
	version 9.2

	syntax , costvar(varname) benvar(varname)

	local stockflowc :  char `costvar'[stockflow]
	if "`stockflowc'"~="flow" {
		di _n as error "Warning: Variable specified by -costvar()- option, " as res "`costvar',"
		di    as error "is not a -flow- variable.  Execution continues."
*		error 198
	}

	local stockflowb :  char `benvar'[stockflow]
	if "`stockflowb'"~="flow" {
		di _n as error "Warning: Variable specified by -benvar()- option, " as res "`benvar',"
		di    as error "is not a -flow- variable.  Execution continues."
*		error 198
	}

end  /* End of subroutine isflow */

program define keepwts, rclass
    version 9.2
*   This subroutine reads the weight expression and determines which variables need to be kept
*     to satisfy it.  Written 1/30/09

    syntax , weights(string) [confirm]
    tokenize `weights', parse(" +-/*()^")

    while "`1'"~="" {  /* Go thru the tokens one at a time */

        local anymath
        foreach X in + - / * ( ) ^ int round abs floor ceil max min { /* Check if it is a math symbol */
            if "`X'"=="`1'" {
                local anymath "yes"
            }
        }

        local isnmbr                 /* Check if it is a number */
        capture confirm number `1'
        if _rc == 0 {
            local isnmbr "yes"
        }
        if "`anymath'`isnmbr'"=="" { /* If not a math symbol or number, strip of aidsproj prefixes */
            local unmodvar = 1
            local oname: subinstr local 1  "bl_" "",  count(local ishere)
            if `ishere' {
                local wtvlst = "`wtvlst'" + " `oname'"
                local unmodvar = 0
            }
            else {
                foreach s of numlist  1/6    {
                    local oname: subinstr local 1  "s`s'_" "",  count(local ishere)
                    if `ishere' {
                        local wtvlst = "`wtvlst'" + " `oname'"
                        local unmodvar = 0
                    }
                }
            }
            if `unmodvar' {
                local wtvlst = "`wtvlst'" + " `1'"
            }
        }
        macro shift
    }
    local wtvlst : list uniq wtvlst

    if "`confirm'" ~= "" {
        foreach X in `wtvlst' {
            capture confirm variable `X'
            if _rc~=0 {
                di _n as error "The variable, " as res "`X'" as error " is used in the weights expression, but does not exist in the baseline data."
                exit 198
            }
        }
    }
    return local wtvlst   `wtvlst'

end /* End of subroutine keepwts */

prog def chkmoresign, rclass
    args difv
    local morewhat : char `difv'[morebadgood]
    return local morebadgood `morewhat'
end /* End of subroutine chkmoresign */

prog def chklength
    args snmbr macro1 macro2
    local len1 : length local macro1
    local len2 : length local macro2
    local newlength = `len1' + `len2'
    if `newlength'>c(maxstrvarlen)-2 {
		local sver = c(stata_version)
		di as err _n "You are currently running version " as res "`sver'" as err " of Stata and you have asked -aidsdif-"  ///
				_n "to compare a baseline scenario to " as res "`snmbr'" as err " or more alternative scenarios."
		di    as err "With this version of Stata, -aidsdif- cannot handle this many scenarios and scenario names."  ///
				_n "If you don't have access to a more recent version of Stata,"  ///
				_n "try shortening or omiting the scenario names.  If you get this same error,"  ///
				_n "you must specify fewer scenarios."
		exit 1000
	}


end   /* End of subroutine chklength */

program define makedyly
//	Create the variables -dy- and -ly- by cumulating annual deaths
//	In Version 2.13, annual deaths are not discounted before cumulation
//	However, death-years averted, defined as the difference between the
//	variable -dy- in two different scenarios, will be discounted in main program.
	syntax , refyear(integer)
	by code: gen dy = sum(d) if year >=`refyear'
		lab var dy "Years of death"
		char define dy[stockflow] flow
		char define dy[scale] 1
		char define dy[morebadgood] moreisbad
		char define dy[unit] death-year
	by code: gen ly = -dy   if year >=`refyear'
		lab var ly "Years of life (Negative of years of death)"
		char define ly[stockflow] flow
		char define ly[scale] 1
		char define ly[morebadgood] moreisgood
		char define ly[unit] life-year
end		/*	End of subroutine makedyly */

program define makevsly
	syntax , refyear(integer)
	by code: gen dy = sum(d) if year >=`refyear'
		lab var dy "Years of death"
		char define dy[stockflow] flow
		char define dy[scale] 1
		char define dy[morebadgood] moreisbad
		char define dy[unit] death-year
	by code: gen ly = -dy   if year >=`refyear'
		lab var ly "Years of life (Negative of years of death)"
		char define ly[stockflow] flow
		char define ly[scale] 1
		char define ly[morebadgood] moreisgood
		char define ly[unit] life-year
	by code: gen vsly = ly*gnipc_hat/1000   if year >=`refyear'
		lab var vsly "Value of years saved ($1,000's)"
		char define vsly[stockflow] flow
		char define vsly[scale] 1000
		char define vsly[morebadgood] moreisgood
		char define vsly[unit] life-year
end		/*	End of subroutine makevsly */

program define pvbyhorizon
*	Compute the present value of the future stream of the variable -var-
*	up through each year of the projection period.  In the last year of the
*	projection period, -horizon-, this variable converges to the value
*	of -var-'s present value over the entire planning horizon.

	syntax , var(name) scenario(integer) refyear(integer) horizon(integer) drate(real)

	local s `scenario'
	local v `var'

	qui gen double pvbyh`s'_`v' = .
	forvalues hor = `horizon'(-1)`refyear'  {
		tempvar pvbyh`hor'_`s'_`v'
		qui by code (year): egen `pvbyh`hor'_`s'_`v'' = total(vrefyr * dif`s'_`v') if year>`refyear'-1 & year<`hor'+1
			lab var `pvbyh`hor'_`s'_`v'' "PV impact of S`s' on `v' through horizon year `hor'"
			qui replace `pvbyh`hor'_`s'_`v'' =. if year > `refyear'

*	The following command should do the functions of the next 5.  Why doesn't it work?
*		by code (year): replace pvbyh`s'_`v' = `pvbyh`hor'_`s'_`v''[1]  if year==`hor'

		qui levelsof code
		foreach c of numlist `r(levels)' {
			sum `pvbyh`hor'_`s'_`v'' if year==`refyear' & code ==`c', meanonly
			qui replace pvbyh`s'_`v' = `r(mean)'  if year==`hor' & code==`c'
		}
	}
	lab var pvbyh`s'_`v' "PV of impact of S`s' on `v' by horizon at r = `drate'"
	local thisscale : char dif`s'_`v'[scale]
		char define pvbyh`s'_`v'[scale] `thisscale'
	local thisbadgood : char dif`s'_`v'[morebadgood]
		char define pvbyh`s'_`v'[morebadgood] `thisbadgood'

end		/* End of subroutine -pvbyhorizon-  */

program define grc1leg_opts, rclass
	syntax [anything] [, LEGendfrom(string)    ///
		POSition(string) ring(integer -1) span     ///
		name(passthru) SAVing(string asis) * ]

	return local legendfrom legendrom(`legendfrom')
	return local position position(`position')
	return local ring ring(`ring')
	return local span `span'
	return local options `options'

end		/* End of subroutine -grc1leg_opts-	*/

program define ylbl_line, rclass
	//  Used only for -graph-, not -gr7- graphs
	//  Version 2.41 15Feb2020
	//  In Version 2.50 26Feb2020, disable -ylblfrmt-'s formatting in favor of version 3.1 of -rescale-
	//  However, since -rescale- sets the same format on all -ncostlist- variables,
	//  consider using format set by ylbl_line to set format for an individual graph.
	//  Still used to set -yline-
	syntax varlist
	local ymax = 0
	local ymin = 1

	foreach v of varlist `varlist' {
		sum `v' if year>`horizon'-6 | year<`refyear'+6, meanonly
		local ymax = max(2*r(max), `ymax', 0.5 )
		local ymin = min(2*r(min), `ymin' )
	}
	if `ymax' > 0 & `ymin' <= 0  local yline yline(0, lcolor(black))
	if abs(`ymax' - `ymin')  <= 2  local ylblfmt format(%9.1f)
	if abs(`ymax' - `ymin')  <= .2  local ylblfmt format(%9.2f)
	if abs(`ymax' - `ymin')  <= .02  local ylblfmt format(%9.3f)
	if abs(`ymax' - `ymin')  <= .002  local ylblfmt format(%9.4f)
	if abs(`ymax' - `ymin')  <= .0002  local ylblfmt format(%9.5f)

	return local ylblfmt `ylblfmt'
	return local yline `yline'
	return local ymax `ymax'
	return local ymin `ymin'
end

*   Version 1.20 adds the -weights- option
*   Version 1.21 replaces the STB -mmerge- command with STATA's -merge- command
*   Version 1.22 uses the -saving- name to leave behind the -aidsdif- data file as well as the graph files
*   Version 1.30 addresses problems with file saving
*   Version 1.40 adds the user-designated scenario names to graph legends
*   Version 1.41 modifies graph subtitles to display "weighted average" notation only where appropriate
*   Version 1.42 replaces calls to tabstat with calls to tabstat2 in order to improve table output
*   Version 1.43 varies line patterns in connected graphs
*   Version 1.44 adds the -noaltshrink- undocumented option, which enlarges text in combined graphs.
*   Version 1.45 removes the diagnostic print of the -altshrink- macro, adds ($1,000's) to label for ce_s*,
*       restricts -altshrink- to 2 or more scenarios and pads margin for CE bar graphs.
*   Version 1.46 fixes bug that prevented combining gr7 graphs and traps mistaken specification of more than one variable
*       in costvar and benvar options and removes abs() from calculation of re?_x and ce?_x variables.
*   Version 1.50 adds the subroutine -chklength- to provide an error trap for long scenario descriptions.
*	Version 2.0 changes:
*		Changed to accept death-years or life-years instead of deaths as the benefit variable.
*		Changed to an -rclass- program so that it can leave behind matrices of costs and benefits.
*		Additional tabstat commands added in -genout2- to compute matrices of accumulated benefits and costs
*		The subroutine -titlev- now accepts death-years (dy) or life-years(ly) as alternative benefit measures
*	Version 2.1 changes:
*		Added lines 395-411 modeled on lines 450-466, so that analysis by -country- also generates ce matrices
*	Version 2.12 changes:
*		Added lines 507-523 modeled on lines 450-466, so that analysis by -world- also generates ce matrices
*	Version 2.13 chnages:
*		Lines 64-77: Change default -benvar- to -ly-.  Lines 124-6, 175-7 create the variables -dy- and -ly-
*		which previously had to be created by user after running -aidsproj- and before running aidsdif.
*	Version 2.14 changes:
*		Removed offending extra "}" at end of program makedyly
*	Version 2.15 21feb2014: Fix saving logic to properly handle the -replace- option in lines 22-24
*		Remove unused -drate()- option from subroutine -makedyly-
*		Improve error trapping on the -refyear- option
*	Version 2.16:
*		Under grcombine option, move graph subtitle from component graphs to the combined graph
*		Add legend options: embed and below (the default).  Also user may specify his own legend options, such as "off"
*		Add graph type option grcebyh to graph the cost-effectiveness ratios as a function of the horizon.
*		Truncate the y-axis of the re and the cebyh graphs to twice the max or min of last 5 or first 5 values,
*		unless the user specifies the noTRUNCATE option.
*	Version 2.20 :
*		Change calculation of cost-effectiveness ratios for -regions- and -world- to define the CE
*		as the ratio of the present values of the raw sum of total costs to the raw sum of total benefits
*		rather than as weighted averages of ratios, which is now the case.  Not relevant for -countries-.
*		Add the -putexcel- option, but it only works for Stata 13 and above.
*	Version 2.21 :
*		Use Vince Wiggins' -grc1leg- instead of Stata's -gr combine-.
*		Add the functionality of the grc1leg options: legendfrom, position, ring, span
*		Make sure yline(0) is not specified when 0 is below the range of y
*	Version 2.22 (1 June 2014):
*		Minor bug fixes for -grc1leg-
*	Version 2.30 (22 Feb 2015): Fix bug in subroutine dyly to set the -morebadgood- characteristic for -dy-
*		Relax constraint that benefit variable must be a -flow- variable, allowing ART person-years to be a benefit.
*		Edit the line "local revmin = min(2*r(min), `revmin' ), changing the 2*r(min) to 0.5*r(min)
*		Relax the requirement that ly or dy be listed as part of the varlist to be kept if one of them is a benvar
*	Version 2.31 (26 Feb 2015): Generalize the -chklength- program to apply to all of the lists in program -makelists-,
*		and to be flexible to the Stata version.
*	Version 2.32 (11 May 2015): Investigate the behavior of dy and ly beginning in line 150.  The program "makedyly"
*		begins in line #1783.  Why do I want to cumulate deaths to make "death years"?  Why not just subtract
*		the population.  Reason is that in this model the population does not decline when the number of AIDS deaths
*		increases. Change size option in local legendb to size(small). Elaborate the default units when -benunit()- is blank.
*		Change the xtitle for the cebyh graphs to "Year of Horizon".
*		Add the margin(right) option to the ytitle(...) specs.  Adds an artifact minus sign.
*	Version 2.33 (25 Mar 2016): Program -grc1leg- was mistakenly required.  Fix this by replacing
*		one remaining occurrence of grc1leg with grc1led_aidscost (line 1430)
*		Pull the -size- option from the grc1leg option.
*	Version 2.34 (12 Apr 2016): Reverse the version 1.44 decision to make altshrink the default.
*		If the option is unspecified, set it to [altshrink/noaltshrink/ if [grcombine/grc1leg] is selected.
*		Replace _ly with _`benvar' so that computations can be made with different benefit variables.
*	Version 2.35 (10 Oct 2017): Change "impact of scenarios on" to "impact on" in panel titles.
*		(Allows full title to display on these panels with options -grc-, grc1leg- and without -altshrink-.)
*	Version 2.36 (15 Feb 2018):  To allow path designation, specify the .gph extension in the macros `svng'.
*	Version 2.37 (6 Feb 2019): Under construction: Add -outpath- & -inpath- options
*	Version 2.38 (6 Feb 2019): Under construction: Allow the -costvar- option to specify a cost other than -costtot-
*	Version 2.39 (4 Dec 2019): Allow a benefit variable which is a multiple of -ly * gnipc- in line 1914
*	Version 2.40 (6Feb2019): Change vertical axis label on CE graph to read cost-benefit ratio when benvar is -vsly-
*	Version 2.41 (15Feb2019): Condition format of dif_* variables on range between minimum & maximum.
*		Add subroutine -ylbl_line- to save repeated code. Fix -makegph- to properly color-code the scenarios
*	Version 2.50 (19Mar2020) Add a plot of the ICER, here it's for the -graph- option only, not yet the -gr7-.
*		Replace -yline(0)- with -yline(0, lcolor(black))-.  Fix titles and decimal places on graphs.
*		Requires version 3.1 of -rescale-.
*	Version 2.51 (8Dec2020) In defning the saving macros, line 1075-82, enclose assembled graph name in double quotes.