prog define flexematrix, rclass
*!	program flexematrix		version 1.4		4May2018
	syntax varlist, [cmdline(string) shares(string) actual verbose USEShareeqs]

*	Syntax 1:
*		flexematrix wrt1 wrt2 , cmdline(string)
*	  		wrt1, wrt2			Two price variables on the rhs of the cost function
*			cmdline(string)		Command line specifying cost function just estimated
*		This syntax is mysterious ...
*	Syntax 2:
*		flexematrix wrt1 wrt2 , cmdline(string) shares(sh1 sh2)
*	  		wrt1, wrt2			Two price variables on the rhs of the cost function
*	  		sh1, sh2			Two share variables in the data that correspond 1:1 with (wrt1, wrt2)
*			cmdline(string)		Command line specifying cost function just estimated
*		This syntax produces the matrix of substitution elasticities
*	Syntax 3:
*		flexematrix wrt1 wrt2 , cmdline(string) shares(sh1 sh2 sh3) useshareeqs
*	  		wrt1, wrt2, wrt3	A mix of K-1 price variables used in estimation with one
*									on which data was normalized (and thus is not used in the regression)
*	  		sh1, sh2, sh3		Share variables in the data that correspond 1:1 with (wrt1, wrt2, wrt3)
*			cmdline(string)		Command line specifying cost function just estimated
*		This syntax produces the matrix of substitution elasticities which can include those derived
*		from a Kth equation added by the -addest- command

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

	if `"`cmdline'"'=="" {
		if `"`e(cmdline)'"'=="" {
			di as err "The option cmdline(string) is omitted and e(cmdline) is empty."  ///
				_n "When the user omits the cmdline(string) option, this program expects"  ///
				_n "to find a command line stored in e(cmdline) by a previous estimation command."
			exit 198
		}
		else {
			local cmdline `e(cmdline)'
		}
	}

	local nprices : word count `varlist'
	mat def elasticities = J(`nprices',`nprices',.)
	mat rownames elasticities = `varlist'
	mat colnames elasticities = `varlist'

	mat def se_elasticities = J(`nprices',`nprices',.)
	mat rownames se_elasticities = `varlist'
	mat colnames se_elasticities = `varlist'

	mat def p_elasticities = J(`nprices',`nprices',.)
	mat rownames p_elasticities = `varlist'
	mat colnames p_elasticities = `varlist'

	if "`shares'"~="" {
		local nshares : word count `shares'
		if `nprices'~=`nshares' {
			di _n as err "If -shares- option is selected, number of shares must equal the number of prices"
			exit 198
		}
		tokenize `shares'
	}

	 if "`useshareeqs'"=="" {
		local restofprices `varlist'
		local sh1 = 1
		foreach p1  of varlist `varlist' {
			local restofprices : list restofprices - p1
			local sh2 = `sh1'
			foreach p2 of varlist `p1' `restofprices' {
				local wrt1 ln_`p1'
				local wrt2 ln_`p2'
				if "`shares'"~="" {
					local matchedshares shares(``sh1'' ``sh2'')
				}

				`shhh' di "flexelasticity `wrt1' `wrt2', cmdline(`cmdline') `matchedshares' `actual'"
				`shhh' flexelasticity `wrt1' `wrt2', cmdline(`cmdline') `matchedshares' `actual'
				`shhh' di _n(2)
				if "`matchedshares'"~="" & "`actual'"~="" {
					`shhh' di as txt "Using actual shares, " _c
				}
				if "`matchedshares'"~="" & "`actual'"=="" {
					`shhh' di as txt "Using fitted shares, " _c
				}
				`shhh' di as txt "elasticity between inputs with prices " as res "`p1'" as txt " and " as res "`p2'" ///
					_n as txt "    Elasticity = " as res `r(sigma)' as txt " with P-value = " as res `r(p)'

				mat def elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(sigma)'
				mat def se_elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(se_sigma)'
				mat def p_elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(p_sigma)'
				mat def elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(sigma)'
				mat def se_elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(se_sigma)'
				mat def p_elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(p_sigma)'
				local sh2 = `sh2'+1
			}
			local sh1 = `sh1'+1
		}
	} // End of code block

	if "`useshareeqs'"=="useshareeqs"  {  //  Code for using coefficients from share equations instead of from the cost function
		local restofprices `varlist'
		local sh1 = 1
		foreach p1  of varlist `varlist' {
			local restofprices : list restofprices - p1
			local sh2 = `sh1'
			foreach p2 of varlist `p1' `restofprices' {
				di "Variables are: `p1' and `p2' and shares are: ``sh1'' and ``sh2''"  /* DEBUG */
				local wrt1 ln_`p1'
				local wrt2 ln_`p2'
				if "`shares'"~="" {
					local matchedshares shares(``sh1'' ``sh2'')
				}
				capture confirm numeric variable `wrt1', exact
				if _rc~=0 {
					local wrt1 ``sh1''
				}
				capture confirm numeric variable `wrt2', exact
				if _rc~=0 {
					local wrt2 ``sh2''
				}

				`shhh' flexelasticity `wrt1' `wrt2', cmdline(`cmdline') `matchedshares' `actual'
				`shhh' di _n(2)
				if "`matchedshares'"~="" & "`actual'"~="" {
					`shhh' di as txt "Using actual shares, " _c
				}
				if "`matchedshares'"~="" & "`actual'"=="" {
					`shhh' di as txt "Using fitted shares, " _c
				}
				`shhh' di as txt "elasticity between inputs with prices " as res "`p1'" as txt " and " as res "`p2'" ///
					_n as txt "    Elasticity = " as res `r(sigma)' as txt " with P-value = " as res `r(p)'

				mat def elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(sigma)'
				mat def se_elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(se_sigma)'
				mat def p_elasticities[colnumb(elasticities,"`p2'"),colnumb(elasticities,"`p1'")] = `r(p_sigma)'
				mat def elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(sigma)'
				mat def se_elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(se_sigma)'
				mat def p_elasticities[colnumb(elasticities,"`p1'"),colnumb(elasticities,"`p2'")] = `r(p_sigma)'
				local sh2 = `sh2'+1
			}
			local sh1 = `sh1'+1
		}
	}

	di _n as txt "Estimated Allen Elasticities of Substitution:"
	matrix list elasticities, format(%12.4f)
	di _n as txt "Standard errors of estimated Allen Elasticities of Substitution:"
	matrix list se_elasticities, format(%12.4f)
	di _n as txt "p-values of tests of Ho: Estimated Allen Elasticities of Substitution = 0:"
	matrix list p_elasticities, format(%12.4f)
	return matrix elasticities = elasticities
	return matrix se_elasticities = se_elasticities
	return matrix p_elasticities = p_elasticities
end

*	Based on glematrixversion 1.0		6/6/2012
*	Version 1.1  Add capability to extract coefficients from share equations
*	Based on quadematrix 1.1		7/17/2012
*	Version 1.2  8/8/2017 Allow e(cmdline) to substitute for the option cmdline(string)
*	Version 1.4  4May2018 Move the command local wrt1 from line 71 to outside the loop to line 69