*! version 1.0.5  09feb2000
program define boxcox2 , eclass
	version 6.0
	if replay() {
		Display `0'
		exit
	}

	syntax varlist [if] [in] [fweight iweight] /*
		*/ [,  NOLOG ITERate( int 50000 ) MODEL(string)/*
		*/ noConstant LRTEST NOLOGLR FROM(string) /*
		*/ NOTRans(varlist) GRaph SAVing(string) MEAN MEDian /*
		*/ Level(integer $S_level) GRVars LStart(real 1.0) /*
		*/ DELta(real 0.0 ) Zero(real 0.0 ) Generate(string)  * ]

/* !! make it so that constant only model works for default lhsonly */

	marksample touse
	markout `touse' `notrans'

	tempname sig bvec

	global T_bvec `bvec'

	if "`nologlr'" != "" {
		local nologlr "nolog"
	}

	local its `iterate'
	local iterate "iterate(`iterate')"

	mlopts mlopts, `options' `iterate'

/* Old options */

	if `zero' != 0.0 {
		di in gr "zero option is no longer supported -- ignored)"
	}

	if `delta' != 0.0 {
		di in gr "delta option is no longer supported -- ignored)"
	}

	if "`generat'" != "" {
		di in gr "generate option is no longer supported -- ignored)"
	}

	if "`grvars'" != "" {
		di in gr "grvars option is no longer supported -- ignored)"
	}
	if "`graph'" != "" {
		di in gr "(graph option is no longer supported -- ignored)"
	}
	if "`saving'" != "" {
		di in gr "(saving option is no longer supported -- ignored)"
	}
	if "`mean'" != "" {
		di in gr "(mean option is no longer supported -- ignored)"
	}
	if "`graph'" != "" {
		di in gr "(median option is no longer supported -- ignored)"
	}



/* Determine Model */
	mparse, `model'
	local model "`s(model)'"

/* Set Globals from options */

	if "`weight'" != "" {
		tempvar wvar
		quietly gen double `wvar' `exp' if `touse'
		global T_wtexp "[`weight'=`wvar']"
	}
	else {
		global T_wtexp
	}

      	global T_nocns "`constan'"


/* make global names of transformed variables */
/* perhaps reset `model' */

	tokenize `varlist'

	local dep "`1'"
	global T_dep="`dep'"
	macro shift

	global T_parm ""
	local T_i 1

	while "`1'" ~= "" {
		if "`model'"=="lhsonly" {
			local notrans `notrans' `1'
		}
		else {
			global T_parm  $T_parm `1'
			local T_i=`T_i'+1
		}
		macro shift
	}

	if `T_i'==1 {
		if "`model'" != "lhsonly" {
			di in red "Model not identified"
			di in red "No transformed variables leaves lambda " /*
			*/ "unidentified"
			di in red "Specify lhsonly model or transformed " /*
			*/ "variables"
			exit 102
		}
	}
	global T_notr `notrans'

/* check for collinear variables */


	_rmcoll $T_parm $T_notr
	global T_notr `r(varlist)'
	tokenize $T_parm
	global T_parm
	local i 1
	while "``i''" != "" {
		global T_notr : subinstr global T_notr "``i''" "" , /*
			*/ word count(local ct)
		if `ct' {
			global T_parm $T_parm ``i''
		}
		local i = `i' + 1
	}


/* would remove collinearity separately -- historical
	if "$T_parm" != "" {
		_rmcoll($T_parm)
		global T_parm `r(varlist)'
	}

	if "$T_notr" != "" {
		_rmcoll($T_notr)
		global T_notr `r(varlist)'
	}
*/



/*	exit with error if an observation is <=0 */

	tokenize `varlist'
	if "`model'"=="lhsonly" {
		ChkVar `1' `touse'
	}
	else {
		if "`model'"=="rhsonly" {
			mac shift
		}
		while "`1'" != "" {
			ChkVar `1' `touse'
			mac shift
		}
	}




