program define contourplot, rclass
*!	program contourplot		version 3.2	by Mead Over	10Jan2019

*	This program assumes that a preceding estimation command
*	has estimated a non-linear function of at least two variables,
*	say -y- and -x-, where the non-linearity is captured by interactions of
*	Stata's continuous factor variables, of the form c.x#c.y.
*	The objective of this program is to construct a graph
*	to trace out level sets of the z-variable at the values
*	defined by the option -ccuts- in order to reveal
*	the shape of a three dimensional surface with the dependent variable
*	on the z-axis (coming out of the page) and the designated y- and
*	x-variables on the vertical (or y) axis and the horizonatal
*	(or x) axis respectivley.

	version 11.1
	syntax [varlist (fv default=none)] [if/] , zvar(string) yvar(string) xvar(string)  ///
			[xvalues(numlist) yvalues(numlist) slope(real -9) ///  Optional from ver 2.2
			chain(passthru) chainlhs(passthru) chainwrt(passthru) ///
			restore(string) cmdline(string)  ///
			intercept(real 0) level(real 95) verbose replace force ///
			yweight(real 1) name(passthru) saving(passthru) keep ///
			xtitle(string) ytitle(string)  note(passthru) ///
			xmin(string) xmax(string) ymin(string) ymax(string)  ///
			ccuts(passthru) xlabel(passthru) text(passthru)   ///
			TItle(passthru)  xlog ylog ///  Options added July 2016
			LTItle(string) LSUBtitle(string) LSIze(string) LSSize(string)  ///
			ylabel(passthru) hilite(string) INTerp(passthru)  ///  Options added June 2017
			timer(real 0) CONTstat(string) FVstat(string) SE   /// Options added Aug & Sep 2017
			MARGOPTions(string) *]  // Added 5May2018

*-----------------------------------------------------------------------------
*  Optional factor variable expression
*	a factor variable expression which reslves to a single factor variable
*	specification can optionally be used to set the values of those
*	factor variables for the execution of contourplot
*-----------------------------------------------------------------------------
*  Required parameters:
*	yvar(string)		Variable on vertical axis (typically untransformed) to be assigned name `yvar'_at
*	xvar(string)		Variable on horizontal axis (typically untransformed) to be assigned name `xvar'_at
*	zvar(string)		Predicted variable on z axis (typically untransformed) to be assigned name `zvar'_hat_at
*-----------------------------------------------------------------------------
*  Optionally unwind variable transformations:
*	chain(string)			o	Unwind all three variables, x, y and z
*		Not yet implemented
*	chainlhs(string)		o	Unwind the lhs or z-variable with this trans
*	chainwrt(string)		o	Unwind the x- and y- variables with this trans
*-----------------------------------------------------------------------------
*  Other optional parameters:
*	xvalues(numlist)	list of values of -xvar- at which y will be estimated
*	yvalues(numlist)	list of values of -yvar- at which y will be estimated
*	yweight(real)		The weight to place on the y-variable in computing
*							total output for the average cost function
*	slope(real)			Slope of an expansion ray to be superimposed on the plot
*							(default is to omit plotting a ray fromthe origin)
*	intercept(real)		Intercept on y axis of the ray (default = 0)
*	replace				Overwrite any existing variables with _at
*							suffixes
*	verbose				Requests output of margins commands
*	force				Forces execution despite the prefix of the dependent
*							variable not being the expected one.
*	restore(string)		Name of previously stored estimates to be recalled.
*						Only necessary if they are not currently in memory.
*	xtitle(string)		Specify the title on the X-axis.  Default is X's label
*	ytitle(string)		Specify the title on the Y-axis.  Default is Y's label
*	note(string)		Note to appear at the bottom of the plot. (Default is no note.)
*	text(y x "string")	A text annotation to the graph positioned at (x, y)
*	xmin(string)		Default would be the minimum or maximum of the
*	xmax(string)		corresponding variable/
*	ymin(string)
*	ymax(string)
*	title(passthru)		Title at top of the contourplot.
*	ltitle(string)		Legend title: (Default: Isocost fromtiers)
*	lsubtitle(string)	Legend subtitle (Default: $1,000's)
*	lsize(string)		Legend size (Default: medium)
//  Options added July 2016
*	xlog				Scale the x-axis with option xscale(log)
*	ylog				Scale the y-axis with option yscale(log)
*	LTItle(string)		Legend title
*	LSUBtitle(string) 	Legend subtitle
*	LSIze(string)		Legend title size
*	LSSize(string)		Legend subtitle size
*	xlabel(passthru)	labels for the x-axis
*	ylabel(passthru)	labels for the y-axis
//  Options added June, Aug 2017
*	hilite(string)
*	INTerp(passthru)	Passed to twoway contourline.  See help twoway_contour##interpolation
*	timer(real 0)
*	CONTstat(string)
*	FVstat(string)
*	SE				Request -margins- to compute standard errors and upper & lower confidence bounds
*	*				Additional options are passed to the countourline command as `options'
//  Options added May, 2018
*	margoptions(string)	Additional -at- options to specify value(s) of rhs variables
*					at which the contour's are constructed.

