prog define wtsreg, eclass
*!	Version 0.3 by Mead Over 31May2024

*	Program -wtsreg-
*	This program estimates a regression of a dependent variable on a set
*	right-hand-side or explanatory variables under the two constraints:
*		1. Each coefficient is within the unit interval
*		2. The sum of the coefficients is equal to 1.0

*	A set of coefficients obeying these properties can be characterized
*	as a set of weights.  In the absent of a constant term,
*	post-multiplication of an NxK theKx1 tranpose of a weight vctor would
*	compute a weighted average of the K columns (i.e. variables) in X.

*	If the -mlogit- model is selected for imposing constraints on the coefficients
*	(the default), this program follows the suggestion in the last example
*	discussed in Isabel Canette's useful blog post.

*	Source: Isabel Canette, "How do I fit a linear regression with interval
*	(inequality) constraints in Stata?"  Accessed 28May2024:
*	"https://www.stata.com/support/faqs/statistics/linear-regression-with-interval-constraints/"

/*
Syntax:
	wtsreg depvar indepvars [aw iw] [if] [in] [, -wtsreg options- -nl options-


wtsreg Options:
	model(string)
	  mlogit            Use -nl- to estimate a multinomial logit model: Each weight in unit interval, wts sum to one

The following alternative models are not yet implemented:
	  regress           Use -nl- to estimate an OLS model: Unconstrained linear regression.
                          For comparison only.  The most efficient way to estimate by OLS is to use -regress-.
	  invlogit			Use -nl- to estimate an inverse logit model: Each weight in unit interval, sum of wts unconstrained
	  quad              Use -nl- to estimate a quadratic model: Each weight positive, sum of wts unconstrained

	delta(real)         Criterion for -nl- convergence. Following Canette (op. cit.), default is 1e-7
	initial(initial_values)
	                    See {help nl##initial} for -nl-'s syntax for specifying starting values.
	nlstart(string)     Optionally how -wtsreg- sets the starting values of the _b coefficients
	                      for the -nl- iterations.  An alternative to
						  and incompatible with -nl-'s option -initial()-.

	  CNSreg            Get the starting values from -cnsreg-. The default.
	  REGress           Get the starting values from -regress-.
	  naive             Set the starting values of all the coefficients to the inverse of the number of donor units.

	toglobals			Set global macros for all of the returned local macros
	verbose             Show intermediate results of -nl- and -nlcom-
	displayopts(string) Options passed to estimates display
	nlstore(estname)    Store the estimation results from the -nl- estimation command with name -estname-
	                        Without this option, only results of the -nlcom- unwinding
							of the coefficient estmates are stored
	iterate(#)          Maximm number of iterations of the -nl- and -nlcom- routines.
	                        Passed through to both -nl- and -nlcom-

	*					Other options passed to -nl-


*/
	version 15