/* IF THIS IS A RHS ONLY MODEL DO IT HERE AND EXIT */
	if "`model'" =="rhsonly"  {
		tempname bvec bvec2 sig_r
		global T_bvec `bvec2'
		global T_sig `sig_r'

/* if there is a constant, do constant only model  */
		if "$T_nocns" == "" {
			tempvar err1  sig1 llft
			tempname LL0
			quietly {

				regress `dep' if `touse' $T_wtexp
				predict double `err1', residual
				replace `err1'=`err1'*`err1'
				summ `err1' $T_wtexp if `touse', meanonly
				scalar `sig1'=r(mean)
				gen double `llft'=-.5*(1+ln(2*_pi)+ln(`sig1'))
				summ `llft' $T_wtexp if `touse', meanonly
				scalar `LL0'=r(sum)
			}
			local  df0 2
		}

		if "`from'" != "" {
			local start "init(`from')"
		}
		else {
			if "`lstart'" != "" {
				local start "init(`lstart', copy)"
			}
		}
		di in gr _n "Estimating full model"

		#delimit ;
		ml model rdu0  "boxco2_l RHS" (lambda: `dep' = )
			if  `touse'  $T_wtexp, technique(nr) search(off)
			`start' `mlopts' `nolog'
			max  ;
		#delimit cr

		tempname llf bf Vf bpar
		scalar `llf' = e(ll)
		local ic=e(ic)
		local rc=e(rc)
		local wt "`e(wtype)'"
		local wex "`e(wexp)'"

		matrix `Vf'=e(V)
		matrix `bpar'=e(b)

/* make call to evaluator with converged lambda */
		global ML_y1 `dep'
		global ML_samp `touse'

		tempname lamc bc
		matrix `bc'=e(b)
		scalar `lamc'=`bc'[1,1]

		tempvar lltmp
		gen double `lltmp' = 0
		boxco2_l RHS `lltmp' `bc' `lltmp'

		matrix `bf'=$T_bvec
		macro drop ML_y1
		macro drop ML_samp


/* Do LR tests lambda=-1,0,1 */
		tempname ll_t0 chit0 p_t0 ll_t1 chit1 p_t1
		tempname ll_tm1 chitm1 p_tm1

		lincompR `ll_t1'  `touse' /* program for computing LL
						when lam=1*/
		scalar `chit1'=2*(`llf'-`ll_t1')
		scalar `p_t1'=chiprob(1, `chit1')


		invcompR `ll_tm1' `touse' /* program for computing LL
						when lam=-1 */
		scalar `chitm1'=2*(`llf'-`ll_tm1')
		scalar `p_tm1'=chiprob(1, `chitm1')

		logcompR  `ll_t0' `touse' /* program for computing LL
						when lam=0 */
		scalar `chit0'=2*(`llf'-`ll_t0')
		scalar `p_t0'=chiprob(1, `chit0')

		if "`lrtest'" != ""	{
			di in gr _n "Estimating comparison models for LR tests"
			local firstlr 1

			local notr "$T_notr"
			local parm "$T_parm"

			tempname lltr pv p chi2 df


			tokenize $T_notr
			while "`1'" != "" {
				global T_notr : subinstr global T_notr /*
					*/ "`1'" "", word

				#delimit ;
				ml model rdu0  "boxco2_l RHS"
					(lambda: `dep' = )
					if `touse'  $T_wtexp, technique(nr)
					search(off) `mlopts'
					`start'  max `nologlr' ;
				#delimit cr

				scalar `lltr'=2*(`llf'-e(ll))
				scalar `pv' = chiprob(1, `lltr' )

				matrix `chi2' = nullmat(`chi2') ,	/*
					   */ cond(`lltr'==., -1, `lltr')
				matrix `p' = nullmat(`p') , 		/*
					*/ cond(`pv'==., -1, `pv')
				matrix `df' = nullmat(`df') , 1

				local pname " `pname' `1'"
				local pcol " `pcol' Notrans "

				local c2name " `c2name' `1'"
				local c2col " `c2col' Notrans "

				global T_notr "`notr'"
				macro shift
			}

			local n : word count $T_parm

			tokenize $T_parm

			while "`1'" != "" {
				global T_parm : subinstr global T_parm /*
					*/  "`1'" "", word
				if `n'==1 {
					di in bl _n /*
					   */ "Comparison model for LR " /*
					   */ "test on `1' is a linear " /*
					   */ "regression"
					di in bl "Lambda is not identified " /*
					   */ "in the restricted  model "

					quietly {
regress `dep' $T_notr if `touse' $T_wtexp, $T_nocns
tempvar err llf2
tempname sigmat llt
predict double `err', residuals
replace `err'=`err'^2
summ `err' if `touse' $T_wtexp, meanonly
scalar `sigmat'=r(mean)
gen double `llf2'=-.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
summ `llf2' $T_wtexp if `touse'
scalar `llt'  =r(sum)
est scalar ll = `llt'

local dfs 2
					}
				}
				else {
					#delimit ;
					ml model rdu0  "boxco2_l RHS"
						(lambda: `dep'= )
						if  `touse'  $T_wtexp,
						technique(nr) search(off)
						`mlopts'
						`start'  max `nologlr' ;
					#delimit cr
					local dfs 1
				}

				scalar `lltr'=2*(`llf'-e(ll))
				scalar `pv' = chiprob(`dfs', `lltr' )

				matrix `chi2' = nullmat(`chi2') ,	/*
					   */ cond(`lltr'==., -1, `lltr')
				matrix `p' = nullmat(`p') , 		/*
					*/ cond(`pv'==., -1, `pv')
				matrix `df' = nullmat(`df') , 1

				local pname " `pname' `1'"
				local pcol " `pcol' Trans "

				local c2name " `c2name' `1'"
				local c2col " `c2col' Trans "


				global T_parm "`parm'"
				macro shift
			}

			matrix colnames `p'= `pname'
			matrix coleq `p' = `pcol'

			matrix colnames `chi2'= `c2name'
			matrix coleq `chi2' = `c2col'

			matrix colnames `df'= `c2name'
			matrix coleq `df' = `c2col'
		}

/* CALL POST PROGRAM */
		BasePost `touse' "`model'" `bf' `Vf' `bpar'

		if "`lrtest'" != "" {
			est matrix chi2m `chi2'
			est matrix pm `p'
			est matrix df `df'
		}

		if "$T_nocns" == "" {
			tempname llt
			scalar `llt'=2*(`llf'-`LL0')
			est scalar ll0=`LL0'
			est scalar df_r=`df0'
			est scalar chi2=`llt'
		}

		est scalar ll=`llf'

		est scalar chi2_t1=`chit1'
		est scalar p_t1 =`p_t1'
		est scalar ll_t1=`ll_t1'

		est scalar chi2_tm1=`chitm1'
		est scalar p_tm1 =`p_tm1'
		est scalar ll_tm1=`ll_tm1'

		est scalar chi2_t0=`chit0'
		est scalar p_t0 =`p_t0'
		est scalar ll_t0=`ll_t0'

		est scalar ic=`ic'
		est scalar rc=`rc'
		est local wtype="`wt'"
		if "`wt'" != "" {
			est local wexp="`exp'"
		}
		est local depvar  "`dep'"
		est local lrtest "`lrtest'"
		est local chi2type "LR"
		est local predict "boxco2_p"
		est local model  "`model'"
		if "`ntrans'"=="" & "$T_nocns" != "" {
			est local ntrans
		}
		else 	est local ntrans "yes"

		est local cmd "boxcox2"
		macro drop T_*

/* Call display programs */
		Display, level(`level')
		exit

	}
/* END RHS ONLY MODEL */


/* ESTIMATE COMPARISON MODEL, If there is a constant in the model  */

	if "$T_nocns" == "" & `its' > 0 {
		di in gr "Estimating comparison model"

		tempname bvec1 LL0
		matrix `bvec1'=1

		ml model rdu0 "boxco2_l ConsOnly" ( `dep' = ) /*
		    */    if  `touse'  $T_wtexp, technique(nr)  search(off) /*
		    */  init(1, copy)   max noheader `nolog'

		scalar `LL0'=e(ll)
		local df0 3
	}


/* IF THIS IS A LEFT HAND SIDE ONLY MODEL DO IT HERE  */

	if "`model'" == "lhsonly" {
		tempname btvec bvec
		global T_bvec "`btvec'"

	 	/* get starts */


		if "`from'" != "" {
			di "from is `from'"
			local start "init(`from')"
		}
		else {
			if "`lstart'" != "" {
				local start "init(`lstart', copy)"
			}
		}

		di in gr _n "Estimating full model"

		#delimit ;
		ml model rdu0 "boxco2_l LHS" (ntrans: `dep' = )
		     if  `touse'  $T_wtexp, technique(nr)
  		    `mlopts' `nolog' search(off)
		     `start'   max  ;
		#delimit cr

/* Save off value of LL and other values of interest */

		tempname llf
		scalar `llf'=e(ll)

		local ic=e(ic)
		local rc=e(rc)
		local wt "`e(wtype)'"
		local wex "`e(wexp)'"

		tempname bf Vf bpar

		matrix `Vf'=e(V)
		matrix `bpar'=e(b)

/* make call to evaluator with converged lambda */
		global ML_y1 `dep'
		global ML_samp `touse'

		tempname lamc bc
		matrix `bc'=e(b)
		scalar `lamc'=`bc'[1,1]

		tempvar lltmp
		gen double `lltmp' = 0
		boxco2_l LHS `lltmp' `bc' `lltmp'

		matrix `bf'=$T_bvec
		macro drop ML_y1
		macro drop ML_samp

/* Do LR tests theta=-1,0,1 in LHS model */
		tempname ll_t1 chit1 p_t1 ll_tm1 chitm1 p_tm1
		tempname ll_t0 chit0 p_t0

		lincompL `ll_t1'  `touse' /* program for computing LL
						when theta=1*/
		scalar `chit1'=2*(`llf'-`ll_t1')
		scalar `p_t1'=chiprob(1, `chit1')


		invcompL `ll_tm1' `touse' /* program for computing LL
						when theta =-1 */
		scalar `chitm1'=2*(`llf'-`ll_tm1')
		scalar `p_tm1'=chiprob(1, `chitm1')

		logcompL  `ll_t0' `touse' /* program for computing LL
						when theta =0 */
		scalar `chit0'=2*(`llf'-`ll_t0')
		scalar `p_t0'=chiprob(1, `chit0')

/* Do LR tests if requested */

		if "`lrtest'" != "" {
			tempname p chi2 df
			di _n in gr "Estimating comparison models for LR tests"

			local firstlr 1

			local notr "$T_notr"
			local parm "$T_parm"

			tempname lltr pv

			tokenize $T_notr
			while "`1'" != "" {
				global T_notr : subinstr global T_notr /*
					*/ "`1'" "", word

				#delimit ;
				ml model rdu0 "boxco2_l LHS"
					(ntrans: `dep' = )
					if  `touse'  $T_wtexp, technique(nr)
					search(off) `mlopts'
					`start'  max `nologlr' ;
				#delimit cr

				scalar `lltr'=2*(`llf'-e(ll))
				scalar `pv' = chiprob(1, `lltr' )

				matrix `chi2' = nullmat(`chi2') ,	/*
					   */ cond(`lltr'==., -1, `lltr')
				matrix `p' = nullmat(`p') , 		/*
					*/ cond(`pv'==., -1, `pv')
				matrix `df' = nullmat(`df') , 1

				local pname `pname' `1'
				local pcol `pcol' Notrans

				local c2name `c2name' `1'
				local c2col `c2col' Notrans

				global T_notr "`notr'"
				macro shift
			}

			tokenize $T_parm
			while "`1'" != "" {
				global T_parm : subinstr global /*
					*/ T_parm "`1'" "", word

				#delimit ;
				ml model rdu0  "boxco2_l LHS"
					(ntrans: `dep' = )
					if  `touse'  $T_wtexp, technique(nr)
					search(off) `mlopts'
					`start'  max `nologlr' ;
				#delimit cr


				scalar `lltr'=2*(`llf'-e(ll))
				scalar `pv' = chiprob(1, `lltr' )

				matrix `chi2' = nullmat(`chi2') ,	/*
					   */ cond(`lltr'==., -1, `lltr')
				matrix `p' = nullmat(`p') , 		/*
					*/ cond(`pv'==., -1, `pv')
				matrix `df' = nullmat(`df') , 1

				local pname `pname' `1'
				local pcol `pcol' Notrans

				local c2name `c2name' `1'
				local c2col `c2col' Notrans

				global T_parm "`parm'"
				macro shift
			}
			matrix colnames `p'= `pname'
			matrix coleq `p' = `pcol'

			matrix colnames `chi2'= `c2name'
			matrix coleq `chi2' = `c2col'

			matrix colnames `df'= `c2name'
			matrix coleq `df' = `c2col'
		}

/* CALL POST PROGRAM */
		BasePost `touse' "`model'" `bf' `Vf' `bpar'

		if "`lrtest'" != "" {
			est matrix chi2m `chi2'
			est matrix pm `p'
			est matrix df `df'
		}


		if "$T_nocns" == "" {
			local llt=2*(`llf'-`LL0')
			est scalar ll0=`LL0'
			est scalar df_r=`df0'
			est scalar chi2=`llt'
		}

		est scalar ll=`llf'

		est scalar chi2_t1=`chit1'
		est scalar p_t1 =`p_t1'
		est scalar ll_t1=`ll_t1'

		est scalar chi2_tm1=`chitm1'
		est scalar p_tm1 =`p_tm1'
		est scalar ll_tm1=`ll_tm1'

		est scalar chi2_t0=`chit0'
		est scalar p_t0 =`p_t0'
		est scalar ll_t0=`ll_t0'

		est scalar ic=`ic'
		est scalar rc=`rc'
		est local wtype="`wt'"
		if "`wt'" != "" {
			est local wexp="`exp'"
		}

		est local depvar "`dep'"
		est local lrtest "`lrtest'"
		est local chi2type "LR"
		est local predict  "boxco2_p"
		est local model  "`model'"

		if "`ntrans'"=="" & "$T_nocns" != "" {
			est local ntrans ""
		}
		else 	est local ntrans "yes"

		est local cmd "boxcox2"
		macro drop T_*

		Display, level(`level')
		exit
	}
/* END LHS ONLY MODEL */

/* THIS SECTION DOES FULL MODEL WITH COMMON LAMBDA OR
		LAMBDA ON THE RHS AND THETA ON THE LHS */

	tempname bvec bvec2 sig
	global T_bvec "`bvec'"
	global T_sig  "`sig'"	 /* use the globals to work on
					the matrix and the scalar in the
					evaluator program */

	if "`model'" == "theta" {
		local lrhs "/theta"
		local evalpgm "Theta"
	}
	else {
		local evalpgm "Lambda"
	}

/* Get starts for full model */

	if `lstart' != 1 {
		di in bl "lstart only valid for LHS or RHS model" /*
		*/ "--option ignored"
	}

	if "`from'" != "" {
		local starts  "init(`from')"
	}
	else {
		if "`model'"=="theta" {
		   	matrix `bvec'= (1, 1)
 		}
    		else {
	   		matrix `bvec'= 1
    		}

		local starts "init(`bvec', copy)"
	}

	di in gr _n "Estimating full model"

 	ml model rdu0 "boxco2_l `evalpgm'" (ntrans: `dep' = ) /*
		*/  `lrhs' if  `touse'  $T_wtexp, technique(nr)  /*
  		*/ search(off) `mlopts' /*
		*/ `starts'   max  `nolog'

	tempname llf bf Vf bpar
	scalar `llf' = e(ll)

	local ic=e(ic)
	local rc=e(rc)
	local wt "`e(wtype)'"
	local wex "`e(wexp)'"

	matrix `Vf'=e(V)
	matrix `bpar'=e(b)

/* make call to evaluator with converged lambda */
		global ML_y1 `dep'
		global ML_samp `touse'

		tempname lamc bc
		matrix `bc'=e(b)

		tempvar lltmp
		gen double `lltmp' = 0
		boxco2_l `evalpgm' `lltmp' `bc' `lltmp'

		matrix `bf'=$T_bvec

		macro drop ML_y1
		macro drop ML_samp

/* Do LR tests lambda=theta=-1,0,1  Note Lambda and Theta
	models have same restricted LL in this case  */

	tempname chit1 p_t1 chitm1 p_tm1 chit0 p_t0 ll_t1 ll_tm1 ll_t0

	lincompT `ll_t1'  `touse' /* program for computing LL
						when theta=1*/
	scalar `chit1'=2*(`llf'-`ll_t1')
	scalar `p_t1'=chiprob(1, `chit1')

	invcompT `ll_tm1' `touse' /* program for computing LL
						when theta=-1 */
	scalar `chitm1'=2*(`llf'-`ll_tm1')
	scalar `p_tm1'=chiprob(1, `chitm1')

	logcompT  `ll_t0' `touse' /* program for computing LL
						when theta =0 */
	scalar `chit0'=2*(`llf'-`ll_t0')
	scalar `p_t0'=chiprob(1, `chit0')

	if "`lrtest'" != ""	{

		di _n in gr "Estimating comparison models for LR tests"
		tempname p chi2 df

		local firstlr 1

		local notr "$T_notr"
		local parm "$T_parm"

		tempname lltr pv

		tokenize $T_notr
		while "`1'" != "" {
			global T_notr : subinstr global T_notr /*
				*/  "`1'" "", word

			#delimit ;
			ml model rdu0 "boxco2_l `evalpgm'"
				(ntrans: `dep' = ) `lrhs'
				if  `touse'  $T_wtexp, technique(nr)
				search(off) `mlopts'
				`start'  max `nologlr' ;
			#delimit cr

			scalar `lltr'=2*(`llf'-e(ll))
			scalar `pv' = chiprob(1, `lltr' )

			matrix `chi2' = nullmat(`chi2') ,	/*
				   */ cond(`lltr'==., -1, `lltr')
			matrix `p' = nullmat(`p') , cond(`pv'==., -1, `pv')
			matrix `df' = nullmat(`df') , 1

			local pname `pname' `1'
			local pcol `pcol' Notrans

			local c2name `c2name' `1'
			local c2col `c2col' Notrans

			global T_notr "`notr'"
			macro shift
		}

		local n : word count $T_parm

		tokenize $T_parm
		while "`1'" != "" {
			global T_parm : subinstr global /*
				*/  T_parm "`1'" "", word

			if `n'==1  & "`model'"=="theta" {
				di in bl _n "Comparison model " /*
					*/ "for `1' is left-hand " /*
					*/ "side only."
				di in bl "Lambda is not " /*
					*/ "identified in the "
				di in bl "restricted model "

				local dfs 2
				#delimit ;
				ml model rdu0 "boxco2_l LHS"
					(ntrans: `dep' = ) if  `touse'
					$T_wtexp, technique(nr)
					search(off) `mlopts'
					`start'  max `nologlr' ;
				#delimit cr
			}
			else {
				local dfs 1
				#delimit ;
				ml model rdu0  "boxco2_l `evalpgm'"
					(ntrans: `dep' = )
					`lrhs' if  `touse'
					$T_wtexp, technique(nr)
					search(off) `mlopts' `start'
					max `nologlr' ;
				#delimit cr
			}

			scalar `lltr'=2*(`llf'-e(ll))
			scalar `pv' = chiprob(`dfs', `lltr' )

			matrix `chi2' = nullmat(`chi2') ,	/*
				   */ cond(`lltr'==., -1, `lltr')
			matrix `p' = nullmat(`p') , cond(`pv'==., -1, `pv')
			matrix `df' = nullmat(`df') , 1

			local pname `pname' `1'
			local pcol `pcol' Notrans

			local c2name `c2name' `1'
			local c2col `c2col' Notrans

			global T_parm "`parm'"
			macro shift
		}
		matrix colnames `p'= `pname'
		matrix coleq `p' = `pcol'

		matrix colnames `chi2'= `c2name'
		matrix coleq `chi2' = `c2col'

		matrix colnames `df'= `c2name'
		matrix coleq `df' = `c2col'

	}

/* CALL POST PROGRAM */
	BasePost `touse' "`model'" `bf' `Vf' `bpar'

	if "`lrtest'" != "" {
		est matrix chi2m `chi2'
		est matrix pm `p'
		est matrix df `df'
	}

	if "$T_nocns" == "" {
		local llt=2*(`llf'-`LL0')
		est scalar df_r=`df0'
		est scalar chi2=`llt'
		est scalar ll0=`LL0'
	}


	est scalar ll=`llf'
	est scalar chi2_t1=`chit1'
	est scalar p_t1 =`p_t1'
	est scalar ll_t1=`ll_t1'

	est scalar chi2_tm1=`chitm1'
	est scalar p_tm1 =`p_tm1'
	est scalar ll_tm1=`ll_tm1'

	est scalar chi2_t0=`chit0'
	est scalar p_t0 =`p_t0'
	est scalar ll_t0=`ll_t0'

	est scalar ic=`ic'
	est scalar rc=`rc'
	est local wtype="`wt'"
	if "`wt'" != "" {
		est local wexp="`exp'"
	}
	est local depvar "`dep'"
	est local lrtest "`lrtest'"
	est local chi2type "LR"
	est local predict  "boxco2_p"

	if "`ntrans'"=="" & "$T_nocns" != "" {
		est local ntrans ""
	}
	else {
		est local ntrans "yes"
	}

	est local model  "`model'"
	est local cmd "boxcox2"

	macro drop T_*

/* Call display programs */
	Display, level(`level')
end


program define mparse, sclass
	sret clear
	syntax, [LHSonly RHSonly LAMbda THETA ]

	if "`lhsonly'`rhsonly'`lambda'`theta'" == "" {
		sret local model "lhsonly"
	}
	else {
		local sum =  ("`lhsonly'" != "") + ("`rhsonly'" != "") /*
			*/ + ("`lambda'" != "") + ("`theta'" != "")
		if `sum' > 1 {
			di in red /*
*/ "model specifications are mutually exlusive:  specify only one model"
			exit 198
		}
		sret local model "`lhsonly'`rhsonly'`lambda'`theta'"
	}
end

program define bhead
	di " "
	tempname p
	local dft= e(df_m)-e(df_r)
	#delimit ;
	di _col(51) in gr "Number of obs"  _col(67) "= " _col(68)
		 %10.0g in ye e(N);
	di _col(51) in gr "LR chi2(" in ye `dft' in gr ")"
		_col(67) "=" _col(69) %10.2f in ye e(chi2);
	scalar `p' = chiprob(`dft',e(chi2) );
	di in gr "Log likelihood = " in ye e(ll) _col(51) in gr "Prob > chi2"
		_col(67) "=" _col(69) in ye  %10.3f `p';
	di " ";
	#delimit cr
end

program define btitle
	#delimit ;
	di _dup(78) in gr "-" ;
	di  %9s in gr "`e(depvar)'"  %1s in gr _col(10) "|"
		%40s in gr _col(12) "Coef.   Std. Err.       z     P>|z|"
		%10s _col(59) "[95% Conf. Interval]";
	di _dup(9) in gr "-" "+" _dup(68) "-";
	#delimit cr
end

program define cdisp
	di _dup(78) in gr "-"
	di " "
	di in gr "Estimates of scale-variant parameters"
	di _dup(24) in gr "-"
	di %11s in gr _col(10) "|      Coef."

	tempname b
	matrix `b'=e(b)
	local eqs : coleq `b'
	local names : colnames `b'
	local cnt 1
	tokenize `eqs'
	if "`1'" == "Notrans" {
		di _dup(9) in gr "-" "+" _dup(14) "-"
		di in ye "Notrans" _col(10) %1s in gr "|"
	}
	while "`1'" == "Notrans" {
		local coef :  word `cnt' of `names'
		di %8s in gr "`coef'" _col(10) %1s "|"  _col(13) /*
			*/ in ye %9.0g `b'[1,`cnt']
		local cnt = `cnt' +1
		macro shift
	}

	if "`1'" == "Trans" {
		di _dup(9) in gr "-" "+" _dup(14) "-"
		di in ye "Trans" _col(10) %1s in gr "|"
	}
	while "`1'" == "Trans" {
		local coef :  word `cnt' of `names'
		di %8s in gr "`coef'" _col(10) %1s  "|" _col(13) /*
			*/ in ye %9.0g `b'[1,`cnt']
		local cnt = `cnt' +1
		macro shift
	}
	while "`1'" != "sigma" {
		local cnt = `cnt' +1
		macro shift
	}
	if "`1'" == "sigma" {
		di _dup(9) in gr "-" "+" _dup(14) "-"
		di %8s in gr "/sigma" _col(10) %1s "|"  _col(13) in ye /*
			*/  %9.0g `b'[1,`cnt']
		di _dup(24) in gr "-"
	}
end

program define cdisplr
	di _dup(78) in gr "-"
	di " "
	di in gr "Estimates of scale-variant parameters"
	di _dup(57) in gr "-"
	di %11s in gr _col(10) "|      Coef." _col(23) %7s /*
		*/ "chi2(" in ye "df" in gr ")" /*
		*/  _col(32) %10s "P>chi2(" in ye "df"  in gr")" /*
		*/ _col(48) %10s "df of chi2"
	tempname b chi p df
	matrix `b'=e(b)
	matrix `chi'=e(chi2m)
	matrix `p'=e(pm)
        matrix `df'=e(df)

	local eqs : coleq `b'
	local names : colnames `b'
	local cnt  1
	local pcnt 1
	tokenize `eqs'
	if "`1'" == "Notrans" {
		di _dup(9) in gr "-" "+" _dup(47) "-"
		di in ye "Notrans" _col(10) %1s in gr "|"
	}
	while "`1'" == "Notrans" {
		local coef :  word `cnt' of `names'
		if "`coef'" == "_cons"   {
			di %8s in gr "`coef'" _col(10) %1s "|"  _col(13) /*
				*/ in ye %9.0g `b'[1,`cnt']
		}
		else {
			di %8s in gr "`coef'" _col(10) %1s "|"  _col(13) /*
				*/ in ye %9.0g `b'[1,`cnt'] _col(25) %7.3f /*
				*/ `chi'[1,`pcnt'] _col(35) %7.3f /*
				*/ `p'[1,`pcnt'] _col(50) %2.0f `df'[1,`pcnt']
			local pcnt = `pcnt' +1
		}
		local cnt = `cnt' +1
		macro shift
	}


	if "`1'" == "Trans" {
		di _dup(9) in gr "-" "+" _dup(47) "-"
		di in ye "Trans" _col(10) %1s in gr "|"
	}

	while "`1'" == "Trans" {
		local coef :  word `cnt' of `names'
		di %8s in gr "`coef'" _col(10) %1s "|"  _col(13) in ye %9.0g /*
			*/ `b'[1,`cnt'] _col(25) %7.3f `chi'[1,`pcnt'] /*
			*/ _col(35) %7.3f `p'[1,`pcnt'] /*
			*/_col(50) %2.0f `df'[1,`pcnt']
		local cnt = `cnt' +1
		local pcnt = `pcnt' +1
		macro shift
	}
	while "`1'" != "sigma" {
		local cnt = `cnt' +1
		macro shift
	}
	if "`1'" == "sigma" {
		di _dup(9) in gr "-" "+" _dup(47) "-"
		di %8s in gr "/sigma" _col(10) %1s "|"  _col(13) in ye %9.0g /* 			*/ `b'[1,`cnt']
		di _dup(57) in gr "-"
	}

end


program define BasePost, eclass
	args touse model bf Vf bpar

	tempname b1 b V V1 newb

	matrix `b1'=`bf'

	local size = colsof(`b1')

	local eqnames ""
	local names ""
	local cnt 1
	local first 1


	tokenize $T_notr
	while "`1'" != "" {
		local eqnames `eqnames' Notrans
		local names `names' `1'
		if `first'==1 {
			matrix `newb'=`b1'[1,`cnt']
			local first 0
		}
		else 	matrix `newb'=`newb',`b1'[1,`cnt']
		macro shift
		local cnt =`cnt' +1
	}

	if "$T_nocns" == "" {
		local eqnames `eqnames' Notrans
		local names `names' _cons
		if `first'==1 {
			matrix `newb'=`b1'[1,`size']
			local first  0
		}
		else 	matrix `newb'=`newb',`b1'[1,`size']
	}

	tokenize $T_parm
	while "`1'" != "" {
		local eqnames `eqnames' Trans
		local names `names' `1'
		if `first'==1 {
			matrix `newb'=`b1'[1,`cnt']
			local first 0
		}
		else 	matrix `newb'=`newb',`b1'[1,`cnt']
		macro shift
		local cnt =`cnt' +1
	}

	matrix colnames `newb' = `names'
	matrix coleq `newb' = `eqnames'

	if "`model'"=="lhsonly" {
		local bparms "theta:_cons"
	}

	if "`model'"=="rhsonly" | "`model'"=="lambda"  {
		local bparms "lambda:_cons"
	}

	if "`model'"=="theta" {
		local bparms "lambda:_cons theta:_cons"
	}

	matrix `b1'=`bpar',sqrt($T_sig)
	matrix colnames `b1'=  `bparms'  sigma:_cons
	matrix `b'=`newb',`b1'

	local nmes : colnames `b'
	local eqnames : coleq `b'

	local size = colsof(`b')

	matrix `V'=J(`size',`size',0)
	matrix `V1'=`Vf'

	if "`model'" != "theta" {
		matrix `V'[`size'-1,`size'-1]=`V1'[1,1]
	}
	else {
		local low = `size'-2
		local high=`size'-1
		matrix `V'[`low',`low']=`V1'[1,1]
		matrix `V'[`low',`high']=`V1'[1,2]
		matrix `V'[`high',`low']=`V1'[2,1]
		matrix `V'[`high',`high']=`V1'[2,2]
	}
	matrix colnames `V'=`nmes'
	matrix coleq `V' = `eqnames'

	matrix rownames `V' = `nmes'
	matrix roweq `V' = `eqnames'

	local df=colsof(`b')
	estimates post `b' `V', dep(`e(depvar)') obs(`e(N)') /*
		*/  esample(`touse')

	est scalar df_m=`df'
end

program define lincompR
	args ll touse
	/* this program computes the value of the log likelikehood for
		the RHSONLY	model when lambda=1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse
		*/

	tokenize $T_parm
	local i 1
	local trans ""

	while "`1'" != "" {
		tempvar t`i'
		gen double `t`i''=`1'-1
		local trans `trans' `t`i''
		macro shift
		local i=`i' +1
	}

	quietly {
		regress $T_dep $T_notr `trans' if `touse' $T_wtexp, $T_nocns
		tempvar err llfv
		tempname sigmat

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llfv' = -.5*(1+ln(2*_pi)+ ln(`sigmat')) if `touse'
		summ `llfv' $T_wtexp, meanonly

		scalar `ll'  =r(sum)
	}
end

program define invcompR
	args ll touse
	/* this program computes the value of the log likelikehood
			for the RHSONLY model when lambda=-1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tokenize $T_parm
	local i 1
	local trans ""

	while "`1'" != "" {
		tempvar t`i'
		gen double `t`i''=1-(1/`1')
		local trans `trans' `t`i''
		macro shift
		local i=`i' +1
	}


	quietly	 {
		regress $T_dep $T_notr `trans' if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
		scalar `sigmat'=r(mean)
		gen double `llf'=-.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll' = r(sum)
	}
end

program define logcompR
	args ll touse
	/* this program computes the value of the log likelikehood
			for the RHSONOLY model when lambda=0
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */


	tokenize $T_parm
	local i 1
	local trans ""

	while "`1'" != "" {
		tempvar t`i'
		gen double `t`i''=ln(`1')
		local trans `trans' `t`i''
		macro shift
		local i=`i'+1
	}

	quietly	 {
		regress $T_dep $T_notr `trans' if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llf'= -.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll'=r(sum)
	}
end




program define lincompL
	args ll touse
	/* this program computes the value of the log likelikehood for
			the LHSONLY	model when theta=1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	quietly	{
		regress $T_dep $T_notr  if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llf'=-.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll'  =r(sum)
	}
end

program define invcompL
	args ll touse
	/* this program computes the value of the log likelikehood
			for the LHSONLY model when theta=-1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tempvar yt

	quietly {
		gen double `yt'=1-(1/$T_dep) if `touse'
		regress `yt' $T_notr  if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llf' = -.5*(1 + ln(2*_pi)+ ln(`sigmat') ) /*
			*/ -2*ln($T_dep) if `touse'
		summ `llf' $T_wtexp , meanonly
		scalar `ll' = r(sum)
	}
end

program define logcompL
	args ll touse
	/* this program computes the value of the log likelikehood
			for the LHSONOLY model when theta=0
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tempvar yt llf err
	tempname sigmat

	quietly	 {
		gen double `yt'=ln($T_dep) if `touse'
		regress `yt' $T_notr if `touse' $T_wtexp, $T_nocns

		predict double `err', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
  		scalar `sigmat'=r(mean)
		gen double  `llf'  = -.5*(1 + ln(2*_pi)+ ln(`sigmat') ) /*
			*/  -ln($T_dep) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll'=r(sum)
	}
end



program define lincompT
	args ll touse
	/* this program computes the value of the log likelikehood for
			the THETA model when theta=1 & lambda=1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tempvar yt

	tokenize $T_parm
	local i 1
	local trans ""

	quietly {
		gen double `yt'=$T_dep-1 if `touse'
		while "`1'" != "" {
			tempvar t`i'
			gen double `t`i''=`1'-1 if `touse'
			local trans "`trans' `t`i''"
			macro shift
			local i=`i' +1
		}
		regress $T_dep $T_notr  `trans' if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err' if `touse', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llf'=-.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll'  =r(sum)
	}
end

program define invcompT
	args ll touse
	/* this program computes the value of the log likelikehood
			for the THETA model when theta=-1 and lambda= -1
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tempvar yt


	tokenize $T_parm
	local i 1
	local trans ""

	quietly {
		gen double `yt'=1-(1/$T_dep) if `touse'
		while "`1'" != "" {
			tempvar t`i'
			gen double `t`i''=1-(1/`1') if `touse'
			local trans " `trans' `t`i'' "
			macro shift
			local i=`i' +1
		}

		regress `yt' $T_notr `trans'  if `touse' $T_wtexp, $T_nocns
		tempvar err llf
		tempname sigmat

		predict double `err' if `touse', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
   		scalar `sigmat'=r(mean)
		gen double `llf' = -.5*(1 + ln(2*_pi)+ /*
			*/ ln(`sigmat') )-2*ln($T_dep) if `touse'
		summ `llf' $T_wtexp , meanonly
		scalar `ll' = r(sum)
	}
end

program define logcompT
	args ll touse
	/* this program computes the value of the log likelikehood
			for the THETA model when theta=0 and lambda=0
		arguments
			1 is name of scalar to hold restricted  LL
			2 is name of touse */

	tempvar yt llf err
	tempname sigmat


	tokenize $T_parm
	local i 1
	local trans ""

	quietly {
		gen double `yt'=ln($T_dep) if `touse'
		while "`1'" != "" {
			tempvar t`i'
			gen double `t`i''=ln(`1') if `touse'
			local trans " `trans' `t`i'' "
			macro shift
			local i=`i' +1
		}

		regress `yt' $T_notr `trans' if `touse' $T_wtexp, $T_nocns

		predict double `err' if `touse', residuals
		replace `err'=`err'^2
		summ `err' if `touse' $T_wtexp, meanonly
  		scalar `sigmat'=r(mean)
		gen double  `llf'  = -.5*(1 + ln(2*_pi)+ ln(`sigmat') )/*
			*/  -ln($T_dep) if `touse'
		summ `llf' $T_wtexp, meanonly
		scalar `ll'=r(sum)
	}
end


program define lrdispR
	di " "
	di _dup(57) in gr "-"
	di in gr _col(4) "Test" _col(17) "Restricted" _col(32) /*
		*/ "LR statistic"  _col(50) "P-Value"
	di  in gr _col(5) "H0:"  _col(15) "log likelihood" _col(35) "X~chi2" /*
		*/ _col(49) "Pr > chi2"
	di _dup(57) in gr "-"
	di in gr "lambda = -1" _col(17)  %10.0g in ye e(ll_tm1) /*
		*/ _col(32)  in ye %8.2f e(chi2_tm1) /*
		*/  _col(51) %5.3f in ye e(p_tm1)

	di in gr "lambda =  0" _col(17) %10.0g in ye e(ll_t0) /*
		*/ _col(32) in ye %8.2f e(chi2_t0) /*
		*/ _col(51)  %5.3f in ye e(p_t0)

	di in gr "lambda =  1" _col(17)  %10.0g in ye e(ll_t1) /*
		*/ _col(32) in ye %8.2f e(chi2_t1) /*
		*/ _col(51)  %5.3f in ye e(p_t1)

	di _dup(57) in gr "-"
end

program define lrdispL
	di " "
	di _dup(57) in gr "-"
	di in gr _col(4) "Test" _col(17) "Restricted" _col(32) /*
		*/ "LR statistic"  _col(50) "P-Value"
	di  in gr _col(5) "H0:"  _col(15) "log likelihood" _col(35) "X~chi2" /*
		*/ _col(49) "Pr > chi2"
	di _dup(57) in gr "-"
	di in gr "theta = -1" _col(17)  %10.0g in ye e(ll_tm1) /*
		*/ _col(32)  in ye %8.2f e(chi2_tm1) /*
		*/  _col(51) %5.3f in ye e(p_tm1)

	di in gr "theta =  0" _col(17) %10.0g in ye e(ll_t0) /*
		*/ _col(32) in ye %8.2f e(chi2_t0) /*
		*/ _col(51)  %5.3f in ye e(p_t0)

	di in gr "theta =  1" _col(17)  %10.0g in ye e(ll_t1) /*
		*/ _col(32) in ye %8.2f e(chi2_t1) /*
		*/ _col(51)  %5.3f in ye e(p_t1)

	di _dup(57) in gr "-"
end

program define lrdispT
	di " "
	if "`e(model)'" == "lambda" {
		di _dup(57) in gr "-"
		di in gr _col(4) "Test" _col(17) "Restricted" _col(32) /*
			*/ "LR statistic" /*
			*/  _col(50) "P-Value"
		di  in gr _col(5) "H0:"  _col(15) "log likelihood" /*
			*/  _col(35) "X~chi2" _col(49) "Pr > chi2"
		di _dup(57) in gr "-"

		di in gr "lambda = -1" _col(17)  %10.0g in ye e(ll_tm1) /*
			*/ _col(32)  in ye %8.2f e(chi2_tm1) /*
			*/  _col(51) %5.3f in ye e(p_tm1)

		di in gr "lambda =  0" _col(17) %10.0g in ye e(ll_t0) /*
			*/ _col(32) in ye %8.2f e(chi2_t0) /*
			*/ _col(51)  %5.3f in ye e(p_t0)

		di in gr "lambda =  1" _col(17)  %10.0g in ye e(ll_t1) /*
			*/ _col(32) in ye %8.2f e(chi2_t1) /*
			*/ _col(51)  %5.3f in ye e(p_t1)

		di _dup(57) in gr "-"
	}
	else {
		di _dup(63) in gr "-"
		di in gr _col(4) "Test" _col(23) "Restricted" _col(38)
			*/ "LR statistic"  _col(56) "P-Value"
		di  in gr _col(5) "H0:"  _col(21) "log likelihood" /*
			*/ _col(41) "X~chi2" _col(55) "Pr > chi2"

		di _dup(63) in gr "-"

		di in gr "theta=lambda = -1" _col(23)  %10.0g in ye /*
			*/ e(ll_tm1) _col(38)  in ye %8.2f e(chi2_tm1) /*
			*/  _col(57) %5.3f in ye e(p_tm1)

		di in gr "theta=lambda =  0" _col(23) %10.0g in ye e(ll_t0) /*
			*/ _col(38) in ye %8.2f e(chi2_t0) /*
			*/ _col(57)  %5.3f in ye e(p_t0)

		di in gr "theta=lambda =  1" _col(23)  %10.0g in ye e(ll_t1)/*
			*/ _col(38) in ye %8.2f e(chi2_t1) /*
			*/ _col(57)  %5.3f in ye e(p_t1)

		di _dup(63) in gr "-"
	}

end


program define ChkVar
	args var touse
	qui count if `var'<=0 & `touse'
	if r(N) {
		di in red /*
		*/ "`1' contains observations that are not strictly positive"
		exit 411
	}
end


program define Display
	syntax [, Level(integer $S_level)]
	local level "level(`level')"

	bhead
	btitle
	if "`e(model)'"=="lambda" |"`e(model)'"=="rhsonly"   {
		_diparm lambda, `level'
	}
	if  "`e(model)'"=="theta" {
		_diparm lambda, `level'
		_diparm theta, `level'
	}
	if  "`e(model)'"=="lhsonly" {
		_diparm theta, `level'
	}


	if "`e(lrtest)'" != "" {
		cdisplr
	}
	else {
		cdisp
	}

	if "`e(model)'"=="rhsonly"   {
		lrdispR
	}
	if  "`e(model)'"=="lhsonly" {
		lrdispL
	}
	if "`e(model)'"=="lambda" |  "`e(model)'"=="theta" {
		lrdispT
	}
end

exit

Use of globals
--------------

	$T_bvec		name of matrix containing coefficient vector
	$T_wtexp 	"[`weight'=`wvar']"
	$T_nocns	"" or "noconstant"
	$T_dep		name of dependent variable
	$T_parm		list of var names:  transformed variables
	$T_notr		list of var names:  nontransformed variables
	$T_sig		name of scalar containing value of sigma