program define margcost, rclass
*!	version 1.7 30Apr2018  by Mead Over, CGD

	syntax varname [if] [in], cmdline(passthru)  ///
		[ chain(passthru) chainlhs(passthru) chainwrt(passthru) lambda(passthru) force  ///
		replace post postname(string) addest verbose colnames(string) KEEPhat ///
		fvset1(varname fv) fvsetto1(real 1) fvset2(varname fv) fvsetto2(real 1) * ]

	di _n as txt "Executing -margcost- with respect to the variable -`varlist'-" ///
		_n "**************************************************************"

	marksample touse

	if "`chainrhs'"~= "" {
		di as err "The option -chainrhs- is not a valid option for the program -margcost-."  ///
			_n "Use the option -chainwrt- instead."
			exit 198
	}

	local wrt `varlist'
	if ("`chain'"~="" | "`chainwrt'"~="" ) {
		if strpos("`wrt'","_")==0 {
			di as err "The variable name {it:`wrt'} does not contain an underscore character."
			exit 198
		}
		else {
			local unwrt = substr("`wrt'",strpos("`wrt'","_")+1,.)
		}
	}
	else {
		local unwrt `wrt'
	}

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

	`shhh' firstpartial `wrt', `cmdline' `chain' `chainlhs' `chainwrt' `lambda' `force'
		`shhh' return list

		if "`r(chainlhs)'"~="" {
			local chainlhs `r(chainlhs)'
			local translhs `r(chainlhs)'_
		}
		if "`r(chainwrt)'"~="" {
			local chainwrt `r(chainwrt)'
			local transwrt `r(chainwrt)'_
		}

		local tokenswowrt `r(tokenswowrt)'
		local wrt `r(wrt)'
		local eqname `r(eqname)'
		local expression `r(expression)'
		local tokenswithwrt `r(tokenswithwrt)'
		local genexpression `r(genexpression)'
		local margexpression `r(margexpression)'
		local depvar `r(depvar)'

	foreach i of numlist 1/2 {
		if "`fvset`i''"~=""  {
			local newfvset`i' = subinstr("`fvset`i''","i.","`fvsetto`i''.",1)
			di as txt "The factor variable " as res "`fvset`i''" as txt " is set to " as res "`newfvset`i''"  ///
				as txt " in all expressions."
			local expression = subinstr(`"`expression'"',"`fvset`i''","`newfvset`i''",.)
			local genexpression = subinstr(`"`genexpression'"',"`fvset`i''","`newfvset`i''",.)
			local margexpression = subinstr(`"`margexpression'"',"`fvset`i''","`newfvset`i''",.)
		}
	}

	//	Joint test of the hypothesis that -unwrt- contributes to explaining
	//	variance of the dependent variable, total cost.
	di as txt _n "For all coefficients involving " as res "`wrt'" as txt ", a joint test of significance yields:"
	testparm `tokenswithwrt'
		local p_`unwrt' = `r(p)'
		if "`addest'"~="" `shhh' addest, name("pval_`unwrt'") value(`p_`unwrt'')
		if "`r(chi2)'"~="" {
			local Chi2_`unwrt' = `r(chi2)'
			if "`addest'"~="" `shhh' addest, name("Chi2_`unwrt'") value(`Chi2_`unwrt'')
		}
		if "`r(F)'"~="" {
			local F_`unwrt' = `r(F)'
			if "`addest'"~="" `shhh' addest, name("F_`unwrt'") value(`F_`unwrt'')
		}

	//	Compute the marginal cost of -unwrt- at each observation
	//	using -firstpartial-'s -r(genexpression)-
		tempvar d_depvar_d_unwrt
		if "`replace'"~="" {
			capture drop `depvar'_hat
			capture drop `d_depvar_d_unwrt'
		}
		if "`e(k_eq)'"~="" & real("`e(k_eq)'")>1 & real("`e(k_eq)'")<. {
			local eq eq(`depvar')
		}
		`shhh' predict `depvar'_hat, xb `eq'

		//	Use -firstpartial-'s returned "genexpression, unless its missing.
		if `"`genexpression'"'~="" {
			capture noi `shhh' gen `d_depvar_d_unwrt'  = `genexpression'
		}
		else {
			capture noi `shhh' gen `d_depvar_d_unwrt'  = `expression'
		}
		if _rc~=0 {
			di as err "Cannot compute a marginal cost at each observation."
			local skip yes
		}
		else {
			lab var `d_depvar_d_unwrt' "Estimated marginal cost at each observation"
		}
		if "`keephat'"=="" {
			capture drop `depvar'_hat
		}

	//	Compute the mean, median, p25 and p75 of the marginal cost over selected observations
	if "`skip'"=="" {
		`shhh'  summ `d_depvar_d_unwrt' if `touse' & e(sample), detail
			local p25_MC_`unwrt' = `r(p25)'
			if "`addest'"~="" `shhh' addest, name("p25_MC_`unwrt'") value(`p25_MC_`unwrt'')
			local med_MC_`unwrt' = `r(p50)'
			if "`addest'"~="" `shhh' addest, name("med_MC_`unwrt'") value(`med_MC_`unwrt'')
			local p75_MC_`unwrt' = `r(p75)'
			if "`addest'"~="" `shhh' addest, name("p75_MC_`unwrt'") value(`p75_MC_`unwrt'')
			local mn_MC_`unwrt' = `r(mean)'
			if "`addest'"~="" `shhh' addest, name("mn_MC_`unwrt'") value(`mn_MC_`unwrt'')
	}

	if "`post'"~="" {
		tempname existing
		est store `existing'
	}

	//	Marginal cost at the user-selected values of all r.h.s. variables
	//		- Exclude observations on which the firstpartial is undefined.
	//		- Use -firstpartial-'s returned "margexpression, unless its missing.
		di _n as txt "Marginal cost of " as res "`unwrt'" as txt " :"  ///
			_n as txt "-----------------------------"
		if "`skip'"=="" {
			local ifgenexp & `d_depvar_d_unwrt'<.
		}
		if `"`margexpression'"'~="" {
			`shhh' margins if `touse' `ifgenexp', vsquish  ///
				expression(`margexpression') `post' `options'
		}
		else {
			`shhh' margins if `touse' `ifgenexp', vsquish  ///
				expression(`expression') `post' `options'
		}
		tempname margb_mat margV_mat margout_mat margp_mat
		matrix define `margb_mat' = r(b)
			local margb = el(`margb_mat',1,1)
		matrix define `margV_mat' = r(V)
		matrix define `margout_mat' = r(table)
		matrix define `margp_mat' = `margout_mat'[4,1...]
			local margp = el(`margp_mat',1,1)
		if "`post'"~="" {
			tempname margeb margeV
			matrix define `margeb' = e(b)
			matrix define `margeV' = e(V)
			if "`postname'"~="" {
				est store `postname'
				di as txt "Margins results stored as: " as res "`postname'"
			}
			est restore `existing'
		}
		if "`addest'"~="" `shhh' addest, name("MC_`unwrt'") value(`margb')
		if "`addest'"~="" `shhh' addest, name("p_MC_`unwrt'") value(`margp')

		if "`shhh'"=="quietly" {
			di 	   as txt "As estimated by Stata's -margins- command, the partial derivative of total cost" ///
				_n as txt "with respect to " as res "`unwrt'" as txt " is " as res `margb' as txt " with p-value " as res `margp'
		}

		local cfnms : colfullnames `margb_mat'

		if "`colnames'"~="" {
			local ncfnms = wordcount(`"`cfnms'"')
			local ncolnames = wordcount(`"`colnames'"')
			if `ncfnms' ~= `ncolnames' {
				di as err "Margins result matrix has " as txt "`ncfnms'" as err " columns."  ///
					_n as err "-colnames()- option has " as txt "`ncolnames'" as err " entries."  ///
					_n as err "The two should match."
					exit 198
			}
			else {
				local cfnms `colnames'
			}
		}
		matrix roweq `margout_mat' = `unwrt'
		matrix colnames `margout_mat' = `cfnms'

	//	Results from -firstpartial-
	return local chain `chain'
	return local chainwrt `chainwrt'
	return local chainlhs `chainlhs'
	return local tokenswowrt `tokenswowrt'
	return local wrt `wrt'
	return local depvar `depvar'
	return local eqname `eqname'
	return local tokenswithwrt `tokenswithwrt'
	return local margexpression `margexpression'
	return local genexpression `genexpression'
	return local expression `expression'

	//	Results of evaluating the derivative
	return matrix b_marg = `margb_mat'
	return matrix V_marg = `margV_mat'
	return matrix p_marg = `margp_mat'
	return matrix margout = `margout_mat'
	if "`post'"~="" {
		return matrix eb_marg = `margeb'
		return matrix eV_marg = `margeV'
	}
	//	Results of the joint test on all coefficients
	return scalar b_marg = `margb'
	return scalar p_marg = `margp'
	return scalar p25_MC_`unwrt' = `p25_MC_`unwrt''
	return scalar med_MC_`unwrt' = `med_MC_`unwrt''
	return scalar p75_MC_`unwrt' = `p75_MC_`unwrt''
	return scalar mn_MC_`unwrt' = `mn_MC_`unwrt''
	if "`Chi2_`unwrt''"~="" {
		return scalar Chi2_`unwrt' = `Chi2_`unwrt''
	}
	if "`F_`unwrt''"~="" {
		return scalar F_`unwrt' = `F_`unwrt''
	}
	return scalar p_`unwrt' = `p_`unwrt''

end
*	version 1.0 6nov2016	Originally called -margcosttest.ado-
*	version 1.1 13jan2017	First published version
*	version 1.2 10Feb2017	Change option boxcox(passthru) to lambda(passthru)
*	version 1.3 20Oct2017	Add shhh to the -margins- command, so display of
*		its full results requires the -verbose- option
*	version 1.4 28Mar2018	When -verbose- is omitted, display (only) the margins
*		estimate and p-value.
*		Add margins output to returned scalars and results stored by -addest-.
*		Fixed bug in option -post-. Add error trap for invalid -chainrhs- option
*	version 1.5 28Mar2018	Add the p25 and p75 estimates of MC to saved output
*	version 1.6 14Apr2018	Since -firstpartial- only returns a "genexpression"
*		when it would be different from "expression", in lines 80-85 allow
*		-margcost- to use "expression" when "genexpression" is blank.
*		Allow the -generate- command with either "expression" or "genexpression"
*		to fail when these expressions contain facyor variables.
*	version 1.7 30Apr2018	Add options -fvset1- and -fvsetto1- to specify
*		replacement "i." in factor variables with a value set by -fvsetto1".
*		Allow the sample to be defined by [if] [in] as well as e(sample)
*	********************************************************************
*	The command : margins i.wave##i.Intervention, dydx(sqrt_srh) expression(exp(xb())/2*sqrt(srh))
*	seems to unwind the reg ln_tc ...........  sqrt_srh  i.wave##i.Intervention#c.sqrt_srh
*	without resorting to the -firstpartial- command.  Trick is to use sqrt(srh) rather than
*	sqrt_srh in order to avoid -margins- taking the derivative of sqrt_srh in the denominator.