if replay() {
  if "`e(cmd)'"!="wtsreg" {
    error 301
  }
  syntax [, COEFlegend verbose * ]
	if "`verbose'"~="verbose" {
		local shhh quietly
	}
  local displayopts `options'
}
else {

	syntax varlist [aw iw] [if] [in] [, model(string) /* INITconst(string) */  ///
		delta(real 1e-7) toglobals verbose displayopts(string) iterate(passthru)  ///
		nlstore(name) initial(passthru) nlstart(passthru) *]

	if "`verbose'"=="verbose" {
		local shhh
		local Xeq Xeq
		local log log
	}
	else {
		local catch cap
		local shhh quietly
		local log nolog
	}

	if "`model'"=="" {
		local model mlogit
	}
	else {
		opts_exclusive "`ols' `invlogit' `mlogit'" model
	}

	gettoken depvar wtsrhs : varlist
		local Nwtsrhs = wordcount("`wtsrhs'")

	`shhh' di _n as txt "Dep. var. is: " as res "`depvar'"
	`shhh' di as txt "RHS variables are: " as res "`wtsrhs'"
	`shhh' di as txt "Number of RHS variables are: " as res "`Nwtsrhs'"
	`shhh' di as txt "The -varlist- macro is unchanged: " as res "`varlist'"
	`shhh' di as txt "The estimation model is: " as res "`model'" _n
	`shhh' di as txt "The user-specified -nlstart()- option is: " as res "`nltart'" _n
	`shhh' di as txt "The user-specified -initial()- option is: " as res "`initial'" _n


	//	Get coefficient starting values for -nl- optimization
	//	Alternative to -nl-'s option -initial()-
	opts_exclusive "`initial' `nlstart'"
	if "`initial'"==""  {

		get_eb_start `wtsrhs' `if' `in' , depvar(`depvar') `nlstart' `verbose'

		if "`r(eb_start)'"~="" {
			tempname eb_start
			mat def `eb_start'= r(eb_start)
			local initial initial(`eb_start')
		}
	}

	if "`model'"=="mlogit"  {

		//  Build the denominator
		local k = 1
		foreach v of varlist `wtsrhs' {
			if `k'==1 {
				local denom 1
			}
			else {
				local denom `denom' + exp({t`k'})
			}
			local k = `k' + 1
		}

		//  Define the -Nwtsrhs- coefficients and construct the non-linear expression for the estimand
		local k = 1
		foreach v of varlist `wtsrhs' {
			if `k'==1 {
				local ma1 1/(`denom')
			}
			else {
				local ma`k' exp({t`k'})/(`denom')
			}
			local rhsvar = word("`wtsrhs'",`k')
			local nlsyntx `nlsyntx' `plus'  `ma`k''*`rhsvar'
			local synline`k' `plus'  `ma`k''*`rhsvar'

			local plus +
			local k = `k' + 1
		}
		local nlsyntx `nlsyntx' + {_cons}

		`shhh' di as txt "Syntax for the RHS of the -nl- comand is:"
		if "`verbose'"=="verbose" {
			foreach i of numlist 1/`Nwtsrhs' {
				di as res "    `synline`i''"
			}
			di as res "    + {_cons}"
		}

		`shhh' di _n as txt "Begin -nl- iterations:" _n

		******************Estimate the model using-nl-********************
		`catch' `Xeq' nl (`depvar' = `nlsyntx') [`weight'`exp'] `if' `in',   ///
			delta(`delta') variables(`wtsrhs') `log' `iterate' `initial' `options'
		******************************************************************

		`shhh' di _n as txt "End of -nl- iterations:"

		//	N.B.: -variables()- enables subsequent -margins, dydx(rhs)-
		if "`e(wexp)'"~="" {
			tempvar wttmp
			qui gen `wttmp' `e(wexp)'
			sum `wttmp', meanonly
			ereturn scalar sumofwts = `r(sum)'
		}
		if "`nlstore'"~="" {
			est store `nlstore'
			di as txt "Results of the -nl- estimation command stored as: " as res "`nlstore'"
		}

		//  Grab the -nl- returned macros
		foreach lcl in cmdline cmd predict marginsnotok marginsok v  ///
		  ce rhs type sexp params depvar properties  {
			local `lcl' `e(`lcl')'
		}

		foreach sclr in df_m df_r mss mms msr rmse ll r2 r2_a dev  ///
		  converged N rss df_t tss cj gm_2 delta k k_eq_model  ///
		  rank log_t ic sumofwts {
			local `sclr' = e(`sclr')
		}

		foreach mat in init b V {
			mat def nl_`mat' = e(`mat')
		}

		//  Build the denominator
		local k = 1
		foreach v of varlist `wtsrhs' {
			if `k'==1 {
				local denom 1
			}
			else {
				local denom `denom' + exp(_b[t`k':_cons])
			}
			local k = `k' + 1
		}

		//  Define the -Nwtsrhs- coefficients and construct the -nlcom- command
		local k = 1
		foreach v of varlist `wtsrhs' {
			if `k'==1 {
				local na1 1/(`denom')
			}
			else {
				local na`k' exp(_b[t`k':_cons])/(`denom')
			}
			local rhsvar = word("`wtsrhs'",`k')
			local nlcomsyntx `nlcomsyntx' (`rhsvar': `na`k'')
			local k = `k' + 1
		}
		local nlcomsyntx `nlcomsyntx' (_cons:_b[_cons:_cons])

		`shhh' di _n as txt "Begin -nlcom- iterations:" _n

		******************************************************************
		//  Constant term (_cons: `na0') added by MO
		`catch' `Xeq' nlcom `nlcomsyntx',  `iterate' post
		******************************************************************

		`shhh' di _n as txt "End of -nlcom- iterations:"

		//	Post the b, V matrices from -nlcom-, replacing those from -nl-
		mat def eb = e(b)
		mat def eV = e(V)
		ereturn repost  b=eb V=eV

		//  Post r(table) from -nlcom-
		mat def table = r(table)
		ereturn matrix table = table

		//  Post the init, b and V matrices from -nl-
		foreach mat in init b V {
			ereturn matrix nl_`mat' = nl_`mat'
		}

		foreach sclr in df_m df_r mss mms msr rmse ll r2 r2_a dev  ///
		  converged N rss df_t tss cj gm_2 delta k k_eq_model  ///
		  rank log_t ic sumofwts {
			ereturn scalar `sclr' = ``sclr''
		}

		foreach lcl in cmdline cmd marginsnotok marginsok v  ///
		  ce rhs type sexp params depvar properties  {
			ereturn local `lcl'  ``lcl''
		}

		//  By changing the returned "predict" macro, the coefficients produced
		//  by -nlcom- can be used to project ou of sample
		ereturn local predict regres_p

		ereturn scalar Nwts = wordcount(e(rhs))
		ereturn local cmd "wtsreg"
		ereturn local origcmd1 "nl"
		ereturn local origcmd2 "nlcom"
		ereturn local title "Lin. reg. by -nl- so _b are positive weights and sum to one"
		ereturn local depname `depvar'

		if "`initial'"~="" {
			local nlstartmat = substr("`initial'", strpos("`initial'","(")+1,.)
			local nlstartmat = subinstr("`nlstartmat'",")","",.)
			confirm matrix `nlstartmat'
			ereturn matrix nlstartmat = `nlstartmat'
		}


		mat def sumofwts = e(b)[1,1..`e(Nwts)']*J(`e(Nwts)',1,1)
		ereturn scalar sumofcoefs = sumofwts[1,1]
	}
	else {
		di as err "Only the -mlogit- model option is avaialble.  Other models are under construction"
		exit 198
	}
}
/*	Preserve this alternative code for the header in case it is needed later
		wtsHeader , mss(`mss') nobs(`N') mms(`mms') r2(`r2') rss(`rss')  ///
			df_r(`df_r') msr(`msr') r2_a(`r2_a')  rmse(`rmse') tss(`tss')  ///
			df_t(`df_t') dev(`dev')
*/
  // Display results
  if `e(sumofwts)'<. {
  	di _n as txt "Sum of weights is: " as res `e(sumofwts)'
  }
  Header
  ereturn display , `coeflegend' `displayopts'

  if "`model'"=="mlogit" & "verbose"~="" {
	`shhh' di as txt "Acknowledgement:"  ///
		_n _col(4) "The code for the -mlogit- model estimation option is based on the Stata blog " ///
		_n _col(4) `""How do I fit a linear regression with interval (inequality) constraints in Stata?""'  ///
		_n _col(4) "by Isabel Canette (nodate) "  ///
		`"{stata `"view browse "https://www.stata.com/support/faqs/statistics/linear-regression-with-interval-constraints/""':Accessed 28May2024.}"'
  }

