*! version 1.1.0  Sep 2021



program quaidsce_estat

	version 12

	if "`e(cmd)'" != "quaidsce" {
		exit 301
	}

	gettoken key 0 : 0, parse(", ")
	local lkey = length(`"`key'"')

	if `"`key'"' == substr("expenditure", 1, max(3, `lkey')) {
		DoExp `0'
	}
	else if `"`key'"' == substr("uncompensated", 1, max(6, `lkey')) {
		DoUncomp `0'
	}
	else if `"`key'"' == substr("compensated", 1, max(4, `lkey')) {
		DoComp `0'
	}
	else {
		di as error "invalid subcommand `key'"
		exit 321
	}

end


/// EXPENDITURE ELASTICITY

program DoExp, rclass

	syntax [if] [in]

	marksample touse

	if "`e(lnprices)'" != "" {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(lnprices)' {
			local lnp`i' `var'
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}
	else {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(prices)' {
			tempvar lnp`i'
			gen double `lnp`i'' = ln(`var')
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}

	if "`e(lnexpenditure)'" != "" {
		local lnexp `e(lnexpenditure)'
	}
	else {
		tempvar exp
		gen double `exp' = ln(`e(expenditure)')
		local lnexp `exp'
	}

	if "`e(lhs)'" != "" {
		local i 1
		foreach var of varlist `e(lhs)' {
			local w`i' `var'
			local `++i'
		}
	}

	//PROCESSING VARIABLES (MARGINS WONT WORK DUE TO PARAMETER NAMES, NLCOM DOES NOT PRODDUCE SE)

	local ndemo = `e(ndemos)'
	quietly {
		if "`e(censor)'" == "censor" {
		foreach x of varlist w* `lnpr' cdfw* pdfw* `lnexp' duw* `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
		else {
		foreach x of varlist w* `lnpr' `lnexp' `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
	}

	//AIDS
	tempname lnpindex
	scalar `lnpindex'= `e(anot)'
	forvalue i=1/`=e(ngoods)' {
		scalar `lnpindex'= `lnpindex' + _b[alpha:alpha_`i']*`lnp`i''m
		forvalue j=1/`e(ngoods)' {
			if `j'>=`i' {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`j'_`i']*(`lnp`i''m*`lnp`j''m))
			}
			 else {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`i'_`j']*(`lnp`i''m*`lnp`j''m))
			}
		}
	}

	forvalue i=1/`=e(ngoods)' {
			tempname gsum`i'
			scalar `gsum`i''= 0
			local j=`i'
			forvalue ii=`j'/`=e(ngoods)' {
				scalar `gsum`i''= `gsum`i'' + _b[gamma: gamma_`ii'_`i']*`lnp`i''m
				}
				}


	//When quadratics
	if "`e(quadratic)'" == "quadratic" {
		tempname bofp
		scalar `bofp'= _b[beta:beta_1]*`lnp1'm
		forvalues i = 2/`=e(ngoods)' {
			 scalar `bofp'= `bofp' + _b[beta:beta_`i']*`lnp`i''m
		}
	}

	//When demographics
	if `ndemo' > 0 {
		tempname cofp mbar
		scalar `cofp'= 1 //It is OK to set 1 because below we set a multiplication
		scalar `mbar'= 1 //It is OK because I need to add a "1"

		forvalue i=1/`=e(ngoods)' {
			tempname betanz`i' cofp`i'
			scalar `betanz`i''=_b[beta:beta_`i']
			scalar `cofp`i''= 0
		}


		foreach var of varlist `e(demographics)' {
			forvalue i=1/`=e(ngoods)' {
				 scalar `cofp`i''= `cofp`i''+(`var'm*_b[eta:eta_`var'_`i'])
				 scalar `betanz`i''=`betanz`i''+(`var'm*_b[eta:eta_`var'_`i'])
									}
		scalar `mbar'= `mbar' + (_b[rho:rho_`var']*`var'm)
		}



			forvalue i=1/`=e(ngoods)' {
				scalar `cofp`i''= `cofp`i''*`lnp`i''m
				scalar `cofp'= `cofp'*`cofp`i''
			}
	}

	//FUNCTION EVALUATOR (PREDICTED SHARE)
	forvalues i = 1/`=e(ngoods)' {
		if `ndemo' == 0 {
		scalar we`i' = _b[alpha:alpha_`i']+_b[beta:beta_`i']*(`lnexp'm-`lnpindex') + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex')^2
			}
		}
		else {
		scalar we`i' = _b[alpha:alpha_`i']+`betanz`i''*(`lnexp'm-`lnpindex'-ln(`mbar')) + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar'))^2
				}
			}
		//When censor
		if "`e(censor)'" == "censor" {
		 scalar we`i' = we`i'*cdfw`i'm + _b[delta:delta_`i']*pdfw`i'm
		}
	}

	//FUNCTION EVALUATOR (ELASTICITY)
	forvalues i = 1/`=e(ngoods)' {
		if `ndemo' == 0 {
			global ie`i' "(1+_b[beta:beta_`i']/w`i'm)"
			if "`e(quadratic)'" == "quadratic" {
				 global ie`i' "(1+1/w`i'm*(_b[beta:beta_`i']+2*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')))"
			}
		}
		else {
			global ie`i' "1+`betanz`i''/w`i'm"
			if "`e(quadratic)'" == "quadratic" {
				 global ie`i' "(1+1/w`i'm*(`betanz`i''+2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))))"
			}
		}
			//When censor
		if "`e(censor)'" == "censor" {
			global ie`i' "(1+1/we`i'*((cdfw`i'm*((${ie`i'}-1)*w`i'm))+_b[tau:M_`i']*pdfw`i'm*(w`i'm-_b[delta:delta_`i']*duw`i'm)))"
				// "${ie`i'}-1" --> to remove the "+1" from original formula
				// "*w`i'm" --> to revert 1/w`i'm from the original formula
		}

	}

	tempname elasi sei
	mat elasi=J(1,`e(ngoods)',0)
	mat sei=J(1,`e(ngoods)',0)
	forvalues i = 1/`=e(ngoods)' {
		qui nlcom(${ie`i'})
		mat elasi[1,`i'] = r(b)
		mat sei[1,`i'] = r(V)
		mat sei[1,`i'] = sqrt(sei[1,`i'])
	}
	ret mat elas_i elasi
	ret mat se_elas_i sei

end


/// UNCOMPENSATED PRICE ELASTICITY

program DoUncomp, rclass

	syntax [if] [in]

	marksample touse

	if "`e(lnprices)'" != "" {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(lnprices)' {
			local lnp`i' `var'
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}
	else {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(prices)' {
			tempvar lnp`i'
			gen double `lnp`i'' = ln(`var')
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}

	if "`e(lnexpenditure)'" != "" {
		local lnexp `e(lnexpenditure)'
	}
	else {
		tempvar exp
		gen double `exp' = ln(`e(expenditure)')
		local lnexp `exp'
	}


	if "`e(lhs)'" != "" { // This doesnt work if there is no censoring (what is it for?)
		local i 1
		foreach var of varlist `e(lhs)' {
			local w`i' `var'
			local `++i'
		}
	}

	//PROCESSING VARIABLES (MARGINS WONT WORK DUE TO PARAMETER NAMES, NLCOM DOES NOT PRODDUCE SE)

	local ndemo = `e(ndemos)'
	quietly {
		if "`e(censor)'" == "censor" {
		foreach x of varlist w* `lnpr' cdfw* pdfw* `lnexp' duw* `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
		else {
		foreach x of varlist w* `lnpr' `lnexp' `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
	}

	//AIDS
	tempname lnpindex
	scalar `lnpindex'= `e(anot)'
	forvalue i=1/`=e(ngoods)' {
		scalar `lnpindex'= `lnpindex' + _b[alpha:alpha_`i']*`lnp`i''m
		forvalue j=1/`e(ngoods)' {
			if `j'>=`i' {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`j'_`i']*(`lnp`i''m*`lnp`j''m))
			}
			 else {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`i'_`j']*(`lnp`i''m*`lnp`j''m))
			}
		}
	}


	forvalue i=1/`=e(ngoods)' {
			tempname gsum`i'
			scalar `gsum`i''= 0
			local j=`i'
			forvalue ii=`j'/`=e(ngoods)' {
				scalar `gsum`i''= `gsum`i'' + _b[gamma: gamma_`ii'_`i']*`lnp`i''m
				}
				}

	//When quadratics
	if "`e(quadratic)'" == "quadratic" {
		tempname bofp
		scalar `bofp'= _b[beta:beta_1]*`lnp1'm
		forvalues i = 2/`=e(ngoods)' {
			 scalar `bofp'= `bofp' + _b[beta:beta_`i']*`lnp`i''m
		}
	}

	//When demographics
	if `ndemo' > 0 {
		tempname cofp mbar
		scalar `cofp'= 1 //It is OK to set 1 because below we set a multiplication
		scalar `mbar'= 1 //It is OK because I need to add a "1"

		forvalue i=1/`=e(ngoods)' {
			tempname betanz`i' cofp`i'
			scalar `betanz`i''=_b[beta:beta_`i']
			scalar `cofp`i''= 0
		}


		foreach var of varlist `e(demographics)' {
			forvalue i=1/`=e(ngoods)' {
				 scalar `cofp`i''= `cofp`i''+(`var'm*_b[eta:eta_`var'_`i'])
				 scalar `betanz`i''=`betanz`i''+(`var'm*_b[eta:eta_`var'_`i'])
									}
		scalar `mbar'= `mbar' + (_b[rho:rho_`var']*`var'm)
		}



			forvalue i=1/`=e(ngoods)' {
				scalar `cofp`i''= `cofp`i''*`lnp`i''m
				scalar `cofp'= `cofp'*`cofp`i''
			}

	}

	//FUNCTION EVALUATOR (PREDICTED SHARE)
	forvalues i = 1/`=e(ngoods)' {
		if `ndemo' == 0 {
		scalar we`i' = _b[alpha:alpha_`i']+_b[beta:beta_`i']*(`lnexp'm-`lnpindex') + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex')^2
			}
		}
		else {
		scalar we`i' = _b[alpha:alpha_`i']+`betanz`i''*(`lnexp'm-`lnpindex'-ln(`mbar')) + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar'))^2
				}
			}
		//When censor
		if "`e(censor)'" == "censor" {
		 scalar we`i' = we`i'*cdfw`i'm + _b[delta:delta_`i']*pdfw`i'm
		}
	}

	//FUNCTION EVALUATOR (ELASTICITY)

	forvalues i = 1/`=e(ngoods)' {
	forvalues j = 1/`=e(ngoods)' {
		local de=cond(`i'==`j',1,0)
		if `ndemo' == 0 {
			if `j'==`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']*(_b[alpha:alpha_`i']+`gsum`i''))))"
			}
			if `j'>`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']*(_b[alpha:alpha_`j']+`gsum`j''))))"
			}
			else {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(_b[beta:beta_`i']*(_b[alpha:alpha_`j']+`gsum`j''))))"
			}

			if "`e(quadratic)'" == "quadratic" {
				if `j'==`i' {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`i']+`gsum`i'')-(_b[beta:beta_`i']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
				if `j'>`i' {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`j']+`gsum`j'')-(_b[beta:beta_`j']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
				else {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`j']+`gsum`j'')-(_b[beta:beta_`j']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
			}
		}
		else {
			if `j'==`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i'')*(_b[alpha:alpha_`i']+`gsum`i'')))"
			}
			if `j'>`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i'')*(_b[alpha:alpha_`j']+`gsum`j'')))"
			}
			else {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(`betanz`i'')*(_b[alpha:alpha_`j']+`gsum`j'')))"
			}

			if "`e(quadratic)'" == "quadratic" {
				if `j'==`i' {
				 global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`i']+`gsum`i'')-(`betanz`i''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}
				if `j'>`i' {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`j']+`gsum`j'')-(`betanz`j''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}

				else {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`j']+`gsum`j'')-(`betanz`j''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}
			}
		}
			//When censor
		if "`e(censor)'" == "censor" {
			global ue`i'`j' "(-`de'+1/we`i'*(cdfw`i'm*((${ue`i'`j'}+`de')*w`i'm) + _b[tau:p`j'_`i']*pdfw`i'm*(w`i'm-_b[delta:delta_`i']*duw`i'm)))"
				// "${ue`i'`j'}+`de'" --> to remove lower delta from original formula
				// "w`i'm" --> to revert 1/w`i'm from the original formula
		}
	}
	}

	tempname elasu seu
	mat elasu=J(`e(ngoods)',`e(ngoods)',0)
	mat seu=J(`e(ngoods)',`e(ngoods)',0)

	forvalues i = 1/`=e(ngoods)' {
		forvalues j = 1/`=e(ngoods)' {
			qui nlcom(${ue`i'`j'})
			mat elasu[`i',`j'] = r(b)
			mat seu[`i',`j'] = r(V)
			mat seu[`i',`j'] = sqrt(seu[`i',`j'])
		}
	}
	ret mat elas_u elasu
	ret mat se_elas_u seu

end


/// COMPENSATED PRICE ELASTICITY

program DoComp, rclass

	syntax [if] [in]

	marksample touse

	if "`e(lnprices)'" != "" {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(lnprices)' {
			local lnp`i' `var'
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}
	else {
		local i 1
		tempname lnpr
		local lnpr ""
		foreach var of varlist `e(prices)' {
			tempvar lnp`i'
			gen double `lnp`i'' = ln(`var')
			local lnpr `lnpr' `lnp`i''
			local `++i'
		}
	}

	if "`e(lnexpenditure)'" != "" {
		local lnexp `e(lnexpenditure)'
	}
	else {
		tempvar exp
		gen double `exp' = ln(`e(expenditure)')
		local lnexp `exp'
	}

	if "`e(lhs)'" != "" {
		local i 1
		foreach var of varlist `e(lhs)' {
			local w`i' `var'
			local `++i'
		}
	}

	//PROCESSING VARIABLES (MARGINS WONT WORK DUE TO PARAMETER NAMES, NLCOM DOES NOT PRODDUCE SE)

	local ndemo = `e(ndemos)'
	quietly {
		if "`e(censor)'" == "censor" {
		foreach x of varlist w* `lnpr' cdfw* pdfw* `lnexp' duw* `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
		else {
		foreach x of varlist w* `lnpr' `lnexp' `e(demographics)' {
		sum `x'
		scalar `x'm=r(mean)
		}
		}
	}

	//AIDS
	tempname lnpindex
	scalar `lnpindex'= `e(anot)'
	forvalue i=1/`=e(ngoods)' {
		scalar `lnpindex'= `lnpindex' + _b[alpha:alpha_`i']*`lnp`i''m
		forvalue j=1/`e(ngoods)' {
			if `j'>=`i' {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`j'_`i']*(`lnp`i''m*`lnp`j''m))
			}
			 else {
				scalar `lnpindex'= `lnpindex' + 0.5*(_b[gamma:gamma_`i'_`j']*(`lnp`i''m*`lnp`j''m))
			}
		}
	}

	forvalue i=1/`=e(ngoods)' {
			tempname gsum`i'
			scalar `gsum`i''= 0
			local j=`i'
			forvalue ii=`j'/`=e(ngoods)' {
				scalar `gsum`i''= `gsum`i'' + _b[gamma: gamma_`ii'_`i']*`lnp`i''m
				}
				}

	//When quadratics
	if "`e(quadratic)'" == "quadratic" {
		tempname bofp
		scalar `bofp'= _b[beta:beta_1]*`lnp1'm
		forvalues i = 2/`=e(ngoods)' {
			 scalar `bofp'= `bofp' + _b[beta:beta_`i']*`lnp`i''m
		}
	}

	//When demographics
	if `ndemo' > 0 {
			tempname cofp mbar
		scalar `cofp'= 1 //It is OK to set 1 because below we set a multiplication
		scalar `mbar'= 1 //It is OK because I need to add a "1"

		forvalue i=1/`=e(ngoods)' {
			tempname betanz`i' cofp`i'
			scalar `betanz`i''=_b[beta:beta_`i']
			scalar `cofp`i''= 0
		}


		foreach var of varlist `e(demographics)' {
			forvalue i=1/`=e(ngoods)' {
				 scalar `cofp`i''= `cofp`i''+(`var'm*_b[eta:eta_`var'_`i'])
				 scalar `betanz`i''=`betanz`i''+(`var'm*_b[eta:eta_`var'_`i'])
									}
		scalar `mbar'= `mbar' + (_b[rho:rho_`var']*`var'm)
		}


			forvalue i=1/`=e(ngoods)' {
				scalar `cofp`i''= `cofp`i''*`lnp`i''m
				scalar `cofp'= `cofp'*`cofp`i''
			}
	}

	//FUNCTION EVALUATOR (PREDICTED SHARE)
	forvalues i = 1/`=e(ngoods)' {
		if `ndemo' == 0 {
		scalar we`i' = _b[alpha:alpha_`i']+_b[beta:beta_`i']*(`lnexp'm-`lnpindex') + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex')^2
			}
		}
		else {
		scalar we`i' = _b[alpha:alpha_`i']+`betanz`i''*(`lnexp'm-`lnpindex'-ln(`mbar')) + `gsum`i''
			if "`e(quadratic)'" == "quadratic" {
				 scalar we`i' = we`i' + (_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar'))^2
				}
			}
		//When censor
		if "`e(censor)'" == "censor" {
		 scalar we`i' = we`i'*cdfw`i'm + _b[delta:delta_`i']*pdfw`i'm
		}


	}


	//FUNCTION EVALUATOR (UNCOMP ELASTICITY)
	forvalues i = 1/`=e(ngoods)' {
	forvalues j = 1/`=e(ngoods)' {
		local de=cond(`i'==`j',1,0)
		if `ndemo' == 0 {
			global ie`i' "(1+_b[beta:beta_`i']/w`i'm)"
			if `j'==`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']*(_b[alpha:alpha_`i']+`gsum`i''))))"
			}
			if `j'>`i' {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']*(_b[alpha:alpha_`j']+`gsum`j''))))"
			}
			else {
			global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(_b[beta:beta_`i']*(_b[alpha:alpha_`j']+`gsum`j''))))"
			}

			if "`e(quadratic)'" == "quadratic" {
				global ie`i' "(1+1/w`i'm*(_b[beta:beta_`i']+2*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')))"
				if `j'==`i' {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`i']+`gsum`i'')-(_b[beta:beta_`i']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
				if `j'>`i' {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`j']+`gsum`j'')-(_b[beta:beta_`j']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
				else {
					global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(_b[beta:beta_`i']+(2*_b[lambda:lambda_`i']/exp(`bofp'))*(`lnexp'm-`lnpindex'))*(_b[alpha:alpha_`j']+`gsum`j'')-(_b[beta:beta_`j']*_b[lambda:lambda_`i']/exp(`bofp')*(`lnexp'm-`lnpindex')^2)))"
				}
			}
		}
		else {
			global ie`i' "(1+`betanz`i''/w`i'm)"
			if `j'==`i' {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i'')*(_b[alpha:alpha_`i']+`gsum`i'')))"
			}
			if `j'>`i' {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i'')*(_b[alpha:alpha_`j']+`gsum`j'')))"
			}
			else {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(`betanz`i'')*(_b[alpha:alpha_`j']+`gsum`j'')))"
			}

			if "`e(quadratic)'" == "quadratic" {
			global ie`i' "(1+1/w`i'm*(`betanz`i''+2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))))"
				if `j'==`i' {
				 global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`i']+`gsum`i'')-(`betanz`i''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}
				if `j'>`i' {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`j'_`i']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`j']+`gsum`j'')-(`betanz`j''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}
				else {
				global ue`i'`j' "(-`de'+1/w`i'm*(_b[gamma:gamma_`i'_`j']-(`betanz`i''+(2*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp'))*(`lnexp'm-`lnpindex'-ln(`mbar')))*(_b[alpha:alpha_`j']+`gsum`j'')-(`betanz`j''*_b[lambda:lambda_`i']/exp(`bofp')/exp(`cofp')*(`lnexp'm-`lnpindex'-ln(`mbar'))^2)))"
				}
			}
		}

		global ce`i'`j' "${ue`i'`j'}+${ie`i'}*w`j'm"

			//When censor
		if "`e(censor)'" == "censor" {
			global ue`i'`j' "(-`de'+1/we`i'*(cdfw`i'm*((${ue`i'`j'}+`de')*w`i'm) + _b[tau:p`j'_`i']*pdfw`i'm*(w`i'm-_b[delta:delta_`i']*duw`i'm)))"
				// "${ue`i'`j'}+`de'" --> to remove lower delta from original formula
				// "w`i'm" --> to revert 1/w`i'm from the original formula

			global ie`i' "(1+1/we`i'*((cdfw`i'm*((${ie`i'}-1)*w`i'm))+_b[tau:M_`i']*pdfw`i'm*(w`i'm-_b[delta:delta_`i']*duw`i'm)))"
				// "${ie`i'}-1" --> to remove the "+1" from original formula
				// "*w`i'm" --> to revert 1/w`i'm from the original formula

			global ce`i'`j' "${ue`i'`j'}+${ie`i'}*we`j'"

		}
	}
	}

	//FUNCTION EVALUATOR (COMP ELASTICITY)

	tempname elasc sec
	mat elasc=J(`e(ngoods)',`e(ngoods)',0)
	mat sec=J(`e(ngoods)',`e(ngoods)',0)
	forvalues i = 1/`=e(ngoods)' {
		forvalues j = 1/`=e(ngoods)' {
			qui nlcom(${ce`i'`j'})
			mat elasc[`i',`j'] = r(b)
			mat sec[`i',`j'] = r(V)
			mat sec[`i',`j'] = sqrt(sec[`i',`j'])
		}
	}
	ret mat elas_c elasc
	ret mat se_elas_c sec

end