asmprobit_lf.ado

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

*!version 1.0.2  23feb2005
program asmprobit_lf
	version 9

	args todo b lnf g negH

	GetCov "`b'"
	if `r(bail)' {
		/* force -ml- to do a step half */
		scalar `lnf' = .
		exit 0
	}
	tempname d R rho
	matrix `R' = r(R)
	matrix `d' = r(d)
	if ($MPROBIT_CHOLESKY==0) matrix `rho' = r(rho)

	/* compute utilities */
	local ncov = $MPROBIT_NVAR + $MPROBIT_NCOR
	tempvar u u1
	local util `"`u'"'
	if $MPROBIT_NATTRIB > 0 {
		matrix score double `u' = `b', eq(#1)		
		local eq = 1
	}
	else {
		gen double `u' = 0		
		local eq = 0
	}
	if $MPROBIT_NCASECOEF > 0 {
		forvalues i=1/$MPROBIT_NALT {
			if (`i'==$MPROBIT_BASE) continue

			matrix score double `u1' = `b' if $MPROBIT_ALTERNATIVES==`i', eq(#`++eq')
			qui replace `u' = `u' + `u1' if $MPROBIT_ALTERNATIVES==`i'
			drop `u1'
		}
	}

	if (`todo' > 1) local todo1 = 1
	else local todo1 = `todo'

	if (`todo' < 0) {
		/* computing probabilties for predict */
		local lf `6'
		local y $MPROBIT_TOUSE
	}
	else {
		tempvar lf
		gen double `lf' = 0.0
		local y $ML_y1
	}

	/* go ahead and set the seed now */
	if ($MPROBIT_METHOD==3) set seed $MPROBIT_SEED
	tempname err
	scalar `err' = 0
	mata: _mprobit_simulator(`"`y'"',`"`util'"',`"$MPROBIT_ID"',`"`lf'"',		///
		`"$MPROBIT_ATTRIBS"',`"$MPROBIT_ASCRS"',`"$MPROBIT_CASEVARS"',		///
		`"$MPROBIT_CSCRS"',$MPROBIT_CONST,`"$MPROBIT_VSCRS"',			///
		`"$MPROBIT_ALTERNATIVES"',`"$MPROBIT_ORDER"',`"`R'"',`todo1',		///
		$MPROBIT_METHOD,$MPROBIT_NPOINTS,$MPROBIT_BURN,$MPROBIT_ANTITHETICS,	///
		$MPROBIT_CHOLESKY,$MPROBIT_BALANCED,`"$MPROBIT_TMAT"',`"`err'"')

	local iff if $ML_y1==1
	if `err' {
		scalar `lnf' = .
		/* force -ml- to do a step half */
		if $MPROBIT_TRACE {
			di as error "covariance is not positive definite"
			if $MPROBIT_CHOLESKY {
				tempname V
				di "R"
				matrix list `R', noheader
				matrix `V' = `R'*`R''
			}
			else {
				local V `R'
			}
			di
			di "covariance"
			matrix list `V', noheader
		}
	}
	else if `todo'>=0 & `todo'<=2 {
		mlsum `lnf' = `lf'
	}
	else if `todo'==3 {
		/* special call from asmprobit for scores to call _robust */  
		qui summarize `lf' $MPROBIT_WTOPT `iff'
		scalar `lnf' = r(sum)
	}
	if (`todo'<=0) exit 0

	GetACGrad "`g'" "`lnf'" `todo' 

	local eq = `r(eq)'
	if (`todo'==3) local ig = `r(ig)'
	tempname g1
	tokenize $MPROBIT_VSCRS
	local k = 0
	forvalues j=1/$MPROBIT_NVAR {
		local k = `k' + 1
		qui replace ``k'' = ``k''*`d'[1,`j'] `iff' 
		if `todo' < 3 {
			mlvecsum `lnf' `g1' = ``k'', eq(`++eq')
			matrix `g' = (`g',`g1')
		}
		else { 
			qui summarize ``k'' $MPROBIT_WTOPT `iff'
			matrix `g'[1,`++ig'] = r(sum)
		}
	}
	forvalues j=1/$MPROBIT_NCOR {
		local k = `k' + 1
		if ($MPROBIT_CHOLESKY==0) qui replace ``k'' = ``k''*(1.0-`rho'[1,`j']^2) `iff'
		if `todo' < 3 {
			mlvecsum `lnf' `g1' = ``k'', eq(`++eq')
			matrix `g' = (`g',`g1')
		}
		else {
			qui summarize ``k'' $MPROBIT_WTOPT `iff'
			matrix `g'[1,`++ig'] = r(sum)
		}
	}
	if (`todo'!=2) exit 0

	mata: _mprobit_outer_prod(`"$MPROBIT_ASCRS"', `"$MPROBIT_CSCRS"',     ///
		`"$MPROBIT_VSCRS"', `"$ML_y1"', $MPROBIT_NALT, $MPROBIT_BASE, ///
		"$MPROBIT_ALTERNATIVES", $MPROBIT_BALANCED, `"$MPROBIT_ID"', "`negH'"')
end

