program define qscope, rclass
*!	program qscope	version 1.31 31Mar2019 by Mead Over, Svetlana Pivovarova and Shunjie Tu

*	Calculates the quasi economy of scope
*	A post-estimation command that requires that estimation results be in memory

	version 11.1
	syntax [if] [in] , OUTputs(varlist) EPSilon(numlist) [DEPvar(varname)  ///
		chain(string) chainlhs(string) chainrhs(string) chainwrt(string) LAMbda(real 0)  ///
		mean median REFerence(string) CONTstat(string) FVstat(string)  ///
		spec1(string) spec2(string) spec3(string) spec4(string) spec5(string) joint(string)  ///
		allset1(varname fv) allsetto1(real 1) allset2(varname fv) allsetto2(real 1)  /// Added 12Oct2018
		allset3(varname fv) allsetto3(real 1) allset4(varname fv) allsetto4(real 1)  ///
		allset5(varname fv) allsetto5(real 1) allset6(varname fv) allsetto6(real 1)  /// Added 31Mar2019
		allset7(varname fv) allsetto7(real 1) allset8(varname fv) allsetto8(real 1)  ///
		verbose debug df(real 9999) estimcheck *]  //  Added 13Nov2018

	// Default is NOT to check for estimability (Added in version 1.24)
	if "`estimcheck'"=="" {
		local noestimcheck noestimcheck
	}
	else {
		local noestimcheck
	}

	local noutputs: word count `outputs'
	if `noutputs' <2 {
		di as err "Economies of scope are only defined when there are at least 2 outputs."  ///
			_n "Here the list of outputs is: " as res "`outputs'"
		exit 198
	}

	if "`verbose'"==""  {
		local shh quietly
	}

	marksample touse
		local andtouse & `touse'

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

	if "`chain'"~="" & ("`chainlhs'"~="" | "`chainrhs'"~="")  {
		di as err "Specifying the option -chain- precludes specifying"  ///
			_n	"either of -chainlhs- or -chainwrt- and has the same effect"  ///
			_n	"as specifying them both with the same transformation choice."
		exit 198
	}
	if "`chain'"~="" & ("`chainlhs'"=="" & "`chainrhs'"=="")  {
		local chainlhs `chain'
		local chainrhs `chain'
	}

	if "`chainrhs'" == "ln"  {
		foreach e of numlist `epsilon'  {
			if `e' <= 0  {
				di as txt "For -ln- transformation on right-hand-side output variables,"   ///
					_n "non-positive values of epsilon produce missing values of -qscope-"
			}
		}
	}
	if ("`chain'" == "boxcox" | "`chainlhs'" == "boxcox" |"`chainrhs'" == "boxcox") & `lambda' == 0  {
		di as err "If any of chain(boxcox), chainrhs(boxcox) or chainlhs(boxcox) is specified,"  ///
			_n "the user must also specify a non-zero lambda."
		di as err "A non-zero lambda must be specified with either chainrhs(boxcox) or chainlhs(boxcox)."
		exit 198
	}

	local _ "_"
	if "`chainrhs'" == "boxcox"  {
		local prfxrhs "bc"
	}
	else if "`chainrhs'" ~= ""  {
		local prfxrhs "`chainrhs'"
	}
	else {
		local prfxrhs
		local _
	}

	if "`contstat'"=="p50" {
		local contstat median
	}
	else if "`contstat'"=="" {
		local contstat mean
	}
	if "`fvstat'"=="" {
		local fvstat mean
	}

	local wrd1 : word 1 of `reference'
	capture confirm number `wrd1'
	if _rc~=0 {
		local reftype stat
		if "`reference'"~="" & "`mean'`median'"~= "" {
			di as err "The reference firm can be specified by either the option -reference(stat)-"  ///
				_n as err "or by specifying one of the two options -mean- or -median-."  ///
				_n as err "Here the user has specified both -reference(" as txt "`reference'" as err ") and "  ///
				_n the option(s) : " as txt "`mean'" "`median'"
			exit 198
		}
		if "`mean'`median'`reference'" == "" | "`mean'" == "mean"  {
			local statsum mean
*			local statmarg mean
		}
		else if "`median'" == "median" | "`reference'"=="p50" | "`reference'"=="median"  {
			local statsum p50
*			local statmarg median
		}
		else {
			local found no
			foreach stat in mean min max p1 p5 p10 p25 p75 p90 p95 p99 {
				if "`reference'"=="`stat'" {
					local statsum `stat'
*					local statmarg `stat'
					local found yes
				}
			}
			if "`found'"=="no" {
				di as err "Error in the option -reference(stat)- or -reference(numlist)-."///
				_n as err "For the reference(stat)- version, the selected statistic, -stat-," ///
				_n as err "must be one of: mean median min max p1 p5 p10 p25 p50 p75 p90 p95 p99."  ///
				_n as err "For the -reference(numlist)- option, the -numlist- must be"  ///
				_n as err "a list of K numerical values, where K is the number of "  ///
				_n as err "output variables listed in the -outputs()- option."  ///
				_n as err "Here the reference option contains: " as res "`reference'"
				exit 198
			}
		}
		display  ///
			_n as txt "For the designated output variables " as res "`outputs'" as txt ","  ///
			_n as txt "joint and specialized costs are computed for the reference output vector"  ///
			_n as txt "composed of the " as res "`statsum'" as txt " values of those variables." ///
			_n as txt "The other continuous r.h.s. variables are set to their " as res "`contstat'" as txt " values"  ///
			_n as txt "and all discrete r.h.s. factor variables set to their " as res "`fvstat'" as txt " values."
	}
	else {
		local reftype nmbrlist
		local nrefwrds : word count `reference'
		if `nrefwrds'~=`noutputs' {
			di as err "When the reference firm is selected by the -reference(numlist)- option,"  ////
			_n as err "the number of values in the -numlist- must equal the number of output variables"  ///
			_n as err "specified in the outputs option."
			exit 198
		}
		local j = 1
		foreach output of local outputs {
			local yval : word `j' of `reference'
			local ref`output' = `yval'
			local j = `j' + 1
		}
		display  ///
			_n as txt "For the designated output variables " as res "`outputs'" as txt ","  ///
			_n as txt "joint and specialized costs are computed for the"  ///
			_n as txt "reference output vector " as res "`reference'"  ///
			_n as txt "The other continuous r.h.s. variables are set to their " as res "`contstat'" as txt " values"  ///
			_n as txt "and all discrete r.h.s. factor variables set to their " as res "`fvstat'" as txt " values."
	}

	// Assemble the at options to fix either factor or continuous variables
	// at other than the default values for all specialized & joint estimates
	// This loop is new to version 1.3
	foreach i of numlist 1/8 {
		if "`allset`i''"~=""  {
			if strpos("`allset`i''",".")>0 {  //  Remove factor variable prefix
				local allset`i' = substr("`allset`i''",strpos("`allset`i''",".")+1,.)
			}
			local atallset `atallset' `allset`i'' = `allsetto`i''
			di as txt "The variable " as res "`allset`i''" as txt " is set to " as res "`allsetto`i''"  ///
				as txt " in all expressions."
		}
	}


	// assemble the -margins- command's -at- syntax for joint production
	local i = 2
	local j = 1
	foreach output of local outputs  {
		if "`reftype'"=="stat" {
			_summraw `prfxrhs'`_'`output' if `touse', stat(`statsum') chainrhs(`chainrhs') lambda(`lambda')  prfxrhs(`prfxrhs') `debug'
			local yrefvalue = `r(rawstat)'
		}
		else if "`reftype'"=="nmbrlist" {
			local yrefvalue = `ref`output''
		}
		else {
			di as err "Unknown type of reference value."
			exit 198
		}

		//	Works whether -chainrhs- specifies a transformation or is blank
		if "`chainrhs'" == "boxcox"  {
			local at = (`yrefvalue'^`lambda' - 1)/`lambda'
		}
		else  {
			local at = `chainrhs'(`yrefvalue')
		}

		local joint_at `joint_at' `prfxrhs'`_'`output'=`at'

		//	syntax to retrieve the total cost of each specialized firm
		local spec_costs `spec_costs'  `comma' _b[`i'._at]
		local msng_costs `msng_costs'  `comma'  .
		local i = `i' + 1
		local j = `j' + 1
		local comma ,
		tempname som_`output'  // Specialized output matrix for this output
	}
	// -at- options for joint costs
	if "`joint'"~=""{
		di _n as txt "Total cost for joint production is estimated with the following additional -at- specification(s): "  ///
		   _n as res "`joint'"
	}
*	local options_joint "at((`contstat') _continuous (`fvstat') _factor `joint_at' `joint')"  /* Version 1.23  */
	local options_joint "at((`contstat') _continuous (`fvstat') _factor `joint_at' `joint' `atallset')"  /* Version 1.3  */

	// -at- options for the specialized costs
	// and the part of qscope formula related to specialized costs
	local firste yes
	tempname qs mat_e mat_qsc mat_se mat_z mat_p mat_tc mat_cs spec_tc rtn_spec_tc mat_som som
	est sto depvar
	foreach e of numlist `epsilon'  {
		capture matrix drop `som'
		qui est res depvar
		local eps = `e'
		local qscope_missing = `eps' <= 0 & ("`chainrhs'" == "ln" |  "`chainrhs'" == "log10" | ("`chainrhs'" == "boxcox" & `lambda' == 0))

		if `qscope_missing'  {
			local sum_spclzd = .
			local varofcost = .
			local tc = .
			local sum_spclzd_p = .
			local qscope = .
			local se = .
			local z_value = .
			local p = .
			mat `spec_tc' = ( `msng_costs' )
			local fill_som = `noutputs' * "`msng_costs' , "
			local fill_som = substr("`fill_som'",1,strlen("`fill_som'")-2)

			mat `som'  = ( `eps' , `fill_som' )
		}
		else  {
			if `eps' < 0 | `eps' > 1/`noutputs'  {
				di as err "Epsilon is required to be between 0 and 1/K, where K = the number of outputs."
				di as err "Here epsilon is " as res "`eps'" as err " and the number of outputs is " as res "`noutputs'."
				exit 198
			}

			local i = 1
			local margoptions `options_joint'
*			local sum_spclzd 0   //  Replace this with the following two initializations
			local sum_spclzd
			local plus
			foreach output1 of local outputs  {
				// assemble the -at- syntax for the margins commands
				// for each of the specialized outputs
				if "`reftype'"=="stat" {
					_summraw `prfxrhs'`_'`output1' if `touse', stat(`statsum') chainrhs(`chainrhs') lambda(`lambda')  prfxrhs(`prfxrhs') `debug'
						local rawoutput1 `r(rawvar)'
						local yrefvalue = `r(rawstat)'
				}
				else if "`reftype'"=="nmbrlist" {
					local rawoutput1 `output1'
					local yrefvalue = `ref`output1''
				}
				else {
					di as err "Unknown type of reference value."
					exit 198
				}

				//	Works whether -chainrhs- specifies a transformation or is blank
				if "`chainrhs'"~="boxcox"  {
					local at1 = `chainrhs'((1 - (`noutputs' - 1)*`eps')*`yrefvalue')
				}
				else  {  //  for boxcox
					local at1 =    ///
						(   ///
							(( (1 - (`noutputs' - 1)*`eps') * `yrefvalue' )^`lambda'  -  1)   ///
										/`lambda'   ///
						)
				}
				local spec_at `prfxrhs'`_'`output1'=`at1'

				local spec_out = (1 - (`noutputs' - 1)*`eps')*`yrefvalue'
				matrix `som_`output1'' = ( `spec_out' )
				local  som_colnames `rawoutput1'
				local  som_eqnames  `rawoutput1'

				//   syntax for the other outputs
				local restofoutputs : list outputs - output1
				foreach output2 of local restofoutputs  {
					if "`reftype'"=="stat" {
						_summraw `prfxrhs'`_'`output2' if `touse', stat(`statsum') chainrhs(`chainrhs') lambda(`lambda')  prfxrhs(`prfxrhs') `debug'
							local rawoutput2 `r(rawvar)'
							local yrefvalue = `r(rawstat)'
					}
					else if "`reftype'"=="nmbrlist" {
						local rawoutput2 `output2'
						local yrefvalue = `ref`output2''
					}
					else {
						di as err "Unknown type of reference value."
						exit 198
					}

					//	Works whether -chainrhs- specifies a transformation or is blank
					if "`chainrhs'"~="boxcox"  {
						local at2 = `chainrhs'(`eps'*`yrefvalue')
					}
					else  {  //  for boxcox
						local at2 =    ///
							(   ///
								(( `eps' * `yrefvalue' )^`lambda'  -  1)   ///
											/`lambda'   ///
							)
					}
					local addon `prfxrhs'`_'`output2'=`at2'
					local spec_at `spec_at' `addon'
					local non_spec_out = `eps'*`yrefvalue'
					matrix `som_`output1'' = (nullmat(`som_`output1'') , `non_spec_out' )
					local  som_colnames `som_colnames' `rawoutput2'
					local  som_eqnames  `som_eqnames' `rawoutput1'
				}
				if "`firste'"=="yes" & "`spec`i''"~=""{
					di _n as txt "The specialized cost " as res "`i'" as txt " is estimated with the additional -at- specification(s): "  ///
						_n as res "`spec`i''"
				}
				local margoptions `margoptions' at((`contstat') _continuous (`fvstat') _factor `spec_at' `spec`i'' `atallset')
				local i = `i'+1
*				local sum_spclzd `sum_spclzd' + _b[`i'._at]  // Replaced with following line to avoid misleading "0 +" at beginning of nlcom expression
				local sum_spclzd `sum_spclzd' `plus' _b[`i'._at]
				local plus +

				matrix colnames `som_`output1'' = `som_colnames'
				matrix coleq  `som_`output1'' = `som_eqnames'

				mat `som'  = (nullmat(`som') , `som_`output1'')
			}  	//  Next output
			local firste no

			// Performing margins estimation:
			`shh' di _n as txt "The following output is from Stata's {help margins} command for epsilon = " as res "`eps'" as txt ":"
			if "`chainlhs'" == ""  {
				`shh' margins if `touse', predict(equation(`depvar')) ///
					 `margoptions' post `options' `noestimcheck'
			}
			if "`chainlhs'" == "sqrt"  {
				`shh' margins if `touse', expression(predict(equation(`depvar'))^2) ///
					 `margoptions' post `options' `noestimcheck'
			}
			if "`chainlhs'" == "ln" | ("`chainlhs'" == "boxcox" & `lambda' == 0)  {
				`shh' margins if `touse', expression(exp(predict(equation(`depvar')))) ///
					 `margoptions' post `options' `noestimcheck'
			}
			if "`chainlhs'" == "log10"  {
				`shh' margins if `touse', expression(10^(predict(equation(`depvar')))) ///
					 `margoptions' post `options' `noestimcheck'
			}
			if "`chainlhs'" == "boxcox" & `lambda' ~= 0  {
				`shh' margins if `touse',   ///
					expression(   ///
						(`lambda'*predict(equation(`depvar')) +1)^(1/`lambda')  ///
					) ///
					 `margoptions' post `options' `noestimcheck'
			}

			//  For this epsilon, total cost of the firm specialized in each output
			local  som_colnames : colnames `som'
			local  som_eqnames  : coleq  `som'
			mat `som' = ( `eps' , `som' )
				matrix colnames `som' = epsilon `som_colnames'
				matrix coleq  `som'  = epsilon `som_eqnames'
