asmprobit.ado

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

*! version 1.0.7  03jun2005
* asmprobit - multinomial probit model for alternative specific
* 		and case specific variables
	
program asmprobit, eclass byable(onecall)
	version 9

	if _by() {
		local BY `"by `_byvars'`_byrc0':"'
	}
	if replay() {
		if `"`BY'"' != "" error 190
		if `"`e(cmd)'"' != "asmprobit" error 301
		Replay `0'
		exit
	}
	qui syntax [anything] [fw pw iw] [if] [in], case(varname) [ ///
		Cluster(varname) VCE(passthru) *]


	if `"`vce'"' != "" {
		tempname id
		local idopt idcluster(`id')
		local case0 "`case'"
		if `"`cluster'"' == "" {
			quietly gen `id' = `case'
			local clopt cluster(`case')
			local caopt case(`id')
		}
		else {
			quietly gen `id' = `cluster'
			local clopt cluster(`cluster')
			local caopt case(`case')
			local bsgropt group(`case')
		}
		if "`weight'" != "" {
			local wgt [`weight'`exp']
		}
		local vceopts jkopts(`clopt') ///
			bootopts(`clopt' `idopt' `bsgropt')
		`by' _vce_parserun asmprobit, `vceopts' : ///
			`anything' `wgt' `if' `in', ///
			`caopt' `vce' `options'
		if "`s(exit)'" != "" {
			ereturn local casevar `case'
			ereturn local clustvar `cluster'
			if "`cluster'" == "" {
				local cmd1 `"`e(command)'"'
				local cmd2 : subinstr local cmd1 "`id'" "`case'"
				ereturn local command `"`cmd2'"'
			}
			exit
		}
		qui syntax [anything] [fw pw iw] [if] [in], case(varname) [ * ]
		local 0 `"`anything' [`weight'`exp'] `if' `in', case(`case0') `options'"'
	}
	cap noi `BY' Estimate `0'
	local rc = _rc
	cap mata : mata drop _mprobit_panel_info
	macro drop MPROBIT_*
	exit `rc'
end

program Estimate, eclass byable(recall)
	syntax varlist [if] [in] [fw pw iw],		/// 
		case(varname) 				///
		ALTernatives(varname) [			///
		noCONstant				///
		CASEVars(varlist)			///
		BASEalternative(string)			///
		SCALEalternative(string)		///
		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)				///
		TECHnique(string)			///
		INITBhhh(integer 0)			///
		ITERate(passthru)			///
		VCE(string)				///
		NOPREserve				///
		Level(cilevel)				///
		DEBug					///
		MLInit					///
		TRace					///
		noLOg					///
		noDROP					///
		* ]

	gettoken choice attribs : varlist

	mlopts mlopts rest, `options'
	PostError , `rest'
	local diopts level(`level')

	cap confirm string variable `case'
	if _rc == 0 {
		di as error "case() variable cannot be a string variable"
		exit 108 
	}
	cap assert float(`case')==float(floor(`case')) 
	if _rc != 0 {
		di as error "case() variable `case' must be coded as integers"
		exit 459 
	}
	local nattrib: word count `attribs'
	local ncasev: word count `casevars'
	if `nattrib'==0 & `ncasev'==0 & "constant"!="" {
		di as error "too few variables specified"
		exit 102
	} 

	marksample touse
	markout `touse' `case' `alternatives' `casevars', strok
	if (`"`cluster'"'!="") {
		markout `touse' `cluster'
		local robust robust
	}
	preserve
	qui keep if `touse'

	cap assert `choice'==0 | `choice'==1 
	if _rc != 0 {
		di as error "{p}depvar `choice' must be a 0-1 variable indicating " ///
		 "which alternative is the outcome{p_end}"
		exit 450
	}

	tempvar altern altlabels
	cap qui egen `altern' = group(`alternatives'), lname(`altlabels')
	if _rc > 0 {
		/* should not happen */
		di as error "failed to renumber alternative levels"
		error _rc
	}
	sort `case' `altern'
	tempname casestat
	local bfwts = 0
	if "`weight'" != "" {
		tempvar wt 
		qui gen `wt' `exp'
		local bfwts = ("`weight'"=="fweight")
	}
	mata: _mprobit_validate_choice(`"`choice'"',`"`altern'"',`"`case'"', ///
		`"`casevars'"', `"`wt'"', `bfwts', `"`touse'"', `"`casestat'"')
	if `casestat'[5,1] == 0 {
		di as error "no cases remain after excluding invalid observations"
		exit 2000
	}
	if `casestat'[5,1] == 1 {
		di as error "only one case remains after excluding invalid observations"
		exit 2001
	}
	if `casestat'[6,1] > 0 {
		/* dropped observations, need to regroup */
		qui keep if `touse'
		if _N == 0 {
			di as error "no cases remain after excluding invalid observations"
			exit 2000
		}
		drop `altern'
		cap qui egen `altern' = group(`alternatives'), lname(`altlabels')
		if _rc > 0 {
			/* should not happen */
			di as error "failed to renumber alternative levels"
			error _rc
		}
		sort `case' `altern'
		/* recall to setup panel info */
		mata: _mprobit_validate_choice(`"`choice'"',`"`altern'"',`"`case'"', ///
			`"`casevars'"', `"`wt'"', `bfwts', `"`touse'"', `"`casestat'"')
	}
	local balanced = `casestat'[4,1]
	if "`drop'" == "" {
		/* remove collinear variables */
		_rmcoll `attribs', noconstant
		local attribs `"`r(varlist)'"'
		local nattrib: word count `attribs'

		_rmcoll `casevars', `constant'
		local casevars `"`r(varlist)'"'
	}

	macro drop MPROBIT_*

	local fullcov 1
	local gencorstr 1
	local cholesky = 1
	local lstd = length("`stddev'")
	if "`stddev'"=="" | "`stddev'"==substr("heteroskedastic",1,max(3,`lstd')) {
		local stddev "heter"
	}
	else if "`stddev'" == substr("homoskedastic",1,max(3,`lstd')) {
		local stddev "homo"
		local cholesky 0
		local fullcov 0
	}
	else {
		MatrixName `stddev', type(stddev)
		local stddev `r(name)'
		local varmattype `r(matspec)'
		local cholesky 0
		local fullcov 0
	}
	local lcor = length("`correlation'")
	if "`correlation'"=="" | "`correlation'"==substr("unstructured",1,max(3,`lcor')) {
		local correlation "unstruct"
		local cholesky = ("`stddev'"=="heter" & "`from'"=="")
	}
	else if "`correlation'" == substr("independent",1,max(3,`lcor')) {
		local cholesky = 0
		local fullcov 0
		local correlation "indep"
	}
	else if "`correlation'" == substr("exchangeable",1,max(3,`lcor')) {
		local cholesky = 0
		local fullcov 0
		local correlation "exch" 
	}
	else {
		MatrixName `correlation', type(correlation)
		local correlation `r(name)'
		local cormattype `r(matspec)'

		local cholesky 0
		local fullcov 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 lmth = length("`intmethod'")

		if ("`intmethod'"==substr("halton",1,max(3,`lmth'))) local imeth = 1
		else if ("`intmethod'"==substr("hammersley",1,max(3,`lmth'))) local imeth = 2
		else if ("`intmethod'"==substr("random",1,max(3,`lmth'))) local imeth = 3
		else {
			di as error "option intmethod() must specify halton, hammersley, " ///
			 "or random" 
			exit 198
		}
	}		

	local iff if `choice'==1
	if ("`cluster'"!="") local clopt2 cluster(`cluster')
       	if ("`robust'"!="" | "`weight'"=="pweight") local pseudo pseudo
	if "`weight'"=="iweight" | "`weight'"=="fweight" {
		local wtoptinit [`weight'`exp']
		local wtopt [`weight'`exp']
	}
	else if "`weight'" == "pweight" {
		local wtopt [iweight`exp']	
	}
	
	tempname corfix
	/* do not use frequencies */
	cap confirm string variable `alternatives'
	if _rc == 0 {
		qui tabulate `altern'
	}
	else {
		tempname altlevels
		qui tabulate `alternatives', matrow(`altlevels')
		local altlabels : value label `alternatives'
	}

	local nalt = r(r)
	if `nalt' > 20 {
		di as error "there are `nalt' alternatives in `alternatives'; " ///
		 "the maximum number of alternatives is 20"
		exit 149
	}
	
	AlternativeIndex "`altlevels'" `nalt' "`altlabels'" "`basealternative'" ///
		"basealternative" 
	local ibase = r(index)
	AlternativeIndex "`altlevels'" `nalt' "`altlabels'" "`scalealternative'" ///
		"scalealternative"
	local iscale = r(index)
	if `iscale'<. & `iscale'==`ibase' {
		di as error "scalealternative() cannot equal the basealternative()"
		exit 459	
	}

	local varmat = ("`varmattype'"!="")
	if `gencorstr' {
		CorStruct `nalt' "`correlation'" "`stddev'" `ibase' `iscale' `varmat' ///
			"`varmattype'" 
		tempname corstr
	}
	else {
		ValidateCorStruct "`corstr'" `nalt' "`cormat'" "`stddev'" `ibase' `iscale' ///
			"`cormattype'" `varmat' "`varmattype'"
	}
	tempname order0
	local ncor = r(ncor)
	local nvar = r(nvar)
	local ibase = `r(ibase)'
	local iscale = `r(ifixed)'
	matrix `corstr' = r(corstr)
	matrix `corfix' = r(corfix)
	matrix `order0' = r(order)

	local ncov = `ncor' + `nvar'
	if (`ncov'==0) local cholesky = 0
	local reparam = (`cholesky' & "`transform'"=="")

	local d1 d1
	ParseTech `"`technique'"' `initbhhh' "`vce'" 

	local technique `"`s(technique)'"'
	local tech `"`s(tech)'"'
	local d1 "`s(d1)'"
	local vce0 "`s(vce0)'"
	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 `nalt' <= 2 {
		local intpoints = 0
		local intburn = 0
		local initbhhh = 0 
	}
	else {
		if `intpoints' < 0 {
			di as error "intpoints() must be greater than 0"
			exit 198
		}
		if `intpoints' == 0 {
			/* number of points in the sequence */
			if (`imeth'==3) local intpoints  = 100*`nalt'
			else local intpoints = 50*`nalt'
		}
		if `intburn' < 0 {
			di as error "intburn() must be greater than 0"
			exit 198
		}
		else if `intburn'>0 & `imeth'==3 {
			di as error "intburn() cannot be used with intmethod(random)"
			exit 198
		}
		local critopt crittype(log simulated-`pseudo'likelihood)
	}

	global MPROBIT_NALT = `nalt'
	global MPROBIT_ALTERNATIVES "`altern'"
	global MPROBIT_ID "`case'"
	if ("`intseed'" == "") global MPROBIT_SEED = "`c(seed)'"
	else if (`imeth' == 3) 	global MPROBIT_SEED = "`intseed'"
	else {
		di as error "{p}intseed() is not allowed when specifying the " ///
		 "Halton or Hammersley sequences{p_end}"
		exit 184
	}
	
	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_FULLCOV = `fullcov'
	global MPROBIT_CORRELATION = "`correlation'"
	global MPROBIT_ANTITHETICS = `antithetics'
	global MPROBIT_CORSTRUCT = "`corstr'"
	global MPROBIT_CORFIXED = "`corfix'"
	global MPROBIT_ORDER = "`order0'"
	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 = `ibase'
	global MPROBIT_WTOPT `"`wtopt'"'
	
	if (`nalt' == 2) local intpoints = 0
	global MPROBIT_NPOINTS = `intpoints'	

	if `ncov' > 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

	GetModel `choice', case(`case') attribs(`attribs') ibase(`ibase') ///
		alternatives(`altern') altname(`alternatives') `constant' ///
		casevars(`casevars') 

	local ncoef = `r(ncoef)'
	if `casestat'[5,1] <= 10*(`ncov'+`ncoef') {
		if `casestat'[5,1] <= `ncov'+`ncoef' {
			di as error "there are only `=`casestat'[5,1]' cases " ///
			 "to estimate `=`ncov'+`ncoef'' parameters"
			exit 2001
		}
		di in gr "{p 0 7}Note: there are only `=`casestat'[5,1]' cases " ///
		 "to estimate `=`ncov'+`ncoef'' parameters{p_end}"
	}
	local altmodel `"`r(altmodel)'"'
	local casemodels `"`r(casemodels)'"'
	local rcfnms `"`r(rcfnms)'"'

	global MPROBIT_NCASECOEF = `r(ncasecf)'
	global MPROBIT_CONST = `r(const)'
	if $MPROBIT_NCASECOEF > 0 {
		local wald = `wald' + `nalt' - 1
		global MPROBIT_ALTEQS `"`r(alteqs)'"'
	}
	tempname b order ic ilog
	if "`from'" == "" {
		local vconstr = 0
		if `"`constraints'"' != "" {
			VParStripe 0
			matrix `b' = J(1,`ncoef'+`ncov',0)
			matrix colnames `b' = `rcfnms' `r(vcfnms)'

			CheckConstraints "`b'" `ncov' `"`constraints'"'
			/* vconstr is the number of constraints on the variance parameters 	*/
			/* bconstr is the number of constraints on the regression parameters 	*/
			local vconstr = `e(vconstr)'
			local bconstr = `e(bconstr)'
			if `bconstr' > 0 {
				/* pass the constraint matrix on to clogit for inital estimates */
				tempname C 
				matrix `C' = e(C)
				local cmat cmat(`C')
			}	
		}
		if `reparam' & `vconstr' > 0 {
			/* constraints on variance/covariance parameters, cannot */
			/* reparameterize 					 */
			local reparam = 0
			global MPROBIT_CHOLESKY = `cholesky'
		}
		if $MPROBIT_CHOLESKY {
			/* pivot base and scale to 1st and second position */
			matrix `order' = J(`nalt',1,2)
			matrix `order'[1,1] = `ibase'
			if `iscale' < . {
				matrix `order'[2,1] = `iscale'
				local j = 3
			}
			else local j = 2

			local k = 1
			local reorder = (`order'[1,1]!=1 | `order'[2,1]!=2)
			forvalues i=`j'/`nalt' {
				while `k'==`ibase' | `k'==`iscale' {
					local `++k'
				}
				matrix `order'[`i',1] = `k++'
				local reorder = (`reorder' | `order'[`i',1]!=`i')
			}
			if (`reorder')	{
				PivotMat "`corstr'" "`order'" 1
				PivotMat "`corfix'" "`order'" 0
			}
			global MPROBIT_ORDER = "`order'"
		}
		if "`mlinit'" == "" {
			/* cannot use ml search(on) since it will try to scale the coefficients */
			InitMProbit `choice' `iff' `wtoptinit', case(`case') attribs(`attribs') /// 
				alternatives(`altern')  ibase(`ibase') casevars(`casevars')     ///
				`constant' `cmat' coefnames(`rcfnms')

			matrix `b' = r(from)
			local seaopt search(off)
			local initopt init(`b')
		}
		else local seaopt search(norescale) 
	}
	else {
		local seaopt search(off)
		local initopt init(`from')
	}
	if `ncov' > 0 {
		VParStripe $MPROBIT_CHOLESKY
		local thetas `"`r(thetas)'"'
		local vcfnms `"`r(vcfnms)'"'
		if "`from'"=="" & "`mlinit'"=="" {
			tempname from
			matrix  `from' = (`b',J(1,`ncov',0))
			matrix colnames `from' = `rcfnms' `vcfnms'
		}

		GenTmat `nvar' `ncor' "`corstr'" 
		tempname T
		matrix `T' = r(T)
		global MPROBIT_TMAT `T'
	}
	forvalues i=1/$MPROBIT_NCASECOEF {
		tempvar c`i'
		gen double `c`i'' = 0.0
		local cscrs `"`cscrs' `c`i''"'
	}
	global MPROBIT_CSCRS `"`cscrs'"'
	if $MPROBIT_NCASECOEF>0 & "`pseudo'"!="" {
		/* check that we will be able to generate scores */
		/*  after estimation				 */