scobit.ado

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

*! version 1.5.6  31jan2005
program define scobit, eclass byable(onecall) prop(ml_score or swml)
	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}
	`BY' _vce_parserun scobit : `0'
	if "`s(exit)'" != "" {
		exit
	}

	version 6, missing
	if replay() {
		if "`e(cmd)'"!="scobit" { error 301 }
		if _by() { error 190 }
		Display `0'
		error `e(rc)'
		exit
	}
	`BY' Estimate `0'
end

program define Estimate, eclass byable(recall)
	syntax varlist(numeric) [if] [in] [fw iw pw] [, OR ASIS FROM(string) /*
	*/ Level(cilevel) noCONstant Robust CLuster(varname) /*
	*/ OFFset(varname numeric) SCore(string) noLOg * ]
	if _by() { 
		_byoptnotallowed score() `"`score'"'
	}

	mlopts mlopts, `options'
	local cns `s(constraints)'

/* Check syntax. */

	if `"`score'"'!="" {
		local nword : word count `score'
		if `nword'==1 & substr("`score'",-1,1)=="*" { 
			local score = /*
			*/ substr("`score'",1,length("`score'")-1)
			local score `score'1 `score'2
			local nword 2
		}
		confirm new variable `score'
		local nword : word count `score'
		if `nword' != 2 {
			di in red "score() must contain the name of " /*
			*/ "two new variables"
			exit 198
		}
		local scname1 : word 1 of `score'
		local scname2 : word 2 of `score'
		tempvar scvar1 scvar2
		local scopt "score(`scvar1' `scvar2')"
	}
	if "`constan'"!="" {
		local nvar : word count `varlist'"
		if `nvar' == 1 {
			di in red "independent variables required with " /*
			*/ "noconstant option"
			exit 102
		}
	}

/* Mark sample. */

	marksample touse

	if `"`cluster'"'!="" {
		markout `touse' `cluster', strok
		local clopt cluster(`cluster')
	}
	if "`offset'"!="" {
		markout `touse' `offset'
		local offopt "offset(`offset')"
	}
	if "`weight'" == "pweight" | "`robust'`cluster'" != "" {
		local robust "robust"
		/* probit does not allow iweights and robust so use crittype
		 * instead of robust option */
		local crtype crittype("log pseudolikelihood")
	}

/* Count obs and check values of `y'. */

	gettoken y xvars : varlist

	qui count if `touse'
	local N `r(N)'

	if `N' == 0 { error 2000 }
	if `N' == 1 { error 2001 }

	qui count if `y'==0 & `touse'
	if r(N)==0 | r(N)==`N' {
		di in red "outcome does not vary; remember:"
		di in red _col(35) "0 = negative outcome,"
		di in red _col(9) /*
		*/ "all other nonmissing values = positive outcome"
                exit 2000
        }


	if "`log'"!="" {
		local qui "quietly"
	}

/* Run logistic to drop any variables and observations. */

	if "`asis'"=="" {
		`qui' logit `y' `xvars' [`weight'`exp'] if `touse', /*
		*/ iter(0) `offopt' nocoef nolog `constan' `crtype'
		qui replace `touse' = e(sample)
		_evlist
		local xvars `s(varlist)'
		tempname lb
	}

/* Estimate logistic model. */

	if "`from'"=="" {
		if "`robust'`cns'`cluster'"!="" | "`weight'"=="pweight" {
			local nowarn  "nowarn"
		}
		else {
			tempname llp
		}
		`qui' di _n in gr "Fitting logistic model:"

		`qui' logit `y' `xvars' [`weight'`exp'] if `touse', /*
		*/ `offopt' nocoef asis `constan' `crtype'

		if "`llp'"!="" { scalar `llp' = e(ll) }

		tempname from
		mat `from' = e(b)
		mat coleq `from' = `y'
	}

/* Estimate full model. */

	`qui' di _n in gr "Fitting full model:"

	ml model d2 scob_lf (`y': `y' = `xvars', `constan' `offopt') /lnalpha /*
	*/ [`weight'`exp'] if `touse', collinear missing max nooutput /*
	*/ nopreserve init(`from') search(off) `robust' `clopt' `scopt' /*
	*/ `log' waldtest(0) `mlopts' title(Skewed logistic regression)

	est local cmd

	qui count if `y'==0 & `touse'
	est scalar N_f = r(N)
	est scalar N_s = e(N) - e(N_f)
	est scalar alpha = exp([lnalpha]_b[_cons])

        if "`score'"!="" {
		label var `scvar1' "Score index for x*b from scobit"
		rename `scvar1' `scname1'
		label var `scvar2' /*
			*/ "Score index for /lnalpha from scobit"
		rename `scvar2' `scname2'
		est local scorevars `scname1' `scname2'
	}

	if "`llp'"!="" {
		est local chi2_ct "LR"
		est scalar ll_c = `llp'
		if e(ll) < e(ll_c) & reldif(e(ll),e(ll_c)) < 1e-5 {
			est scalar chi2_c = 0
				/* otherwise, let it be negative when
				   it does not converge
				*/
		}
		else	est scalar chi2_c = 2*(e(ll)-e(ll_c))
	}

	est local offset1
	est local offset "`offset'"

	est local predict "scob_p"
        est local cmd     "scobit"

	Display, `or' level(`level') `nowarn'
	error `e(rc)'
end

program define Display
	syntax [, Level(cilevel) OR NOwarn]

	local crtype = upper(substr(`"`e(crittype)'"',1,1)) + /*
		*/ substr(`"`e(crittype)'"',2,.)

        di in gr _n "`e(title)'" _col(49) "Number of obs     =" /*
	*/ in ye _col(70) %9.0g e(N) _n /*
	*/ in gr _col(49) "Zero outcomes     =" /*
        */ in ye _col(70) %9.0g e(N_f) _n /*
        */ in gr "`crtype' = " in ye %9.0g e(ll) /*
        */ in gr _col(49) "Nonzero outcomes  =" /*
        */ in ye _col(70) %9.0g e(N_s) _n

	ml di, noheader level(`level') `or' first plus
	_diparm lnalpha, level(`level')
        di in smcl in gr "{hline 13}{c +}{hline 64}"
	_diparm lnalpha, level(`level') exp label("alpha")
        di in smcl in gr "{hline 13}{c BT}{hline 64}"

	if "`e(chi2_ct)'"=="LR" {
		if e(chi2_c) < 1e5 {
			local fmt "%9.2f"
		}
		else 	local fmt "%9.2e"
		di in gr "Likelihood-ratio test of alpha=1:   " /*
		*/ in gr "chi2(" in ye "1" in gr ") =" in ye `fmt' /*
		*/ e(chi2_c) in gr "    Prob > chi2 = " in ye %6.4f /*
		*/ chiprob(1, e(chi2_c))
	}

	if `"`nowarn'"'=="" {
		di _n in gr /*
		*/ "note: Likelihood-ratio tests are recommended for " /*
		*/ "inference with scobit models."
	}
	_prefix_footnote
end