end  /* End of main program wtsreg */

cap prog drop Xeq
program define Xeq
*      Version 1.0             Copied from a Stata -*_examples.ado- file.
*       Utility to echo a command to the results window (and thus to the log file)
*       and then execute it. Used e.g. in programs or loops where the command would be hidden.
        version 14

        di as txt _n
        di as txt "-> " as inp _asis `"`0'"'
*       di as txt _n

        `0'

end   /*  End of subroutine Xeq  */

program def Header
*	Copied from Stata's -nl- version 2.7.3  12may2021
        #delimit ;
        di;
        di in gr "      Source {c |}      SS            df       MS";
        di in gr "{hline 13}{c +}{hline 34}" _c;
        di in gr _col(53) "Number of obs = " in ye %10.0fc e(N) ;
        di in gr "       Model {c |} " in ye %10.0g e(mss)
                " " %10.0f e(df_m) " " %11.0g e(mms) _c;
        di in gr _col(53) "R-squared     = " in ye %10.4f e(r2);
        di in gr "    Residual {c |} " in ye %10.0g e(rss) " "
                %10.0f e(df_r) " " %11.0g e(msr) _c;
        di in gr _col(53) "Adj R-squared = " in ye %10.4f e(r2_a) ;
        di in gr "{hline 13}{c +}{hline 34}" _c;
        di in gr _col(53) "Root MSE      =  " in ye %9.0g e(rmse) ;
        di in gr "       Total {c |} " in ye %10.0g e(tss)
                " " %10.0f e(df_t) " " %11.0g e(tss)/e(df_t) _c;
        di in gr _col(53) "Res. dev.     =  " in ye %9.0g e(dev) ;
        #delimit cr
end  /* end of subroutine Header */

program def wtsHeader
	syntax [, mss(real 0) nobs(int 0) mms(real 0) r2(real 0) rss(real 0)   ///
		df_r(int 0) msr(real 0) r2_a(real 0)  rmse(real 0) tss(real 0)  ///
		df_t(int 0) dev(real 0)  ]
*	Copied from Stata's -nl- version 2.7.3  12may2021
        #delimit ;
        di;
        di in gr "      Source {c |}      SS            df       MS";
        di in gr "{hline 13}{c +}{hline 34}" _c;
        di in gr _col(53) "Number of obs = " in ye %10.0fc `N' ;
        di in gr "       Model {c |} " in ye %10.0g `mss'
                " " %10.0f `df_m' " " %11.0g `mms' _c;
        di in gr _col(53) "R-squared     = " in ye %10.4f `r2';
        di in gr "    Residual {c |} " in ye %10.0g `rss' " "
                %10.0f `df_r' " " %11.0g `msr' _c;
        di in gr _col(53) "Adj R-squared = " in ye %10.4f `r2_a' ;
        di in gr "{hline 13}{c +}{hline 34}" _c;
        di in gr _col(53) "Root MSE      =  " in ye %9.0g `rmse' ;
        di in gr "       Total {c |} " in ye %10.0g `tss'
                " " %10.0f `df_t' " " %11.0g `tss'/`df_t' _c;
        di in gr _col(53) "Res. dev.     =  " in ye %9.0g `dev' ;
        #delimit cr