if `timer'>0 {
	timer clear `timer'
	timer on `timer'
}

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

*	Optionally restore previous estimation results
if "`restore'"~="" {
	est restore `restore'  // name of previous estimation result
}
*	Retrieve information from preceding estimation command
local depvar: word 1 of `e(depvar)'
capture confirm numeric variable `depvar'
if _rc~= 0 {
	di as err "Estimation results must be restored"
	exit 198
}
local cmdline `e(cmdline)'  //  Only used to label saved data

*	Define the sample based on the previous estimation with
*	the possibility of further exclusions based on the [if/]
*	in this program's command.
marksample touse
	tempvar esmpl
	gen byte `esmpl' = e(sample)
	markout `touse' `esmpl'

local marglist `varlist'
/*
if "`marglist'"=="" {
	qui replace `touse' = 1
}
else {
*/
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 specifications"  ///
			_n as err "provided as arguments may not contain the general term -i.-."  ///
			_n as err "Here the arguments provided are: `marglist'"
		exit 198
	}
	drop `r(varlist)'
}

*	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'

*	User specified legend title, subtitle and size
	if `"`ltitle'"'=="" {
		local ltitle Isocost frontiers
	}
	if `"`lsubtitle'"'=="" {
		local lsubtitle ($1,000's)
	}
	if `"`lsize'"'~="" {
		local lsize , size(`lsize')
	}
	if "`lssize'"~=""{
		local lssize , size(`lssize')
	}

*	User specified  X- and Y-titles:
if "`xtitle'"=="" {
	local xtitle : variable label `xvar'
}
if "`ytitle'"=="" {
	local ytitle : variable label `yvar'
}

*	Suppress -margins- computation of the standard errors and upper and lower bounds
if "`se'" == "se" {
	local nose
}
else {
	local nose nose
}

*	User specified minimima and maxima for the plot axes.
if "`xmin'"=="" {
	sum `xvar', meanonly
		local xmin = `r(min)'
}
else {
	local xmin = real("`xmin'")
}
local andxmin  & `xvar' > `xmin'

if "`xmax'"=="" {
	sum `xvar', meanonly
		local xmax = `r(max)'
}
else {
	local xmax = real("`xmax'")
}
local andxmax  & `xvar' <= `xmax'

if "`ymin'"=="" {
	sum `yvar', meanonly
		local ymin = `r(min)'
}
else {
	local ymin = real("`ymin'")
}
local andymin  & `yvar' > `ymin'
if "`ymax'"=="" {
	sum `yvar', meanonly
		local ymax = `r(max)'
}
else {
	local ymax = real("`ymax'")
}
local andymax  & `yvar' <= `ymax'

*	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 fvstat mean  // default
}

*	User specified log scale for x- and/or y-axis
if "`xlog'"~="" {
	local xlog xscale(log)
	local xmin = max(`xmin',.1)
}
if "`ylog'"~="" {
	local ylog yscale(log)
	local ymin = max(`ymin',.1)
}

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

capture confirm numeric variable `xvar'
if _rc==0 {
	qui sum `xvar' if `touse', d
	local median = r(p50)
	local elastmax = `median'
	local p5 = r(p5)
	local p25 = r(p25)
	local p75 = r(p75)
	local p95 = r(p95)
	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)) ///
		text(\`y4text' `p95' "95th Percentile", orient(vertical) size(medium)	///
			placement(w))
}
else {
	local elastmax = 1000
}

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 combinations of the x- and y-variables"

*		First, generate the set of at() options that define the grid
*			(For rayplot, predicted values must be generated along a ray from the origin.)
*			(For contourplot, in contrast, the predicted values should be generated
*			across a grid which is evenly spaced in the x-y plane.)

if "`xvalues'"=="" {
	tempvar f_of_x autogenx
	if "`transx'"=="none" {
		qui gen `f_of_x' =  `xvar'
		local andfofxmin = "`andxmin'"
		local andfofxmax = "`andxmax'"
	}
	else {
		qui gen `f_of_x' =  `transx'(`xvar')
		local andfofxmin = "& `f_of_x' > " + strofreal(`transx'(real(substr("`andxmin'",strpos("`andxmin'",">")+1,.))))
		local andfofxmax = "& `f_of_x' <= " + strofreal(`transx'(real(substr("`andxmax'",strpos("`andxmax'","<=")+1,.))))
	}

	sum `f_of_x' if `touse' `andfofxmin' `andfofxmax' `andif', meanonly
	qui gen `autogenx' if `touse' = autocode(`f_of_x',20,r(min),r(max))
	qui levelsof `autogenx', local(atspecsx)
}
else {
	foreach x of numlist `xvalues' {
		if "`transx'"=="none" {
			local transxfunc
		}
		else {
			local transxfunc `transx'
		}
		local tx = round(`transxfunc'(`x'), 0.001)
		local atspecsx `atspecsx'  `tx'
	}
}

if "`yvalues'"=="" {
	tempvar f_of_y autogeny

	if "`transy'"=="none" {
		qui gen `f_of_y' =  `yvar'
		local andfofymin = "`andymin'"
		local andfofymax = "`andymax'"
	}
	else {
		qui gen `f_of_y' =  `transy'(`yvar')
		local andfofymin = "& `f_of_y' > " + strofreal(`transy'(real(substr("`andymin'",strpos("`andymin'",">")+1,.))))
		local andfofymax = "& `f_of_y' <= " + strofreal(`transy'(real(substr("`andymax'",strpos("`andymax'","<=")+1,.))))
	}

	sum `f_of_y' if `touse' `andfofymin' `andfofymax' `andif', meanonly
	qui gen `autogeny' if `touse' = autocode(`f_of_y',20,r(min),r(max))
	qui levelsof `autogeny', local(atspecsy)
}
else {
	foreach y of numlist `yvalues' {
		if "`transy'"=="none" {
			local transyfunc
		}
		else {
			local transyfunc `transy'
		}
		local ty = round(`transyfunc'(`y'), 0.001)
		local atspecsy `atspecsy'  `ty'
	}
}

