_mprobitestimator.ado

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

*! version 1.0.3  01apr2005
* _mprobitestimator - multinomial probit 
	
program _mprobitestimator, eclass 
	syntax varlist [if] [fw pw iw], [	 	///
		case(varname) 				///
		ALTernatives(varname)			///
		ONErpc					///
		noCONSTant				///
		casevars(varlist)			///
		BASEalternative(string)			///
		SCALEalternative(string)		///
		ITERate(integer 50)			///
		noTRANSform				///
		CORRelation(string)			///
		STDdev(string)				///
		INTMethod(string)			///
		INTPoints(integer 0)			///
		INTBurn(integer 0)			///
		INTSeed(string) 			///
		ANTIthetics				///
		CONSTraints(numlist > 0 integer)	///
		Robust					///
		CLuster(varname)			///
		FROM(string)				///
		Level(cilevel) 				///
		TECHnique(string)			///
		INITBhhh(integer 0)			///
		VCE(string)				///
		NOPREserve				///
		DEBug					///
		TRace					///
		noLOg					///
		* ]

	gettoken choice attribs : varlist
	mlopts mlopts rest, `options'

	local onerpc = ("`onerpc'" != "")
	if ("`rest'" != "") PostError , onerpc(`onerpc') `rest'

	local diopts level(`level')

	macro drop MPROBIT_*

	local fullcov 1
	local gencorstr 1
	local cholesky = 1

	if "`stddev'"=="" | substr("`stddev'",1,3)=="het" {
		local stddev "heter"
	}
	else if substr("`stddev'",1,3) == "hom" {
		local stddev "homo"
		local cholesky 0
		local fullcov 0
	}
	else {
		MatrixName "`stddev'" "stddev"
		local stddev `r(name)'
		local varmattype `r(matspec)'
		local varcopy `r(copy)'
		local fullcov 0
		local cholesky 0
	}
	if "`correlation'"=="" | substr("`correlation'",1,3)=="uns" {
		local correlation "unstruct"
		local cholesky = ("`stddev'" == "heter")
	}
	else if substr("`correlation'",1,3)=="ind" {
		local cholesky = 0
		local fullcov = 0
		local correlation "indep"
	}
	else if substr("`correlation'",1,3)=="exc" {
		local cholesky = 0
		local fullcov = 0
		local correlation "exch" 
	}
	else {
		MatrixName "`correlation'" "correlation"
		local correlation `r(name)'
		local cormattype `r(matspec)'
		local corcopy = `r(copy)'

		local fullcov 0
		local cholesky 0
		local gencorstr 0
		tempname corstr
		matrix `corstr' = `correlation'
		/* save the name of the original matrix */
		local cormat = "`correlation'"
		local correlation user
	}

	local antithetics = ("`antithetics'"!="") 

	if ("`intmethod'"=="") local imeth = 2
	else {
		local meth = substr("`intmethod'",1,3)

		if ("`meth'"=="hal" | "`meth'"=="Hal")	local imeth = 1
		else if ("`meth'"=="ham" | "`meth'"=="Ham") local imeth = 2
		else if (substr(`"`intmethod'"',1,1)=="r") local imeth = 3
		else {
			di as error "option intmethod() must specify halton, hammersley " ///
			 "or random" 
			exit 198
		}
	}		
	if (`"`cluster'"'!="") local robust robust
       	if ("`robust'"!="") local pseudo pseudo
	
	tempname N mchoices choices corfix

	if ("`weight'" == "fweight") local tabfrq `weight'`exp'

	qui tabulate `alternatives' [`tabfrq'], matcell(`N') matrow(`mchoices')

	/* save the original order of alternatives */
	matrix `choices' = `mchoices'

	local nchoice = r(r)
	local balanced = 1
	if `onerpc' == 0 {
		local ni = `N'[1,1]
		local m = rowsof(`mchoices')
		local i = 1
		while `balanced' & `++i'<=`m' {
			local balanced = (`N'[`i',1] == `ni')
		}
		local  noscvars  noscvars
		/* -asmprobit- has already marked the sample and kept the 	*/
		/*  valid observations. Use -ml ... if `choice'- 		*/
		local iff if `choice'==1
		local clopt2 cluster(`cluster')
		if `nchoice' > 2 {
			local critopt crittype(log simulated-`pseudo'likelihood)
			local amax = `mchoices'[1,1]
			forvalues i=2/`nchoice' {
				if (`mchoices'[`i',1]>`amax') local amax = `mchoices'[`i',1]
			}
			if (`amax'>9) global MPROBIT_DELIMIT _
		}
	}
	else {
		local robopt `robust'
		local clopt  `"cluster(`cluster')"'
		local iff `if'

		if "`basealternative'" == "" {
			/* mimic -mlogit-: use maximum frequency as default base alternative */
			local ni = `N'[1,1]
			local basealternative = `mchoices'[1,1]
			forvalues i=2/`nchoice' {
				if `N'[`i',1] > `ni' {
					local ni = `N'[`i',1]
					local basealternative = `mchoices'[`i',1]
				}
			}
		}
	}
	local fbase = ("`basealternative'" == "")
	local fscale = ("`scalealternative'" == "")

	forvalues i=1/`nchoice' {
		if (`fbase'==0)	local fbase = (`basealternative'==`mchoices'[`i',1])
		if (`fscale'==0) local fscale = (`scalealternative'==`mchoices'[`i',1])
	}
	if `fbase' == 0 {
di as error "{p}base alternative `basealternative' is not one of the " ///
"alternatives in `alternatives';" ///
" use {help tabulate##|_new:tabulate} for a list of values{p_end}"
		exit 111
	}
	if `fscale' == 0 {
di as error "{p}scale alternative `scalealternative' is not one of " ///
"the alternatives in `alternatives'; " ///
"use {help tabulate##|_new:tabulate} for a list of values{p_end}"
		exit 111
	}

	local varmat = ("`varmattype'"!="")
	if `gencorstr' {
		CorStruct  "`mchoices'" "`correlation'" "`stddev'" "`basealternative'" ///
			"`scalealternative'" `varmat' `varcopy' "`varmattype'"
		tempname corstr
		matrix `corstr' = r(corstr)
		matrix `corfix' = r(corfix)
		local ncor = r(ncor)
		local nvar = r(nvar)
		local stddev "`r(stddev)'"
		local basealternative = `r(base)'
		local scalealternative = `r(fixed)'
	}
	else {
		ValidateCorStruct "`corstr'" "`mchoices'" "`cormat'" "`stddev'" ///
			"`basealternative'" "`scalealternative'" `corcopy'      ///
			"`cormattype'" `varmat' `varcopy' "`varmattype'"
		local ncor = r(ncor)
		local nvar = r(nvar)
		local basealternative = `r(base)'
		local scalealternative = `r(fixed)'
		matrix `corstr' = r(corstr)
		matrix `corfix' = r(corfix)
	}
	local ncov = `ncor' + `nvar'
	if (`ncov'==0) local cholesky = 0
	local reparam = (`cholesky' & "`transform'"=="")

	local d1 d1
	local tech `technique'
	if `onerpc' == 0 {
		ParseTech `"`technique'"' `initbhhh' "`vce'" 

		local technique `"`s(technique)'"'
		local tech `"`s(tech)'"'
		local d1 "`s(d1)'"
		local vce "`s(vce)'"
		local tech2 "`s(tech2)'"
		local vce2 "`s(vce2)'"
		local d21 "`s(d21)'"
		local initbhhh = `s(initbhhh)'
		local techmeth "`s(techmeth)'"
		
		if `nchoice' <= 2 {
			local intpoints = 0
			local initbhhh = 0 
		}
		else if `intpoints' <= 0 {
			/* number of points in the sequence */
			if (`imeth'==3) local intpoints  = 100*`nchoice'
			else local intpoints = 50*`nchoice'
		}
	}
	else {
		/* number of quadrature points */
		if (`nchoice' <= 2) local intpoints = 0
		else if (`intpoints' <= 0) local intpoints = 15

		if ("`technique'"=="") local tech nr 
	}

	local nattrib: word count `attribs'

	/* choices vector must be in the same order as the 	*/
	/* row stripe of corstruct				*/
	global MPROBIT_ALTLEVELS `mchoices'
	global MPROBIT_ALTERNATIVES "`alternatives'"
	global MPROBIT_ID "`case'"

	if ("`intseed'" == "") global MPROBIT_SEED = "`c(seed)'"
	else global MPROBIT_SEED = "`intseed'"
	
	global MPROBIT_METHOD = `imeth'
	global MPROBIT_BURN = `intburn'
	global MPROBIT_TRACE = ("`trace'"!="")
	if (`reparam') {
		global MPROBIT_CHOLESKY 1
		local cholesky 0
	}
	else global MPROBIT_CHOLESKY = `cholesky'
	global MPROBIT_CORRELATION = "`correlation'"
	global MPROBIT_ANTITHETICS = `antithetics'
	global MPROBIT_FULLCOV = `fullcov'
	global MPROBIT_CORSTRUCT = "`corstr'"
	global MPROBIT_CORFIXED = "`corfix'"
	global MPROBIT_NVAR = `nvar'
	global MPROBIT_NCOR = `ncor'
	global MPROBIT_ATTRIBS `"`attribs'"'
	global MPROBIT_NATTRIB = `nattrib'
	global MPROBIT_CASEVARS  `"`casevars'"'
	global MPROBIT_BALANCED = `balanced'
	global MPROBIT_BASE = `basealternative'
	
	if (`nchoice' == 2) local intpoints = 0
	global MPROBIT_NPOINTS = `intpoints'	

	if `ncov' > 0 {
		GenTmat `nvar' `ncor' "`corstr'" `fullcov'
	
		tempname T
		matrix `T' = r(T)
		global MPROBIT_TMAT `T'
		if `onerpc' == 0 {
			forvalues i=1/`ncov' {
				tempvar v`i'
				gen double `v`i'' = 0.0
				local vscrs `"`vscrs' `v`i''"'
			}
			global MPROBIT_VSCRS `vscrs'
		}
	}
	if `nattrib' > 0 {	
		forvalues i=1/`nattrib' {
			tempvar a`i'
			gen double `a`i'' = 0.0
			local ascrs `"`ascrs' `a`i''"'
		}
		global MPROBIT_ASCRS `ascrs'
		local wald = -1
	}
	else local wald = 0

	if `onerpc' {
		/* need one for the base -mprobit_lf- will pass it on to	*/
		/* mata routine _mprobit 					*/
		tempvar cb
		gen double `cb' = 0
		global MPROBIT_CSCRS `cb'
	}
	if ("`constraints'" != "") local conopt constraints(`constraints')

	tempname from1 caselvls
	if "`from'" != "" {
		/* note altlevels uses the original order of alternatives */
		CoefFrom `"`from'"', alternatives(`alternatives') altlevels(`choices') 	///
			attribs(`attribs') casevars(`casevars') base(`basealternative') ///
			`constant' cholesky(`cholesky') choice(`choice')

		if `reparam' & `"`constraints'"'!="" {
			CheckConstraints "`from1'" `ncov' `"`constraints'"'
			/* constraints on variance/covariance parameters, cannot */
			/* reparameterize 					 */
			if (`e(vconstr)' > 0) {
				local reparam = 0
				global MPROBIT_CHOLESKY = `cholesky'
			}
		}
	}
	else {
		/* note altlevels uses the original order of alternatives */
		InitMProbit `choice' `iff' [`weight'`exp'], case(`case') attribs(`attribs')       ///
			alternatives(`alternatives') altlevels(`choices') base(`basealternative') ///
			`constant' casevars(`casevars') cholesky(`cholesky') `conopt'

		local vconstr = `r(vconstr)'
		if `reparam' & `vconstr' > 0 {
			local reparam = 0
			global MPROBIT_CHOLESKY = `cholesky'
		}
	}
	
	matrix `from1' = r(from)
	local altmodel `"`r(altmodel)'"'
	local casemodels `"`r(casemodels)'"'
	local thetas `"`r(thetas)'"'
	local rcfnms `"`r(rcfnms)'"'
	local baselab `"`r(baselab)'"'
	global MPROBIT_CASEEQS `"`r(caseeqs)'"'
	global MPROBIT_NCASECOEF = `r(ncasecf)'
	global MPROBIT_CONST = `r(const)'
	if $MPROBIT_NCASECOEF > 0 {
		local wald = `wald' - `nchoice' + 1
		matrix `caselvls' = r(caselevels)
		global MPROBIT_CASELEVELS `caselvls'
	}
	if `reparam' {
		/* put parameter estimates into cholesky metric for optimization */
		Reparameterize "`from1'" 0 `"`rcfnms'"'
		local thetas = `"`r(thetas)'"'
	}
	if `onerpc' == 0 { 
		forvalues i=1/$MPROBIT_NCASECOEF {
			tempvar c`i'`j'
			gen double `c`i'`j'' = 0.0
			local cscrs `"`cscrs' `c`i'`j''"'
		}
		global MPROBIT_CSCRS `"`cscrs'"'
	}
	if "`tech2'"=="" & `reparam' {
		local tech2 nr
		local vce2  oim
		if ("`vce'"=="opg") local d21 d2
		else local d21 d1
	}
	if "`tech2'" != "" {
		local novce  novce
		local waldopt2 wald(`wald')
	}
	else {
		if ("`vce'"!="") local vceopt  vce(`vce')
		local waldopt wald(`wald')
	}

	local converged = 0
	local rc = 0
	if `initbhhh' > 0 {
		/* `onercp'!=0, we use d2 and nr but the algorithm is actually 	*/
		/*  bhhh technique and opg for negative hessian 		*/
		ml model d2 mprobit_lf `altmodel' `casemodels' `thetas' `iff' 	///
			[`weight'`exp'], technique(nr) collinear `nopreserve' 	///
			init(`from1') `conopt' `noscvars' `critopt'

		if "`log'" == "" {
			di _n in smcl in gr "(setting optimization to BHHH)"
			if ("`techmeth'" != "") local log1 log
		}

		ml maximize, iterate(`initbhhh') search(off) nooutput `trace' ///
			`dopts' `log' novce nowarning

		matrix `from1' = e(b)
		local converged = e(converged)
		local rc = e(rc)
		if `converged' & "`tech2'"=="" {
			local tech2 nr
			local vce2  oim
			if ("`vce'"=="opg") local d21 d2
			else local d21 d1
			local waldopt2 wald(`wald')
		}
	}
	if `converged' == 0 {
		if "`debug'" != "" {
			global MPROBIT_TRACE = 0

			ml model d1debug mprobit_lf `altmodel' `casemodels' `thetas' `iff' ///
				[`weight'`exp'], `nopreserve' collinear init(`from1') `noscvars'
			local dopts gradient showstep trace
		}
		else {
			/* let -ml- approximate the negative hessian (unless tech=bhhh) */
			ml model `d1' mprobit_lf `altmodel' `casemodels' `thetas' `iff'     ///
				[`weight'`exp'], technique(`tech') `nopreserve' collinear   ///
				`conopt' init(`from1') `noscvars' `vceopt' `robopt' `clopt' ///
				`mlopts' `critopt' `waldopt'

			local nooutput nooutput
			if ("`log1'" != "") di _n in smcl in gr "(switching optimization to " ///
				"`techmeth')"
		}
		ml maximize, iterate(`iterate') search(off) `nooutput' `trace' `novce' ///
			`dopts' `log'

		local converged = e(converged)
		local rc = e(rc)
	}
	if "`tech2'" != "" {
		matrix `from1' = e(b)
		if `reparam' {
			/* put parameter estimates into desired metric */
			if ("`log'" == "") di _n in gr "Reparameterizing to standard " ///
				"deviation/correlation metric" 

			Reparameterize "`from1'" 1 `"`rcfnms'"' "`choices'"

			local thetas = `"`r(thetas)'"'
			global MPROBIT_CHOLESKY = 0
		}
		ml model `d21' mprobit_lf `altmodel' `casemodels' `thetas' [`weight'`exp']     ///
			`iff', technique(`tech2') vce(`vce2') collinear `conopt' init(`from1') ///
			`nopreserve' `mlopts' `critopt' `waldopt2' 

		qui ml maximize, search(off) nooutput iteronly1 iterate(0) nowarning

		ereturn scalar converged = `converged'
		ereturn scalar rc = `rc'
	}

	if "`robust'"!=""  & `onerpc'==0 {
		tempname V b g negH lnf
		global ML_y1 = "`choice'"
		scalar `lnf' = 0
		matrix `V' = e(V)
		matrix `from1' = e(b)

		mprobit_lf 3 "`from1'" "`lnf'" "`g'" "`negH'"  
		/* temporary varaibles for covariance scores */
		if `nattrib' > 0 {
			/* dummy names for coefficients */
			foreach  scr of global MPROBIT_ASCRS {
				local robnms `"`robnms' `scr':_cons"'
			}
		}
		if $MPROBIT_NCASECOEF > 0 {
			/* dummy names for coefficients */