xtmixed.ado

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

*! version 1.1.7  23may2005
program xtmixed, eclass byable(onecall)
	version 9
	
	if _by() {
		local by "by `_byvars'`_byrc0':"
	}
	if replay() {
		if "`e(cmd)'" != "xtmixed" {
			error 301
		}
		if _by() {
			error 190
		}
		Display `0'
		exit
	}

	capture noisily `by' Estimate `0'
	local rc = _rc
	capture mata: _xtm_cleanup()
	exit `rc'
end

program Estimate, sortpreserve eclass byable(recall)
	local globopts `"Level(cilevel) PPARSEPOST"'
	local globopts `"`globopts' REML MLe VARiance EMITERate(integer 20)"'
	local globopts `"`globopts' EMTOLerance(real 1e-10) EMLOG EMDOTs"'
	local globopts `"`globopts' noGRoup noHEADer ESTMetric"'
	local globopts `"`globopts' noSTDerr noLRtest noLOg EMONLY NOEMSHOW"'
	local globopts `"`globopts' noFETable noRETable"'

	// undocumented, for use by predict and estat

	local globopts `"`globopts' grouponly getblups(string)"'

	// mlopts

	local mlopts `"NONRTOLerance NRTOLerance(string)"'
	local mlopts `"`mlopts' TRace GRADient HESSian showstep"'
	local mlopts `"`mlopts' technique(string) SHOWNRtolerance"'
	local mlopts `"`mlopts' ITERate(string) TOLerance(string)"'
	local mlopts `"`mlopts' LTOLerance(string) GTOLerance(string)"'
	local mlopts `"`mlopts' DIFficult"'

	local globallow `"`globopts' `mlopts'"'

	// parse

	_parse expand cmd glob : 0 , common(`globallow')

	forvalues k=1/`cmd_n' {				// allow if/in anywhere
		local 0 `"`cmd_`k''"'
		syntax [anything] [if] [in] [, * ]

		if `"`if'"' != `""' {
			if `"`glob_if'"' != `""' {
				di as error "multple if conditions not allowed"
				exit 198
			}
			local glob_if `"`if'"'
		}
		if `"`in'"' != `""' {
			if `"`glob_in'"' != `""' {
				di as error "multiple in ranges not allowed"
				exit 198
			}
			local glob_in `"`in'"'
		}
	}

	forvalues k=0/`=`cmd_n'-1' {			// Parse subcmds
		ParseCmd `k' `cmd_`=`k'+1''
		local allnms `allnms' `levnm_`k'' `varnames_`k''
	}
	if "`depname'" == "" {
		di as err "dependent variable not specified"
		exit 198
	}
	local allnms `allnms' `depname'

	local 0 `"`glob_if' `glob_in'"'			// set estimation sample
	syntax [if] [in]
	marksample touse
	local allnms : subinstr local allnms "_all" "" , word all
	markout `touse' `allnms' , strok

	tempvar one					// processing macs
	qui gen byte `one' = 1 if `touse'
	local lev_beg 0
	local last_levnm "______"
	tempname lev_nvars sublevels
	if (`cmd_n' > 1)  matrix `sublevels' = J(1,`=`cmd_n'-1', 0)

	forvalues k=0/`=`cmd_n'-1' {			// reconcile levels
		local k1 = `k' - 1			// and sublevels

		if "`levnm_`k''" != "`last_levnm'" {

			if `k' > 0 {			// level var repeated
				if 0`:list levnm_`k' in levnms' {	
					if `"`levnm_`k''"' != `""' {
						RepeatedLevelError `levnm_`k''
						exit 198
					}
				}
				local levnms `levnms' `levnm_`k''

				FixupConstant `lev_beg' `k1' `last_const'
				SetSublevels  `lev_beg' `k1' `sublevels'
				RMcolLevel `lev_beg' `k1' `depname'	    ///
					   "`lev_constant'" "`lev_isfctrs'" ///
					    `touse' `lev_varnames'
			}

			local lev_varnames `""`varnames_`k''""'
			if "`constant_`k''" != "" {
				local lev_constant 0
			}
			else {
				local lev_constant 1
			}
			local last_const   = cond("`constant_`k''"=="", `k', -1)
			local lev_isfctrs  `isfctr_`k''
			local lev_beg      `k'
		}
		else {
			local lev_varnames `"`lev_varnames' "`varnames_`k''""'
			local lev_isfctrs `lev_isfctrs' `isfctr_`k''

			if "`constant_`k''" != "" {
				local lev_constant `lev_constant' 0
			}
			else {
				local lev_constant `lev_constant' 1
			}
			if "`constant_`k''" == "" {
				local last_const = `k'
			}
		}

		local last_levnm `levnm_`k''
		local kh `k'
	}

	FixupConstant `lev_beg' `kh' `last_const'
	SetSublevels  `lev_beg' `kh' `sublevels'
	RMcolLevel `lev_beg' `kh' `depname' "`lev_constant'"		///
		   "`lev_isfctrs'" `touse' `lev_varnames'

	tempname fctrlevs

	forvalues k=0/`=`cmd_n'-1' {
		if `isfctr_`k'' {			// create factors
			tempvar ifctr
			qui egen long `ifctr' = group(`varnames_`k'') if `touse'
			label variable `ifctr'				///
			      `"R.`:subinstr local varlist_`k' " " " R."'"'
			local varlist_`k' `ifctr'

			sum `ifctr' , meanonly
			matrix `fctrlevs' = (nullmat(`fctrlevs'), r(max))
		}
		else {					// handle ts ops
			tsrevar `varnames_`k'' if `touse'
			local varlist_`k' `r(varlist)'
		}

		if "`constant_`k''" == "" {		// handle constants
			local varlist_`k' `varlist_`k'' `one'
			local varnames_`k' `varnames_`k'' _cons
		}

		if `k' > 0 & !`isfctr_`k'' {
			foreach var of local varlist_`k' {
				matrix `fctrlevs' = (nullmat(`fctrlevs'), 1)
			}
							// Fixup covariance
			if (`:list sizeof varlist_`k'' + 		///
			    ("`constant_`k''"=="")) == 1 {
				local cov_`k' = "identity"
			}

		}

	}
	tsrevar `depname' if `touse'			// handle ts depvar
	local depvar `r(varlist)'

	foreach nm of local levnms {			// encode levels
		tempvar level 
		if ("`nm'" == "_all")  local nm `one'
		local nms `nms' `nm'
		qui egen long `level' = group(`nms') if `touse'
		local levelvars `levelvars' `level'
	}

	forvalues k=1/`=`cmd_n'-1' {		// check for empty levels
		if "`varlist_`k''" == "" {
			if !`k' { 
				local lab fixed effects
			}
			else {
				local lab group `levnm_`k'':
			}
			di as error "`lab' equation empty"
			exit 102
		}
	}

	// global options are in glob_op  -- your parsing goes here
	local 0 `", `glob_op'"'
	syntax [ , `globopts' *]

	mlopts mlopts, `options'
	local 0 `", `mlopts'"'
	syntax [, technique(string) *]
	if `:list posof "bhhh" in technique' {
		di as err "option technique(bhhh) not allowed"
		exit 198
	}

// error checking options goes here

	if "`mle'" != "" {
		local ml ml
	}
	if "`ml'" != "" {
		if "`reml'" != "" {
			di as err "ml and reml may not be specified at " _c
			di as err "the same time"
			exit 198
		}
		local reml 0
		local method ML
	}
	else {
		local reml 1
		local method REML
	}
	if "`emonly'" != "" & "`emdots'" == "" & "`noemshow'" == "" {   
		local emlog emlog
	}
	if "`emlog'" != "" & "`emdots'" != "" {
		di as err "emlog and emdots may not be specified together"
		exit 198
	}
	if "`variance'" != "" & "`estmetric'" != "" {
		di as err "variance and estmetric may not be specified together"
		exit 198
	}

	if "`pparsepost'" != "" {
		ereturn clear
		foreach mac in depvar depname cmd_n {
			ereturn local `mac' ``mac''
		}
		forvalues l = 0/`=`cmd_n'-1' {
			foreach mac in levnm isfctr cov	k sublevels	///
				       constant varlist varnames {
				ereturn local `mac'_`l' ``mac'_`l''
			}
		}
		exit
	}
							// set up _xtm_model
	tempvar nn
	qui gen long `nn' = _n if `touse'
	gsort -`touse' `levelvars' `nn'

	tempname obs re_n
	qui count if `touse'
	scalar `obs' = r(N)
	scalar `re_n' = `cmd_n' - 1
	
	mata: _xtm_setup()

	// obtain blups if specified and exit
	if `"`getblups'"' != "" {
		CheckBlups `getblups'
		SetBlupMats `touse'
		GetBlups `getblups'
		exit	
	}
		
	if `reml' {
		local crittype log restricted-likelihood
	}
	else {
		local crittype log likelihood
	}
	if "`log'" != "" {
		local q qui
	}

	tempname delta deltaf theta thetaf gamma s h err beta Vbeta V