if `timer'>0 {
	timer off `timer'
	`shhh' timer list `timer'
	local timeb4marg `r(t`timer')'
	`shhh' di as txt "Seconds from start of contourplot until execution of margins:" `timeb4marg'
	timer on `timer'
}

di _n as txt "Use Stata's -margins- command to generate the predicted dependent variables" ///
	_n " over a grid of the x- and y-variables.  The command syntax is:"
di as txt "margins `marglist' if `yvar'>`ymin' & `xvar'>`xmin' `andtouse', `nose'"  ///
	_n "expression(`predictz')"  ///
	_n "at((`contstat') _continuous (`fvstat') _factor `txvar'=(`atspecsx') `tyvar'=(`atspecsy') `margoptions' )"

`shhh' margins `marglist' if `yvar'>`ymin' & `xvar'>`xmin' `andtouse', `nose' expression(`predictz')  ///
	at((`contstat') _continuous (`fvstat') _factor `txvar'=(`atspecsx') `tyvar'=(`atspecsy')  ///
	`margoptions' )

if `timer'>0 {
	timer off `timer'
	`shhh' timer list `timer'
	local timeofmarg = `r(t`timer')' - `timeb4marg'
	`shhh' di as txt "Seconds to execute margins: " `timeofmarg'
	timer on `timer'
}

*	Keep only two columns of the -at- matrix to generate the _at variables
	tempname at atx aty table
	mat define `at' = r(at)
	mat define `atx' = `at'[1...,"`txvar'"]
	mat define `aty' = `at'[1...,"`tyvar'"]

svmat `atx'
	rename `atx'1 `txvar'_at // e.g. sqrt or ln of the number of X variable
svmat `aty'
	rename `aty'1 `tyvar'_at // e.g. sqrt or ln of the number of Y variable
	genatvar `xvar', trans(`transx')
	genatvar `yvar', trans(`transy')
		format `txvar'_at `tyvar'_at `xvar'_at `yvar'_at %12.0fc
	matchvarlbl `xvar' `yvar' `txvar' `tyvar',  ///
		newvar(`xvar'_at `yvar'_at `txvar'_at `tyvar'_at)
	local newvrbls `xvar'_at `yvar'_at `txvar'_at `tyvar'_at

	mat define `table' = r(table)'
	if "`replace'"=="replace" {
		if `"`nose'"' == "" {
			dropany b se z pvalue ll ul df crit eform
		}
		else {
			dropany b
		}
	}

*	Store margins command's predicted values of the -zvar- as "_hat_at" variables
*		For simplicity of notation, the "_at" suffix is omitted from the other
*		variables created by -margins-, which are the se_, ll_, and ul_ variables.
svmat `table', names(col)
	rename b `zvar'_hat_at
		lab var `zvar'_hat_at "Predicted value of `zvar' at (`xvar'_at, `yvar'_at)"
	if `"`nose'"' == "" {
		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)"
		dropany z pvalue df crit eform
		format `zvar'_hat_at se_`zvar'_hat ll_`zvar'_hat ul_`zvar'_hat %12.0fc
		local newvrbls `newvrbls' `zvar'_hat_at se_`zvar'_hat ll_`zvar'_hat ul_`zvar'_hat
	}
	else {
		format `zvar'_hat_at %12.0fc
		local newvrbls `newvrbls' `zvar'_hat_at
	}

