xtmixed.ado

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

		   _col(73) as res %6.4f e(p_c)
		return local conserve conserve
	}
end

program DiVarComp, rclass
	syntax [, level(cilevel) VARiance noLRtest]

	local depvar `e(depvar)'
	if strpos("`depvar'",".") {
		gettoken ts rest : depvar, parse(".")
		gettoken dot depvar : rest, parse(".")
	}
	
	local dimx `e(k_f)'

	// display header 

	di
	di as txt "{hline 29}{c TT}{hline 48}
	local k = length("`level'")
	di as txt _col(3) "Random-effects Parameters" _col(30) "{c |}" ///
		_col(34) "Estimate" _col(45) "Std. Err." _col(`=61-`k'') ///
		`"[`=strsubdp("`level'")'% Conf. Interval]"'
	di as txt "{hline 29}{c +}{hline 48}

	local zvars `e(revars)'
	local dimz `e(redim)'

	// loop over levels

	local foot = 1
	local pos `dimx'
	local levs : word count `e(ivars)'
	forvalues k = 1/`levs' {
		local lev : word `k' of `e(ivars)'
		local vartype : word `k' of `e(vartypes)'
		di as res abbrev("`lev'",12) as txt ": `vartype'" ///
			_col (30)"{c |}"

		GetNames "`vartype'" "`zvars'" "`dimz'" `foot'
		local zvars `s(zvars)'     // collapsed lists
		local dimz `s(dimz)'
		local names `"`s(names)'"'
		if `"`s(footnote)'"' != "" {
			local footnotes `"`footnotes' `s(footnote)'"'
			local ++foot
		}

		local nbeta : word count `names'
	
		DiParms `pos' `nbeta' `"`names'"' `level' `variance'

		local pos = `pos' + `nbeta'
		di as txt "{hline 29}{c +}{hline 48}
	}

	// Residual Variance 

	if "`variance'" == "" {
		local parm exp
		local label sd
		local p 17
	}
	else {
		local parm f(exp(2*@)) d(2*exp(2*@))
		local label var
		local p 16
	}
	qui _diparm lnsig_e, `parm' level(`level') notab
	di as txt _col(`p') "`label'(Residual)" _col(30) "{c |}" ///
	   as res _col(33) %9.0g r(est) ///
	   as res _col(44) %9.0g r(se)  ///
	   as res _col(58) %9.0g cond(missing(r(se)),.,r(lb))  ///
	   as res _col(70) %9.0g cond(missing(r(se)),.,r(ub))
	di as txt "{hline 29}{c BT}{hline 48}
	
	// Footnotes
	if "`e(chi2_c)'" != "" & "`lrtest'" == "" {
		DiLRTest
		local conserve `r(conserve)'
	}
	
	forvalues k = 1/`=`foot'-1' {
		local note: word `k' of `footnotes'
		local indent = 3 + length("`k'")
		di as txt `"{p 0 `indent' 4} `note'{p_end}"'
	}
	if "`conserve'" != "" {
		di as txt _n "{p 0 6 4}Note: {help j_xtmixedlr##|_new:LR test is conservative} and provided only for reference{p_end}"
		return local note note
	}
end

program DiParms
	args pos nb names cilev var

	local stripes : coleq e(b), quoted
	forvalues k = 1/`nb' {
		local label : word `k' of `names'	
		local eq : word `=`pos'+`k'' of `stripes'
		GetParmEqType `eq' `var'
		local parm `s(parm)'
		local label `"`s(type)'`label'"'	
		local diparmeq `"`s(diparmeq)'"'

		local p = 29 - length("`label'")
		qui _diparm `diparmeq', `parm' level(`cilev') notab
		di as txt _col(`p') "`label'" _col(30) "{c |}" ///
		   as res _col(33) %9.0g r(est) ///
		   as res _col(44) %9.0g r(se)  ///
		   as res _col(58) %9.0g cond(missing(r(se)),.,r(lb))  ///
		   as res _col(70) %9.0g cond(missing(r(se)),.,r(ub))
	}