// fit reference regression model
	
	if `cmd_n' == 1 | "`lrtest'" == "" {
		MyReg `depvar' `varlist_0' if `touse', reml(`reml')
		local ll_0 `r(ll0)'
		local gam `r(lnsig)'
	}

	if `cmd_n' == 1 { 	// model is linear regression 
		local ll_2 `ll_0'
		mat `beta' = e(b)    // from regress
		mat `Vbeta' = e(V)  
		local c = colsof(`Vbeta')
		if !`reml' & `c'{         // ml uses different sigma^2
			mat `Vbeta' = `Vbeta' * e(df_r)/e(N)
		}
		mat `gamma' = `gam'
		mat colnames `gamma' = lnsig_e:_cons
		mat `V' = 0
		if "`stderr'" == "" {
			mat `V' = 1/(2*e(N))
		}
		local lrtest nolrtest
		local emconv 1
	}
	else {
// starting values

		mata: _xtm_start("`delta'",0.375)   

// EM iterations
		if `emiterate' {
		   `q' IterateEM `delta', crit(`crittype') iter(`emiterate') ///
			tol(`emtolerance') `emlog' `emdots' 
		    local emconv = r(converged)
		}
		mat `deltaf' = `delta'

// Transform to NR parameterization
		mata: _xtm_delta_to_theta("`theta'","`delta'","`err'")
		if scalar(`err') {
			di as err "initial values/EM solution not feasible"
			exit 430 
		}
		mat `thetaf' = `theta'
		mata: _xtm_ml_eqlist("eqlist")

		local df_c : word count `eqlist'
		mat `h' = J(1,`df_c',1)
	
		if "`emonly'" == "" {
			`q' di
			`q' di as txt "Performing gradient-based optimization: " 
			cap noi `q' ml model d0 mixed_ll `eqlist', ///
				max init(`theta', copy) ///
				search(off) crittype(`crittype') dzeros(`h') ///
				collinear missing nopreserve `mlopts' 
			if _rc {	
				if _rc == 1 {
					exit 1
				}
				local noconv noconv
				mat `delta' = `deltaf'
				mat `theta' = `thetaf'
			}
			else {
				mat `theta' = e(b)
				mata: _xtm_theta_to_delta("`delta'","`theta'","`err'")	
				local opt `e(opt)'	
				local ml_method `e(ml_method)'
				local technique `e(technique)'
			}
		}
		scalar `s' = -1 
		mata: _xtm_mixed_ll("`delta'","`s'",0)
		local ll_1 = r(ll)
		scalar `s' = r(mse)
		local ll_2 = `ll_1'

// Reparameterize
		mata: _xtm_theta_to_gamma("`gamma'","`theta'","`s'")
		if "`noconv'" == "" & "`emonly'" == "" & "`stderr'" == "" {
			mat `h' = e(ml_d0_s), 1
			local c = colsof(`theta')
			forvalues i = 1/`c' {
			    mat `h'[1,`i'] = `h'[1,`i'] * ///
				abs(`gamma'[1,`i']/`theta'[1,`i'])
			}
			`q' di 
			`q' di as txt "Computing standard errors:" 
			cap ml model d0 mixed_ll_reparm ///
			    `eqlist' /lnsig_e, ///
			    max init(`gamma',copy) search(off) iter(0) ///
			    dzeros(`h') collinear missing nopreserve `mlopts'
			if _rc {
				if _rc == 1 {
					exit 1
				}
				local nostd nostd
				local ll_2 = `ll_1'
			}
			else {
				local ll_2 = e(ll)
				if reldif(`ll_1',`ll_2') > 1e-10 {
					local nostd nostd
				}
			}
			if "`nostd'" != "" {
				di as err "standard error calculation failed"
			}
		}

// Get final beta
        	mata: _xtm_mixed_ll("`delta'","`s'",1) 
		mata: _xtm_beta("`beta'","`Vbeta'","`s'")
		if "`noconv'"=="" & "`emonly'"=="" & ///
			"`nostd'" == "" & "`stderr'"=="" {
			mat `gamma' = e(b)
			mat `V' = e(V)	
		}
		else {
			local eqnames : subinstr local eqlist "/" "", all
			mat coleq `gamma' = `eqnames' lnsig_e
			mat colnames `gamma' = _cons
			local c = colsof(`gamma')
			mat `V' = J(`c',`c',0)
		}
	}

