hetprob.ado

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

*! version 3.10.1  20dec2004
program define hetprob, eclass byable(onecall) prop(ml_score)
	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}
	`BY' _vce_parserun hetprob, jkopts(eclass) : `0'
	if "`s(exit)'" != "" {
		exit
	}

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

program define Estimate, eclass byable(recall)
	syntax varlist(ts) [if] [in] [fweight pweight iweight] , /*
		*/ het(string)  [Robust CLuster(varname numeric) /*
		*/ SCore(string) OFFset(varname numeric) /*
		*/ FROM(string) noCONstant noLOg noSKIP noLRtest /*
		*/ Level(cilevel) MLMethod(string) ASIS DIFficult * ]
	if _by() {
		_byoptnotallowed score() `"`score'"'
	}


	/* force difficult option in full model */
	local diff "difficult"
	local p_vars `varlist'
	mlopts mlopt, `options'
	local cns `s(constraints)'
	local p_off `offset'
	local p_cons `constan' 
	if "`mlmethod'" == "" {
		local mlmetho d2
	}

	marksample touse

	if "`weight'" != "" { local wgt [`weight'`exp'] }
	if "`cluster'" != "" { local clopt "cluster(`cluster')" }
	if "`score'" != "" {
		local wcount : word count `score'
		if `wcount'==1 & substr("`score'",-1,1)=="*" { 
			local score = /*
			*/ substr("`score'",1,length("`score'")-1)
			local score `score'1 `score'2 
			local wcount 2
		}
		if `wcount' != 2 {
			di in red "score() requires two new variable names"
			exit 198
		}
		local s1 : word 1 of `score'
		local s2 : word 2 of `score'
		confirm new var `s1' `s2'
		tempvar sc1 sc2
		local scopt "score(`sc1' `sc2')"
	}

	local 0 `het'
	syntax varlist(numeric min=1 ts) [, OFFset(varname numeric) /*
			*/	noCONstant *]
	local heti `varlist'
	local hetf `offset'
	local offset ""

	markout `touse' `cluster' `offset' `heti' `hetf', strok
	_rmcoll `heti' if `touse', constant
	local heti "`r(varlist)'"

/* Check hetero arguments */
	if "`constan'`options'" != "" {
		noi di in red "Heteroskedastic option(s) not permitted: " _c
		noi di in red "`constan' `options'"
		exit 101
	}

/* process hetero offset */
	if "`hetf'" != "" {
		local hetoff "offset(`hetf')"
		markout `touse' `hetf'
	}

/* Check predictor equation */
	tokenize `p_vars'
	local dep `1'
	tsunab dep : `dep'
	local depn : subinstr local dep "." "_"
	mac shift
	local ind `*'
	if "`p_cons'" == "noconstant" & "`ind'" == "" {
		noi di in red "independent varlist required with " _c
		noi di in red "noconstant option."
		exit 100
	}

/* process predictor offset */
	if "`p_off'" != "" {
		local offset "offset(`p_off')"
		markout `touse' `p_off'
	}


/* check independent variables for collinearity */
	noi _rmcoll `ind' if `touse', `p_cons'
	local ind "`r(varlist)'"
 
	if "`log'" != "" { 
		local log quietly 
	}
	else {
		local log noisily
	}

	if "`asis'" == "" {
		`log' di
		`log' probit `dep' `ind' `wgt' if `touse', iter(0) nocoef nolog
		qui replace `touse' = e(sample)
		tempname lb
		mat `lb' = e(b)
		local ind : colnames(`lb')
		tokenize `ind'
		local last : word count `ind'
		if "``last''" == "_cons" { local `last' }
		local ind "`*'"
	}
	else {
		`log' di
	}

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

/* output options */
	if "`robust'`cluster'" != "" | "`weight'" == "pweight" {
		local robust "robust"
	}

	if "`ind'" == "" | "`robust'`p_cons'" != "" {
		local skip ""
	}

	if "`robust'" != "" {
		local lrtest nolrtest
	}


	if "`from'" != "" {
		local init "init(`init')"
	}

/* code for LR test of full model with hetero vs. full model no hetero */
	local lr_OK 0
	if "`lrtest'" == "" {
		tempname ll_c chi2_c
		* compute full probit without heteroskedastic equation

		`log' di in gr "Fitting probit model:"
		`log' capture `log' probit `dep' `ind' `wgt' if `touse', /*
				*/ `p_cons' `offset' nocoef asis
		if _rc == 0 {
			scalar `ll_c' = e(ll)        
			scalar `chi2_c' = e(chi2)    
			if "`init'" == "" & `ll_c' < . { 
				local lr_OK 1
				tempname myb 
				mat `myb' = get(_b)
				mat coleq `myb' = `depn'
				local init "init(`myb')"
			}
		}
	}

/* constant only model */
	if "`skip'" != "" {
		local continu "continue"

		/* starting values from nonhetero probit */
		tempname p_1 p_2 Ival
		quietly probit `dep' if `touse' `wgt', `offset'
		local fulll = e(ll)
		capture matrix `p_1' = e(b)
		if _rc {
			local confrom ""
		}
		else {
			matrix coleq `p_1' = `dep':
			local rcount : word count `heti'
			matrix `p_2' = J(1,`rcount',0)
			mat colnames `p_2' = `heti'
			mat coleq `p_2' = lnsigma2
			mat `Ival' = `p_1', `p_2'
			local confrom "init(`Ival', copy)"
		}
        
		`log' di _n in gr "Fitting constant-only model:"
		#delimit ;
		`log' ml model `mlmethod' hetpr_lf 
			(`depn': `dep' = , `offset') 
			(lnsigma2: = `heti', nocons `hetoff') 
			if `touse' `wgt', 
			`confrom' max missing nopreserve wald(0) `mlopt'
			nooutput search(off) `diff' ;
		#delimit cr

		/* use info for good starting values */
		if "`init'" == "" { 
			tempname myh
			mat `myh' = e(b)
			local init "init(`myh', copy)"
		}

	}
	else {
		local continu "wald(1)"
	}

	if "`from'" != "" {
		local init "init(`from')"
	}


	`log' display _n in gr "Fitting full model:"
	#delimit ;
	`log' ml model `mlmethod' hetpr_lf 
		(`depn': `dep' = `ind', `p_cons' `offset') 
		(lnsigma2: `heti', nocons `hetoff') if `touse' `wgt', 
		max missing nopreserve `mlopt'
		`continu' search(off)
		`clopt' `scopt' `robust' nooutput `init'
		title(Heteroskedastic probit model) `diff' ;
	#delimit cr
	* if "`log'" == "quietly" { di _n }

	if "`scopt'" != "" {
		rename `sc1' `s1'
		if "`s2'" != "" {
			rename `sc2' `s2'
		}
		est local scorevars `s1' `s2'
	}

	if "`lrtest'`cns'" == "" {
		est scalar ll_c =  `ll_c'       
		est scalar chi2_c = 2*(e(ll) - e(ll_c))       
		local fulldof : word count `heti'
		est scalar df_m_c = `fulldof'
		est scalar p_c = chiprob(e(df_m_c),e(chi2_c))
		est local chi2_ct "LR"
	}
	else {
		qui test [lnsigma2]
		est scalar chi2_c = r(chi2)
		est scalar df_m_c = r(df)
		est scalar p_c    = r(p)
		est local chi2_ct "Wald"
	}

	est scalar N_s = `npos'
	est scalar N_f = `nneg'

	est local method "ml"
	est local predict "hetpr_p"
	est local cmd "hetprob"

	Display, level(`level') `lrtest'
end

program define Display
	syntax [, Level(cilevel) noLRtest ]

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

	di in gr _n "Heteroskedastic probit model" _col(49) /*
		*/ "Number of obs     =" in ye _col(70) %9.0g e(N)
	di in gr _col(49) "Zero outcomes     =" /*
		*/ in ye _col(70) %9.0g e(N_f) 
	di in gr _col(49) "Nonzero outcomes  =" /*
		*/ in ye _col(70) %9.0g e(N_s) _n
	if !missing(e(df_r)) {
		di in gr _col(49) "F(" %4.0f in ye `e(df_m)' /*
			*/ in gr "," in ye %7.0f e(df_r) /*
			*/ in gr ")" _col(67) "=" in ye _col(70) %9.2f e(F) 
		di in gr "`crtype' = " in ye %9.0g e(ll) /*
			*/ in gr _col(49) "Prob > F" _col(67) "=" /*
			*/ in ye _col(70) /*
			*/ in ye %9.4f Ftail(e(df_m),e(df_r),e(F)) _n
	}
	else {
		di in gr _col(49) "`e(chi2type)' chi2(" in ye `e(df_m)' /*
			*/ in gr ")" _col(67) "=" in ye _col(70) %9.2f e(chi2) 
		di in gr "`crtype' = " in ye %9.0g e(ll) /*
			*/ in gr _col(49) "Prob > chi2" _col(67) "=" /*
			*/ in ye _col(70) /*
			*/ in ye %9.4f chiprob(e(df_m),e(chi2)) _n
	}

	ml di, noheader level(`level')
	if "`e(chi2_ct)'" == "LR" {
                di in green "Likelihood-ratio test of lnsigma2=0: " _c
	}
	else {
                di in green "Wald test of lnsigma2=0:             " _c
	}
	di in green "chi2(" in ye "`e(df_m_c)'" in gr ") = " /*
		*/ in ye %8.2f e(chi2_c) _c
	di in green "   Prob > chi2 = " in ye %6.4f /*
		*/ chiprob(e(df_m_c),e(chi2_c))
end
exit

Note:  The above could be coupled with the _crcshdr.ado code if we save some
       additional words -- right now that program uses "Censored" and 
       "Uncensored" instead of what we want here "Positive" and "Negative".