end

program GetNames, sclass
	args type zvars dimz foot

	gettoken dim dimz : dimz
	forvalues k = 1/`dim' {
		gettoken tok1 zvars : zvars
		local fullvarnames `fullvarnames' `tok1'
		local len = length("`tok1'")
		if substr("`tok1'",1,2) == "R." {
			local w = substr("`tok1'",3,`len') 
			local tok1 = "R." + abbrev("`w'",8)
		}
		else {
			local tok1 = abbrev("`tok1'",8)
		}
		local varnames `varnames' `tok1'
	}
	if ("`type'" == "Unstructured") {
		forvalues j = 1/`dim' {
			local w : word `j' of `varnames'
			local names `"`names' "(`w')""'
		}
		forvalues j = 1/`dim' {
			forvalues k = `=`j'+1'/`dim' {
				local w1 : word `j' of `varnames'
				local w2 : word `k' of `varnames'
				local names `"`names' "(`w1',`w2')""'	
			}
		}
	}
	else if ("`type'" == "Independent") {    
		forvalues j = 1/`dim' {
			local w : word `j' of `varnames'
			local names `"`names' "(`w')""'
		}
	}
	else {    // Identity or Exchangeable
		local ex = ("`type'" == "Exchangeable")
		if (`dim' == 1) {		// check for factor variable
			local w `varnames'
			local names `""(`w')""'
			if `ex' {
				local names `"`names' "(`w')""'
			}
		}
		else if (`dim' == 2) {
			local w1 : word 1 of `varnames'
			local w2 : word 2 of `varnames'
			local names `""(`w1' `w2')""'
			if `ex' {
local names `"`names' "(`w1',`w2')""'
			}
		}
		else {
			local k : length local varnames
			if `k' > 20 {		// too long
				local w1 : word 1 of `varnames'	
				local w2 : word `dim' of `varnames'
				local names `""(`w1'..`w2')(`foot')""'
				if `ex' {
local names `"`names' "(`w1'..`w2')(`foot')""'
				}
				local footnote `""(`foot') `fullvarnames'""'
			}
			else {
				local names `""(`varnames')""'
				if `ex' {
local names `"`names' "(`varnames')""'
				}
			}
		}
	}

	sreturn local zvars "`zvars'"
	sreturn local dimz "`dimz'"
	sreturn local names `"`names'"'
	sreturn local footnote `"`footnote'"'
end

program GetParmEqType, sclass
	args eq var
	
	if substr("`eq'",1,1) == "l" {	  // log standard deviation
		if "`var'" == "" {	  // se/corr metric
			local parm exp
			local type sd
		}			
		else {			  // var/cov metric
			local parm f(exp(2*@)) d(2*exp(2*@))
			local type var
		}
		local deq `eq'
	}
	else { // substr("`eq'",1,1) == "a"  atanh correlation
		if "`var'" == "" {        // se/corr metric
			local parm tanh
			local type corr
			local deq `eq'
		}
		else {			  // var/cov metric
			ParseEq `eq'
			local eq2 lns`r(n1)'_`r(n2)'_`r(n3)'
			local eq3 lns`r(n1)'_`r(n2)'_`r(n4)'
			local deq `eq' `eq2' `eq3'	
			local parm f(tanh(@1)*exp(@2+@3))
			local parm `parm' d((1-(tanh(@1)^2))*exp(@2+@3)
			local parm `parm' tanh(@1)*exp(@2+@3)
			local parm `parm' tanh(@1)*exp(@2+@3)) 
			local type cov
		}
	}

	sreturn local parm `"`parm'"'
	sreturn local type `"`type'"'
	sreturn local diparmeq `"`deq'"'
end

program ParseEq, rclass
	args eq

	// I've got "eq" == "atr#_#_#_#", and I need the four #'s
	// returned as r(n1), r(n2), r(n3), and r(n4)

	local len = length("`eq'")
	local eq = substr("`eq'",4,`len')
	forvalues k = 1/4 {
		gettoken n`k' eq : eq, parse(" _")
		return local n`k' `n`k''
		gettoken unscore eq : eq, parse(" _")
	}
end

program SaveGroupInfo, eclass
	args levnames

	qui count if e(sample)
	ereturn scalar N = `r(N)'
	local levels : list uniq levnames
	local k : word count `levels'
	if `k' == 0 {
		exit
	}
		
	tempname gmin gmax gavg Ng
	mat `Ng' = J(1,`k',0)
	mat `gmin' = J(1,`k',0)
	mat `gavg' = J(1,`k',0)
	mat `gmax' = J(1,`k',0)

	forvalues i = 1/`k' {
		GroupStats `:word `i' of `levels'' if e(sample)	
		mat `Ng'[1,`i'] = r(ng)
		mat `gmin'[1,`i'] = r(min)
		mat `gavg'[1,`i'] = r(avg)
		mat `gmax'[1,`i'] = r(max)
		local ++i
	}
	ereturn matrix N_g `Ng'
	ereturn matrix g_min `gmin'
	ereturn matrix g_avg `gavg'
	ereturn matrix g_max `gmax'
end

program GroupStats, rclass
	syntax name [if]
	marksample touse

	if "`namelist'" == "_all" {
		tempvar one
		qui gen byte `one' = 1 if `touse'
		local namelist `one'
	}
	tempname T
	qui {
		bysort `touse' `namelist': /// 
			gen long `T' = cond(_n==_N,_N,.) if `touse'
		sum `T' if `touse'
	}
	return scalar ng = r(N)
	return scalar min = r(min)
	return scalar max = r(max)
	return scalar avg = r(mean)
end

program IterateEM, rclass
	syntax name(name=theta) [, iter(integer 1) tol(real 1.0) ///
			crit(string) emlog emdots ]

	tempname thetaold s lold lnew delta err
	mat `thetaold' = `theta'

	if "`emlog'" == "" {
		local st *
	}
	di as txt _n "Performing EM optimization: " _c
	if "`emdots'" != "" {
		local didot `"di as txt "." _c"'
	}
	else {
		di
		`st' di
	}
	scalar `lold' = -1
	scalar `lnew' = 1
	
	local i 0	
	while `i' <= `iter'-1 & ///
	   reldif(scalar(`lold'),scalar(`lnew')) > `tol' {
		mata: _xtm_em_iter("`thetaold'","`theta'","`s'")
		if missing(r(ll)) {
			if "`emdots'" != "" {
				di _n
			}
			di as err "likelihood evaluates to missing"
			exit 430
		}
`st'di as txt "Iteration `i':" _col(16) "`crit' = " as res %10.0g `r(ll)'
		`didot'
		mat `thetaold' = `theta'
		scalar `lold' = scalar(`lnew')
		scalar `lnew' = r(ll)
		local ++i
	}
	scalar `s' = -1
	mata: _xtm_mixed_ll("`theta'","`s'",0)
`st'di as txt "Iteration `i':" _col(16) "`crit' = " as res %10.0g `r(ll)'
	if "`emdots'" != "" {
		di 
	}
	return scalar converged = (`i'<`iter')
end

program MyReg, rclass
	syntax anything(name=model) [if] [, reml(integer 0)]
	
	gettoken depvar xvars : model

	// strip _cons off the end of xvars
	local k: word count `xvars'
	if `k' {
		if "`:word `k' of `xvars''" == "_cons" {
			local xvars : subinstr local xvars "_cons" ""
		}
		else {
			local nocons nocons
		}
	}
	else {
		local nocons nocons
	}
	qui regress `depvar' `xvars' `if', `nocons'
	tempname lnsig ll0 logdet fact
	if `reml' {	
		mata: _xtm_logdetr00("`logdet'")
		scalar `lnsig' = log(e(rmse))
		scalar `ll0' = -e(df_r)/2*(log(2*_pi) + 1 + 2*`lnsig') ///
			       -`logdet'
	}
	else {		
		scalar `lnsig' = log(e(rss)/e(N))/2
		scalar `ll0' = e(ll)
		
	}
	return scalar lnsig = `lnsig'
	return scalar ll0 = `ll0'
end

program SetBlupMats 
	args touse
	tempname gamma theta delta err s2

	if `"`e(cmd)'"' != "xtmixed" {
di as err "BLUP calculations require a previously fitted {cmd:xtmixed} model"
		exit 301
	}

	qui count if `touse'
	if r(N) != e(N) {
	    di as err "BLUP calculation failed; estimation data have changed"
	    exit 459
	}

	mat `gamma' = e(b)
	local p = colsof(`gamma')
	scalar `s2' = exp(2*`gamma'[1,`p'])
	mat `gamma' = `gamma'[1,`=e(k_f)+1'...]

	local k 0
	capture {
		 mata: _xtm_gamma_to_theta("`theta'","`gamma'","`err'")
		 local k = scalar(`err')
		 mata: _xtm_theta_to_delta("`delta'","`theta'","`err'")
		 local k = `k' + scalar(`err')
	}
	if `k' | _rc {
		di as err "BLUP calculation failed; e(b) unsuitable"
		exit 459
	}
	local ll 0
	capture {
		mata: _xtm_mixed_ll("`delta'","`s2'",2)
		local ll = r(ll)
	}
	if _rc | reldif(`ll',e(ll)) > 1e-12 {
	    di as err "BLUP calculation failed; estimation data have changed"
	    exit 459
	}
end

program CheckBlups
	local k : word count `0'
	if mod(`k',2) {
		di as err "getblups() incorrectly specified"
		exit 198
	}
	forvalues i = 1/`=`k'/2' {
		gettoken level 0 : 0
		gettoken stub 0 : 0
		if `:list stub in stublist' {
			di as err "stubs must be unique"
			exit 198
		}
		local stublist `stublist' `stub'
		if `"`level'"' != "_all" {
			unab level : `level'
		}
		local ivars `e(ivars)'
		local p : list posof "`level'" in ivars
		if !`p' {
			di as err `"`level' is not a group variable"'
			exit 198
		}
		local m 0
		forvalues j = `p'/`:word count `ivars'' {
			if `"`:word `j' of `ivars''"' != "`level'" {
				continue, break
			}
			local np : word `j' of `e(redim)'	
			forvalues k = 1/`np' {
				confirm new variable `stub'`=`m'+`k''
			}
			local m = `m' + `np'
		}
	}
end

program GetBlups
	mata: _xtm_bmat = _xtm_blup(0) // yes, _xtm_bmat is external

	local k : word count `0'
	forvalues i = 1/`=`k'/2' {
		gettoken level 0 : 0
		gettoken stub 0 : 0
		if `"`level'"' != "_all" {
			unab level : `level'
		}
		GetBlupsLevel `level' `stub'		
	}
end

program GetBlupsLevel
	args lvar stub

	local ivars `e(ivars)'
	local zvars `e(revars)'

	// get z variables
	local w : word count `ivars'
	forvalues i = 1/`w' {
		local lev : word `i' of `ivars'
		local dim : word `i' of `e(redim)'
		if "`lev'" != "`lvar'" {
			local dont *
		}
		forvalues j = 1/`dim' {
			gettoken z zvars : zvars 
			`dont'local zs `zs' `z'
		}
		local dont
	}

	local uivars : list uniq ivars
	local p : list posof "`lvar'" in uivars 	// p is the true level
	local nz : word count `zs'

	forvalues i = 1/`nz' {
		qui gen double `stub'`i' = . 
		label var `stub'`i' `"BLUP r.e. for `lvar': `:word `i' of `zs''"'
	}
	mata: _xtm_blup_save("`stub'",`p')
end

exit