program GetCov, rclass
	args b 

	/* variances and correlations */
	tempname D L rho var
	matrix `D' = I($MPROBIT_NALT) 
	if $MPROBIT_CHOLESKY {
		/* cholesky factor of covariances */
		local rjj lnl
		local rij l
		local L `D'
		local R `D'
	}
	else {
		/* variances and correlations */
		tempname R
		matrix `L' = I($MPROBIT_NALT)
		local rjj lnsigma
		local rij atanhr
	}
	if $MPROBIT_NCOR > 0 {
		matrix `rho' = J(1,$MPROBIT_NCOR,.)
	}
	else {
		matrix `rho' = J(1,1,.)
	}
	if $MPROBIT_NVAR > 0 {
		matrix `var' = J(1,$MPROBIT_NVAR,.)
	}
	else {
		matrix `var' = J(1,1,.)
	}
	
	local rc = 0
	local bail = 0
	forvalues j=1/$MPROBIT_NALT {
		local cj = $MPROBIT_CORSTRUCT[`j',`j']
		if `cj' < . {
			if `var'[1,`cj'] >= . {
				if $MPROBIT_FULLCOV {
					if ($MPROBIT_CHOLESKY) local jj `j'_`j'
					else local jj `j'
				}
				else local jj P`cj'
				local k = colnumb(`b',"`rjj'`jj':")
				if `k' >= . {
					di as error "could not find coefficient associated " ///
					 "with parameter `rjj'`jj':_cons"
					local rc = 497
					local bail = 1
					continue, break
				}
				matrix `var'[1,`cj'] = exp(`b'[1,`k']) 
				if `var'[1,`cj'] >= 1.0e5 {
					if $MPROBIT_TRACE {
						di as error "invalid standard deviation " ///
						 "estimate `rjj'`jj' = " `b'[1,`k']
					}
					local bail = 1
					continue, break
				}
			}
			matrix `D'[`j',`j'] = `var'[1,`cj']
		}
		else {
			matrix `D'[`j',`j'] = $MPROBIT_CORFIXED[`j',`j']
		}
		local j1 = `j'+1
		forvalues i=`j1'/$MPROBIT_NALT {
			local ci = $MPROBIT_CORSTRUCT[`i',`j']
			if `ci' < . {
				if `rho'[1,`ci'] >= . {
					if ($MPROBIT_FULLCOV) local ij `i'_`j'
					else local ij P`ci'

					local k = colnumb(`b',"`rij'`ij':")
					if `k' >= . {
						local rc = 497
						local bail = 1
						di as error "could not find coefficient " ///
						 "associated with parameter `rij'`ij':_cons"
						continue, break
					}
					if $MPROBIT_CHOLESKY {
						matrix `rho'[1,`ci'] = `b'[1,`k']
					}
					else {
						matrix `rho'[1,`ci'] = tanh(`b'[1,`k'])
					}
				}
				matrix `L'[`i',`j'] = `rho'[1,`ci']
			}
			else {	
				matrix `L'[`i',`j'] = $MPROBIT_CORFIXED[`i',`j']
			}
			
			if ($MPROBIT_CHOLESKY==0) matrix `L'[`j',`i'] = `L'[`i',`j']
		}
		if (`bail') continue, break
	}

	if (`rc'>0) exit `rc'
	if `bail' == 0 {
		if ($MPROBIT_CHOLESKY == 0) matrix `R' = `D'*`L'*`D'
		return matrix R = `R'
		return matrix d = `var'
		return matrix rho = `rho'
	}
	return local bail `bail'
end

program GetACGrad, rclass
	args g lnf todo

	tempname g1
	local ncov = $MPROBIT_NVAR + $MPROBIT_NCOR

	if `todo' < 3 {
		/* need to call -mlvecsum- once for eq(1) even though the derivatives	*/
		/* will be bogus							*/
		if ($MPROBIT_NATTRIB > 0) local scr : word 1 of $MPROBIT_ASCRS
		else local scr : word 1 of $MPROBIT_CSCRS
	
		mlvecsum `lnf' `g' = `scr', eq(1)
	}
	else {
		loca ns = $MPROBIT_NATTRIB+$MPROBIT_NCASECOEF*($MPROBIT_NALT-1)+`ncov'
		matrix `g' = J(1,`ns',0)
	}
	/* now get the real gradient						*/
	local ig = 0
	local eq = 0
	if $MPROBIT_NATTRIB > 0 {
		foreach scr of global MPROBIT_ASCRS {
			qui summarize `scr' $MPROBIT_WTOPT if $ML_y1==1
			matrix `g'[1,`++ig'] = r(sum)
		}
		local `++eq'
	}
	if $MPROBIT_NCASECOEF > 0 {
		if $MPROBIT_NATTRIB == 0 {
			if ($MPROBIT_BASE == 1) local i1 = 2
			else local i1 = 1
			foreach scr of global MPROBIT_CSCRS {
				qui summarize `scr' $MPROBIT_WTOPT if $MPROBIT_ALTERNATIVES==`i1'
				matrix `g'[1,`++ig'] = r(sum)
			}
			local `++i1'
		}
		else local i1 = 1

		forvalues i=`i1'/$MPROBIT_NALT {
			if ($MPROBIT_BASE == `i') continue
			foreach scr of global MPROBIT_CSCRS {
				qui summarize `scr' $MPROBIT_WTOPT if $MPROBIT_ALTERNATIVES==`i'
				if (`todo' < 3) {
					matrix `g' = (`g',r(sum))
				}
				else {
					matrix `g'[1,`++ig'] = r(sum)
				}
			}
		}
		local eq = `eq' + $MPROBIT_NALT - 1
	}
	return local eq = `eq'
	return local ig = `ig'
end