*!Version 1.2 14Dec2016 by Mead Over, Center for Global Devleopment
*!	Companion program to be executed by the help file: firstpartial.sthlp

program define firstpartial_example_boxcox
	version 11.1

	capture which flexmake
	if _rc~=0 {
		di as err "Running this example requires that the user-written program -{search flexmake}- already be installed."  ///
		_n "User can install this program by clicking on the above link and installing the flexmake program"  ///
		_n "from CGD's Stata library at " as res "net from http://www.cgdev.org/doc/stata."
		exit 198
	}

	set more off
	`0'

	capture prog drop firstpartial
end

program Msg
	di as txt
	di as txt "-> " as res `"`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 res _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 setup

	*	Load the data from Bill Greene's website
	capture import delimited using "http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF4-4.txt", delimiters(" ", collapse) varn(2) clear
	if _rc~=0 {
		Msg import delimited using "http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF4-4.txt", delimiters(" ", collapse) varn(2) clear
		di as err "no; data in memory would be lost"
		di as err "First {stata clear:clear} your memory before running these examples"
		exit
	}

	Xeq set more off
	Xeq import delimited using "http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF4-4.txt", delimiters(" ", collapse) varn(2) clear

	Xeq drop v11
	Xeq gen plbypf = pl/pf
	Xeq 	lab var plbypf "pl/pf"
	Xeq gen pkbypf = pk/pf
	Xeq 	lab var pkbypf "pk/pf"
	Xeq gen costbypf = 1000000*cost/pf
	Xeq 	lab var costbypf "1000000*cost/pf"

	Xeq *	Example using the boxcox option in flexmake.
	Xeq *	This option replaces ln_q with bc_q

	Xeq *	Here are the revised set of cross-equation constraints
	Xeq constraint 1  _b[ln_costbypf:c.ln_pkbypf] = _b[sk:_cons]
	Xeq constraint 2  _b[ln_costbypf:c.ln_plbypf] = _b[sl:_cons]
	Xeq constraint 3  _b[ln_costbypf:c.ln_pkbypf#c.ln_pkbypf#c.__half] = _b[sk:ln_pkbypf]
	Xeq constraint 4  _b[ln_costbypf:c.ln_plbypf#c.ln_plbypf#c.__half] = _b[sl:ln_plbypf]
	Xeq constraint 5 _b[ln_costbypf:c.ln_pkbypf#c.ln_plbypf] = _b[sl:ln_pkbypf]
	Xeq constraint 6 _b[ln_costbypf:c.ln_pkbypf#c.ln_plbypf] = _b[sk:ln_plbypf]
	Xeq constraint 7 _b[ln_costbypf:c.bc_q#c.ln_pkbypf] = _b[sk:bc_q]
	Xeq constraint 8 _b[ln_costbypf:c.bc_q#c.ln_plbypf] = _b[sl:bc_q]

	*	Define the variables to hold the results of iterated execution of sureg
	Xeq gen rss_1 = .
	Xeq 	lab var rss_1 "Sum of squared residuals for the cost equation"
	gen iter = .
	Xeq 	lab var iter "Iteration number of search for best lambda"
	Xeq gen ll = .
	Xeq 	lab var ll "Log-likelihood conditional on lambda"
	Xeq gen lambda = .
	Xeq 	lab var lambda "Boxcox parameter at which equations are estimated"


end

program define iterate_boxcox
	local iter = 1
	quietly {
		foreach lambda of numlist .01(.01).2 .2(.1).9  {

			replace lambda = `lambda' if _n == `iter'

			noi di _n(2) as txt "----------------------------------------------------------"
			noi di _n as txt "Execute -flexmake- with lambda = " as res "`lambda'"  ///
			_n as inp "flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) ln boxcox(`lambda') replace globals half nomake"  ///
			_n as txt "(Output omitted)"

			flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) ln boxcox(`lambda') replace globals half nomake

			noi di _n as txt "Execute -sureg- with lambda = " as res "`lambda'"
			Ms2 sureg (_char(36)depvar _char(36)ownprices _char(36)ownprices2 _char(36)crossprices _char(36)ownoutputs _char(36)ownoutputs2 _char(36)crossoutputs _char(36)outputsxprices ) (sk _char(36)ownprices _char(36)ownoutputs ) (sl _char(36)ownprices _char(36)ownoutputs ) , constraints(1/8)
			sureg ($depvar $ownprices $ownprices2 $crossprices $ownoutputs $ownoutputs2 $crossoutputs $outputsxprices ) (sk $ownprices $ownoutputs ) (sl $ownprices $ownoutputs ) , constraints(1/8)
			noi di as txt "(Output omitted)"


			replace rss_1 = `e(rss_1)'  if _n==`iter'
			replace iter = `iter'       if _n==`iter'
			replace ll   = `e(ll)'      if _n==`iter'
			local iter = `iter' + 1
		}
	}

	Msg *	Here are the fit criteria across the iterations
	Xeq list iter lambda rss_1 ll if iter<.

end

program define findbestfit
	NewName loglkhd
	Xeq preserve
	Xeq 	gsort -ll
	Xeq 	global maxllit = iter[1]
	Xeq 	global maxlllambda = round(lambda[1],.01)
	Xeq capture noi restore

	Xeq di as txt _n "The best fitting boxcox parameter lambda is: " as res "${maxlllambda}" as txt " achieved on iteration " as res "${maxllit}"

	Xeq *	We can graph the log-likelihood criterion against the values of lambda
	Xeq twoway line ll lambda, sort name(loglkhd, replace)

end

program define reestimate_best_fits
	quietly {
		local boxcostbest
		foreach lambda of numlist ${maxlllambda} .12  {
			noi di _n(2) as txt "----------------------------------------------------------"
			noi di _n as txt "Execute -flexmake- with lambda = " as res "`lambda'"  ///
			_n as inp "flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) ln boxcox(`lambda') replace globals half nomake"  ///
			_n as txt "(Output omitted)"

			flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) ln boxcox(`lambda') replace globals half nomake

			noi di _n as txt "Execute -sureg- with lambda = " as res "`lambda'"
			Ms2 sureg (_char(36)depvar _char(36)ownprices _char(36)ownprices2 _char(36)crossprices _char(36)ownoutputs _char(36)ownoutputs2 _char(36)crossoutputs _char(36)outputsxprices ) (sk _char(36)ownprices _char(36)ownoutputs ) (sl _char(36)ownprices _char(36)ownoutputs ) , constraints(1/8)
			noi di as txt "(Output omitted)"

			sureg ($depvar $ownprices $ownprices2 $crossprices $ownoutputs $ownoutputs2 $crossoutputs $outputsxprices ) (sk $ownprices $ownoutputs ) (sl $ownprices $ownoutputs ) , constraints(1/8)
			if "`boxcostbest'"=="" {
				est store BoxcoxBest
				local boxcostbest full
			}
			else {
				est store Boxcox2ndBest
			}
		}
	}
end

program show_table
	Xeq	*	Display the translog, the best boxcox, and the 2nd best boxcox
	Xeq estimates table Translog BoxcoxBest Boxcox2ndBest, t stat(N r2_1 r2_2 r2_3 ll)
end
*	Version 1.0 15 Aug2016 by Mead Over, Center for Global Devleopment
*	Version 1.2 14Dec2016	Add the -ln- option to flexmake when boxcox is specified.
*		Change the prefix on boxcox-transformed variables from "boxcox_" to "bc_".


/*
*	Setup: use data from Christensen and Greene(1976)
import delimited using "http://people.stern.nyu.edu/wgreene/Text/Edition7/TableF4-4.txt", delimiters(" ", collapse) varn(2) clear
drop v11
genl plbypf = pl/pf
genl pkbypf = pk/pf
genl costbypf = 1000000*cost/pf

*	Example using the boxcox option in flexmake.
*	This option replaces ln_q with bc_q

*	Here are the revised set of cross-equation constraints
constraint 1  _b[ln_costbypf:c.ln_pkbypf] = _b[sk:_cons]
constraint 2  _b[ln_costbypf:c.ln_plbypf] = _b[sl:_cons]
constraint 3  _b[ln_costbypf:c.ln_pkbypf#c.ln_pkbypf#c.__half] = _b[sk:ln_pkbypf]
constraint 4  _b[ln_costbypf:c.ln_plbypf#c.ln_plbypf#c.__half] = _b[sl:ln_plbypf]
constraint 5 _b[ln_costbypf:c.ln_pkbypf#c.ln_plbypf] = _b[sl:ln_pkbypf]
constraint 6 _b[ln_costbypf:c.ln_pkbypf#c.ln_plbypf] = _b[sk:ln_plbypf]
constraint 7 _b[ln_costbypf:c.bc_q#c.ln_pkbypf] = _b[sk:bc_q]
constraint 8 _b[ln_costbypf:c.bc_q#c.ln_plbypf] = _b[sl:bc_q]

*	We iterate across values of the boxcox parameter (here called lambda)
*	We then choose the estimation result with the highest log-likelihood
gen rss_1 = .
	lab var rss_1 "Sum of squared residuals for the cost equation"
gen iter = .
	lab var iter "Iteration number of search for best lambda"
gen ll = .
	lab var ll "Log-likelihood conditional on lambda"
gen lambda = .
	lab var lambda "Boxcox parameter at which equations are estimated"
local iter = 1
quietly {
	foreach lambda of numlist .01(.01).2 .2(.1).9  {

		replace lambda = `lambda' if _n == `iter'
		flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) ln boxcox(`lambda') replace globals half

		sureg ($depvar $ownprices $ownprices2 $crossprices $ownoutputs $ownoutputs2 $crossoutputs $outputsxprices )   ///
			(sk $ownprices $ownoutputs ) (sl $ownprices $ownoutputs ) , constraints(1/8)

		replace rss_1 = `e(rss_1)'  if _n==`iter'
		replace iter = `iter'       if _n==`iter'
		replace ll   = `e(ll)'      if _n==`iter'
		local iter = `iter' + 1
	}
}

*	Here are the fit criteria across the iterations
list iter lambda rss_1 ll if iter<.

preserve
	gsort -ll
	local maxllit = iter[1]
	local maxlllambda = round(lambda[1],.01)
restore

di as txt _n "The best fitting boxcox parameter lambda is: " as res "`maxlllambda'"  ///
	as txt " achieved on iteration " as res "`maxllit'"

*	We can graph the log-likelihood criterion against the values of lambda
twoway line ll lambda, name(loglkhd, replace)

*	Now re-estimate the cost function conditional on the best fitting value of lambda
flexmake costbypf, outputs(q) prices(pkbypf plbypf) shares(sk sl) boxcox(`maxlllambda') replace globals half

sureg ($depvar $ownprices $ownprices2 $crossprices $ownoutputs $ownoutputs2 $crossoutputs $outputsxprices )   ///
	(sk $ownprices $ownoutputs ) (sl $ownprices $ownoutputs ) , constraints(1/8)

est store Boxcox

firstpartial bc_q , cmdline(`e(cmdline)') chainwrt(boxcox) chainlhs(ln) lambda(`maxlllambda')
	local margexpression `r(margexpression)'
	local genexpression `r(genexpression)'

sum bc_q, det
margins , expression(`margexpression') atmeans at(bc_q=(`r(mean)' `r(p10)' `r(p25)' `r(p50)' `r(p75)' `r(p90)'))

marginsplot, xlabel(,angle(45) format(%9.3f)) ytitle("Marginal cost of output") name(boxcox, replace)

est table Translog Boxcox, t stat(N aic r2_1 rss_1)

end
*/