end  /* end of subroutine wtsHeader */

prog define get_eb_start , rclass
*	Define the matrix of starting values for -nl-

	syntax varlist [if] [in] , depvar(varname) [nlstart(string) nlstore(string) verbose]

	local wtsrhs `varlist'
	local Nwtsrhs = wordcount("`varlist'")


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

	if "`nlstart'"=="" & "`verbose'"~="" {
		di as txt "Without user-specified starting values for -nl-, -wtsreg- defaults to using -cnsreg- starting values."
	}

	//  Use here to apply listwise deletion to -regress- or -cnsreg-
	marksample touse

	//  Vector of starting coefficients to be returned
	//  Coefficients mst be in the same order that they appear
	//  on -nl-'s command line.
	tempname eb_start

	// Default is -eb_start- is from -cnsreg-
	if "`nlstart'"=="" | "`nlstart'"=="cnsreg" {
		//	Define the sum of the named estimated coefficients from -regress- or -cnsreg-
		foreach v of varlist `wtsrhs' {
			local sumoflincoefs `sumoflincoefs' `plus' _b[`v']
			local plus +
		}

		constraint 1  `sumoflincoefs' == 1
		cap cnsreg `depvar' `wtsrhs' if `touse' , constraint(1)

		if _rc>0  {
			di as err "At -wtsreg- initiation, -cnsreg- fails to provide starting values."  ///
				_n "Try a different -nlstart()- or -initial()- option."
			error
		}
		`shhh' cnsreg `depvar' `wtsrhs' if `touse' , constraint(1)

		if "`nlstore'"~="" {
			est store cnsreg_`nlstore'
		}
		mat def `eb_start' = e(b)[1,1..`Nwtsrhs']
	}
	else if "`nlstart'"=="regress" | "`nlstart'"=="reg" | "`nlstart'"=="ols" {
		cap reg `depvar' `wtsrhs' if `touse' , constraint(1)
		if _rc>0 {
			di as err "At -wtsreg- initiation, -regress- fails. Try a different -nl_init- option."
			error
		}
		`shhh' reg `depvar' `wtsrhs' if `touse'

		if "`nlstore'"~="" {
			est store reg_`nlstore'
		}
		mat def `eb_start' = e(b)[1,1..`Nwtsrhs']
	}
	else if "`nlstart'"=="naive" | "`nlstart'"=="equal" {
		foreach v of varlist `wtsrhs' {
			mat def `eb_start' = nullmat(`eb_start'), 1/`Nwtsrhs'
		}
	}
	else {
		di as err "The -nlstart()- option is misspecified"
		error 198
	}

	return matrix eb_start = `eb_start'

