heckman.ado

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

			*/ _col(59) in gr "Prob > chi2 = " in ye %6.4f e(p_c)
		di in smcl in gr "{hline 78}"
		_prefix_footnote
		exit e(rc)
	}
end

program define DiLambda
	args level

	tempname z lb ub

	scalar `z' = invnorm((100 + `level') / 200)
	scalar `lb' = e(lambda) - `z'*e(selambda)
	scalar `ub' = e(lambda) + `z'*e(selambda)

	local llb : di %9.0g `lb'
	if length("`llb'") > 9 {
		local lbfmt "%8.0g"
	}
	else	local lbfmt "%9.0g"

	local uub : di %9.0g `ub'
	if length("`uub'") > 9 {
		local ubfmt "%8.0g"
	}
	else	local ubfmt "%9.0g"

	di in smcl in gr _col(3) "    lambda {c |}  " /*
		*/  in ye %9.0g e(lambda) "  " %9.0g e(selambda)    /*
		*/  _col(58) `lbfmt' `lb' "   "  `ubfmt' `ub'

end


program define Display2
	syntax [, Level(cilevel) ]

	_crcshdr
	est display, level(`level') plus
        di in smcl in gr "         rho {c |} " in ye %10.5f e(rho)
        di in smcl in gr "       sigma {c |} " in ye %10.0g e(sigma)
	di in smcl in gr "      lambda {c |} " in ye %10.0g e(lambda) /*
						*/ "  " %9.0g e(selambda) 
	di in smcl in gr "{hline 13}{c BT}{hline 64}"
	_prefix_footnote
end

program define Reparm, eclass

			/* somewhat superseded by _diparm, but kept for
			 * lambda and so sigma and rho can be saved in e() */

	tempname b v d tau lns rho sig tmp lambda lamse
	mat `b' = get(_b)
	mat `v' = get(VCE)

	mat `tmp' = `b'[1,"athrho:_cons"]
	scalar `tau' = `tmp'[1,1]
	mat `tmp' = `b'[1,"lnsigma:_cons"]
	scalar `lns' = `tmp'[1,1]

	scalar `rho' = (exp(2*`tau')-1) / (exp(2*`tau')+1)
	scalar `sig' = exp(`lns')
	scalar `lambda' = `rho'*`sig'

	mat `d' =  ( `sig'*4*exp(2*`tau')/(exp(2*`tau')+1)^2 , `lambda' )
	mat `v' =  (`v'["athrho:_cons".."lnsigma:_cons",      /*
		*/   "athrho:_cons".."lnsigma:_cons"]  )
	mat `v' = `d'*`v'*`d''
	scalar `lamse' = sqrt(`v'[1,1])

	est scalar rho    = `rho'
	est scalar sigma  = `sig'
	est scalar lambda = `lambda'
	est scalar selambda  = `lamse'
	/*  Double saves for backward compatibility */
	global S_1 = `rho'
	global S_2 = `sig'
	global S_3 = `lambda'
	global S_4 = `lamse'
end