// Post Parameters

	if `"`varnames_0'"' == "" {
		PostParms `depname', re(`gamma') vre(`V') touse(`touse')
	}
	else {
		PostParms `depname' `varnames_0', beta(`beta') ///
			vbeta(`Vbeta') re(`gamma') vre(`V') touse(`touse')
	}

// Out of the variance estimates, count how many sigmas and rhos 

	mat `gamma' = e(b)
	local eqs : coleq `gamma'
	local k_rs 0
	forval i = `=e(k_f)+1'/`e(k)' {
		local eq : word `i' of `eqs'
		if substr("`eq'",1,3) == "lns" {
			local ++k_rs
		}
	}
	ereturn scalar k_rs = `k_rs'
	ereturn scalar k_rc = e(k_r) - e(k_rs)

	if "`emonly'" != "" {
		ereturn local emonly emonly
		ereturn scalar converged = `emconv'
	}
	else {
		ereturn scalar converged = ("`noconv'"=="")
	}
	ereturn local depvar `depname'
	ereturn scalar ll = `ll_2'

// Comparison test
	
	if `cmd_n' > 1 & "`lrtest'" == "" {
		ereturn scalar ll_c = `ll_0'
		ereturn scalar df_c = `df_c'
		ereturn scalar chi2_c = 2*(e(ll) - e(ll_c))
		if e(chi2_c) < 0 {
			ereturn scalar chi2_c = 0
		}
		ereturn scalar p_c = chi2tail(e(df_c),e(chi2_c))
		if e(df_c) == 1 & e(chi2_c) > 1e-5 {
			ereturn scalar p_c = 0.5*e(p_c)
		}
	}