*			mat list `som' , title(som)
			mat `spec_tc' = ( `spec_costs' )
			// qscope = (sum of specialized costs - total cost)/(total cost)
			// or qscope = [(sum of specialized costs)/(total cost)] - 1
			local tc = _b[1._at]
			local sum_spclzd_p = `sum_spclzd'
			`shh' nlcom(qscope: (`sum_spclzd')/_b[1._at]-1), df(`df') post
			local qscope = _b[qscope]
			local se = _se[qscope]
			local z_value = `qscope'/`se'
			local p = 2*(1-normal(abs(`qscope'/`se')))  //  Assuming a z-dist

		}

		mat `mat_som'  = (nullmat(`mat_som') \ `som')
		mat `mat_e'  = (nullmat(`mat_e') \ `e')
		mat `mat_qsc'  = (nullmat(`mat_qsc') \ `qscope')
		mat `mat_se' = (nullmat(`mat_se') \ `se')
		mat `mat_z' = (nullmat(`mat_z') \ `z_value')
		mat `mat_p' = (nullmat(`mat_p') \ `p')
		mat `mat_tc' = (nullmat(`mat_tc') \ `tc')
		mat `mat_cs' = (nullmat(`mat_cs') \ `sum_spclzd_p')
		mat `rtn_spec_tc' = (nullmat(`rtn_spec_tc') \ `spec_tc')

	}

	mat `qs' = `mat_e' , `mat_cs' , `mat_tc' , `mat_qsc' , `mat_se' , `mat_z' , `mat_p'
	matrix colnames `qs' = epsilon sum_spclzd  tc qscope se z p
	matrix colnames `rtn_spec_tc' = `outputs'
	matrix define `rtn_spec_tc' = `qs'[1...,1 ], `rtn_spec_tc', `qs'[1...,2 ]
	local somfullnames : colfullnames `som'
	matrix colnames `mat_som' = `somfullnames'

	`shh' {
*	Display table of output combinations for specialized and non-specialized firms for each value of -epsilon-
		di _n as txt "The following table gives the amount of each output"  ///
			_n "produced by each specialized and non-specialized firm."  ///
			_n "If K outputs are specified, there are K sets of K columns each."  ///
			_n "After the first, each column title consists of an identifer "  ///
			_n "for the specialization of the firm (on top) and the name"  ///
			_n "of the output appearing in the column (on bottom). "  ///
			_n "This matrix can be retrieved as " as res "r(spec_outputs)"
		mat list `mat_som', title("Output combinations for specialized & non-specialized firms")

