program define flexfit, eclass
*!	Version 1.5	29Dec2017 by Mead Over
*	This program is a post-estimation command for Stata estimation programs.
*	The program takes as a required option -chainlhs(string) which specifies
*	one of four eligible monotonic transformations that were applied to the
*	dependent variable prior to estimation.

*	Then this program:
*		a.	Computes the predicted values over the estimation sample
*			from the estimated equation
*			(1)	Computes the squared correlation of -depvar- and -dephat- as e(r2_hat)
*		b.	Passes the predicted values through the inverse of the
*			transformation function
*		c.	Computes the square root of the mean of the sum of squared deviations
*			between the actual values of the dependent variable and
*			the unwound predicted values and appends it to the -ereturn-
*			space as the scaler
*				-e(flex_rmsd)-
*		d.	Computes the correlation coefficient between the actual and
*			predicted values of the untransformed dependent variable
*			and its square and adds them to the ereturn space as scalars:
*				-e(flex_r)-
*				-e(flex_r2)-
*		e.	The -atspec- option is not yet implemented.
*	Any other options are passed through to the predict command

	syntax [name] [if] [in], chainlhs(string) [ lambda(real 0)  ///
			double replace depvar(varname) /* atspec(string) */  * ]

	if "`e(cmd)'" =="" {
		di _n as err "No estimation results are currently in memory."
		exit 198
	}

	marksample touse
		tempvar esample
		gen byte `esample' = 1.0 if e(sample) < .
		markout `touse'  `esample'
		local andtouse & `touse'

	tempvar tmpdepvar dephat undepvar undephat
	if "`depvar'"=="" {
		local depvar = word("`e(depvar)'",1)
	}
	clonevar `tmpdepvar' = `depvar'

	*	Generated fitted values from the preceding estimation
	predict `double' `dephat'  if `touse', xb `options'

	*	Unwinding the lhs transformation:
	if "`chainlhs'" == "sqrt"  {
		qui gen `double' `undepvar' = `tmpdepvar'^2  if `touse'
		qui gen `double' `undephat' = `dephat'^2  if `touse'
	}
	else if "`chainlhs'" == "ln" {
		qui gen `double' `undepvar' =  exp(`tmpdepvar')  if `touse'
		qui gen `double' `undephat' =  exp(`dephat')  if `touse'
	}
	else if "`chainlhs'" == "log10" {
		qui gen `double' `undepvar' =  10^(`tmpdepvar')  if `touse'
		qui gen `double' `undephat' =  10^(`dephat')  if `touse'
	}
	else if "`chainlhs'" == "boxcox" | "`chainlhs'" == "bc" {
		if `lambda' > 0 & `lambda' < .  {
			qui gen `double' `undepvar' =  (`lambda'*`tmpdepvar' +1)^(1/`lambda')  if `touse'
			qui gen `double' `undephat' =  (`lambda'*`dephat' +1)^(1/`lambda')  if `touse'
		}
		else {
			di _n as err "When the -chainlhs- option specifies the boxcox transformation,"  ///
				_n as err "the option lambda must be specified and given a positive value."  ///
				_n as err "Here the value of lambda is " as res "`lambda'"
			exit 198
		}
	}
	else {
		di _n as err "The option -chain(trans)- is required and the specified transformation"  ///
			_n as err "must be one of the following: sqrt, ln, log10, boxcox."
		exit 198
	}

	*	Computing the summary statistics
	qui corr `depvar' `dephat'  if `touse'
		ereturn scalar r2_hat = r(rho)^2
	di _n as txt "Over " as res %-9.0fc `r(N)' as txt "observations, the correlation between the predicted"  ///
		_n as txt "and actual dependent variable is: " as res %-15.4fc `r(rho)'  ///
		_n as txt "and its squared value , -r2_hat-, is: " as res as res %-15.4fc `r(rho)'^2

	tempvar ssdev
	qui gen `double' `ssdev' = (`undepvar' - `undephat')^2  if `touse'
		sum `ssdev'  if `touse', meanonly
		ereturn scalar flex_rmsd = sqrt(r(mean))
		ereturn scalar flex_N = r(N)
	di _n as txt "Over " as res %-9.0fc `r(N)' as txt "observations, the root mean squared deviation"  ///
		_n as txt "of the unwound predicted from the unwound actual dependent variable is: " as res %-15.4fc sqrt(`r(mean)')

	qui corr `undepvar' `undephat'  if `touse'
		ereturn scalar flex_r = r(rho)
		ereturn scalar flex_r2 = r(rho)^2
	di _n as txt "Over " as res %-9.0fc `r(N)' as txt "observations, the correlation between the unwound predicted"  ///
		_n as txt "and the unwound actual dependent variable, -flex_r-, is: " as res %-15.4fc `r(rho)'  ///
		_n as txt "and its squared value, -flex_r2-, is: " as res as res %-15.4fc `r(rho)'^2

	//	When user spcifies [name], flexfit creates a new variable to contain the unwound version
	//  of -depvar-.  With the -replace- optionk this new variable replaces
	//	any preexisting variable with this same name
	if "`namelist'"~="" & "`atspec'"=="" {
		if "`replace'"=="replace" {
			capture drop `namelist'
		}
		qui gen `double' `namelist' = `undephat' if `touse'
			lab var `namelist' "Unwound version of predicted `depvar': Generated by flexfit"
		ereturn local flexname `namelist'
	}

	/*	Under construction
	if "`namelist'"~="" & "`atspec'"~="" {
		if "`replace'"=="replace" {
			capture drop `namelist'
		}
		di as txt "Unwound predicted variable generated for `atspec'"
		qui gen `double' `namelist' = `undephat'
			lab var `namelist' "Unwound version of predicted `depvar': Generated by flexfit"
		ereturn local flexname `namelist'
	}
	*/

end
*	Version 1.0	13Feb2017 by Mead Over
*	Version 1.1	17Feb2017 the macro e(depvar) contains multiple dependent variables,
*			choose the first one.
*	Version 1.2	20Feb2017 Fix bug in above change.
*	Version 1.3	21Oct2017 Add output on the squared correlation between -depvar- and -dephat-: e(r2_hat)
*		Same as r2 for -regress-, but not computed by many other estimators
*	Version 1.4 22Nov2017 Small edits and bug fixes.
*	Version 1.5 29Dec2017 Allow -bc- as a substitute for the -boxcox- option
*	Ideas for the future:
*		Allow -atspec- so that predicted value is for specified values
*		of some right-hand-side variables.  Perhaps replace -predict- with
*		-margins-, at() with asobserved on all but the named variables,
*		which is set to the -at- specific values.
*		Adapt to perform same calculations on an equation other than the first
*		of a multi-equation system estimator like -sureg-.