*!	Version 1.4 23Jan2017 by Mead Over to accompany version 1.5 of -isoplot-

program define isoplot_example
	version 11.1

	set more off
	`0'

	capture prog drop isoplot
end

program Msg
	di as txt
	di as txt "-> " as inp `"`0'"'
end

program Ms2
*	Use this version for messages containing macros
*	Use _char(36) for $ (the global macro delimiter)
*	Use _char(96) for ` (the backtick)
*	Use _char(39) for ' (the foretick or single quote)
	di as txt
	local intro `"as txt "-> ""'
	while `"`0'"'~="" {
		local nextchar = strpos(`"`0'"', "_char(")
		if `nextchar' >0 {
			local stuff = substr(`"`0'"',1,`nextchar'-1)
			local macdelim = substr(`"`0'"',`nextchar',9)
			di `intro' as inp `"`space'`stuff'"' `macdelim' _c
			local 0 = substr(`"`0'"',`nextchar'+9,.)
			if "`macdelim'"=="_char(39)" {
				local space " "
			}
			else {
				local space
			}
			local intro
		}
		else {
			di `intro' as inp `"`0'"'
			local 0
		}
	}
end

program Xeq
	di as txt
	di as txt `"-> "' as inp _asis `"`0'"'
	`0'
end

program define NewFile
	args fn
	capture confirm new file `fn'
	if _rc {
		di as txt "{p 0 2 2}"
		di "Example cannot be run because you already have a"
		di "file named {res:`fn'}"
		di "{p_end}"
		error 602
	}
end

program define NewName
	args fn
	capture graph rename `fn' `fn'
	if _rc==0 {
		di as txt "{p 0 2 2}"
		di "Example cannot be run because you already have a"
		di "graph in memory named {res:`fn'}."
		di "Either drop this graph or simply click {stata graph drop _all}."
		di "{p_end}"
		di as err _n "graph `fn' already exists"
		exit 110
	}
end

program define cost_dgp, rclass
*	Data generating process for simulating cost function data with Cost = f(y2,y1; w)
*	Generate the two output variables such that y2 has a uniform distribution
*	and y1 has a specified correlation with y2.

*	Procedure recommended by Bill Gould here:
*	http://www.stata.com/support/faqs/statistics/multivariate-normal-observations/

syntax [, seed(real 54321) n(real 10000) NORMal INVNormal rho(real 0)  ///
	alpha(real 100) beta(real 3) scale(real 1) gamma(real 1) sigma(real 0.2)]

*	seed  : Seed for random number generation
*	n     : Size of simulated data set
*	invnormal: Simulated data is joint normal distributed (default is uniform data)
*	rho   : Correlation coefficient between y1 and y2
*	alpha : Intercept parameter
*	beta  : Complementarity parameter
*	scale : Scale elasticity parameter
*	gamma : Elasticty of cost w.r.t. the factor price, w
*	sigma : Elasticty of cost w.r.t. the factor price, w