end  /* End of subroutine -get_eb_start-  */
exit

*   Ver.	Date		Modifications
*	0.1		30May2024	Added replay() code, option -nlstore()-,
*						  aw and iw weights, display of "sumofwts"
*						  Edited output under the -verbose- option
*	0.2		31May2024	Add option -iterate(#)- to both -nl- and -nlcom-
*						Add `if' & `in' to -nl-, but avoid listwise deletion,
*						  because synthetic control does not use it.
*						Code the -depvar- for -nl-
*	0.3		31May2024	Add option -nstart(cnsreg)- to provide default
*						  starting values for -nl-'s coefficient estimates

*TODO:
*	1.	One or two other estimation approaches (??):
*			a.	-invlogit- which constrains the coefficients to the unit intrval,
*				but not the summation to 1.0
*			b.	-sqrt()-, -exp()- or -ln()- transformations which constrain the
*				coefficients to be non-negative or positive, but not the summation to 1.0
*			a.	-regress- or -cnsreg- for debugging comparison only.
*				Clearly less efficient than simply using Stata's -regress- or -cnsreg-.
*				but not the signs or magnitudes of the coefficients
*	2.	Test -wtsreg- performance with factor variables and alternative
*		clustering options.
*	3.	Add option to use -nl- instead of -nl- following Canette's companion blog post.
*	4.	Accommodate the group nature of weights typical of synthetic control
*		For sythetic control, the observations are nested within G groups
*		defined by the distinct treatment and control units.
*		Typically each distinct unit (i.e. group of observations)
*		is assigned an "importance" weight such that the sum of the G
*		distinct weights equals one. -wtsreg- could recognize this
*		possible structure of the weights by:
*			a.	Providing descriptive statistics on the weights
*				within and between groups.
*	5.	Could a program be designed to call -wtsreg- repeatedly by -ml-
*		in order to find the set of group weights that best fit the
*		sythetic control training data?