program define rayplot, rclass
	version 11.1
*!	program rayplot		version 2.4 7May2019		by Mead Over

	syntax [varlist (fv default=none)] [if] [in],  ///
		ZVar(string) YVar(string) XVARs(varlist) XVALues(numlist) SLopes(numlist) ///
		[chain(passthru) chainlhs(passthru) chainwrt(passthru) ///
		restore(string) cmdline(string)  note(passthru) ///
		intercept(real 0) level(real 95) verbose replace force ///
		yweight(real 1) GPHPREfix(string) saving(string) keep  ///
		eqname(passthru)  ///  Added July, 2016
		xtitle1(string) xtitle2(string)  ///  Added Feb, 2017
		TCYTitle(string) RESYTitle(string) ACYTitle(string) MCYTitle(string) ACMCYTitle(string) ///
		ytitle2(string) SUBTitle(string) SE  ///  Add SE option 13May2018
		boxcox(real 0.5) DEStination(string asis)  /// Added Oct, 2017 to enable -boxcox- transformation of outputs
		CONTstat(string) FVstat(string) MARGOPTions(string)  /// Added Jan2018 along with -marglist-
		yrescale  old resnote ]

if "`old'"~="" {
	// Unnecessarily lengthened gph filenames
	local _along_ray_  _along_ray_
	local _along_ray   _along_ray
}

marksample touse
	tempvar esmpl
	gen byte `esmpl' = e(sample)
	markout `touse' `esmpl'

local marglist `varlist'

if "`marglist'"~="" {
	fvrevar `varlist'
	if wordcount("`r(varlist)'") > 1 {
		di as err "Unlike the -margins- command, this program, -rayplot-, can only"   ///
			_n as err "compute and graph a single projection.  Thus the margin specification"  ///
			_n as err "provided as an argument may not contain the general term -i.-."  ///
			_n as err "Here the argument provided is: `marglist'"
		exit 198
	}
	drop `r(varlist)'
	// Enable replacing multi-level factor variables in the command liine
	// with specific levels of those variables from -marglist-
	tokenize `marglist' , parse(" #")
	while "`1'"~="" {
		if "`1'"~="#" {
			if strpos("`1'",".")==0 {
				di as err "The margin specification can contain only level-specific factor variables."   ///
					_n as err "Here the specification is: `marglist'"
				exit 198
			}
			local op = substr("`1'",1,strpos("`1'",".")-1)
			local v = substr("`1'",strpos("`1'",".")+1,.)
			local ops `ops' `op'  // List of specific levels
			local vs `vs'   `v'  // List of associated variable names
		}
		macro shift
	}
}

if "`verbose'"~="verbose" {
	local shhh quietly
}

*	Suppress -margins- computation of the standard errors and upper and lower bounds
if "`se'" == "se" {
	di _n as txt "Specifying the -se- option increases -rayplot-'s execution time."
	local nose
}
else {
	di as txt "Omitting the -se- option suppresses -rayplot-'s calculation and display of forecast/confidence intervals."
	local nose nose
}

*	Optionally restore previous estimation results
if "`restore'"~="" & "`cmdline'"~= "" {
	di as err _n "The options -restore- and -cmdline- are substitutes for one another." ///
		_n "Use one or another or, if the appropriate information is in -ereturn-, neither."
	exit 198
}

tempvar esample
if "`restore'`cmdline'"== "" {
	di as txt _n "Using retained estimation results"
	if "`e(cmd)'"=="" {
		di as err _n "No estimation results are retained in memory." ///
			_n "Use either the -restore()- or the -cmdline()- option."
		exit 198
	}

	cleancmdline, cmdline("`e(cmdline)'")
		local depvar `r(depvar)'
		local estcmd `r(estcmd)'
		if "`eqname'"=="" {
			local eqname eqname(`r(eqname)')
		}
		local cmdline `r(cmdline)'
		local origcmdline `r(origcmdline)'
		gen byte `esample' = e(sample)
		markout `esample'
}
else if "`restore'"~="" & "`cmdline'"== "" {
*	Retrieve information from preceding estimation command
		est restore `restore'  // e.g. ols_sqrt
		local cmdline `e(cmdline)'
		local depvar `e(depvar)'
		gen byte `esample' = e(sample)
		markout `esample'
}
else  { // If we arrive here, only the option -cmdline()- is non-blank
	cleancmdline, cmdline(`"`cmdline'"')
		local depvar `r(depvar)'
		local estcmd `r(estcmd)'
		if "`eqname'"=="" {
			local eqname eqname(`r(eqname)')
		}
		local cmdline `r(cmdline)'
		local origcmdline `r(origcmdline)'
		gen byte `esample' = e(sample)
		markout `esample'
}

*	For each multi-leveled discrete factor variable i.x,
*	substitute the specific level of the variable given in -marglist-

if "`marglist'"~="" {
	tokenize `"`origcmdline'"', parse(" #")
	while "`1'"~="" {
*		di as inp "`1'"
		if strpos("`1'",".") & substr("`1'",1,1)=="i" {
			local k = 1
			local replaced
			foreach v in `vs' {
				if "`v'"==substr("`1'",strpos("`1'",".")+1,.) {
					local op_i = word("`ops'",`k')
					local new1  `op_i'.`v'
					local newcmdline `newcmdline' `new1'
					local replaced replaced
				}
				local k = `k' + 1
			}
			if "`replaced'"=="" {
				local newcmdline `newcmdline' `1'
			}
		}
		else {
			local newcmdline `newcmdline' `1'
		}
		macro shift
	}
	local newcmdline = subinstr("`newcmdline'"," # ","#",.)
	local origcmdline `newcmdline'
}

*	Define the local -xvar- as the first variable in the option -xvars()-
local xvar = word("`xvars'",1)

*	Define -slope- as the first value in the option -slopes(numlist)-
*		********* Need to generalize to allow three or more outputs **********
local slope = real(word("`slopes'",1))

*	Set the transformed versions of these names
transid  `xvar', `chain' `chainwrt'
	local transx `r(trans)'
	local txvar  `r(transvar)'
transid  `yvar', `chain' `chainwrt'
	local transy `r(trans)'
	local tyvar  `r(transvar)'
transid  `zvar', `chain' `chainlhs'
	local transz `r(trans)'
	local tzvar  `r(transvar)'

*	Check to see that the dependent variable in estimates in memory
*	has a prefix indicating that it was transformed as specified here.
	chktrans, varname(`depvar') trans(`transz') `force'

if "`transy'"=="ln" & `slope'==0 & `intercept'==0 {
	di as err "With the logarithmic chain function, -ln-, either the ray intercept" ///
		_n "or the ray slope or both must be positive."
	exit 198
}

*	Generate predicted values of the dependent variable (e.g. total cost) along the ray
*		First, generate the set of at() options that define
*		the ray or expansion path

*	Stata permits a wide range of possible statistics for continuous variables
*	See help margins##atspec for details. Stata will catch illegitimate user input
if "`contstat'"=="" {
        local contstat mean  // default
}
*	Stata permits the class of variables _factor to be set to one of
*	mean, asobserved, base, asbalanced. Stata will catch illegitimate user input
/*
if "`fvstat'"~="" {
	local fvspec (`fvstat') _factor
}
*/

local atspecs
local n_atspecs = 0   // Count the number of observations to be created
foreach x of numlist `xvalues' {
	// Transformed values of the x-variable
	if "`transx'"=="none" {
		local transxfunc `x'
	}
	else {
		if "`transx'"=="bc" {
			local transxfunc ((`boxcox'*(`x') - 1)/`boxcox')
			local bc_option boxcox(`boxcox')
		}
		else {
			local transxfunc `transx'(`x')
		}
	}
	local tx = round(`transxfunc', 0.001)
	// Transformed values of the y-variable
	if "`transy'"=="none" {
		local transyfunc `intercept' + `slope'*`x'
	}
	else {  //  Assumes the same transformation is used on both the x- and y-variables
		if "`transx'"=="bc" {
			local transyfunc ((`boxcox'*(`intercept' + `slope'*`x') - 1)/`boxcox')
			local bc_option boxcox(`boxcox')
		}
		else {
			local transyfunc `transx'(`intercept' + `slope'*`x')
		}
	}
	local ty = round(`transyfunc',0.001)
	// Assemble the atspec for this pair of x- and y-values
	local atspecs `atspecs' at((`contstat') _continuous `fvspec' `margoptions' `txvar'=`tx' `tyvar'=`ty' )
*	local atspecs `atspecs'  at(`margoptions' `txvar'=`tx' `tyvar'=`ty' )
	local n_atspecs = `n_atspecs' + 1
}

*	Compose the title on the x-axis to describe the ray
if `"`xtitle1'"'~="" {
	local xtitle1 `"`xtitle1'"'
	local xtitle2 `"`xtitle2'"'
}
else if `slope'>0 {
	local xlbl : var label `xvar'
	local xtitle1 "`xlbl'"
	if `intercept'==0 {
		local y2xslope = round(`slope',.1)
		local xtitle2 "Assuming a ratio of `y2xslope' units of `yvar' per `xvar'"
	}
	else {
		local intcpt = round(`intercept',.01)
		local slp = round(`slope',.01)
		local xtitle2 "along a ray with `xvar' = `intcpt' + `slp' * `yvar'"
	}
}
else {
	local xtitle1 "Quantity of `xvar'"
	if `intercept'==0 {
		local xtitle2 "assuming no `yvar' services are delivered"
	}
	else {
		local intcpt = round(`intercept',.01)
		local slp = round(`slope',.01)
		local xtitle2 "assuming all facilities deliver `intercept' units of `yvar'"
	}
}

