gamma.ado

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

*! version 1.4.12  21jun2005
program define gamma, eclass byable(recall) prop(ml_score)
	version 7.0, missing
	if replay() {
		if `"`e(cmd)'"' != "gamma" { error 301 } 
		if _by() { error 190 } 
		syntax [, Level(cilevel) noCOEF noHEADer TR ]
		if "`e(sigma)'"=="" {
			local ancilla 1
		}
	}
	else {
		syntax [varlist] [if] [in] [fweight pweight iweight] /*
			*/ [, noCOEF CLuster(string) Dead(string)/*
			*/ DEBUG FROM(string) noHEADer MLMethod(string) /*
			*/ MLOpt(string) SHApe(real 1 ) noCONstant /*
			*/ Level(cilevel) noLOg OFFset(string) /*
			*/ Robust SCore(string) T0(string) SKIP TR /*
			*/ ANCillary(varlist) ANC2(varlist) STrata(varname) *]

		if _by() {
			_byoptnotallowed score() `"`score'"'
		}


		tokenize `varlist'
		local t `1'
		mac shift 
		local rhs `*'
		if "`strata'"~="" {
			if "`ancillary'"~="" | "`anc2'"~="" {
				noi di as err /*
				*/ "options strata() and ancillary()" /*
				*/ " may not be specified together"
				exit 198
			}
			qui xi, prefix(_S) i.`strata'
			local rhs `rhs' _S*
			local ancillary _S*
			local anc2 _S*
		}

		unab t0 : `t0', min(0) max(1) name(t0())
		unab dead: `dead', min(0) max(1) name(dead())
		unab offset: `offset', min(0) max(1) name(offset())

		if "`cluster'"!="" {
			unab cluster:  `cluster', max(1) name(cluster())
			local cluopt cluster(`cluster')
		}
		if "`from'" != "" { local iniopt init(`from') }
		if "`mlmethod'" == "" { local mlmetho = "d2" }
		if "`offset'" !="" { local offopt = "offset(`offset')" }
		mlopts options, `options'

		if "`score'" != "" { 
			local n : word count `score'
			if `n'==1 & substr("`score'",-1,1)=="*" { 
				local score = /*
				*/ substr("`score'",1,length("`score'")-1)
				local score `score'1 `score'2 `score'3 
				local n 3
			}
			if `n' != 3 { 
				di as err /*
			*/ "score() invalid:  three new variable names required"
				exit 198 
			}
			confirm new var `score'
			local scopt "score(`score')"
		}

		if "`weight'" != "" { 
			tempvar wv
			qui gen double `wv' `exp'
			local w [`weight'=`wv']
		}

		tempvar touse 
		mark `touse' `w' `if' `in'
		markout `touse' `t' `rhs' `dead' `t0' `offset'
		markout `touse' `cluster', strok

/*
		preserve
		quietly keep if `touse'
*/

		if "`dead'" != "" {
			local sdead "`dead'"
			capture assert `dead'==0 | `dead'==1 if `touse'
			if _rc { 
				tempvar mydead 
				qui gen byte `mydead' = `dead'!=0 if `touse'
				local dead "`mydead'"
			}
		}
		else {
			tempvar dead 
			qui gen byte `dead'=1
			local sdead 1
		}

		if "`t0'" == "" {
			local t0 0
		}
		capture assert `t0' < `t' if `touse'
		if _rc {
			di as err "`t0' >= `t' in some obs."
			exit 498
		}


		_rmcoll `rhs' `w' if `touse', `constant'
		local rhs "$S_1"
		global S_1

		global EREGd `dead'
		global EREGt0 `t0'
		*global S_K = `shape'  /* Keep for reference if 2 parameters */

		if "`log'"!="" { local nlog="*" }
		tempvar mysamp num den
		tempname b f V  g b0 bc0
		quietly { 
			tempvar swt
			if "`weight'"=="aweight" | "`weight'"=="pweight" {
				tempvar wvn
				summ `wv' if `touse', meanonly 
				gen double `wvn' = `wv'/r(mean)
				local wvngen 1
				qui sum `wvn', meanonly 
				gen double `swt'=r(sum)
			}
			else if "`weight'"!="" {
				local wvn `wv' 
				qui sum `wvn', meanonly 
				gen double `swt'=r(sum)
			}
			else {
				local wvn 1
				count if `touse'
				gen `swt'=r(N)
			}
			gen double `num' = `wvn'*ln(`t') if `touse' & `dead'
			replace `num' = sum(`num')
			global EREGa = `num'[_N]/`swt'
			drop `num' `sw'
			if "`wvngen'"!=""  { drop `wvn' }
		}
		if "`constant'"!="" {
			local skip = "skip"
			`nlog' di as txt "Fitting full model:"
		}
		if "`rhs'" != "" & "`skip'"=="" { 
			summ `dead' if `touse', meanonly
			local mysumd = r(sum)
			summ `t' if `touse', meanonly
			local mysumt = r(sum)
			local cval = -ln(`mysumd'/`mysumt')
			if "`offset'"!="" {
				summ `offset' if `touse', meanonly
				local cval = `cval' - r(mean)
			}
			`nlog' di ""
			`nlog' di as txt "Fitting constant-only model:"
			ml model d2 gamma_d2 (`t': `t'=, `offopt') /*
			*/ (ln_sig: `ancillary') (kappa: `anc2') `w' /*
			*/ if `touse', /*
			*/ init(_cons=`cval' /ln_sig=0 /kappa=1) /*
			*/ missing collin nopreserve wald(0) `mlopt' /*
			*/ max noout `log' `options' search(off) `robust'
			local cont continue
			`nlog' di ""
			`nlog' di as txt "Fitting full model:"
		}
		else {
			
			local cont wald(1)
			if "`iniopt'"=="" {
				if "`constant'"=="" { 
					summ `dead' if `touse', meanonly
					local mysumd = r(sum)
					summ `t' if `touse', meanonly
					local mysumt = r(sum)
					local cval = -ln(`mysumd'/`mysumt')
					local iniopt= /*
				      */ "init(_cons=`cval' /ln_sig=1 /kappa=1)"
				}
				else {
					local iniopt= /*
				      */ "init( /ln_sig=1 /kappa=1)"
				}
			}
		}
		ml model d2  gamma_d2 /*
			*/ (`t': `t'=`rhs' , `offopt' `constant' )  /*
			*/ (ln_sig: `ancillary') (kappa: `anc2') `w' /*
			*/ if `touse', `cont' noout /*
			*/ `robust' `cluopt' `scopt' `iniopt' `mlopt' /*
			*/ missing collin nopreserve /*
			*/ max search(off) `log' `options' 
		if "`e(wtype)'" != "" { 
			est local wexp "`exp'"
		}
		est local title2 "accelerated failure-time form"
		est local predict gamma_p
		est local cmd gamma 
		est local t0 `t0'
		est local dead `sdead'

		global S_E_cmd gamma
	}
	global EREGw
	global EREGd
	global EREGt
	global EREGt0
	global EREGa
	if "`ancillary'" =="" & "`ancilla'" =="" {
		local bg "[ln_sig]_b[_cons]"
		local bb "[kappa]_b[_cons]"
		est scalar sigma = exp(`bg')
		est scalar kappa = `bb'
		est local stcurve="stcurve"
	}

	if "`coef'"=="" {
		if `"`tr'"'!=`""' {
			local hr `"eform(Tm. Ratio)"'
                }

		if "`ancillary'" =="" & "`anc2'"=="" & "`ancilla'" ==""{
			if "`header'" == "" {
		     	    di _n as txt "Gamma regression -- entry time `e(t0)'"
			}
			version 9: ///
			ml di, `header' `hr' level(`level') first plus ///
				title(`e(title2)')
			local bg "[ln_sig]_b[_cons]"
			local sg "[ln_sig]_se[_cons]"
			local iz = invnorm(1-(1-`level'/100)/2)
			local ll = `bg'-`sg'*`iz'
			local uu = `bg'+`sg'*`iz'
			local zz = `bg'/`sg'
			local se = sqrt(exp(`bg')*exp(`bg')*`sg'*`sg')
			local bb "[kappa]_b[_cons]"
			local sb "[kappa]_se[_cons]"
			local bl = `bb'-`sb'*`iz'
			local bu = `bb'+`sb'*`iz'
			local bz = `bb'/`sb'
			local be = sqrt(exp(`bb')*exp(`bb')*`sb'*`sb')
			local ig = 1/exp(`bg')
			local ib = 1/exp(`bb')
			
			est scalar sigma = exp(`bg')
			est scalar kappa = `bb'
		
		di as txt in smcl  "     /ln_sig {c |}  "  as res %9.0g `bg' /*
			*/ "  " %9.0g (`sg') " " %8.2f (`zz') /*
			*/ " " %7.3f (normprob(-abs(`zz'))*2) /*
			*/ _col(58) %9.0g (`ll') "   " %9.0g (`uu')
		di as txt in smcl  "      /kappa {c |}  "  as res %9.0g `bb' /*
			*/ "  " %9.0g (`sb') " " %8.2f (`bz') /*
			*/ " " %7.3f (normprob(-abs(`bz'))*2) /*
			*/ _col(58) %9.0g (`bl') "   " %9.0g (`bu')
		di as txt in smcl "{hline 13}{c +}{hline 64}"
		di as txt in smcl  "       sigma {c |}  "  /*
			*/ as res %9.0g exp(`bg') /*
			*/ "  " %9.0g (`se') " "  /*
			*/ _col(58) %9.0g exp(`ll') "   " %9.0g exp(`uu')
		di as txt in smcl "{hline 13}{c BT}{hline 64}"
			ml_footnote
		}
		else {
			if "`header'" == "" {
		     	    di _n as txt "Gamma regression -- entry time `e(t0)'"
			}
			version 9: ml di, `header' `hr' level(`level') 

		}
	}
end