program define GetIval, rclass 
	args 	    dep		/*
		*/  ind		/*
		*/  nc		/*
		*/  selind	/*
		*/  selnc	/*
		*/  b0		/*
		*/  wgt		/*
		*/  wtype	/*
		*/  wtexp	/*
		*/  touse	/*
		*/  seldep	/* dependent variable is observed
		*/  off		/* regression offset as -offvar- 
		*/  seloffo	/* selection offset as -offset(offvar)-
		*/  rhometh	/* method of computing or handling rho
		*/  twostep	/*
		*/  first	/*
		*/  lvmills

	if "`first'" != "" { local first noisily }
	if "`twostep'" == "" { 
		local do2 "*" 
	}
	
	tempname bprb breg sigma lnsigma rho rho1 athrho Vprb deltbar wtobs
	tempvar pxb delthat mills
	qui {
					/* First-step -- probit */

		`first' probit `seldep' `selind' `wgt' /*
			*/ if `touse', asis `selnc' `seloffo'
		`do2' matrix `Vprb' = get(VCE)
		predict double `pxb', xb, if `touse'
		mat `bprb' = get(_b)
		gen double `mills' = normd(`pxb') / normprob(`pxb') 
		return scalar prbll = e(ll)

					/* Regression strictly to test indep */
		if "`off'" != ""  { 
			local origdep `dep'
			tempname dep
			gen double `dep' = `origdep' - `off'
		}
		if "`wtype'" == "aweight" {
			tempvar wt
			gen double `wt' `wtexp'
			sum `wt' if `touse', meanonly
			replace `wt' = `wt' * r(N) / r(sum)
			sum `wt'  if `touse' & `seldep'
			scalar `wtobs' = r(sum)         
			                /* required, e(N) is rounded */
			regress `dep' `ind' [iweight=`wt'], `nc', /*
				*/ if `touse' & `seldep'
		}
		else { 
			regress `dep' `ind' `wgt', `nc', if `touse' & `seldep' 
			scalar `wtobs' = e(N)
		}
		return scalar regll = -0.5 * `wtobs' * /*
			*/ ( ln(2*_pi) + ln(e(rss)/`wtobs') + 1 )

					/* Compute delta-bar */
		gen double `delthat' = `mills' * (`mills' + `pxb')
		sum `delthat' if `touse' & `seldep', meanonly
		scalar `deltbar' = r(mean)

					/* Second-step -- regression w/ Mills */
		if "`nc'" == "" {
			tempvar one
			gen byte `one' = 1
		}
    		regress `dep' `ind' `one' `mills' `wgt', nocons, /*
			*/ if `touse' & `seldep'
		mat `breg' = get(_b)

					/* Heckman's two-step consistent 
		 			 * estimates of sigma and rho */
		tempname rss ebar adj ratio
		if e(rss) >= . {
			tempvar e
			predict double `e' if `touse', resid
			replace `e' = sum(`e'^2)
			scalar `rss' = `e'[_N]
			if `rss' >= . {
				error 1400
			}
		}
		else    scalar `rss' = e(rss)
		
		scalar `ebar' = `rss' / e(N)
		scalar `sigma' = sqrt(`ebar' + `deltbar'*_b[`mills']^2)

		scalar `rho' = _b[`mills'] / `sigma'
		scalar `rho1' = `rho'
		if abs(`rho') > 1 & `rhometh' >= 1 & `rhometh' <= 2 {
			scalar `rho' = `rho'/abs(`rho')
			if "`twostep'" != "" {
				noi di in blue "note: two-step estimate " /*
				*/ "of rho = "  `rho1' " is being "        /*
				*/ "truncated to " `rho'
			}
			if `rhometh' == 2 {
				scalar `sigma' = _b[`mills']*`rho'
			}
		}
		if `sigma' == 0 {scalar `sigma' = .001}
		if `rho' >= . { scalar `rho' = 0 }

		if "`twostep'" != "" {
			if "`lvmills'" != "" { 
				rename `mills' `lvmills'
				local mills `lvmills'
				label var `mills' "nonselection hazard"
			}
			noi Heck2 `breg' `rho' `sigma' `bprb' `Vprb' `dep' /*
				*/ "`ind'" "`nc'" "`seldep'" "`selind'"    /*
				*/ "`selnc'" "`mills'" "`delthat'" 	   /*
				*/ "`touse'" "`rhometh'" "`rho1'"
		}
		else {
			scalar `athrho' = max(min(`rho',.85), -.85)
			scalar `athrho' = 0.5 * log((1+`athrho') / /*
				*/ (1-`athrho'))
			scalar `lnsigma' = ln(`sigma')
			mat `b0' = `breg'[1,1..(colsof(`breg')-1)] ,   /*
				*/ `bprb' , `athrho' , `lnsigma'
		}


	}
end


/*  Note:  Heck2 assumes that the regression w/ a Mills ratio is the
 *  current set of estimates */

