prog define isoplot
*!	version 1.6		by Mead Over	6Oct2016
	version 11.1
	syntax varlist (max=2 min=2) [if/][in],  ISOVar(varname) CCUTs(string) ///
		[CENTiles values PLOTJoint(string) PLOTSpec(string)   ///
		EPSilon(real 0) TOLerance(real 1) function(string) function2(string) ///
		noSCAtter CONTour LOWess  BWidth(real 0.1) LWidth(passthru)   ///
		LGNDTitle(string) LGNDSubtitle(string) LEGend(passthru) *]

*	varlist				Two existing variables, typically two outputs which
*						appear on the right-hand-side of a cost function.
*						The first variable will be plotted on the y-axis
*						and the second on the x-axis.
*
*  REQUIRED OPTIONS
*	isovar(varname)		Name of the existing variable for which level-set scatter
*						plots are constructed.  Typically this variable is
*						the total cost (or estimated total cost) of producing
*						the multiple outputs that appear on the right-hand-side
*						of a multiple output cost function.
*
*	centiles|values		One of these two options must be specified.
*
*	ccuts(numlist)		List of numbers between 0 and 100, which determines
*						the values at which level sets (or contours)
*						will be plotted.  With the -centiles- option,
*						these are percentiles of the -isovar- istribution.
*						With the option values, these are the values of
*						-isovar-.
*
*  OPTIONAL OPTIONS
*	noScatter			Suppress the scatterplots
*
*	lowess				Superimpose nonparametricly fitted lines onto the
*						contours.  The lines are fitted with Stata's
*						twoway lowess graph command.
*
*	bwidth(0.1)			Bandwidth for the lowess line graphs. Default is 0.1.
*
*	lwidth(linestyle)	Width of the line used in lowess plots.
*
*	contour				Superimpose a contourline constructed by Stata's
*						-twoway contourline- command. If this option is
*						selected, the user is confined to specifying only
*						a single value of ccut().  To draw lines at multiple
*						values, choose the -lowess- option.
*
*	plotjoint(mean|median|p25|p75)
*						Plot the output combination representing joint
*						production of the mean, median 75th or 25th
*						percentile of the two outputs.
*
*	plotspec(stat)	Plot the two output combinations representing the
*	epsilon(e)			(quasi-) specialized production of the two outputs.
*						Combinations plotted will be [e*Y1_m, (1-e)*Y2_m]
*						and at [(1-e)*Y1_m, e*Y2_m], where -m- denotes
*						the -stat- selected by the plot option
*						and e = epsilon in the P&B definition of -qscope-.
*						The default is e = 0.
*
*	function(string)	Overlay a graph of the output appearing on the vertical
*						axis as a function of the output on the horizontal axis.
*						Specify the function by writing y = f(x), where
*						the syntax is the same as in graph twoway function.
*
*	function2(string)	Overlay a second function graph.
*
*	tolerance(tol)
*		If -centiles-
*						Bands of plotted points lie within a range of each
*						of the designated percentiles of -isovar-, where the
*						range is a multiple -tol- of the difference
*						between the upper and lower bounds of the estimated
*						centiles. Default is 1.0
*		If -values-
*						Bands of plotted points lie within a range of each
*						of the designated values of -isovar-, where the
*						range is a percentage -tol- plus or minus that level.
*						Default is 1.0%.
*
*	lgndtitle(string)	Specify the legend's title
*	lgndsubtitle(string)Specify the legend's subtitle.
*
*	*					Any other -graph twoway- options.

	marksample touse
	markout `touse' `isovar'

	if "`if'"~="" {
		local ifopt & `touse'
	}

	if ("`centiles'"~="" & "`values'" ~= "") | ("`centiles'"=="" & "`values'" == "") {
		di as err "User must specify either option -centiles-"  ///
			_n "or option -values-, but not both."
		exit 198
	}

	if strpos("`ccuts'",",")>0  {
		if strpos(substr("`ccuts'",strpos("`ccuts'",","),.),"asis")>0 {
			local ccuts = substr("`ccuts'",1,strpos("`ccuts'",",")-1)
			local asis asis
		}
		else {
			di as err "The only legitimate suboption for the -ccuts()- option is -asis-."  ///
					_n "Here the option -ccuts()- is: " as res "`ccuts'"
			exit 198
		}
	}
	//	Unless the option asis is selected, sort the ccuts in ascending order
	if "`asis'"=="" {
		numlist "`ccuts'", sort
		local ccuts `r(numlist)'
	}
	local nccuts : word count `ccuts'
	tokenize `ccuts'
	if "`centiles'"~="" {
		centile `isovar', centile( `ccuts' )
	}
	else {
		di _n as txt "Isocost bands defined by the following ccuts"  ///
			_n "and their lower and upper tolerance bounds,"  ///
			_n "defined as " as res `tolerance'   ///
			as txt "% above and below the contour value"  ///

		di _n as txt "             {c |}                                    Lower           Upper"  ///
				  _n "    Variable {c |}     Order          Contour         Bound           Bound"  ///
				  _n "{hline 13}{c +}{hline 61}"  ///
				  _n " `isovar'" _col(14) "{c |}" _c
	}
	local nextline
	foreach c of numlist 1/`nccuts' {
		tempvar `isovar'_cut`c'

		if "`centiles'"~="" {
			local lb = `r(c_`c')' - `tolerance'*0.5*(`r(c_`c')' - `r(lb_`c')')
			local ub = `r(c_`c')' + `tolerance'*0.5*(`r(ub_`c')' - `r(c_`c')')
			local disptile`c' : di %-9.1fc `r(c_`c')'
		}
		else {
			local lb = (1-`tolerance'/100)*``c''
			local ub = (1+`tolerance'/100)*``c''
			local disptile`c' : di %-9.1fc ``c''
		}

		//	Generate and plot a separate variable for each value of ccuts.
		//	Each variable is defined as the observations of -isovar- that
		//	lie within the specified band around a value in the list of ccuts.
		qui gen ``isovar'_cut`c'' = `isovar' if `isovar' < `ub' & `isovar' >= `lb' `ifopt'

		if "`scatter'"~="noscatter" {
			local scattercmds `scattercmds' ///
				(scatter `varlist' if ``isovar'_cut`c'' < .  `ifopt')
		}

		if "`contour'"~="" {
			local contourcmds `contourcmds' ///
				(contourline `isovar' `varlist' if ``isovar'_cut`c'' < .  `ifopt', levels(2))
		}
		if "`lowess'"~="" {
			local lowesscmds `lowesscmds' ///
				(lowess `varlist' if ``isovar'_cut`c'' < .  `ifopt', bwidth(`bwidth') `lwidth' )
		}

		if "`centiles'"~="" {
			local lgndcmds `lgndcmds'	`c' `"``c''th %ile: $`disptile`c''"'
		}
		else {
			local lgndcmds `lgndcmds'	`c' `"Within `tolerance'% of $`disptile`c''"'
		}
		local lgndordr = `c'
		if `"`lgndtitle'"' == "" {
			local lgndtitle "Isocost frontiers"
		}
		if `"`lgndsubtitle'"' == "" {
			local lgndsubtitle "of `isovar'"
		}
		if "`centiles'"=="" {
			di as txt `nextline' _col(23) as res "`c'" _col(30) %12.3fc ``c'' _col(44) %12.3fc `lb'  _col(60) %12.3fc `ub'
			local nextline _col(14) "{c |}"
		}
	}

	if "`plotjoint'"~=""  {
		local plotjnt (scatter
		foreach var of varlist `varlist' {
			tempvar `var'jnt
			qui summ `var', detail
			if substr("`plotjoint'",1,4)=="mean" {
				local stat = `r(mean)'
			}
			else if substr("`plotjoint'",1,6)=="median" | "`plotjoint'"=="p50" {
				local stat = `r(p50)'
			}
			else if "`plotjoint'"=="p25" {
				local stat = `r(p25)'
			}
			else if "`plotjoint'"=="p75" {
				local stat = `r(p75)'
			}
			else {
				di as err "The only legitimate values for -stat- in the option -plotjoint(stat)-"  ///
						_n "are -mean-, -median-, -p25- and -p75-. Here the option has the value: "  ///
						as res "`plotjoint'"
				exit 198
			}
			qui gen ``var'jnt' = `stat' if `var' < . & `var'[_n-1]==.
			local plotjnt `plotjnt' ``var'jnt'
			local plotjntvars `plotjntvars' ``var'jnt'
		}
		local lgndcmds `lgndcmds' - "-------------------"

		local plotjnt `plotjnt'  , msize(huge))
		local lgndordr = `lgndordr' + 1
		local lgndcmds `lgndcmds' `lgndordr' "Joint production" "of `plotjoint' outputs" "in one firm"
		local otherplots `plotjnt'
	}

	if "`plotspec'"~=""  {
		local plotspcl (scatter
		local eps = `epsilon'
		local di_eps : di %-9.7g `epsilon'
		foreach var of varlist `varlist' {
			tempvar `var'spcl
			qui summ `var', detail
			if substr("`plotspec'",1,4)=="mean" {
				local stat = `r(mean)'
			}
			else if substr("`plotspec'",1,6)=="median" | "`plotspec'"=="p50" {
				local stat = `r(p50)'
			}
			else if "`plotspec'"=="p25" {
				local stat = `r(p25)'
			}
			else if "`plotspec'"=="p75" {
				local stat = `r(p75)'
			}
			else {
				di as err "The only legitimate values for -stat- in the option -plotspec(stat)-"  ///
						_n "are -mean-, -median-, -p25- and -p75-. Here the option has the value: "  ///
						as res "`plotspec'"
				exit 198
			}
			qui gen ``var'spcl' = `eps'*`stat' if `var' < . & `var'[_n-1]==.
				qui replace ``var'spcl' = (1-`eps')*`stat' if ``var'spcl' >= . & ``var'jnt'[_n-1]<.
				local eps = 1 - `eps'
			local plotspcl `plotspcl' ``var'spcl'
			local plotspclvars `plotspclvars' ``var'spcl'
		}
		if "`otherplots'"==""{
			local lgndcmds `lgndcmds' - "-------------------"
		}
		local plotspcl `plotspcl'  , msize(huge))
		local lgndordr = `lgndordr' + 1
		local lgndcmds `lgndcmds' `lgndordr' "Specialized production" "of `plotjoint' outputs" "in two firms" "with e = `di_eps'"
		local otherplots `otherplots' `plotspcl'
	}

	if "`otherplots'`contourcmds' `lowesscmds'"~="" | {
		local axis y
		foreach var of varlist `varlist' {
			local `axis'vlbl : variable label `var'
			if "``axis'vlbl'"=="" {
				local `axis'vlbl `var'
			}
			local axistitles `axistitles' `axis'title("``axis'vlbl'")
			local axis x
		}
	}

	if "`function'"~=""  {
		local functioncmds (function `function')
	}

	if "`function2'"~=""  {
		local functioncmds `functioncmds' (function `function2')
	}

	if `"`legend'"'=="" {
		local legend  ///
			legend( size(vsmall)  ///
				order(  `lgndcmds' )  ///
				cols(1) pos(3)  ///
				title(`"`lgndtitle'"', size(medium))  ///
				subtitle(`"`lgndsubtitle'"')  ///
			)
	}

	twoway `scattercmds' `otherplots' `contourcmds' `lowesscmds' `functioncmds',  ///
		`axistitles' `legend' `options'

end  /* End of program -isoplot-  */
*	Version 1.2	1Oct2016 Based on version 1.1 of isodraw.ado-.  Add the -values- option
*				and the lgndtitle and lgndsubtitle options.  Add the p75 statistic.
*	Version 1.3	2Oct2016 Separate the plotjoint() option from the specialized option.
*				Order the scatter plot isocost lines in ascending order
*	Version 1.4 3Oct1916 Add the -lowess- and the -contour- options.
*				Allow option -noscatter- to suppress the points.
*				Give user optional control over the legend.
*	Version 1.5 3Oct1916 Add option to plot one or two functions.
*				Add the [in] syntax and add the marksample command
*	Version 1.6 6Oct1916 Change option -specialized- to option -plotspec-