*	Default is no expansion path plotted at all
if `slope' > 0 {
	tempvar ray_ratio
	genl `ray_ratio' = `intercept' + `slope'*`xvar' if `yvar'>`ymin' & `yvar'<`xmax'
		/// Constrain the line's y dimension to be smaller than xmax in order
		/// to prevent the graph from extending the height of the scatter plot.
	local expansionpath (line `ray_ratio' `xvar' if `ray_ratio'>`ymin' & `xvar'>`xmin' `andtouse', sort)
}

	if `"`nose'"' == "" {
		rescale ul_`zvar'_hat ll_`zvar'_hat `zvar'_hat_at, thousands
	}
	else {
		rescale `zvar'_hat_at, thousands
	}

if "`hilite'" ~= "" {
	local hilite 	(scatter \`yvar' \`xvar' if `hilite' , msymbol(Oh) msize(huge))
}

if `timer'>0 {
	timer off `timer'
	`shhh' timer list `timer'
	local timeb4gph `r(t`timer')'
	`shhh' di as txt "Seconds until start of twoway contourline: " `timeb4gph' - `timeofmarg'
*	timer clear `timer'
	timer on `timer'
}

*	Construct the contour plot
sum `xvar'_at, meanonly
	local xmin_at `r(min)'
	local xmax_at `r(max)'
sum `yvar'_at, meanonly
	local ymin_at `r(min)'
	local ymax_at `r(max)'

twoway  ///
	(contourline `zvar'_hat_at `yvar'_at `xvar'_at  ///
		if `yvar'_at>`ymin_at' & `xvar'_at>`xmin_at'  ///
		 & `yvar'_at<`ymax_at' & `xvar'_at<`xmax_at',  ///
		`ccuts' colorlines `interp' `options' )  ///
	(scatter `yvar' `xvar' if `yvar'>`ymin' & `xvar'>`xmin' `andtouse')  ///
	`expansionpath' `hilite' ,  ///
	plegend(on   ///
		title(`"`ltitle'"', `lsize' )  ///
		subtitle(`"`lsubtitle'"' `lssize' )  ///
		pos(2))  ///
	`xlabel' `xlog' xtitle(`"`xtitle'"') ///
	`ylabel' `ylog' ytitle(`"`ytitle'"') ///
	`title'  ///
	legend(off)  `note'  `text' `name' `saving'

if `timer'>0 {
	timer off `timer'
	`shhh' timer list `timer'
	local timeofgph `r(t`timer')'
	`shhh' di as txt "Seconds to execute twoway (contourline) (scatter): " `timeofgph' - `timeb4gph'
*	timer clear `timer'
	timer on `timer'
}