//  data resorted in what follows

	SaveGroupInfo `"`levelvars'"' 

	ereturn local crittype `crittype'
	ereturn local method `method'
	ereturn local opt `opt'
	ereturn local ml_method `ml_method'
	ereturn local technique `technique'
	ereturn local title Mixed-effects `e(method)' regression

	forvalues k=1/`=`cmd_n'-1' {
		if `isfctr_`k'' {
			local varnames_`k'				///
			      `"R.`:subinstr local varnames_`k' " " " R."'"'
		}
		local zvars     `"`zvars' `varnames_`k''"'
		local vartypes  `"`vartypes' `=proper("`cov_`k''")'"'
		local lvnms     `"`lvnms' `levnm_`k''"'
		local dimz	`"`dimz' `:list sizeof varnames_`k''"'
	}
	ereturn local revars `zvars'
	ereturn local vartypes `vartypes'
	ereturn local ivars `lvnms'
	ereturn local redim `dimz'
	ereturn local cmd xtmixed
	ereturn local predict xtmixed_p
	ereturn local estat_cmd xtmixed_estat

	if "`varnames_0'" == "_cons"  | "`varnames_0'" == "" {
		ereturn scalar chi2 = .
		ereturn scalar p = .
		ereturn scalar df_m = 0
	}
	else {