*	Select the transformation to unwind the dependent variable
unwindz , transz(`transz') `eqname'
	local predictz `r(predictz)'

*	Form an index of the x and y variables to serve as the denominator
*	of average cost.  Assume they are transformed the same way
unwindxy , transxy(`transx') xvar(`xvar') yvar(`yvar') yweight(`yweight') ///
	intercept(`intercept') slope(`slope') `bc_option'
	local outputindex `r(oi_x)'  //  Use the x-axis variable as the numeraire for AC
	local xunwound `r(xunwound)'
	local yunwound `r(yunwound)'

confirm numeric variable `xvar'
`shhh' di as txt _n "Description of the y-variable observations used in previous estimation"
`shhh' sum `yvar' if `touse', d
`shhh' di as txt _n "Description of the x-variable observations used in previous estimation"
`shhh' sum `xvar' if `touse', d
local median = r(p50)
local elastmax = `median'
local p5 = r(p5)
local p10 = r(p10)
local p25 = r(p25)
local p50 = r(p50)
local p75 = r(p75)
*local p95 = r(p95)  //  Jan., 2019: Suppress the 95th percentile xline
local xlines xlines(`p25' `median' `p75' `p95')  ///
	text(\`y4text' `p25' "25th Percentile", orient(vertical) size(medium) ///
		placement(w)) ///
	text(\`y4text' `median' "Median", orient(vertical) size(medium)  ///
		placement(w)) ///
	text(\`y4text' `p75' "75th Percentile", orient(vertical) size(medium)	///
		placement(w))

/*  //  Jan., 2019: Suppress the 95th percentile xline (Option?)
///
	text(\`y4text' `p95' "95th Percentile", orient(vertical) size(medium)	///
		placement(w))
*/

if "`replace'"=="replace" {
	dropany at1 at2  ///
		`xvar'_at `yvar'_at  ///
		`txvar'_at `tyvar'_at /* `tzvar'_hat_at  */ ///
		`zvar'_hat_at se_`zvar'_hat ll_`zvar'_hat ul_`zvar'_hat ///
		Avg_`zvar'_hat se_Avg_`zvar'_hat ll_Avg_`zvar'_hat ul_Avg_`zvar'_hat ///
		Mrg_`zvar'_hat se_Mrg_`zvar'_hat ll_Mrg_`zvar'_hat ul_Mrg_`zvar'_hat
}

**********************************************************
*	Generate new variables `xvar'_at, `yvar'_at and `zvar'_hat_at
*	and estimate the value of the dependent variable (e.g. Total Cost)
*	over a grid of x y values
`shhh' di as txt _n "Predict the total cost at selected output values" _n

`shhh' margins `marglist' `if' `in', `atspecs' expression(`predictz') level(`level') `nose'

*	Keep only two columns of the -at- matrix to generate the _at variables
	tempname atx aty at table
	mat define `at' = r(at)
	mat define `atx' = `at'[1...,"`txvar'"]
	mat define `aty' = `at'[1...,"`tyvar'"]
	mat define `at' = `atx',`aty'
svmat `at'
	rename `at'1 `txvar'_at // e.g. sqrt or ln of the number of ART patient-years
	rename `at'2 `tyvar'_at // e.g. sqrt or ln of the number of HIV patient-years
	genatvar `xvar', trans(`transx')
	genatvar `yvar', trans(`transy')
		format `txvar'_at `tyvar'_at `xvar'_at `yvar'_at %12.2fc
	matchvarlbl `xvar' `yvar' `txvar' `tyvar',  ///
		newvar(`xvar'_at `yvar'_at `txvar'_at `tyvar'_at)
	local newvrbls `xvar'_at `yvar'_at `txvar'_at `tyvar'_at
	local atvrbls `xvar'_at `yvar'_at `txvar'_at `tyvar'_at

	mat define `table' = r(table)'
	dropany b se z pvalue ll ul df crit eform

*	Store the predicted values from the margins command as "_hat" or "_hat_at" variables
*		The predicted value of the z-variable is named "`zvar'_hat_at" so that
*		the user can use "`zvar'_hat" for predictions over the entire sample.
*		For simplicity of notation, the "_at" suffix is omitted from the other
*		variables created by -margins-
svmat `table', names(col)
	rename b `zvar'_hat_at
		lab var `zvar'_hat_at "Predicted value of `zvar' at (`xvar'_at, `yvar'_at)"
	local depvrbls `zvar'_hat_at
	local newvrbls `newvrbls' `zvar'_hat_at
	if "`se'"=="se" {
		rename se se_`zvar'_hat
			lab var se_`zvar'_hat "Standard error of `zvar'_hat_at at (`xvar'_at, `yvar'_at)"
		rename ll ll_`zvar'_hat
			lab var ll_`zvar'_hat "Lower limit of `zvar'_hat_at at (`xvar'_at, `yvar'_at)"
		rename ul ul_`zvar'_hat
			lab var ul_`zvar'_hat "Upper limit of `zvar'_hat_at at (`xvar'_at, `yvar'_at)"
		local depvrbls `depvrbls' se_`zvar'_hat ll_`zvar'_hat ul_`zvar'_hat
		local newvrbls `newvrbls' se_`zvar'_hat ll_`zvar'_hat ul_`zvar'_hat
	}
	dropany z pvalue df crit eform
	format `newvrbls'  %12.2fc

*	Compute the formal definition of the economies of scale for a cost function
*	as total cost divided by the sum across outputs i of lambda_i*Yi*MCi
*	where lambda_x = 1/(slope + 1) and lambda_y = slope/(slope + 1)
*	and keep the partial derivative expressions for use below.
if "`transx'"=="none" {
	local chainwrtarg
}
else {
	if "`transx'"=="bc" {
		local chainwrtarg chainwrt(`transx') boxcox(`boxcox')
	}
	else {
		local chainwrtarg chainwrt(`transx')
	}
}

*		Here's the partial along the y dimension
`shhh' firstpartial `tyvar', cmdline("`origcmdline'") chainlhs(`transz') `chainwrtarg'
	local marg_d_TC_d_`yvar' `r(margexpression)'

*		Here's the partial along the x dimension
`shhh' firstpartial `txvar', cmdline("`origcmdline'") chainlhs(`transz') `chainwrtarg'
	local marg_d_TC_d_`xvar' `r(margexpression)'

`shhh' di as txt _n "Estimating the Ray Economies of Scale"
 *  Critical to use -xunwound- and -yunwound- here...
 `shhh' margins `marglist' `if' `in' /*if `xvar'>0 & `yvar'>0*/,  ///
	`atspecs' level(`level') `nose' ///
	expression(  ///
		`predictz'  ///
		    /       ///
		( `xunwound'*`marg_d_TC_d_`xvar'' + `yunwound'*`marg_d_TC_d_`yvar'' )  ///
	)

	mat define `table' = r(table)'
	dropany b se z pvalue ll ul df crit eform
	if "`replace'" ~= "" dropany RES_`zvar'_hat se_RES_`zvar'_hat ll_RES_`zvar'_hat ul_RES_`zvar'_hat
svmat `table', names(col)
	rename b RES_`zvar'_hat
		lab var RES_`zvar'_hat "Ray economies of scale along an expansion path"
	local newvrbls `newvrbls' RES_`zvar'_hat
	local resvrbls RES_`zvar'_hat
	if "`se'"=="se" {
		rename se se_RES_`zvar'_hat
			lab var se_RES_`zvar'_hat "Standard error of RES_`zvar'_hat"
		rename ll ll_RES_`zvar'_hat
			lab var ll_RES_`zvar'_hat "Lower limit of RES_`zvar'_hat"
		rename ul ul_RES_`zvar'_hat
			lab var ul_RES_`zvar'_hat "Upper limit of RES_`zvar'_hat"
		local newvrbls `newvrbls' RES_`zvar'_hat se_RES_`zvar'_hat ll_RES_`zvar'_hat ul_RES_`zvar'_hat
		local resvrbls `resvrbls' se_RES_`zvar'_hat ll_RES_`zvar'_hat ul_RES_`zvar'_hat
	}
	dropany z pvalue df crit eform

	format `newvrbls' %12.2fc

di _n as txt  "The Ray Economies of Scale at each projected point is:"
list `xvar'_at `yvar'_at `zvar'_hat_at `costperxvar' `resvrbls' in 1/`n_atspecs'

*	Index of the ray economies of scale along the ray from the 5th to the 50th percentile of x
numlist "`xvalues'", sort
local nbtwn5and50 = 0
foreach xv of numlist `r(numlist)' {
	if `xv' >= `p5' & `xv' <= `p50' {
		local nbtwn5and50 = `nbtwn5and50' + 1
	}
}

if `nbtwn5and50' < 2 {
	di as err "In order to estimate the Ray Economies of Scale,"  ///
		_n "-rayplot- requires the user to specify the option -xvalues(numlist)-"  ///
		_n "such that at least two of its values are between the 5th and the 50th"  ///
		_n "percentile of the -x- variable.  Here the 5th and 50th percentile are"  ///
		_n "5th percentile = `p5' and 50th percentile = `p50'"  ///
		_n "The option -xvalues(numlist)- is: `xvalues'"

	exit 198
}

`shhh' di as txt "sum RES_`zvar'_hat if `xvar'_at >= `p5' & `xvar'_at <= `p50'"
`shhh' sum RES_`zvar'_hat if `xvar'_at >= `p5' & `xvar'_at <= `p50'
	local RES_elast : di %7.1fc 100*`r(mean)'
	local tc_elast : di %7.1fc 100*(1/`r(mean)')
	local elastlbl "between percentiles 5 (`p5') and 50 (`p50')"

if "`resnote'"~= "" {
	local resnote note("Between the 5th & 50th percentile of output, doubling output is associated with a`tc_elast'% cost increase", span)
}

if `"`tcytitle'"'=="" {
	local ytitle1 "Total cost"
}
else {
	local ytitle1 `"`tcytitle'"'
}

******************PRODUCE GRAPHS OF TC & RES ****************************
*	Display graph of total cost or output function
if "`yrescale'" ~="" {
	rescale `zvar'_hat_at, thousands
		lab var `zvar'_hat_at "Predicted TC"
}

if "`se'"=="se" {
	if "`yrescale'" ~="" {
		rescale ul_`zvar'_hat ll_`zvar'_hat , thousands
	}
		sum ul_`zvar'_hat, meanonly
			local ul_max `r(max)'
		sum ll_`zvar'_hat, meanonly
			local ll_max `r(max)'
		local y4text = `ll_max' + 2*(`ul_max' - `ll_max')/3

	twoway (rarea ul_`zvar'_hat ll_`zvar'_hat `xvar'_at, pstyle(ci)) ///
		(line `zvar'_hat_at `xvar'_at, lwidth(thick)),  ///
		xlabels(#7, format(%9.0fc)) `xlines' ///
		ylabels(#7, angle(horizontal) `ylblfrmt') ///
		xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
		ytitle(`"`ytitle1'"' "Thousands of dollars") ///
		title("Total cost per accounting period" "along an expansion path")  ///
		subtitle(`"`subtitle'"', size(small) span) ///
		`resnote' `note'  ///
		legend(order( 2 1 "`level'% confidence interval") col(1) ring(0) pos(4))  ///
		name(`gphprefix'TC`_along_ray_'ci, replace) saving(`"`destination'/`gphprefix'TC`_along_ray_'ci.gph"', `replace')
	if "`yrescale'" ~="" {
		rescale ul_`zvar'_hat ll_`zvar'_hat `zvar'_hat_at, undo
	}
}
else {
	sum `zvar'_hat, meanonly
		local y4text = `r(min)' + 2*(`r(max)' - `r(min)')/3
	twoway ///
		(line `zvar'_hat_at `xvar'_at, lwidth(thick)),  ///
		xlabels(#7, format(%9.0fc)) `xlines' ///
		ylabels(#7, angle(horizontal) `ylblfrmt') ///
		xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
		ytitle(`"`ytitle1'"' "Thousands of dollars") ///
		title("Total cost per accounting period" "along an expansion path")  ///
		subtitle(`"`subtitle'"', size(small) span) ///
		`resnote' `note'  ///
		legend(off)  ///
		name(`gphprefix'TC`_along_ray', replace) saving(`"`destination'/`gphprefix'TC`_along_ray'.gph"', `replace')
		if "`yrescale'" ~="" {
			rescale `zvar'_hat_at, undo
		}
}

*	Display graph of Ray Economies of Scale and its confidence interval
qui summ `xvar', det
local p5 `r(p5)'
sum RES_`zvar'_hat if `xvar'_at > `p5', meanonly
local y4text = `r(min)' + 3*(`r(max)'-`r(min)')/4

if `"`resytitle'"'=="" {
	local ytitle1 "Returns to scale along an expansion path"
}
else {
	local ytitle1 `"`tcytitle'"'
}

maximorum RES_`zvar'_hat if `xvar'_at > `p5' & `xvar'_at < `p75'
	local toplabel = `r(toplabel)'
	local interval =  `r(interval)'
	local y4text = `r(y4text)'
	local ylblfrmt `r(ylblfrmt)'
if `r(maximorum)'< -1 | `r(maximorum)'>. {
	if `r(maximorum)'< -10^9 | `r(maximorum)'>. {
		local maxofres "undefined at values of `xvar' greater than `p5'"
	}
	else {
		local maxofres `r(maximorum)'
	}
	di as txt "Rayplot only constructs and displays a graph of the RES"  ///
		_n as txt "when its maximum estimated value exceeds -1 and is non-missing." ///
		_n as txt "Here the maximum estimated value of the RES is " as res "`maxofres'".
}
else {

	twoway (line RES_`zvar'_hat `xvar'_at if `xvar'_at > `p5' & `xvar'_at < `p75' & RES_`zvar'_hat > -1, lwidth(thick)) , ///
		xlabels(#7, format(%9.0fc)) `xlines' ///
		ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
		ylines(1,lwidth(medthick)) ///
		xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
		ytitle(`"`ytitle1'"' ) ///
		title("Ray economies of scale along an expansion path")  ///
		subtitle("Ratio of the average to the marginal cost" ///
			"Values above 1.00 denote economies of scale")  ///
		`resnote' `note'  ///
		legend(order( 2 "Ray economies of scale" 1 "`level'% confidence interval") col(1) ring(0) pos(2)) `note'  ///
		name(`gphprefix'RES`_along_ray', replace) saving(`"`destination'/`gphprefix'RES`_along_ray'.gph"', `replace')

	if "`se'"=="se" {
		tempname ll_RES ul_RES
		qui clonevar `ll_RES' = ll_RES_`zvar'_hat
			qui replace `ll_RES' = -1 if `ll_RES' < -1
		qui clonevar `ul_RES' = ul_RES_`zvar'_hat
			qui replace `ul_RES' = 3*RES_`zvar'_hat if `ul_RES' > 3*RES_`zvar'_hat
		maximorum ul_RES_`zvar'_hat ll_RES_`zvar'_hat if `xvar'_at > `p5' & `ll_RES' > -1 & `ul_RES' < 10 & RES_`zvar'_hat > -1
			local toplabel = `r(toplabel)'
			local interval =  `r(interval)'
			local y4text = `r(y4text)'
			local ylblfrmt `r(ylblfrmt)'

		twoway (rarea ul_RES_`zvar'_hat `ll_RES' `xvar'_at if `xvar'_at > `p5' & RES_`zvar'_hat > -1, pstyle(ci)) ///
			(line RES_`zvar'_hat `xvar'_at if `xvar'_at > `p5' & RES_`zvar'_hat > -1, lwidth(thick)) , ///
			xlabels(#7, format(%9.0fc)) `xlines' ///
			ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
			ylines(1,lwidth(medthick)) ///
			xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
			ytitle(`"`ytitle1'"' ) ///
			title("Ray economies of scale along an expansion path")  ///
			subtitle("Ratio of the average to the marginal cost" ///
				"Values above 1.00 denote economies of scale")  ///
			`resnote' `note'  ///
			legend(order( 2 "Ray economies of scale" 1 "`level'% confidence interval") col(1) ring(0) pos(2)) `note'  ///
			name(`gphprefix'RES`_along_ray_'ci, replace) saving(`"`destination'/`gphprefix'RES`_along_ray_'ci.gph"', `replace')
	}
}

**********************************************************
*	Estimate average cost along the ray
`shhh' di as txt _n "Predict the average cost at selected output values" _n

*	Any definition of average cost must be evaluated at a point in the xy plane.
*	The following alternative definitions of average cost all evaluate
*	average cost along the line in the xy plane defined by the parameters
*	-intercept- and -slope-.

/*
*	Weighted average of the average costs, where at any point (x,y)
*	the angle of a ray from the origin to that point, as a proportion
*	of 90 degrees, is the weight on the average cost of x, with
*	one minus this proportion being the weight on the average cost of y
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') expression(   ///
		(1-180*atan(`yunwound'/`xunwound')/_pi/90)*(`predictz'/`yunwound')  ///
				+		///
		(  180*atan(`yunwound'/`xunwound')/_pi/90)*(`predictz'/`xunwound')  ///
	)
*/


/*
*	Average with respect to an average output, where average output
*	is defined as above.
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') expression(   ///
		(`predictz') 		///
			/				///
		(					///
			(1-180*atan(`yunwound'/`xunwound')/_pi/90)*(`xunwound')  ///
					+		///
			(  180*atan(`yunwound'/`xunwound')/_pi/90)*(`yunwound')  ///
		)  ///
	)
*/

/*
*	Average cost of a one unit movement along the ray from the origin
*	where the unit is defined by the Pythagorean theorem
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') ///
	expression( (`predictz') / sqrt(`xunwound'^2+`yunwound'^2) )
*/

/*
*	Average cost of a one unit movement along the ray from the origin
*	(Same as above expression with the `outputindex' in the denominator.)
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') ///
	expression( (`predictz') / (`outputindex') )
*/


*	Average cost per person-year of services
*	(Recommended by Baumol, Panzar & Willig, 1982, p. 49, Def. 3.4.2)
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') ///
	expression( (`predictz') / (`xunwound'+`yunwound') )


/*
*	Average cost of another person-year of service y assuming 1/slope years of x per y
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') ///
	expression( (`predictz') / (`yunwound') )
*/

*	Average cost of another person-year of service x assuming -slope- years of y per x
*		This definition amounts to assuming that the y output is essential to
*		the production of x in the proportion specified along the ray.
*		It attributes all cost to the x output and zero cost to the y output
/*
`shhh' margins `marglist' `if' `in' /*if `xvar'>0 & `yvar'>0*/, `nose'  ///
	`atspecs' level(`level') ///
	expression( (`predictz') / (`xunwound') )
*/

mat define `table' = r(table)'

	dropany b se z pvalue ll ul df crit eform

svmat `table', names(col)
	rename b Avg_`zvar'_hat
		lab var Avg_`zvar'_hat "Average cost of `xvar': TC/`xvar'"
	local newvrbls `newvrbls' Avg_`zvar'_hat
	local avgvrbls Avg_`zvar'_hat
	if "`se'"=="se" {
		rename se se_Avg_`zvar'_hat
			lab var se_Avg_`zvar'_hat "SE of Average cost of `xvar': TC/`xvar'"
		rename ll ll_Avg_`zvar'_hat
			lab var ll_Avg_`zvar'_hat "Lower limit of `level' confidence interval of Avg_`zvar'_hat"
		rename ul ul_Avg_`zvar'_hat
			lab var ul_Avg_`zvar'_hat "Upper limit of `level' confidence interval of Avg_`zvar'_hat"
		dropany z pvalue df crit eform
		local newvrbls `newvrbls' se_Avg_`zvar'_hat ll_Avg_`zvar'_hat ul_Avg_`zvar'_hat
		local avgvrbls `avgvrbls' se_Avg_`zvar'_hat ll_Avg_`zvar'_hat ul_Avg_`zvar'_hat
	}
	if "`yrescale'" ~="" {
		rescale `avgvrbls', thousands
	}
	format `avgvrbls' %12.2fc

	if `"`acytitle'"'=="" {
		local ytitle1 "Average cost of composite output"
	}
	else {
		local ytitle1 `"`acytitle'"'
	}

	if "`se'"=="se" {
		maximorum ul_Avg_`zvar'_hat Avg_`zvar'_hat ll_Avg_`zvar'_hat
			local toplabel = `r(toplabel)'
			local interval =  `r(interval)'
			local y4text = `r(y4text)'
			local ylblfrmt `r(ylblfrmt)'

		twoway (rarea ul_Avg_`zvar'_hat ll_Avg_`zvar'_hat `xvar'_at, pstyle(ci)) ///
			(line Avg_`zvar'_hat `xvar'_at, lwidth(thick)) , ///
			xlabels(#7, format(%9.0fc))  `xlines' ///
			ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
			xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
			ytitle(`"`ytitle1'"' `"`ytitle2'"') ///
			title("Average cost along an expansion path", span)  ///
			subtitle(`"`subtitle'"', size(small) span) ///
			`resnote' `note'  ///
			legend(order( 2 "Average" 1 "`level'% confidence interval") col(1) ring(0) pos(2)) ///
			name(`gphprefix'AC`_along_ray_'ci, replace)  ///
			saving(`"`destination'/`gphprefix'AC`_along_ray_'ci.gph"', `replace')
	}
	else {
		maximorum Avg_`zvar'_hat
			local toplabel = `r(toplabel)'
			local interval =  `r(interval)'
			local y4text = `r(y4text)'
			local ylblfrmt `r(ylblfrmt)'

		twoway  ///
			(line Avg_`zvar'_hat `xvar'_at, lwidth(thick)) , ///
			xlabels(#7, format(%9.0fc))  `xlines' ///
			ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
			xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
			ytitle(`"`ytitle1'"' `"`ytitle2'"') ///
			title("Average cost along an expansion path", span)  ///
			subtitle(`"`subtitle'"', size(small) span) ///
			`resnote' `note'  ///
			legend(off)  ///
			name(`gphprefix'AC`_along_ray', replace)  ///
			saving(`"`destination'/`gphprefix'AC`_along_ray_'ci.gph"', `replace')
	}

**********************************************************
*	The marginal cost of ART at these same points is:

`shhh' di as txt _n "Predict the marginal cost at selected output values" _n

/*
*	Along a ray following the equation y = intercept + slope*x
*	This approach blows up if expansion path is close to either axis
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') expression(  ///
		sqrt( ///
			`xunwound'^2 + (`intercept' + `slope'*`xunwound')^2 ///
		)  ///
		*  ///
		(  ///
			(`marg_d_TC_d_`xvar''/`xunwound')+(`marg_d_TC_d_`yvar''/`yunwound')  ///
		)	///
	)
*/
/*
*	Average of the marginal costs according to the angle of the expansion path
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') expression(  ///
		(1-180*atan(`yunwound'/`xunwound')/_pi/90)*(`marg_d_TC_d_`xvar'')  ///
					+				///
		(180*atan(`yunwound'/`xunwound')/_pi/90)*(`marg_d_TC_d_`yvar'')  ///
	)
*/
/*
*	Marginal cost of another person-year of service y assuming 1/slope years of x per y
`shhh' margins `marglist' `if' `in', `atspecs' level(`level') expression(  ///
		`marg_d_TC_d_`yvar''  ///
					+				///
		(1/`slope')*(`marg_d_TC_d_`xvar'')  ///
	)
*/

*	Marginal cost of another person-year of service x assuming slope years of y per x
`shhh' di _n as txt "    Marginal cost of another person-year of service x assuming slope years of y per x"
`shhh' margins `marglist' `if' `in' , `nose'  ///
	`atspecs' level(`level') expression(  ///
		`marg_d_TC_d_`xvar''  ///
					+				///
		(`slope')*(`marg_d_TC_d_`yvar'')  ///
	)


mat define table = r(table)'

dropany b se z pvalue ll ul df crit eform
svmat table, names(col)
	rename b Mrg_`zvar'_hat
		lab var Mrg_`zvar'_hat "Marginal cost of `xvar': d TC/d `xvar'"
local newvrbls `newvrbls' Mrg_`zvar'_hat
local mrgvrbls Mrg_`zvar'_hat
if "`se'"=="se" {
	rename se se_Mrg_`zvar'_hat
		lab var se_Mrg_`zvar'_hat "SE of Marginal cost of `xvar': d TC/d `xvar'"
	rename ll ll_Mrg_`zvar'_hat
		lab var ll_Mrg_`zvar'_hat "Lower limit of Marginal cost of `xvar': d TC/d `xvar'"
	rename ul ul_Mrg_`zvar'_hat
		lab var ul_Mrg_`zvar'_hat "Upper limit of Marginal cost of `xvar': d TC/d `xvar'"
	drop z pvalue df crit eform
	local newvrbls `newvrbls' se_Mrg_`zvar'_hat ll_Mrg_`zvar'_hat ul_Mrg_`zvar'_hat
	local mrgvrbls `mrgvrbls' se_Mrg_`zvar'_hat ll_Mrg_`zvar'_hat ul_Mrg_`zvar'_hat
}
if "`yrescale'" ~="" {
	rescale `mrgvrbls', thousands
}
format `mrgvrbls'  %12.2fc

*	For each point in output space, list the total, marginal and average cost
*	First the average cost
di _n as txt  "The Average Cost at each projected point is:"
list `xvar'_at `yvar'_at `zvar'_hat_at `avgvrbls'  in 1/`n_atspecs'

*	Then the marginal cost
di _n as txt  "The Marginal Cost at each projected point is:"
list `xvar'_at `yvar'_at `zvar'_hat_at `mrgvrbls' in 1/`n_atspecs'

*	List average and marginal together with their ratio
tempvar ratioAvg2Mrg
qui gen `ratioAvg2Mrg' = Avg_`zvar'_hat/Mrg_`zvar'_hat
	lab var `ratioAvg2Mrg' "Quotient of Average divided by Marginal"
	char define `ratioAvg2Mrg'[varname] ratioAvg2Mrg
di _n as txt  "The ratio of Average to Marginal Cost at each projected point is:"
list `xvar'_at `zvar'_hat_at Avg_`zvar'_hat Mrg_`zvar'_hat `ratioAvg2Mrg' in 1/`n_atspecs', subvarname

if `"`mcytitle'"'=="" {
	local ytitle1 "Marginal cost of composite output"
}
else {
	local ytitle1 `"`mcytitle'"'
}

******************PRODUCE GRAPHS OF MC & AC ****************************
*	Graph the marginal cost function with its confidence interval
if "`se'"=="se" {
	maximorum ul_Mrg_`zvar'_hat Mrg_`zvar'_hat ll_Mrg_`zvar'_hat
		local toplabel = `r(toplabel)'
		local interval =  `r(interval)'
		local y4text = `toplabel'/4
		local ylblfrmt `r(ylblfrmt)'

	twoway (rarea ul_Mrg_`zvar'_hat ll_Mrg_`zvar'_hat `xvar'_at, pstyle(ci)) ///
		(line Mrg_`zvar'_hat `xvar'_at, lwidth(thick)) , ///
		xlabels(#7, format(%9.0fc))  `xlines' ///
		ylines(0,lwidth(medthick)) ///
		ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
		xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
		ytitle(`"`ytitle1'"' `"`ytitle2'"') ///
		title("Marginal cost along an expansion path", span)  ///
		subtitle(`"`subtitle'"', size(small) span) ///
		legend(order( 2 "Marginal" 1 "`level'% confidence interval") col(1) ring(0) pos(2)) `note'  ///
		name(`gphprefix'MC`_along_ray_'ci, replace) saving(`"`destination'/`gphprefix'MC`_along_ray_'ci.gph"', `replace')
}

*	Overlay the MC on the AC without the CI's
if `"`acmcytitle'"'=="" {
	local ytitle1 "Average and marginal cost of composite output"
}
else {
	local ytitle1 `"`acmcytitle'"'
}

maximorum Avg_`zvar'_hat Mrg_`zvar'_hat
	local toplabel = `r(toplabel)'
	local interval =  `r(interval)'
	local y4text = `r(y4text)'
	local ylblfrmt `r(ylblfrmt)'

twoway (line Avg_`zvar'_hat `xvar'_at, lwidth(thick))  ///
	(line Mrg_`zvar'_hat `xvar'_at, lwidth(thick)) , ///
	xlabels(#7, format(%9.0fc))  `xlines' ///
	ylines(0,lwidth(medthick)) ///
	ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
	xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
	ytitle(`"`ytitle1'"' `"`ytitle2'"') ///
	title("Average and marginal cost along an expansion path", span)  ///
	subtitle(`"`subtitle'"', size(small) span) ///
	`resnote' `note'  ///
	legend(order( 1 "Average" 2 "Marginal") col(1) ring(0) pos(2)) ///
	name(`gphprefix'AC_MC`_along_ray', replace) saving(`"`destination'/`gphprefix'AC_MC`_along_ray'.gph"', `replace')

*	Overlay the MC on the AC with the CI's
if "`se'"=="se" {
	maximorum ul_Avg_`zvar'_hat ul_Mrg_`zvar'_hat
		local toplabel = `r(toplabel)'
		local interval =  `r(interval)'
		local y4text = `r(y4text)'
		local ylblfrmt `r(ylblfrmt)'

	twoway (rarea ul_Avg_`zvar'_hat ll_Avg_`zvar'_hat `xvar'_at, pstyle(ci)) ///
		(rarea ul_Mrg_`zvar'_hat ll_Mrg_`zvar'_hat `xvar'_at, pstyle(ci)) ///
		(line Avg_`zvar'_hat `xvar'_at, lwidth(thick))  ///
		(line Mrg_`zvar'_hat `xvar'_at,lwidth(thick)),  ///
		xlabels(#7, format(%9.0fc)) `xlines'  ///
		ylines(0,lwidth(medthick)) ///
		ylabels(0(`interval')`toplabel', angle(horizontal) `ylblfrmt') ///
		xtitle(`"`xtitle1'"' `"`xtitle2'"')  ///
		ytitle(`"`ytitle1'"' `"`ytitle2'"') ///
		yscale(titlegap(1)) ///
		title("Average and Marginal cost along an expansion path", span)  ///
		subtitle(`"`subtitle'"', size(small) span) ///
		`resnote' `note'  ///
		legend(order(2 "Average" 3 "Marginal"  ///
			1 "`level'% confidence intervals")  ///
			col(1) ring(0) pos(2)) ///
		name(`gphprefix'AC_MC`_along_ray_'ci, replace) saving(`"`destination'/`gphprefix'AC_MC`_along_ray_'ci.gph"', `replace')
}
if "`yrescale'" ~="" {
	rescale `avgvrbls' `mrgvrbls', undo
}

**********************WRAP UP, SAVE AND EXIT **************************
di _n as txt "By executing Stata's -margins- command, this program"  ///
	_n as txt "has generated the following new variables."
des `newvrbls'
if "`verbose'"=="verbose" {
	sum `newvrbls'
}

*	Since the number of xvalues at which these are evaluated will typically
*	be small, return these results as matrices.
tempname atmat depmat resmat avgmat mrgmat
mkmat `atvrbls', matrix(`atmat') nomissing
foreach vlist in depvrbls resvrbls avgvrbls mrgvrbls {
	local prfx = substr("`vlist'",1,strpos("`vlist'","vrbls")-1)
	mkmat `xvar'_at ``vlist'', matrix(``prfx'mat') nomissing
}

if "`saving'"~="" {
	tokenize `"`saving'"', parse(" ,")
	local filename `1'
	local comma `2'
	local savereplace `3'
	preserve
	quietly {
		keep `newvrbls'
		gen long n = _n
			lab var n "Order of simulated observations"
		sort n
		keep if `xvar'_at~=.
		order n `newvrbls'
		compress
		notes drop _all
		notes _dta: The command line of the regression was:
		notes _dta: `cmdline'
		label data "Simulated observations generated by -rayplot-"
	}
	save "`destination'/`filename'" `comma' `savereplace'
	di _n as txt "New variables have been saved to the file " as res "`destination'/`filename'"
	restore
}
if "`keep'"=="" {
	di _n as txt "By default the new variables are dropped at program completion." _n ///
		"Specify -keep- to retain them in memory or -saving(filename)- to save them to disk."
	dropany `newvrbls'
}
return local gphprefix `gphprefix'
return local newvrbls `newvrbls'
foreach matname in atmat depmat resmat avgmat mrgmat  {
	return matrix `matname' = ``matname''
}

end	//	End of program rayplot.ado

program chktrans	// This version differs from that in -firstpartial.ado-
*	Check to see if a variable has the expected prefix
	syntax , varname(string) [trans(string) force ]
	local var_prefix = substr("`varname'",1,strpos("`varname'","_")-1)
	if "`var_prefix'"~="`trans'" {
		di as err "A -chain- option has been selected to unwind a " as res "`trans'" as err " function"  ///
			_n as err "on the variable " as res "`varname'" as err ".  But the prefix on the variable" ///
			_n as err "does not match the trans function designated in the -chain- option." _n
		if "`force'"=="force" {
			di as err "Because the " as res "`force'" as err " option has been specified, execution proceeds."
		}
		else {
			di as err "Execution terminates.  To force execution despite this incongruence," ///
				_n	"specify the " as txt "force" as err " option on the call to -rayplot-"
			exit 198
		}
	}
end	// End of program chcktrans

program define getelast, rclass
	version 11.1
	syntax varlist (max=2), [min(real 1) max(real 100000000)]
	tempname currentestimates
	est store `currentestimates'
	tokenize `varlist'
	tempvar ln_depvar ln_rayvar
	gen `ln_depvar' = ln(`1')
	gen	`ln_rayvar' = ln(`2')
	di _n as txt " Over the range of " as res "`rayvar'"  ///
		as txt " from " as res `min' as txt " to " as res `max' as txt "," ///
		as txt _n "the log-log regression yields the following elasticity estimate:"
	reg `ln_depvar' `ln_rayvar' if `2'>`min' & `2'<`max'
		local tc_elast: di  %8.1fc 100* (1 - _b[`ln_rayvar'])

	return local elasticity `tc_elast'
	if `min'>1 | `max'<100000000 {
			sum `2', meanonly
			local min = round(`r(min)',1)
			return local elastlbl " in range `min' to `max'"
	}
	est restore `currentestimates'
end // End of program getelast

program define maximorum, rclass
	version 11.1
	syntax varlist (numeric) [if]

	marksample touse

	local maximorum = -10^10
	foreach var in `varlist' {
		sum `var' if `touse', meanonly
		if "`r(max)'"~="" {
			if `r(max)' > `maximorum' & `r(max)'<=. {
				local varwithmax `var'
				local maximorum = `r(max)'
			}
		}
	}

	local toplabel =  ///
		cond(`maximorum'>50,100*ceil(`maximorum'/100),  ///
		cond(`maximorum'>20,50*ceil(`maximorum'/50),  ///
		cond(`maximorum'>10,20*ceil(`maximorum'/20),  ///
		cond(`maximorum'>5,10*ceil(`maximorum'/10),  ///
		cond(`maximorum'>1,5*ceil(`maximorum'/5),1)))))

	local interval =  ///
		cond(`toplabel'>500000,100000,  ///
		cond(`toplabel'>100000,25000,  ///
		cond(`toplabel'>20000,5000,  ///
		cond(`toplabel'>10000,2500,  ///
		cond(`toplabel'>5000,2000,  ///
		cond(`toplabel'>3000,1000,  ///
		cond(`toplabel'>1000,250,  ///
		cond(`toplabel'>500,200,  ///
		cond(`toplabel'>100,50,  ///
		cond(`toplabel'>50,20,  ///
		cond(`toplabel'>20,10,  ///
		cond(`toplabel'>10,5,  ///
		cond(`toplabel'>1,2,0.2)))))))))))))

	local fmt =  ///
		cond(`toplabel'>500000,"%12.0fc",  ///
		cond(`toplabel'>100000,"%10.0fc",  ///
		cond(`toplabel'>20000,"%9.0fc",  ///
		cond(`toplabel'>5000,"%8.0fc",  ///
		cond(`toplabel'>2000,"%7.0fc",  ///
		cond(`toplabel'>400,"%6.0fc",  ///
		cond(`toplabel'>100,"%6.0fc",  ///
		cond(`toplabel'>50,"%6.0fc","%6.0fc"))))))))
	local y4text = 3*`toplabel'/4

	return local ylblfrmt format(`fmt')
	return local y4text `y4text'
	return local interval `interval'
	return local toplabel `toplabel'
	return local varwithmax `varwithmax'
	return local maximorum = `maximorum'
end	// End of program maximorum

program define matchvarlbl
	version 11.1
	syntax varlist, newvar(varlist)
	tokenize `varlist'
	foreach var in `newvar' {
		local existinglbl : variable label `1'
		label variable `var' `"`existinglbl'"'
		macro shift
	}
end  // End of program matchvarlbl

program define dropany
	version 11.1
	syntax namelist
	novarabbrev {
		foreach var in namelist `namelist' {
			capture drop `var'
		}
	}
end  // End of program dropany

program define transid, rclass
*	Identify the transformation for the named variable
*	and name its transformed version
	syntax name, [chain(string) chainlhs(string) chainwrt(string) ]

	local nwords: word count `chain' `chainwrt' `chainlhs'
	if `nwords'>1 {
		di as err "Only one of options -chain-, -chainwrt- and -chainlhs-" //
			_n "is permitted on a call to this sub-routine."
		exit 198
	}

	if "`chain'`chainlhs'`chainwrt'" == "" {
		local trans none		//  The "none" transformation type added in version 1.8
		local transvar `namelist'
	}
	else {
			local trans `chain'`chainlhs'`chainwrt'
			local transvar `trans'_`namelist'
	}

	//	This is the tranformattion that must be unwound
	return local trans `trans'
	//	The following variable name was included in the estimation
	//	and thus will be recognized by the margins computations.
	return local transvar `transvar'

end	//  End of the program transid

program define unwindz, rclass
	syntax , [transz(string) eqname(string)]
	if "`eqname'"=="" {
		local pred_opt xb
	}
	else {
		local pred_opt equation(`eqname')
	}
	if "`transz'"== "" | "`transz'"== "none" {
		local predictz predict(`pred_opt')
	}
	else {
		if "`transz'"~="sqrt" & "`transz'"~="ln" & "`transz'"~="log10" {
			di as err "Option -transz- must equal -ln-, -sqrt- or -log10-"
			exit 198
		}
		if "`transz'"== "sqrt" {
			local predictz predict(`pred_opt')^2
		}
		if "`transz'"== "ln" {
			local predictz exp(predict(`pred_opt'))
		}
		if "`transz'"== "log10" {
			local predictz 10^predict(`pred_opt')
			di as err "Transformation -log10- no yet implemented."
			exit 198
		}
	}
	return local predictz `predictz'
end	 //	End of program unwindz

program define unwindxy, rclass
*	This subroutine specifies the denominator of the average cost calculation.
*	Possibilities include:
*		1.	The Pythagorean distance along the ray in output space
*		2.	The sum of the outputs on horizontal and vertical axes
*		3.	Disregard one of the outputs, attributing all costs to the other
*		4.	A weighted sum of the outputs, with weights representing the
*			relative "utility" of the two outputs: (e.g. the marginal
*			rate of substitution in a welfare function with these two outputs)
	version 11.1
	syntax, transxy(string) intercept(real) slope(real) ///
		xvar(varlist) yvar(varlist) yweight(real) [ boxcox(real 0.5) ]

	if "`transxy'"~="ln" & "`transxy'"~="log10" & "`transxy'"~="sqrt" & "`transxy'"~="bc" & "`transxy'"~="none" {
		di as err "Transformation on x and y variables must be -ln- or -log10- or -sqrt- or -bc- or -none-"
		exit 198
	}
	if `slope' < 0 {
		di as err "Slope must be non-negative"
		exit 198
	}

	if "`transxy'"=="bc" {
		local xunwound (`boxcox'*c.`transxy'_`xvar' + 1)^(1/`boxcox')
		local yunwound (`boxcox'*c.`transxy'_`yvar' + 1)^(1/`boxcox')
	}
	if "`transxy'"=="sqrt" {
		local xunwound c.`transxy'_`xvar'*c.`transxy'_`xvar'
		local yunwound c.`transxy'_`yvar'*c.`transxy'_`yvar'
	}
	if "`transxy'"=="ln" {
		local xunwound exp(c.`transxy'_`xvar')
		local yunwound exp(c.`transxy'_`yvar')
	}
	if "`transxy'"=="log10" {
		local xunwound 10^(c.`transxy'_`xvar')
		local yunwound 10^(c.`transxy'_`yvar')
	}
	if "`transxy'"=="none" {
		local xunwound `xvar'
		local yunwound `yvar'
	}

	//	If the ray is virtually parallel to the x-axis, ignore the y dimension
	if `slope'< .05 {
		local oi_x `xunwound'
	}
	//	Otherwise, measure output
	else {
*	1.	The Pythagorean distance along the ray in output space
		local oi_pythag  sqrt((`xunwound')^2 + (`yunwound')^2)
*	2.	Use x alone as the output index
		local oi_x `xunwound'
*	2.	Use the sum of x and y along the expansion ray as the output index
		local oi_sum  `xunwound' + `yunwound'
	}
	return local oi_pythag
	return local oi_x
	return local oi_sum
	return local xunwound `xunwound'
	return local yunwound `yunwound'
	return local outputindex `outputindex'

end	//	End of program unwindxy

program define genatvar
	syntax name, trans(string) [ boxcox(real 0.5) ]
	if "`trans'"=="bc" {
		genl `namelist'_at = (`boxcox'*`trans'_`namelist'_at - 1)^(1/`boxcox')
	}
	if "`trans'"=="sqrt" {
		genl `namelist'_at = `trans'_`namelist'_at^2
	}
	if "`trans'"=="ln" {
		genl `namelist'_at = exp(`trans'_`namelist'_at)
	}
	if "`trans'"=="log10" {
		genl `namelist'_at = 10^(`trans'_`namelist'_at)
	}
end	//	End of program genatvar
/*
program cleancmdline, rclass
*	Version 1.2 by Mead Over 19Aug2016
*	This programs accepts the command line from either a single or a
*	multiequation estimator and returns a simplified command line which
*	begins with the dependent variable.  For a multiequation command like
*	-sureg- or -reg3-, this program strips e(cmdline) of the estimation command,
*	the equation name (identified because it is followed by a colon) and all
*	parentheses.

	syntax , cmdline(string)

	//	Break the command line up into chunks
	//	Assumes the command line can be parsed on spaces and a colon
	//	and that the first two tokens are the name of the equation
	//	and the dependent variable.

	local origcmdline `cmdline'
	tokenize `"`cmdline'"', parse(" (")

	// If the first or second token in the command line is a left parenthesis,
	// strip everything up to the -(- character from the beginning of the -cmdline-
	if "`1'"=="(" | "`2'"=="(" {
		local cmdline = substr(`"`cmdline'"',strpos(`"`cmdline'"',"(")+1,.)

		//	If the second token was a -(-, then the first is an estimation command
		if "`2'"=="(" {
			local estcmd `1'
		}
	}
	else {  //  Is the first token an estimation command?
		capture confirm numeric variable `1'
		if _rc==0 {
			di as err "First word in the command line should be one of Stata's"  ///
			_n "estimation commands that is linear in the parameters."  ///
			_n "Examples include {help reg}, {help sureg}, {help frontier}, etc.."  ///
			_n "Here the first word in the command line is the variable: " as res "`1'"
			exit 198
		}
		else {
			local estcmd `1'
			macro shift
			local cmdline `*'
		}
	}

	//	Now consider the command line without the estimation command or the
	//	equation binding prefixes.
	tokenize `"`cmdline'"', parse(" :")

	if "`2'"==":" {
		local eqname `1'
		local depvar `3'
		local cmdline = substr(`"`cmdline'"',strpos(`"`cmdline'"',":")+1,.)
	}
	else {
*		local eqname = `1'   //  replacing this line by the following one in 1.21
		local eqname = word("`e(eqnames)'",1)
		local depvar `1'
	}

	capture confirm numeric variable `depvar'
		if _rc~=0 {
			di as err _n "Unable to find the dependent variable in the command line" ///
				_n "Here it is: -" as res "`depvar'" as err "-"
			di as err _n "which has been extracted from the command line:" ///
				_n as res `"`origcmdline'"'
			exit 198
		}

	//	truncate the command line before the first parenthesis
	local firstparen = min(strpos(`"`cmdline'"',"("),strpos(`"`cmdline'"',")"))
	if `firstparen' >0 & `firstparen' < . {
		local cmdline = substr(`"`cmdline'"',1,`firstparen'-1)
	}

	return local cmdline `cmdline'
	return local origcmdline `origcmdline'
	return local estcmd `estcmd'
	return local eqname `eqname'
	return local depvar `depvar'

*	Version information to program -cleancmdline-
*	Version 1.0 13July2016
*	Version 1.1 13July2016		Only change is to edit the comment in line 20
*	Version 1.2 19Aug2016		Clean simple command lines like: "reg y x1 x2 x1*x2"
end /* End of program cleancmdline */
*/
*	Wrapper program executes -margins- four times to compute
*		total cost along a ray: line 238
*		ray economies of scale: line 306
*		average cost along a ray: line 553
*		marginal cost along a ray: line 662
*
*	Version history:
*	1.1 Generalizes to accept a transformation of the lhs different
*		from that applied to the rhs variables.
*	1.2 Generalizes the divisor in the average cost calculation
*	1.3 Changes name of predicted dependent variable from *_hat to _hat_at
*		and adds variable labels
*	1.4 Add the -cmdline()- option to allow user to supply a command line,
*		which then overrides the command line from the most recent estimation
*		and is an alternative to the restore option.
*	1.5 Cosmetic changes to the graph titles
*	1.6 Replace the hard-coded "artpy- with `xvar'
*	1.7 10/2012 Allow user to specify the -note- option for bottom of graph
*	1.9 11Jul2016 Add the -none- transformation type to transid and unwindxy
*			Generalize the mat define at commands
*	1.10 23Feb2017 Call firstpartial with the original commandline
*					include code for cleancmdline as part of this program
*					Enable xtitle1, 2 options. Replace %12.0fc with %12.2fc
*	1.11 2Mar2017 Allow abreviations of options
*	1.12 6Mar2017 Fix bug in -maximoroum-
*	1.13 3May2017 Add subtitle to the Total Cost curve
*					Suppress printing of ll_RES_tc_hat and ul_RES_tc_hat
*					Fixed `tc_elast' to refer to p50 instead of p75
*	1.14 6Jun2017 Generalize -xtitle1- default, reposition xline titles on TC graph
*					Suppress display of RES graphs when maximum estimated RES < -1.
*					Correct computation of RES average to between 5th & 50th %ile
*	2.00 22Oct2017 Change the SLope(option) to SLopes(numlist)
*					Enable -chainwrt- option to allow -bc- for boxcox
*					When bc is specified, the boxcox parameter is by default 0.5
*					Specifying the option -boxcox(#)- allows user to choose another
*	2.01 22Oct2017  Fix bug in line 59
*	2.02 5Feb2018   When listing RES_tc_hat in line 311, also list ll and ul
*					 change label to read "Ray Economies" not "Ray Elasticity"
*	2.03 31Mar2018  Add trap for the number of xvalues between p5 and p50 < 2
*	2.10 30Apr2018  Add -fvset- options from -margcost- to set factor variables
*					If -replace-, drop the RES_hat variables
*	2.11 11May2018  Add option -gphprefix- to the returned macros
*	2.12 11May2018  Add saving function and -destination- option
*	2.20 13May2018  Make computation of standard errors, confidence intervals, optional
*	2.21 18May2018  Add error trap for fvset`i' containing an empty macro.
*	2.22 21May2018  Enable -destination- option to specify location for -keep- data.
*	2.23 26Dec2018  When there is no transformation define transxfunc, transyfunc
*		in lines 114, 127
*	2.24 8Jan2019	Correct line 127 to properly read: local transyfunc `intercept' + `slope'*`x'
*	2.30 10Jan2019	Specify values of factor variables directly, by allowing a -varlist-
*		and conveying it to the -margins- commands. Remove the -fvset- options.
*		Specify -marksample- and change former occurrences of marksample to markout.
*		Add -contstat-, -fvstat- and -margoptions- options to specify value(s)
*		of rhs variables at which the projections and graphs are constructed.
*  		 Add the "yrescale" option.  Shorten the default filenames unless option -old-
*	2.31 31Jan2019	Comment out internal veresion of -cleancmdline- in favor of ver. 1.5
*		Replace multi-level factor variables in the commmand line (beginning with -i.- )
*		with with the levels specified for those variables in -marglist-
*	Future development:  Add an option e.g. RESXvalues(numlist) to restrict EOS measures
*		by a default to a smaller range, such as less than the median of the output data.
*		Also align code on factor-variable replacemment between this program annd -contourplot-
*	2.32  12Mar2019  Add the -resnote- option to suppress the note on RES between 5th & 50th %ile
*	2.40  7May2019  Change \*`gphprefix' to /`gphprefix to avoid \ escape interpretation
*		Change legacy local macro ART2preART to y2xslope
*		Fix bug in option -saving- which saves the generated "_at" variables