program define Heck2, eclass   /* BetaReg rho sigma Cov(probit) 
				  RegVars RecCons PrbVars Mills touse */
	args	    b		/*  regression Beta
		*/  rho		/*  
		*/  sigma	/*
		*/  bprb	/*  probit Beta
		*/  Vprb	/*  probit VCE
		*/  dep		/*  dependent variable
		*/  ind		/*  independent varlist
		*/  nc		/*  "" or "noconstant"
		*/  seldep	/*  selection eqn dependent variable
		*/  selind	/*  selection eqn indep varlist
		*/  selnc	/*  selection eqn : "" or "noconstant"
		*/  mills	/*  Mills' ratio
		*/  delthat	/*
		*/  touse	/*  To use variable
		*/  rhometh	/*  method of using rho
		*/  rho1	/*  original two-step rho if outide -1,1 */

	tempname XpX1 F Q rho2 VBm 

	/* Compute Heckman's two-step consistent estimates of Cov */

					/* Get X'X inverse from the 
					 * second-step regression */

	if e(rmse) != 0 {
		matrix `XpX1' = get(VCE) / e(rmse)^2
	}
	else {
di in red "insufficient information in sample to estimate heckman model"
		exit 2000
	}

					/*  Heckman's adjustment to Cov  */
	local fulnam `ind'
	if "`nc'" == "" {
		tempvar one
		g byte `one' = 1
		local fulnam "`fulnam' _cons"
	}
	local fulnam "`fulnam' lambda"
	local kreg : word count `fulnam'
	local kreg1 = `kreg' + 1
	qui mat accum `F' = `ind' `one' `mills' `selind' [iw=`delthat'], /*
		*/ `selnc', if `touse' & `seldep'
	mat `F' = `F'[1..`kreg', `kreg1'...]
	mat rowname `F' = `fulnam'
	scalar `rho2' = `rho' * `rho'
	mat `Q' = `rho2' * `F'*`Vprb'*`F''

					/* Finish the variance computation */
	tempvar rho2del
	if (`rho' < -1 | `rho' > 1) & `rhometh' == 3 { scalar `rho2' = 1 }
	qui gen double `rho2del' = 1 - `rho2' * `delthat'  if `touse'
	qui mat accum `VBm' = `ind' `one' `mills' [iw=`rho2del'], /*
		*/ nocons,  if `touse' & `seldep'
	mat `VBm' = `sigma'^2 * `XpX1' * (`VBm' + `Q') * `XpX1'

					/* Build the set of full eqn names */
	local depn : subinstr local dep "." "_"
	local eqnam
	tokenize `fulnam'
	local i 1
	while `i' < `kreg' {
		local eqnam `eqnam' `depn':``i''
		local i = `i' + 1
	}
	if substr("`seldep'",1,2) != "__" {
		local selname `seldep'
		local selname : subinstr local selname "." "_"
	}
	else    local selname select
	local kprb : word count `selind'
	local i 1
	tokenize `selind'
	while `i' <= `kprb' {
		local eqnam `eqnam' `selname':``i''
		local i = `i' + 1
	}
	if "`selnc'" == "" { 
		local eqnam `eqnam' `selname':_cons 
		local kprb = `kprb' + 1
	}
	local eqnam `eqnam' mills:lambda

					/*  Put model and selection 
					 *  equation together 
					 *  Quicker w/ matrix expr, but would
					 *  be hard to understand */

	tempname bfull Vfull zeros Vmills CovMill zeroPrb b0 bmills
	local kreg0 = `kreg' - 1
	mat `b0' = `b'[1,1..`kreg0'] 
	mat `bmills' = `b'[1,`kreg'] 
	mat `bfull' = `b0', `bprb' , `bmills'
	mat `Vmills' = `VBm'[`kreg', `kreg']
	mat `CovMill' = `VBm'[1..`kreg0', `kreg']
	mat `VBm' = `VBm'[1..`kreg0', 1..`kreg0']
	mat `zeros' = J(`kreg0', `kprb', 0)
	mat `zeroPrb' = J(`kprb', 1, 0)
	mat `VBm' = /*
		*/  ( `VBm'       , `zeros'    ,  `CovMills' \      /*
		*/    `zeros''    , `Vprb'     ,  `zeroPrb'  \      /*
		*/    `CovMills'' , `zeroPrb'' ,  `Vmills' )
	mat `VBm' = (`VBm' + `VBm'') / 2
	mat colnames `bfull' = `eqnam'
	mat rownames `VBm' = `eqnam'
	mat colnames `VBm' = `eqnam'

	qui count if `touse'
	capture est post `bfull' `VBm', obs(`r(N)')
	if _rc {
		if _rc == 506 {
			di in red "estimate of rho outside the interval "  /*
				*/ "[-1, 1] has led to a covariance matrix"
			di in red "that is not positive definite"
			* mat `VBm' = 0*`VBm'
			mat `VBm'[1, 1] = 0*`VBm'[1..`kreg0', 1..`kreg0']
		}
		est post `bfull' `VBm', obs(`r(N)')
	}

	est scalar N = `r(N)'
	capture qui test `ind', min
	if _rc == 0 {
		est scalar chi2 = `r(chi2)'
		est scalar df_m = `r(df)'
		est scalar p    = `r(p)'
	}
	qui count if `seldep' == 0 & `touse'
	est scalar N_cens = r(N)	/* early, sic */
	est local chi2type "Wald"
	est local title2 "(regression model with sample selection)"
	est local title "Heckman selection model -- two-step estimates"
	if substr("`mills'", 1, 2) != "__" { est local mills `mills' }
	est scalar rho = `rho'
	if `rho1' != `rho' { est scalar rho1 = `rho1' }
	est scalar sigma = `sigma'
	est scalar lambda = [mills]_b[lambda]
	est scalar selambda = [mills]_se[lambda]
	est local predict "heckma_p"
end

program define Parse50
	gettoken regeq  0 : 0 , parse(" ,")
	gettoken probit 0 : 0 , parse(" ,")

	syntax [if] [in] [fweight aweight] [, Level(cilevel) /*
		*/ noConstant noLOg * ]

	eq ? `regeq'
	local regeq `r(eqname)'
	tokenize "`r(eq)'"
	local dv `1'
	if "`dv'"=="" { local dv `r(eqname)' }
	mac shift
	local vreg `*'

	eq ? `probit'
	local vprobit `r(eq)'

					/* Send things back to caller in 
					 * names used by -syntax- and caller */
	c_local dep `dv'
	c_local depn `dv'
	c_local ind `vreg'
	c_local selind `vprobit'
	c_local nc `constan'
	c_local level `level'

	c_local log `log'
	if "`log'" == "" {
		local showlog "noisily"
	} 
	else    local showlog "quietly"
	c_local showlog `showlog'
	if "`weight'" != "" {
		c_local wtexp	`"`exp'"'
		c_local	wtype	`"`weight'"'
		c_local	wgt	`"[`weight'`exp']"'
	}
	c_local if0 `"`if'"'
	c_local in0 `in'
	c_local mlmetho d2
	c_local option0 `options'
end

exit