di as txt _n "Beginning of the data generating process."  ///
	"  See Gould (19??) in the {help isoplot##references:isoplot help file} for a description of the data generation methods."

	qui describe
	if `r(N)'~=0 {
		di as err "Running these examples requires that you first {stata clear} Stata's memory"
		exit 198
	}
	local rhotxt = subinstr("`rho'",".","",.)
	NewName corr_y2y1`rhotxt'

*	Specify the correlation coefficient between the two outputs
*	Xeq	clear
	Xeq	set seed `seed'
	Xeq	global rho = `rho'
	Msg *	Specify the correlation matrix between the two output variables, y1 and y2.
	Ms2 matrix P = (1,_char(36){rho} \ _char(36){rho}, 1)
		matrix P = (1,${rho} \ ${rho}, 1)
	Msg *	Compute the Cholesky factorization of the correlation matrix
	Xeq matrix A = cholesky(P)
	Msg *	Store the elements of A as local macros:
	Msg local coef_c1 = A[2,2]
		local coef_c1 = A[2,2]
		local di_coef_c1 : di %9.5f `coef_c1'
	Msg local coef_c2 = A[2,1]
		local coef_c2 = A[2,1]

	Xeq	set obs `n'

	if "`normal'"~="" {
		local invnormal invnormal
	}
	Msg *	Generate the output data for the y2 & y1 variables
	Msg *		First generate two intermediate utility variables with zero correlation, c1 and c2.
	Xeq	gen c1 = 10*abs(`invnormal'(uniform()))
	Xeq gen c2 = 10*abs(`invnormal'(uniform()))
	Msg *		Then use the Cholesky decomposition of the correlation matrix to generate the two output variables.
	Xeq gen y2 = c2
	Xeq lab var y2 "y2 = 10*abs(`invnormal'(uniform()))"
	Ms2 gen y1= 1.2*(_char(96)coef_c1_char(39)*c1 + _char(96)coef_c2_char(39)*c2)
		gen y1= 1.2*(`coef_c1'*c1 + `coef_c2'*c2)
	Xeq sum y1, meanonly
	Msg *	By rounding the randomly generated values of the two output variables, we assure some purely specialized firms.
	Xeq qui replace y1 = round(y1 - r(min) + .0001, .01)
	Xeq qui replace y2 = round(y2, .01)
	local addtxt : di %-9.3f r(min)
	local addtxt " - `addtxt'"
	Msg lab var y1 "y1 = 1.2 * (`di_coef_c1'*c1 + `coef_c2'*c2)`addtxt'"
		lab var y1 "y1 = 1.2 * (`di_coef_c1'*c1 + `coef_c2'*c2)`addtxt'"
	Xeq drop c1 c2

	Msg *	The correlation between the two output variables is:
	Xeq corr y2 y1
		local di_rhohat : di %9.3f r(rho)
	Msg So the correlation coefficient between the two outputs is: `di_rhohat'

	Xeq scatter y2 y1, title("Scatter plot of two output variables")  ///
		subtitle("with correlation coefficient `di_rhohat'") ///
		name(corr_y2y1`rhotxt', replace)

	Msg *	Generate the input price variable, w.  (Think "wage".)
	Xeq gen w=20*abs(uniform())
		Xeq lab var w "input price"
	Xeq gen lnw = ln(w)  // allows flexmake to construct ln_w
		Xeq lab var lnw "ln(input price)"

	Msg *	Specify the parameters which define the cost function
	Xeq global alpha = `alpha'
	Xeq global beta = `beta'  //  Complementarity parameter
	Xeq global scale = `scale'  //  Scale parameter
	Xeq global gamma = `gamma'  // Coefficient of ln(w)
	Xeq global sigma = `sigma'  // Size of disturbance term

	Msg *	Generate the true cost of joint production according to a CES structure

	Ms2 gen true = (_char(36)alpha*(y2^_char(36){beta} + y1^_char(36){beta})^(_char(36){scale}/_char(36){beta})) + _char(36)gamma * lnw
		gen true = (${alpha}*(y2^${beta} + y1^${beta})^(${scale}/${beta})) + $gamma * lnw
	Ms2 lab var true "True cost = (_char(36)alpha*(y2^_char(36)beta + y1^_char(36)beta)^(_char(36)s/_char(36)beta)) "
		lab var true "True cost = ($alpha*(y2^$beta + y1^$beta)^($s/$beta)) "

	Msg *	Add stochastic disturbance term
	di as txt _n "Generate the observed cost of joint production by appending a disturbance"
	Ms2 gen observed = true *  exp(_char(36){sigma}*`invnormal'(uniform()))
		gen observed = true *  exp(${sigma}*`invnormal'(uniform()))
	Ms2 lab var observed "observed = true *  exp(_char(36)sigma*`invnormal'(uniform()))"
		lab var observed "observed = true *  exp($sigma*`invnormal'(uniform()))"

	Xeq label data "Simulated cost function data with rho = ${rho}, beta = `beta', scale = `scale'"

di as txt "End of the data generating process"_n(2) "Here are details of the generated data"

	Xeq describe
	Xeq summarize

	return local rho `rho'  //  Specified
	return local rhohat `rhohat'  //  Estimated on simulated data

end

program define show_cent_true
	local rhotxt = subinstr("${rho}",".","",.)
	NewName true_rho`rhotxt'
	di as inp  ///
		_n(3) "*	Construct the isocost frontiers  at four centiles of the true cost function"  ///
		_n(1) "isoplot y2 y1, isovar(true) centiles ccuts(10 20 50 90) name(true_rho`rhotxt', replace)"
	isoplot y2 y1, isovar(true) centiles ccuts(10 20 50 90) name(true_rho`rhotxt', replace)  ///
		title("Example 1: Display four isocost frontiers", span)
end

program define show_cent_obs
	local rhotxt = subinstr("${rho}",".","",.)
	NewName obs_rho`rhotxt'
	di as inp  ///
		_n(3) "*	In real data, the true cost is perturbed by a stochastic disturbance"  ///
		_n(1) "*	So isocost frontiers for observed data would look more like this"  ///
		_n(1) "isoplot y2 y1, isovar(observed) centiles ccuts(10 20 50 90) name(obs_rho`rhotxt', replace)"
	isoplot y2 y1, isovar(observed) centiles ccuts(10 20 50 90) name(obs_rho`rhotxt')  ///
		title("Example 2: Isocost -bands- are wider"  "with a disturbance term", span)
end

prog def show_jntplot
	local rhotxt = subinstr("${rho}",".","",.)
	NewName true_med_ys_rho`rhotxt'
	di as inp  ///
		_n(3) "*	Using the ccuts() option, we can plot isocost frontiers for any value(s) of cost."  ///
		_n(1) "*	For example, the cost of jointly producing the median amounts of the two outputs"  ///
		_n(1) "*	at median factor prices in a single firm can be found by first calculating"  ///
		_n(1) "*	the medians of the two outputs and of the factor price(s)."
		Xeq centile y1
		global med_y1 : display %6.3f r(c_1)
		Msg global med_y1 = r(c_1)
		Xeq centile y2
		global med_y2 : display %6.3f r(c_1)
		Msg global med_y2 = r(c_1)
		Xeq centile lnw
		global med_lnw : display %6.3f r(c_1)
		Msg global med_lnw  = r(c_1)

	di as inp _n "*	Then we substitute these median values into the cost function:" _n
	Ms2 global cost_of_jnt = (_char(36)alpha*(_char(36)med_y1^_char(36)beta + _char(36)med_y2^_char(36)beta)^(_char(36)s/_char(36)beta)) + _char(36)gamma * _char(36)med_lnw
		global cost_of_jnt = (${alpha}*(${med_y1}^${beta} + ${med_y2}^${beta})^(${scale}/${beta})) + $gamma * ${med_lnw}
		global cost_of_jnt : display %9.4f $cost_of_jnt

	di as inp _n "*	The true cost of jointly producing the two outputs is:"  ///
				_n "di \$cost_of_jnt"
	di as res ${cost_of_jnt}

	di as inp  ///
		_n(2) "*	Using the -plotjoint(median)- option, we can superimpose a point"  ///
		_n(1) "*	representing a firm producing median levels of the two outputs"  ///
		_n(1) "*	onto the isocost frontier for the cost of that joint production."
	Ms2 isoplot y2 y1, isovar(true) plotjoint(median) values ccuts(_char(36)cost_of_jnt) name(true_med_ys_rho`rhotxt', replace)
	isoplot y2 y1, isovar(true) plotjoint(median) values ccuts(${cost_of_jnt}) name(true_med_ys_rho`rhotxt', replace)  ///
		title("Example 3: Joint production of median outputs", span)
end

prog def show_specplot
	local rhotxt = subinstr("${rho}",".","",.)
	NewName true_spec_ys_rho`rhotxt'
	di as inp  ///
		_n(2) "*	By specifying the -plotspec(median)- option, we plot two points"  ///
		_n(1) "*	representing two firms each of which produces the median level"  ///
		_n(1) "*	of one output and zero of the other output.  Thus, the sum of the"  ///
		_n(1) "*	outputs of the two specialized firms equals the total output "  ///
		_n(1) "*	of the single two-product firm plotted in the graph -true_med_ys_rho`rhotxt'-."  ///
		_n(1) "*	The cost of this specialized production will be the sum of the costs"  ///
		_n(1) "*	in the two firms.  For the firm specialized in y1, the cost is:"

	Ms2 "global cost_of_y1 = (_char(36)alpha*(_char(36)med_y1^_char(36)beta + 0^$beta)^(_char(36)s/_char(36)beta)) + _char(36)gamma * _char(36)med_lnw}
		global cost_of_y1 = (${alpha}*(${med_y1}^${beta} + 0^${beta})^(${scale}/${beta})) + $gamma * ${med_lnw}
	Ms2 di _char(36)cost_of_y1
	di as res ${cost_of_y1} _n

	Msg *	The cost for the firm specialized in y2 is:
	Ms2 "global cost_of_y2 = (_char(36)alpha*(0^_char(36)beta + _char(36)med_y2^_char(36)beta)^(_char(36)s/_char(36)beta)) + _char(36)gamma * _char(36)med_lnw"
		global cost_of_y2 = (${alpha}*(0^${beta} + ${med_y2}^${beta})^(${scale}/${beta})) + $gamma * ${med_lnw}
	Ms2 di _char(36)cost_of_y2
	di as res ${cost_of_y2}

	di as inp  ///
		_n(2) "*	Using the -plotspec(median)- option, as well as the"  ///
		_n(1) "*	-plotjoint(median)- option, we plot allthree firms"  ///
		_n(1) "*	superimposed on the isocost frontier for each firm"

	Ms2 isoplot y2 y1, isovar(true)  plotspec(median) plotjoint(median) values ccuts(_char(36)cost_of_y1  _char(36)cost_of_y2 _char(36)cost_of_jnt ) name(true_spec_ys_rho`rhotxt', replace)
	isoplot y2 y1, isovar(true)  plotspec(median) plotjoint(median) values ccuts($cost_of_y1  $cost_of_y2 $cost_of_jnt ) name(true_spec_ys_rho`rhotxt', replace)  ///
		title("Example 4: Isocost frontiers" "for joint and specialized production", span)

end

prog def show_speccost
	local rhotxt = subinstr("${rho}",".","",.)
	NewName true_all_rho`rhotxt'
	di _n(2) "*	The total cost of producing the median amounts of both outputs is the sum of the costs of the two firms:"
	di _n(1) as inp "global cost_of_spec = \$cost_of_y1  +  \$cost_of_y2"
		global cost_of_spec = $cost_of_y1  +  $cost_of_y2

	Ms2 di _char(36)cost_of_spec
	di as res ${cost_of_spec}

	di as inp  ///
		_n(2) "*	Adding an isocost frontier at this higher level of cost"  ///
		_n(1) "*	shows graphically the far larger set of output combinations"  ///
		_n(1) "*	that can be produced in a single firm producing both outputs"

	di _n(2) "isoplot y2 y1, isovar(true) plotspec(median) plotjoint(median) values ccuts(\$cost_of_y1  \$cost_of_y2 \$cost_of_jnt \$cost_of_spec) name(true_all_rho`rhotxt', replace)"
	isoplot y2 y1, isovar(true) plotspec(median) plotjoint(median) values ccuts($cost_of_y1  $cost_of_y2 ${cost_of_jnt} ${cost_of_spec}) name(true_all_rho`rhotxt', replace)  ///
		title("Example 5: Cost savings from joint production", span)  ///
		subtitle("Outer frontier shows output combinations achievable"  ///
			"through joint production which cost the same"  ///
			"as producing both outputs in specialized firms", span)

end
*	Version 1.0  30Sep2016 for version 1.2 of isoplot.ado
*	Version 1.1  1Oct2016 Switch roles of y2 and y1
*	Version 1.3  15Oct2016 to accompany version 1.5 of -isoplot-