*	Display table the distribution of the cost of specialized production across the specialized firms for each value of -epsilon-
		di _n as txt "The following table gives the breakdown of specialized cost"  ///
			_n "across the K specialized firms and, in the last column,"  ///
			_n "the total cost of producing all the outputs summed"  ///
			_n "across all K specialized firms."  ///
			_n "This matrix can be retrieved as " as res "r(spec_tc)"
		mat list `rtn_spec_tc', title("Cost of production in each specialized firm and total specialized cost")
	}

*	Display table of -qscope-results for each value of -epsilon-
	local firste yes
	local r = 1
	foreach e of numlist `epsilon' {
*		local z_value = `qs'[`r',2]/`qs'[`r',3]
		if "`firste'"=="yes"  {
			di as txt  ///
		_n _col(1) "Values of" 	_col(12) "Specialzd" _col(24) "  Joint  " _col(36) "      " _col(45) "Standard"  _col(57) "Ratio of   " _col(70) "p-value of"   ///
		_n _col(1) "epsilon"   	_col(12) "  Cost   " _col(24) "  Cost   " _col(36) "qscope" _col(45) "  Error "  _col(57) "qscope/s.e." _col(70) "Ho: q = 0 "  ///
		_n _col(1) "=========" 	_col(12) "=========" _col(24) "=========" _col(36) "======" _col(45) "========"  _col(57) "===========" _col(70) "=========="  ///
				as res _n ///
		_col(1) `qs'[`r',1]	_col(12) `qs'[`r',2] _col(24) `qs'[`r',3] _col(36) %6.2f `qs'[`r',4] _col(45) %8.3f `qs'[`r',5]	_col(57) %9.3f `qs'[`r',6] _col(70) %9.3f `qs'[`r',7]
			local firste no
		}
		else  {
			 di as res  ///
		_col(1) `qs'[`r',1]	_col(12) `qs'[`r',2] _col(24) `qs'[`r',3] _col(36) %6.2f `qs'[`r',4] _col(45) %8.3f `qs'[`r',5]	_col(57) %9.3f `qs'[`r',6] _col(70) %9.3f `qs'[`r',7]
		}
		local r = `r' + 1
	}

	return matrix spec_outputs = `mat_som'
	return matrix spec_tc = `rtn_spec_tc'
	return matrix qscope = `qs'
	return local  epsilon `epsilon'
	return local  statsum  `statsum'
	return local  statmarg `statmarg'
	return scalar noutputs = `noutputs'
	return scalar df = `df'
*	return scalar p = 2*ttail(`df',abs(`qscope'/`se'))  //  Assuming a t-dist with `df' degrees of freedom
	return scalar p = 2*(1-normal(abs(`qscope'/`se')))  //  Assuming a z-dist
	return scalar last_eps = real(word("`epsilon'",wordcount("`epsilon'")))
	return scalar se = `se'
	return scalar lambda = `lambda'
	return scalar sum_spclzd = `sum_spclzd_p'
	return scalar totalcost = `tc'
	return scalar last_qscope = `qscope'

	qui est res depvar
end //  End of program qscope

program define _summraw, rclass
*	Parse a possibly transformed variable to give the name of the raw variable.
*		Possible transformations are -ln-, -log10-, -sqrt-, -bc-.
*	Compute the chosen statistic for the raw variable.
*	Apply the transformation to the computed statistic.
	syntax varname [if] [in], stat(string) [CHAINrhs(string) lambda(real 0) debug prfxrhs(string)]
	tempname rawstat chainstat nobs

	if "`debug'"==""  {
		local shh quietly
	}

	local var `varlist'
	//	If -chainrhs- is present, variable -varname- is assumed transformed
	//  coded by the prefix -prfxrhs- with the transformation
	if "`chainrhs'"~="" {
		if "`chainrhs'"~="`prfxrhs'" & ~("`chainrhs'"=="boxcox" & "`prfxrhs'"=="bc") {
			di as err "If specified, the option -chain- or the option -chainrhs- must correspond"  ///
				_n as err "to one of the allowed transformations. These are: " as res "ln, log10, sqrt, bc"
				exit 198
		}
		if "`chainrhs'"=="ln" {
			local rawvar = substr("`var'", strpos("`var'","ln_")+3,.)
		}
		else if "`chainrhs'"=="log10" {
			local rawvar = substr("`var'", strpos("`var'","log10_")+6,.)
		}
		else if "`chainrhs'"=="sqrt" {
			local rawvar = substr("`var'", strpos("`var'","sqrt_")+5,.)
		}
		else if "`chainrhs'"=="boxcox"|"`chainrhs'"=="bc" {
			local rawvar = substr("`var'", strpos("`var'","bc_")+3,.)
		}
		else {
			di as err "The chain or chainrhs option must be one of the following:"  ///
				_n "    ln, log10, sqrt, boxcox."
			exit 198
		}
	}
	else {
		local rawvar `var'
	}

	//	Compute chosen statistic
	_legit `stat'
		local stat `r(legitstat)'
	`shh' summarize `rawvar' `if' `in', det
			scalar `rawstat' = r(`stat')
			scalar `nobs' = r(N)

	//	Transform the computed statistic
	if "`chainrhs'"~="" {
		if "`chainrhs'"=="ln" | "`chainrhs'"=="log10" | "`chainrhs'"=="sqrt" {
				scalar `chainstat' = `chainrhs'(`rawstat')
		}
		else if "`chainrhs'"=="boxcox"|"`chainrhs'"=="bc" {
			scalar `chainstat' =    ///
				(   ///
					(( `rawstat' )^`lambda'  -  1)   ///
											/`lambda'   ///
				)
		}
		else {
			di as err "The chain or chainrhs option must be one of the following:"  ///
				_n "    ln, log10, sqrt, boxcox."
			exit 198
		}
	}
	else {
		scalar `chainstat' = `rawstat'
	}

	return local argvar `var'
	return scalar N = `nobs'
	return scalar chainstat = `chainstat'
	return scalar rawstat = `rawstat'
	return local stat `stat'
	return local rawvar `rawvar'
	return local chain `chainrhs'
end // End of program _summraw

program define _legit, rclass
	args stat
	if "`stat'"=="median" {
		local stat p50
	}
	local stat_is bad
	foreach legitstat in p1 p5 p10 p25 p50 p75 p90 p95 p99 mean min max {
		if "`stat'"== "`legitstat'" {
			local stat_is good
		}
	}
	if "`stat_is'"=="bad" {
		di as err "The statistic -" as txt "`stat'" as err "- cannot be computed by {help summarize}"
		exit 198
	}
	return local legitstat `stat'
end // End of program _legit


*	Version History: (Versions 1.0 and 1.1 by Shunjie Tu)
*	1.0	24Aug2016	By Shunjie Tu
*	1.1	26Aug2016	add -if-, -in- and options -mean- and -median-,
*					but p50 (median) statistic not allowed with factor variables.
*	1.2	8Sep2016    replace "_all" with "_continuous (mean) _factor"
*					(see the undocumented help for _ms_at_parse now referenced in the help file)
*					Also some cosmetic changes:
*						Change returned values to scalars
*						Widened the columns of the output
*						Add specialised & joint cost to returned value and output
*	1.3	14Sep2016   Limit epsilon to < 1/K
*	1.4	16Sep2016   Add a p-value and a degrees of freedom option, -df()-
*	1.5	16Sep2016   Change the unwind option for cost from chainrhs to chainlhs
*	1.6 20Sep2016	Incorporate Shunjie's version of boxcox transformations, but retain my unwinding
*					Add return locals for statsum & statmarg
*					Should unwind lhs transformation in the -margins- computation
*	1.7 24Sep2016	Allow list of epsilon values (not just one value)
*					Return results for all values of epsilon in a matrix -QSCOPE-
*					Use margins, expression- instead of Delta Method to compute variances of predicted values
*					Accommodate zero among the epsilon values
*	1.8 29Sep2016	Form the matrix QSCOPE as a series of column vectors
*					Add -saving-, -model- and -append- options
*	1.9 9Dec2016	Add -spec1- through -spec5-, and -joint- options
*					Reduce the number of calls of -margins- to one
*					Compute standard error for qscope in -margins- postesimation instead of delta method
*	1.10 4Jan2017	Suppress the saving option and the associated append and model options
*					Adapt to using "bc" as the prefix for the -boxcox- trans on rhs
*	1.11 16Jan2017	When epsilon is zero, set qscope to missing when the rhs transformation
*					is ln or log10 or boxcox(0), regardless of the lhs transformation.
*					Add subroutines _summraw and _legit.
*	1.12 23Jan2017	Move down the printing of the table of -qscope- results to
*					below the computations in order to collect the detailed results
*					on specialized production for transformed variables.
*					Present the results for the breakdown of output and cost
*					across the specialized firms in the matrices r(spec_tc) & r(spec_outputs)
*	1.13 11Feb2017	Add the p-value of the test that qscope = 0 to the standard output
*					Fix bug in unwinding boxcox transformation on lhs
*					Add the -*- to the syntax command for options to be apssed to -margins-
*	1.20 15Feb2017	Add the -reference(stat)- and the -reference(numlist)- options
*					to enable -qscope- to specify that the level of the outputs
*					of a reference firm be other than the sample mean or median
*					of the output vector.  Eligible stats are those in the
*					-summarize, detail- return space.
*					Replace "(`statmarg') _continuous (mean) _factor" with
*					"(`contstat') _continuous (`fvstat') _factor"
*	1.21 23Feb2017	Fix bug in -log10- option in line 362
*	1.22 28Mar2018	Add error trap to catch the invalid "chainwrt" option
*	1.23 5May2018	Regardless of -verbose-, always display options -speck- & -joint-
*	1.24 13Nov2018	Often there are no observations which specialize in only one output.
*					Since this creates "empty cells", we make -noestimcheck- the default.
*	1.30 13Nov2018	Add the -allset`i'- and the -allsetto`i'- options: 1, ..., 4
*	1.31 31Mar2019	Add the -allset`i'- and the -allsetto`i'- options: 5, ..., 8
*					Not documented yet in -qscope.sthlp-