if `"`nose'"' == "" {
	rescale ul_`zvar'_hat ll_`zvar'_hat `zvar'_hat_at, undo
}
else {
	rescale `zvar'_hat_at, undo
}

*******************************************************************
di _n as txt "By executing Stata's -margins- command, this program has generated the following new variables."
des `newvrbls'
if "`verbose'"=="verbose" {
	sum `newvrbls'
}
if "`saving'"~="" {
	tokenize `"`saving'"', parse("( ,)")
	if "`1'"=="saving" & "`2'"=="(" {
		macro shift
		macro shift
	}
	else {
		di as err _n "Option -saving- misspecified.  It reads:" as txt "`saving'"
		exit 198
	}
	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 -contourplot-"
	}
	save `filename' `comma' `savereplace'
	di _n as txt "Simulated data used to construct the contour plot has been saved to the file " as res "`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'
}

if `timer'>0 {
	timer off `timer'
	`shhh' timer list `timer'
	local timetotal `r(t`timer')'
	`shhh' di as txt "Time to execute contourplot: " `timetotal' " seconds equals " `timetotal'/60 " minutes."
	return local timetotal `timetotal'
	return local timeofgph `timeofgph'
	return local timeb4gph `timeb4gph'
	return local timeofmarg `timeofmarg'
	return local timeb4marg `timeb4marg'
}

return add
return local newvrbls `newvrbls'

end	//	End of program contourplot.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 -contourplot-"
			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)

	local maximorum = -10^10
	foreach var in `varlist' {
		sum `var', meanonly
		if `r(max)' > `maximorum' {
			local varwithmax `var'
			local maximorum = `r(max)'
		}
	}
	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
	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 2.1
		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)]
	if "`transz'"== "" | "`transz'"== "none" {
		local predictz predict(xb)
	}
	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(xb)^2
		}
		if "`transz'"== "ln" {
			local predictz exp(predict(xb))
		}
		if "`transz'"== "log10" {
			local predictz 10^predict(xb)
			di as err "Transformation -log10- no yet implemented."
			exit 198
		}
	}
	return local predictz `predictz'
end	 //	End of program unwindz

program define genatvar
	syntax name, trans(string)
	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

*	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 Require that estimation results either already be in memory
*		or be restored by the -restore- option.
*		(The -cmdline- option is irrelevant  to contourplot.ado)
*	1.5 Replace hard coded reference to -artpy-
*	1.6 Add -xtitle- and -ytitle- options
*	1.7 11/5/12 Allow user to specify a note at the bottom of the contourplot with a -note- option
*	2.0 5/1/16 Begin to generalize the program for use on data other than from the MATCH study
*			First step: eliminate the requirement for a variable called "nation"
*			Enable the -if- syntax.  Suppress the -in- syntax
*			Add a -hilite- option expressed as a logical condition
*	2.1 14July2016 Allow transid to be "none"
*	2.2 18July2016 Fix sample selection in the -countourline- command.
*			Make optional the -xvalues- and -yvalues- options.
*			Default is to use 20 equally spaced values within their ranges.
*			Delete -rayplot-'s -unwindxy- subroutine, which is not used here
*			Change the -gphprefix- option to a conventional -name()- option
*	2.3 10Aug2016 Repair bug in -saving- option, which now saves
*			both the gph file and the simulated data
*			Implement -lsize- and -lssize- options
*	2.4 6Jun2017 Replace the word -sample- with `touse'
*			Add the interp(passthru) option.  See help twoway_contour##interpolation
*			Add -hilite- to the syntax.  But it needs to be tested
*	2.5 25Aug2017 Add the timer(real 0) option to track time use by program segment
*			Attempt to enable contourplot to find its own xvalues and yvalues.
*	2.6 25Sep2017 Add options to fix continuous and factor values at values other than their means
*	2.7 3Oct2017 Jeff Pitblado's addition of the NOSE option passed to -margins-
*	2.8 20Oct2017 Reverse Jeff P.'s option so that NOSE is the default, now the
*				option -se- produces s.e.'s, ll's and ul's. Also enhance timer report
*	2.9 1Apr2018 Suppress artifactual horz on vert lines in contourplot
*				by constraining it to lie within the max and min of the
*				xvar_at and yvar_at variables created by -margins-
*	3.0 26Apr2018 Add the sort option to the expansion path in row 404
*				and change the vertical limit to `ray_ratio'>`ymin'
*	3.1 5May2018 Add the -margoptions- option to allow additional -at- specifications
*	3.2 10Jan2019 Add -marglist- to pass categorical variables to the margins commands
*		Apply dropany to the `tzvar'_hat_at variable