boxcox.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 1.2.3  17mar2005
program define boxcox , eclass byable(recall)
	version 7.0, missing
	if _caller()<7 {
		if _by() { error 190 }
		boxcox_6 `0'
		exit
	}
						
	if replay() {
		if `"`e(cmd)'"' != "boxcox"  { error 301 }
		if _by() { error 190 }
		Display `0'
		exit
	}

	syntax varlist(ts) [if] [in] [fweight iweight] /*
		*/ [,  noLOg ITERate( int 50000 ) MODEL(string)/*
		*/ noConstant LRTEST noLOGLR FROM(string) /*
		*/ NOTRans(varlist ts) GRaph SAVing(string) MEAN MEDian /*
		*/ Level(cilevel) GRVars LStart(real 1.0) /*
		*/ DELta(real 0.0 ) Zero(real 0.0 ) Generate(string)  * ] 

	tokenize `varlist'
	macro shift

	if "`*'`notrans'" == "" & "`constant'" != "" {
		di as err "no covariates specified"
		exit 198
	}	
	

	marksample touse 
	markout `touse' `notrans'

	tempname sig bvec

	global T_bvec `bvec'
	
	if "`loglr'" != "" {
		local loglr "nolog"
	}

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

	mlopts mlopts, `options' `iterate'
	local cns `s(constraints)'
	if "`cns'" != "" {
		local lrtest
	}

/* Old options */

	if `zero' != 0.0 {
		di as txt "zero() option is no longer supported --" /*
			*/ " set version < 7"
	}

	if `delta' != 0.0 {
		di as txt "delta() option is no longer supported -- " /*
			*/ "set version < 7"
	}
	
	if "`generate'" != "" {
		di as txt "generate() option is no longer supported -- " /*
			*/ "set version < 7"
	}
	
	if "`grvars'" != "" {
		di as txt "grvars option is no longer supported -- " /*
			*/ "set version < 7"
	}
	if "`graph'" != "" {
		di as txt "graph option is no longer supported -- " /*
			*/ "set version < 7"
	}
	if "`saving'" != "" {
		di as txt "saving() option is no longer supported -- " /*
			*/ "set version < 7"
	}
	if "`mean'" != "" {
		di as txt "mean option is no longer supported -- " /*
			*/ "set version < 7"
	}
	if "`median'" != "" {
		di as txt "median option is no longer supported -- " /*
			*/ "set version < 7"
	}



/* 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']"  */
		global T_wtexp "[`weight'=`wvar']"
	}
	else {
		global T_wtexp 
	}

      	global T_nocns "`constant'"


/* 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 as err "Model not identified"
			di as err "No transformed variables leaves lambda " /*
			*/ "unidentified"
			di as err "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)
						/* -summarize- does not allow
						   iweights, so use -regress-
						   with constant-only model 
						*/
				reg `err1' $T_wtexp if `touse'
				scalar `sig1' = _b[_cons]

				gen double `llft'=-.5*(1+ln(2*_pi)+ln(`sig1'))
				*summ `llft' $T_wtexp if `touse', meanonly
				*scalar `LL0'=r(sum)
				reg `llft' $T_wtexp if `touse'
				scalar `LL0' = _b[_cons]*e(N)
			}
			local  df0 2
		}

		if "`from'" != "" {
			local start "init(`from')"
		}
		else {
			if "`lstart'" != "" {
				local start "init(`lstart', copy)"
			}
		}
		di as txt _n "Fitting full model"
	
		#delimit ;
		ml model rdu0  "boxco_l RHS" (lambda: `dep' = ) 
			if  `touse'  $T_wtexp, technique(nr) search(off)
			`start' `mlopts' `log'
			max  ;
		#delimit cr
		tempname llf bf Vf bpar lsig
		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
		boxco_l RHS `lltmp' `bc' `lltmp'
	
		matrix `bf'=$T_bvec
		scalar `lsig'=$T_sig
		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'=1-chi2(1, `chit1')

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

		if "`lrtest'" != ""	{
			di as txt _n "Fitting 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  "boxco_l RHS" 
					(lambda: `dep' = ) 
					if `touse'  $T_wtexp, technique(nr) 
					search(off) `mlopts' 
					`start'  max `loglr' ;
				#delimit cr
		
				scalar `lltr'=2*(`llf'-e(ll))	
				scalar `pv' = 1-chi2(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 as txt _n /*
					   */ "Comparison model for LR " /*
					   */ "test on `1' is a linear " /*
					   */ "regression"	
					di as txt "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)
reg `err' if `touse' $T_wtexp
scalar `sigmat' = _b[_cons]

gen double `llf2'=-.5*(1 + ln(2*_pi)+ ln(`sigmat') ) if `touse'
*summ `llf2' $T_wtexp if `touse'
*scalar `llt'  =r(sum)
reg `llf2' $T_wtexp if `touse'
scalar `llt'  = _b[_cons]*e(N)
est scalar ll = `llt'

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

				scalar `lltr'=2*(`llf'-e(ll))

				scalar `pv' = 1-chi2(`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' `lsig'

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

		if "$T_nocns" == "" & "`cns'" == "" {
			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'" */
		est local wtype="`weight'"
		if "`wt'" != "" {
			est local wexp="`exp'"
		}
		est local depvar  "`dep'"
		est local lrtest "`lrtest'"
		est local chi2type "LR"
		est local predict "boxco_p"
		est local model  "`model'"
		if "`ntrans'"=="" & "$T_nocns" != "" {
			est local ntrans 
		}
		else 	est local ntrans "yes"

		est local cmd "boxcox"
		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 as txt "Fitting comparison model"

		tempname bvec1 LL0
		matrix `bvec1'=1
	
		ml model rdu0 "boxco_l ConsOnly" ( `dep' = ) /*
		    */    if  `touse'  $T_wtexp, technique(nr)  search(off) /*
		    */  init(1, copy)   max noheader `log'

		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'" != "" {
			local start "init(`from')"
		}
		else {
			if "`lstart'" != "" {
				local start "init(`lstart', copy)"
			}
		}
	
		di as txt _n "Fitting full model"

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

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

		tempname llf lsig
		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
		boxco_l LHS `lltmp' `bc' `lltmp'
	
		matrix `bf'=$T_bvec
		macro drop ML_y1 
		macro drop ML_samp
		scalar `lsig' = $T_sig
		
/* 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'=1-chi2(1, `chit1')

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

/* Do LR tests if requested */

		if "`lrtest'" != "" {
			tempname p chi2 df
			di _n as txt "Fitting comparison models for LR tests"

			local firstlr 1

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

			tempname lltr pv
	
			tokenize $T_notr